[SciPy-User] the skellam distribution

[Ernest Adrogué]

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',
                      shapes="mu1,mu2", extradoc="")


Bye.

-- 
Ernest