Re: request for Poisson CDF inverse

[Jason Detwiler <jasondet@gmail.com>]

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
>