[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
_______________________________________________
SciPy-User mailing list
SciPy...@scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-user

[josef...@gmail.com]

2009/11/8 Ernest Adrogué <eadr...@gmx.net>:

Thanks, I think the distribution of the difference of two poisson
distributed random variables could be useful.

Would you please open an enhancement ticket for this at
http://projects.scipy.org/scipy/report/1

I had only a brief look at it so far, I had never looked at the
Skellam distribution before, and just read a few references.

The "if x < 0 .. else ..." will have to be replace with a
"numpy.where" assignment, since the methods are supposed to work with
arrays of x (as far as I remember)

_rvs could be implemented directly instead of generically (I don't
find the reference, where I saw it, right now).

Documentation will be necessary, a brief description in the
(currently) extradocs, and a listing of the properties for the
description of the distributions currently in the stats tutorial.

I have some background questions, which address the limitation of the
implementation (but are not really necessary for inclusion into
scipy).

The description in R mentions several implementation of Skellam. Do
you have a rough idea what the range of parameters are for which the
implementation using ncx produces good results? Do you know if any
other special functions would produce good results over a larger
range, e.g. using Bessel function?

Wikipedia, http://en.wikipedia.org/wiki/Skellam_distribution , also
mentions (but doesn't describe) the case of Skellam distribution with
correlated Poisson distributions. Do you know what the difference to
your implementation would be?

Tests for a new distribution will be picked up by the generic tests,
but it would be useful to have some extra tests for extreme/uncommon
parameter ranges. Do you have any comparisons with R, since you
already looked it?

Thanks again, I'm always looking out for new useful distributions,
(but I have to find the time to do the testing and actual
implementation).

Josef

[Bruce Southey]

Generally any R code can not be used in numpy because R is GPL. Usually
R code is also licensed under GPL so translation from R to Python/numpy
still maintains the original license. So the code not used by numpy
unless that code is licensed under a BSD compatible license.

You *must* show that you implementation is from a BSD-compatible source
not from the R package. I can see that your code is very simple so there
should be an viable alternative source.

Also, in the _stats function why do you do not re-use the mean and var
variables in computing the g1 and g2 variables?

What are 'x, mu1, mu2' ?
This looks like a scalar implementation so you need to either check that
or allow for array-like inputs.

Bruce

[Johann Cohen-Tanugi]

From what I understand of the initial statement from Ernest:

"
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
"

he used the spec, as defined in this pdf, and did not look at the code
itself. If my interpretation of the small preamble above is correct, I
believe his implementation is not GPL-tainted, right?

Johann

[josef...@gmail.com]

We only need to read the R documentation and not the R code:

pskellam(x,lambda1,lambda2) returns pchisq(2*lambda2, -2*x, 2*lambda1)
for x <= 0 and 1 - pchisq(2*lambda1, 2*(x+1), 2*lambda2) for x >= 0. When
pchisq incorrectly returns 0, a saddlepoint approximation is
substituted, which typically gives at
least 2-figure accuracy.
The quantile is defined as the smallest value x such that F(x) p,
where F is the distribution
function. For lower.tail=FALSE, the quantile is defined as the largest
value x such that
F(x;lower.tail=FALSE) p.
rskellam is calculated as rpois(n,lambda1)-rpois(n,lambda2)
dskellam.

and

The relation of dgamma to the modified Bessel function of the first
kind was given by Skellam
(1946). The relation of pgamma to the noncentral chi-square was given
by Johnson (1959). Tables
are given by Strackee and van der Gon (1962), which can be used to
verify this implementation (cf.
direct calculation in the examples below).

the rest follows from the Wikipedia page (which is also in the list of
references in R docs), there is no copyright on the definition of a
distribution.

>
> Also, in the _stats function why do you do not re-use the mean and var
> variables in computing the g1 and g2 variables?
>
> What are 'x, mu1, mu2' ?

x is the integer at which cdf pr pmf are calculated, mu1,mu2 are the
parameters of the poisson distributions, following wikipedia.

Josef

> This looks like a scalar implementation so you need to either check that
> or allow for array-like inputs.
>
> Bruce
>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy...@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>

[Bruce Southey]

I am not a lawyer! But I do not see that any reference to not seeing the
code. Furthermore, there is insufficient information in the cited
reference for this implementation (but I have not seen the actual code
and would rather not have to see it). But, as Josef pointed out, there
is a Wikipedia source so it should be trivial to show that this code is
independent of the R implementation.

[josef...@gmail.com]

same form different interpretation

"The distributions of the di erence between two independent and two
bivariate (correlated)
Poisson variates are of the same form. However, the interpretation of
the parameters is different.
Assuming that the bivariate Poisson distribution is the correct
distribution, then the
marginal means x and y will be unbiased estimates of 1 + 3 and 2
+ 3, respectively,
instead of the parameters of interest 1 and 2. Therefore, the
parameters of the PD distribution
are not directly connected to the marginal means of the actual Poisson
distributions."

from
Bayesian analysis of the di erences of count data
D. Karlis and I. Ntzoufras
STATISTICS IN MEDICINE
Statist. Med. 2006; 25:1885–1905

they have some funny application to soccer scores
http://stat-athens.aueb.gr/~jbn/publications.htm

Josef
Published online 26 October 2005 in Wiley InterScience
(www.interscience.wiley.com). DOI: 10.1002/sim.2382

[Ernest Adrogué]

9/11/09 @ 09:40 (-0500), thus spake josef...@gmail.com:

> Thanks, I think the distribution of the difference of two poisson
> distributed random variables could be useful.
>
> Would you please open an enhancement ticket for this at
> http://projects.scipy.org/scipy/report/1

Okay, I will.

> I had only a brief look at it so far, I had never looked at the
> Skellam distribution before, and just read a few references.
>
> The "if x < 0 .. else ..." will have to be replace with a
> "numpy.where" assignment, since the methods are supposed to work with
> arrays of x (as far as I remember)

Done.

> _rvs could be implemented directly instead of generically (I don't
> find the reference, where I saw it, right now).

I suppose it could be done, but I can't figure out how.
As far as I understand, you can't derive the poisson parameters
from the skellam parameters alone, you need the correlation
coefficient (between the two poisson rv) too, is that right?

> Documentation will be necessary, a brief description in the
> (currently) extradocs, and a listing of the properties for the
> description of the distributions currently in the stats tutorial.
>
> I have some background questions, which address the limitation of the
> implementation (but are not really necessary for inclusion into
> scipy).
>
> The description in R mentions several implementation of Skellam. Do
> you have a rough idea what the range of parameters are for which the
> implementation using ncx produces good results? Do you know if any
> other special functions would produce good results over a larger
> range, e.g. using Bessel function?

No, sorry, I have no idea. The R paper says that in R the Bessel
function is more accurate, but I don't think this means that the
Bessel function implementation is more accurate in general.
I compared results using the Bessel function in scipy and using
the ncx2 implementation, and the differences were minimal, although
I admit my testing wasn't extensive. Another problem is that
I haven't got a table of values for this particular distribution,
so in case results differ I can't tell which one is more accurate.

> Wikipedia, http://en.wikipedia.org/wiki/Skellam_distribution , also
> mentions (but doesn't describe) the case of Skellam distribution with
> correlated Poisson distributions. Do you know what the difference to
> your implementation would be?

As far as I know, it makes no difference if the two Poisson
variates are correlated or not. The Skellam parameters are
defined as

mu1 = lam1 - rho * sqrt(lam1*lam2)
mu2 = lam2 - rho * sqrt(lam1*lam2)

where lam1 and lam2 are the Poisson means and rho is the
correlation coefficient. So, in case there is correlation it
is implicit in the parameters mu1 and mu2, and it doesn't
make any difference in terms of calculating the pmf or cdf.
Again, I am no statician or mathematician.

> Tests for a new distribution will be picked up by the generic tests,
> but it would be useful to have some extra tests for extreme/uncommon
> parameter ranges. Do you have any comparisons with R, since you
> already looked it?

I have this visual test, although it's usefulness is
questionable. It shows the deviation between observed and
expected frequencies for a Skellam random variable. Ideally
errors should be random and centered around 0. To me
it looks like errors tend to increase around the mean
and are smaller along the tails, but I don't know if this
means something or not.

import numpy
import scipy.stats.distributions

poisson = scipy.stats.distributions.poisson
ncx2 = scipy.stats.distributions.ncx2

# Skellam distribution

class skellam_gen(scipy.stats.distributions.rv_discrete):
def _pmf(self, x, mu1, mu2):

px = numpy.where(x < 0, ncx2.pdf(2*mu2, 2*(1-x), 2*mu1)*2,
ncx2.pdf(2*mu1, 2*(x+1), 2*mu2)*2)

return px
def _cdf(self, x, mu1, mu2):
x = numpy.floor(x)

px = numpy.where(x < 0, ncx2.cdf(2*mu2, -2*x, 2*mu1),
1-ncx2.cdf(2*mu1, 2*(x+1), 2*mu2))

return px
def _stats(self, mu1, mu2):
mean = mu1 - mu2
var = mu1 + mu2

g1 = (mean) / numpy.sqrt((var)**3)
g2 = 1 / var

return mean, var, g1, g2
skellam = skellam_gen(a=-numpy.inf, name="skellam", longname='A Skellam',
shapes="mu1,mu2", extradoc="")

if __name__ == '__main__':

lam1 = 3.4
lam2 = 6.1
n = 5000

poisson_var1 = numpy.random.poisson(lam1, n)
poisson_var2 = numpy.random.poisson(lam2, n)
skellam_var = poisson_var1-poisson_var2

low = min(skellam_var)
high = max(skellam_var)
obs_freq = numpy.histogram(
skellam_var, numpy.arange(low, high+2))[0] / float(n)

rho = numpy.corrcoef(poisson_var1, poisson_var2)[1,0]
mu1 = lam1 - rho * numpy.sqrt(lam1*lam2)
mu2 = lam2 - rho * numpy.sqrt(lam1*lam2)
exp_freq = skellam.pmf(numpy.arange(low, high+1), mu1, mu2)
print obs_freq-exp_freq

# plot
import matplotlib.pyplot as plt
plt.figure().add_subplot(1,1,1).plot(obs_freq-exp_freq)
plt.show()

Regards.

[Ernest Adrogué]

9/11/09 @ 17:01 (+0100), thus spake Johann Cohen-Tanugi:

> From what I understand of the initial statement from Ernest:
> "
> 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
> "
>
> he used the spec, as defined in this pdf, and did not look at the code
> itself. If my interpretation of the small preamble above is correct, I
> believe his implementation is not GPL-tainted, right?

Your interpretation of my preamble is correct.
I didn't look at the source code, I only read the PDF doc
where they explain in plain English how the pmf and cdf are
calculated.

[Bruce Southey]

On 11/09/2009 12:20 PM, Ernest Adrogué wrote:
> 9/11/09 @ 17:01 (+0100), thus spake Johann Cohen-Tanugi:
>
>> From what I understand of the initial statement from Ernest:
>> "
>> 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
>> "
>>
>> he used the spec, as defined in this pdf, and did not look at the code
>> itself. If my interpretation of the small preamble above is correct, I
>> believe his implementation is not GPL-tainted, right?
>>
> Your interpretation of my preamble is correct.
> I didn't look at the source code, I only read the PDF doc
> where they explain in plain English how the pmf and cdf are
> calculated.
>
>

Okay,
Then just provide a suitable reference outside the R package where
someone can derive it independently or get more details. This also
allows people to go back and check the implementation when things are
not as expected. Perhaps the reference that Josef provided provides the
same formulation?

Bruce

[josef...@gmail.com]

On Mon, Nov 9, 2009 at 1:49 PM, Bruce Southey <bsou...@gmail.com> wrote:
> On 11/09/2009 12:20 PM, Ernest Adrogué wrote:
>>   9/11/09 @ 17:01 (+0100), thus spake Johann Cohen-Tanugi:
>>
>>>    From what I understand of the initial statement from Ernest:
>>> "
>>> 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
>>> "
>>>
>>> he used the spec, as defined in this pdf, and did not look at the code
>>> itself. If my interpretation of the small preamble above is correct, I
>>> believe his implementation is not GPL-tainted, right?
>>>
>> Your interpretation of my preamble is correct.
>> I didn't look at the source code, I only read the PDF doc
>> where they explain in plain English how the pmf and cdf are
>> calculated.
>>
>>
> Okay,
> Then just provide a suitable reference outside the R package where
> someone can derive it independently or get more details. This also
> allows people to go back and check the implementation when things are
> not as expected. Perhaps the reference that Josef provided provides the
> same formulation?

If someone wants to check the relationship of chisquare and difference
of two independent poisson, but that's an implementation detail and
not really informative in a doc string.

On an Extension of the Connexion Between Poisson and χ2 Distributions
Author(s): N. L. Johnson
Source: Biometrika, Vol. 46, No. 3/4 (Dec., 1959), pp. 352-363

Josef

[josef...@gmail.com]

2009/11/9 Ernest Adrogué <eadr...@gmx.net>:

>  9/11/09 @ 09:40 (-0500), thus spake josef...@gmail.com:
>> Thanks, I think the distribution of the difference of two poisson
>> distributed random variables could be useful.
>>
>> Would you please open an enhancement ticket for this at
>> http://projects.scipy.org/scipy/report/1
>
> Okay, I will.
>
>> I had only a brief look at it so far, I had never looked at the
>> Skellam distribution before, and just read a few references.
>>
>> The "if x < 0 .. else ..." will have to be replace with a
>> "numpy.where" assignment, since the methods are supposed to work with
>> arrays of x (as far as I remember)
>
> Done.
>
>> _rvs could be implemented directly instead of generically (I don't
>> find the reference, where I saw it, right now).
>
> I suppose it could be done, but I can't figure out how.
> As far as I understand, you can't derive the poisson parameters
> from the skellam parameters alone, you need the correlation
> coefficient (between the two poisson rv) too, is that right?

it's can be generated as difference between two independent
poisson, (see R docs). Correlation only matters for the
interpretation, as you explain below and I found in a reference.

>
>> Documentation will be necessary,  a brief description in the
>> (currently) extradocs, and a listing of the properties for the
>> description of the distributions currently in the stats tutorial.
>>
>> I have some background questions, which address the limitation of the
>> implementation (but are not really necessary for inclusion into
>> scipy).
>>
>> The description in R mentions several implementation of Skellam. Do
>> you have a rough idea what the range of parameters are for which the
>> implementation using ncx produces good results? Do you know if any
>> other special functions would produce good results over a larger
>> range, e.g. using Bessel function?
>
> No, sorry, I have no idea. The R paper says that in R the Bessel
> function is more accurate, but I don't think this means that the
> Bessel function implementation is more accurate in general.
> I compared results using the Bessel function in scipy and using
> the ncx2 implementation, and the differences were minimal, although
> I admit my testing wasn't extensive. Another problem is that
> I haven't got a table of values for this particular distribution,
> so in case results differ I can't tell which one is more accurate

I will look at it, but small numerical inaccuracies, 1e-??, can also be
fixed later, if someone finds a better implementation.,

>
>> Wikipedia, http://en.wikipedia.org/wiki/Skellam_distribution , also
>> mentions (but doesn't describe) the case of Skellam distribution with
>> correlated Poisson distributions. Do you know what the difference to
>> your implementation would be?
>
> As far as I know, it makes no difference if the two Poisson
> variates are correlated or not. The Skellam parameters are
> defined as
>
> mu1 = lam1 - rho * sqrt(lam1*lam2)
> mu2 = lam2 - rho * sqrt(lam1*lam2)
>
> where lam1 and lam2 are the Poisson means and rho is the
> correlation coefficient. So, in case there is correlation it
> is implicit in the parameters mu1 and mu2, and it doesn't
> make any difference in terms of calculating the pmf or cdf.
> Again, I am no statician or mathematician.

I agree, that's what I have seen in the references.

I'm not sure the example is correct. You are simulating two
independent poisson variables, so the difference skellam_var should be
distributed as skellam with
mu1 = lam1
mu2 = lam2
and theoretical rho should be zero. Or am I missing something?

thanks, the numpy.where looks good, I still have to actually run the examples.

Josef

[Ernest Adrogué]

9/11/09 @ 14:36 (-0500), thus spake josef...@gmail.com:

> I'm not sure the example is correct. You are simulating two
> independent poisson variables, so the difference skellam_var should be
> distributed as skellam with
> mu1 = lam1
> mu2 = lam2
> and theoretical rho should be zero. Or am I missing something?

Yes, rho should be zero, but in practice, there may be a
certain amount of correlation even though the variables are
theoretically independent.

If you set rho=0 when it's not actually zero, this will
result in an artificially worse fit, which sort of defeats
the purpose of this test. It would make sense to set rho=0
if we were testing whether the to Poisson variates are
independent, though.

Bye the way, I have opened a ticket here:
http://projects.scipy.org/scipy/ticket/1050

Bye.

[Ernest Adrogué]

9/11/09 @ 11:58 (-0500), thus spake josef...@gmail.com:

> from
> Bayesian analysis of the di erences of count data
> D. Karlis and I. Ntzoufras
> STATISTICS IN MEDICINE
> Statist. Med. 2006; 25:1885–1905
>
> they have some funny application to soccer scores
> http://stat-athens.aueb.gr/~jbn/publications.htm

Yes, I have working on modelling soccer scores for some
time now, and have come to the conclusion that Poisson models
are doomed, because in essence football scores are not
Poisson distributed. If you are interested in these things
let me suggest you to concentrate your efforts on the negative
binomial model, no matter how tempting the Poisson model may be.

Cheers.