*request for Poisson CDF inverse@ 2007-07-06 7:30 Jason Detwiler2007-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 inverse2007-07-06 7:30 request for Poisson CDF inverse Jason Detwiler@ 2007-07-06 9:31 ` Heikki Orsila2007-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 inverse2007-07-06 7:30 request for Poisson CDF inverse Jason Detwiler 2007-07-06 9:31 ` Heikki Orsila@ 2007-07-08 21:54 ` Brian Gough2007-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 inverse2007-07-08 21:54 ` Brian Gough@ 2007-07-08 23:39 ` Heikki Orsila2007-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 inverse2007-07-08 23:39 ` Heikki Orsila@ 2007-07-10 14:49 ` Brian Gough2011-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

*2007-07-10 14:49 ` Brian GoughRe: request for Poisson CDF inverse@ 2011-08-31 10:51 ` Jason Detwiler0 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).