public inbox for gsl-discuss@sourceware.org
 help / color / mirror / Atom feed
* request for Poisson CDF inverse
@ 2007-07-06  7:30 Jason Detwiler
  2007-07-06  9:31 ` Heikki Orsila
  2007-07-08 21:54 ` Brian Gough
  0 siblings, 2 replies; 6+ messages in thread
From: Jason Detwiler @ 2007-07-06  7:30 UTC (permalink / raw)
  To: gsl-discuss

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

^ permalink raw reply	[flat|nested] 6+ messages in thread

* Re: request for Poisson CDF inverse
  2007-07-06  7:30 request for Poisson CDF inverse Jason Detwiler
@ 2007-07-06  9:31 ` Heikki Orsila
  2007-07-08 21:54 ` Brian Gough
  1 sibling, 0 replies; 6+ messages in thread
From: Heikki Orsila @ 2007-07-06  9:31 UTC (permalink / raw)
  To: Jason Detwiler; +Cc: gsl-discuss

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

^ permalink raw reply	[flat|nested] 6+ messages in thread

* Re: request for Poisson CDF inverse
  2007-07-06  7:30 request for Poisson CDF inverse Jason Detwiler
  2007-07-06  9:31 ` Heikki Orsila
@ 2007-07-08 21:54 ` Brian Gough
  2007-07-08 23:39   ` Heikki Orsila
  1 sibling, 1 reply; 6+ messages in thread
From: Brian Gough @ 2007-07-08 21:54 UTC (permalink / raw)
  To: Jason Detwiler; +Cc: gsl-discuss

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

^ permalink raw reply	[flat|nested] 6+ messages in thread

* Re: request for Poisson CDF inverse
  2007-07-08 21:54 ` Brian Gough
@ 2007-07-08 23:39   ` Heikki Orsila
  2007-07-10 14:49     ` Brian Gough
  0 siblings, 1 reply; 6+ messages in thread
From: Heikki Orsila @ 2007-07-08 23:39 UTC (permalink / raw)
  To: Brian Gough; +Cc: Jason Detwiler, gsl-discuss

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

^ permalink raw reply	[flat|nested] 6+ messages in thread

* Re: request for Poisson CDF inverse
  2007-07-08 23:39   ` Heikki Orsila
@ 2007-07-10 14:49     ` Brian Gough
  2011-08-31 10:51       ` Jason Detwiler
  0 siblings, 1 reply; 6+ messages in thread
From: Brian Gough @ 2007-07-10 14:49 UTC (permalink / raw)
  To: Heikki Orsila; +Cc: Jason Detwiler, gsl-discuss

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

^ permalink raw reply	[flat|nested] 6+ messages in thread

* Re: request for Poisson CDF inverse
  2007-07-10 14:49     ` Brian Gough
@ 2011-08-31 10:51       ` Jason Detwiler
  0 siblings, 0 replies; 6+ messages in thread
From: Jason Detwiler @ 2011-08-31 10:51 UTC (permalink / raw)
  To: Brian Gough; +Cc: Heikki Orsila, gsl-discuss

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
>

^ permalink raw reply	[flat|nested] 6+ messages in thread

end of thread, other threads:[~2011-08-31 10:51 UTC | newest]

Thread overview: 6+ messages (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2007-07-06  7:30 request for Poisson CDF inverse Jason Detwiler
2007-07-06  9:31 ` Heikki Orsila
2007-07-08 21:54 ` Brian Gough
2007-07-08 23:39   ` Heikki Orsila
2007-07-10 14:49     ` Brian Gough
2011-08-31 10:51       ` Jason Detwiler

This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).