# [SciPy-User] the skellam distribution

Sun Nov 8 09:16:25 CST 2009

```Hi,

In case somebody is interested, or you want to include it
in scipy. I used these specs here from the R package:
cran.r-project.org/web/packages/skellam/skellam.pdf

Note that I am no statician, somebody who knows what he's
doing (as opposed to me ;) should verify it's correct.

import numpy
import scipy.stats.distributions

# Skellam distribution

ncx2 = scipy.stats.distributions.ncx2

class skellam_gen(scipy.stats.distributions.rv_discrete):
def _pmf(self, x, mu1, mu2):
if x < 0:
px = ncx2.pdf(2*mu2, 2*(1-x), 2*mu1)*2
else:
px = ncx2.pdf(2*mu1, 2*(x+1), 2*mu2)*2
return px
def _cdf(self, x, mu1, mu2):
x = numpy.floor(x)
if x < 0:
px = ncx2.cdf(2*mu2, x*(-2), 2*mu1)
else:
px = 1-ncx2.cdf(2*mu1, 2*(x+1), 2*mu2)
return px
def _stats(self, mu1, mu2):
mean = mu1 - mu2
var = mu1 + mu2
g1 = (mu1 - mu2) / numpy.sqrt((mu1 + mu2)**3)
g2 = 1 / (mu1 + mu2)
return mean, var, g1, g2
skellam = skellam_gen(a=-numpy.inf, name="skellam", longname='A Skellam',