request for Poisson CDF inverse

request for Poisson CDF inverse

[Jason Detwiler]

I am looking to calculate the inverse of the CDF for the Poisson
distribution. I noticed that many of the CDFs coded in gsl have
inverse functions also coded, but this is not the case for the Poisson
distribution (unfortunately for me). I would like to "request" these
functions. I would try to contribute them myself, but I'm not
confident I could code them to the standards of gsl. What are the
chances of my wishes being fulfilled by another developer?

thanks,
Jason Detwiler

Re: request for Poisson CDF inverse

[Heikki Orsila]

On Thu, Jul 05, 2007 at 05:16:38PM -0700, Jason Detwiler wrote:
> I am looking to calculate the inverse of the CDF for the Poisson
> distribution. I noticed that many of the CDFs coded in gsl have
> inverse functions also coded, but this is not the case for the Poisson
> distribution (unfortunately for me). I would like to "request" these
> functions. I would try to contribute them myself, but I'm not
> confident I could code them to the standards of gsl. What are the
> chances of my wishes being fulfilled by another developer?

Are you looking for something like this?

	http://zakalwe.fi/~shd/code/poisson.py

This code maps a probability from range [0, 1) into a number of 
random events happening at a given invocation. This is done in 
Poisson_Process.random().

Example, 1E6 invocations, probability p = 0.05 of an event. This will 
give frequencies for each N (the number of random events happening):

$ ./poisson.py |sort -n |uniq -c |sort -n
...
     75 3
   2482 2
  47754 1
 949688 0

As a comparison:
	1E6*0.05^1 = 50000
	1E6*0.05^2 = 2500
	1E6*0.05^3 = 125

So these values are pretty close to what they should be. This random 
process is only an approximation. I put an intentional limit into 
probabilities so that events that rarer than 1E-9 do not happen unless 
the Poisson_Process constructor is suplied with a value.

-- 
Heikki Orsila			Barbie's law:
heikki.orsila@iki.fi		"Math is hard, let's go shopping!"
http://www.iki.fi/shd

Re: request for Poisson CDF inverse

[Brian Gough]

At Thu, 5 Jul 2007 17:16:38 -0700,
Jason Detwiler wrote:
> I am looking to calculate the inverse of the CDF for the Poisson
> distribution. I noticed that many of the CDFs coded in gsl have
> inverse functions also coded, but this is not the case for the Poisson
> distribution (unfortunately for me). I would like to "request" these
> functions. I would try to contribute them myself, but I'm not
> confident I could code them to the standards of gsl. What are the
> chances of my wishes being fulfilled by another developer?

Hello,

Currently there aren't any inverses for the discrete distributions.

What application do you need it for? (commercial or research?)

-- 
Brian Gough

Re: request for Poisson CDF inverse

[Heikki Orsila]

On Sun, Jul 08, 2007 at 10:54:48PM +0100, Brian Gough wrote:
> At Thu, 5 Jul 2007 17:16:38 -0700,
> Jason Detwiler wrote:
> > I am looking to calculate the inverse of the CDF for the Poisson
> > distribution. I noticed that many of the CDFs coded in gsl have
> > inverse functions also coded, but this is not the case for the Poisson
> > distribution (unfortunately for me). I would like to "request" these
> > functions. I would try to contribute them myself, but I'm not
> > confident I could code them to the standards of gsl. What are the
> > chances of my wishes being fulfilled by another developer?
> 
> Hello,
> 
> Currently there aren't any inverses for the discrete distributions.

Woud you consider a patch? I might convert the python code I posted 
into C, and make that the first discrete distribution having an inverse.

-- 
Heikki Orsila			Barbie's law:
heikki.orsila@iki.fi		"Math is hard, let's go shopping!"
http://www.iki.fi/shd

Re: request for Poisson CDF inverse

[Brian Gough]

At Mon, 9 Jul 2007 02:39:30 +0300,
Heikki Orsila wrote:
> Woud you consider a patch? I might convert the python code I posted 
> into C, and make that the first discrete distribution having an inverse.

The basic code is already there, in the inverse chi-squared
distribution, since the poisson distribution is related to that.  It's
mainly a question of testing it in cdf/test.c.

-- 
Brian Gough

Re: request for Poisson CDF inverse

[Jason Detwiler]

It's been a long time since this was discussed, but I came upon this
problem again in my research, and this time I was able to find a
reasonably fast implementation for this function. To be clear, what I
am after is the inverse of gsl_cdf_poisson_P(k, mu) with respect to
mu, i.e. the value of mu such that gsl_cdf_poisson_P(k, mu) = P for
some k, P.  Note that although the Poisson distribution is discrete,
this inverse of its CDF with respect to mu is a continuous function.

The best way I found to solve this is to root-solve f(mu) =
gsl_cdf_poisson_P(k,mu)-P using Newton's method (out to second order,
since the Poisson distribution has nice derivatives that are easy to
calculate). Here is some pseudocode for anyone interested in
implementing it in gsl (I'm still not confident I could do it up to
gsl standards):

// first guess based on variance = mu for Poisson dist:
erf((mu-k)/sqrt(2*mu))/2+0.5 ~= P
double nSigma = erfInverse(2.*P-1.)*sqrt(2.);
// nSigma*sqrt(mu) = k - mu: square both sides and solve quadratic
function for mu
// mu = [2k + nSigma^2 +/- sqrt([2k + nSigma^2]^2 - 4k^2)]/2: choose + for
double b = k + nSigma*nSigma/2.;
double mu = b + sqrt(b*b - k*k);

// Now root-solve f(mu) using 2nd order Newton's method (converges cubicly):
// mu_next = mu - f(mu)/f'(mu) + f''(mu)*f(mu)^2 / f'(mu)^3 / 2
double delta = mu;
double eps = 1.e-8;
while(fabs(delta/mu) > eps) {
  double f1 = poisson(k, mu);
  double f = (poissonCDF(k, mu) - P)/f1;
  double f2 = 1.-k/mu;
  delta = f*(1. + 0.5*f2*f);
  mu += delta;
}
return mu;

Thanks,
Jason

On Tue, Jul 10, 2007 at 7:49 AM, Brian Gough <bjg@network-theory.co.uk> wrote:
> At Mon, 9 Jul 2007 02:39:30 +0300,
> Heikki Orsila wrote:
>> Woud you consider a patch? I might convert the python code I posted
>> into C, and make that the first discrete distribution having an inverse.
>
> The basic code is already there, in the inverse chi-squared
> distribution, since the poisson distribution is related to that.  It's
> mainly a question of testing it in cdf/test.c.
>
> --
> Brian Gough
>