From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 17658 invoked by alias); 31 Aug 2011 10:51:28 -0000 Received: (qmail 2319 invoked by uid 22791); 31 Aug 2011 00:54:48 -0000 X-SWARE-Spam-Status: No, hits=-0.8 required=5.0 tests=BAYES_20,DKIM_SIGNED,DKIM_VALID,DKIM_VALID_AU,FREEMAIL_FROM,RCVD_IN_DNSWL_LOW X-Spam-Check-By: sourceware.org MIME-Version: 1.0 In-Reply-To: <87hcoc9uvj.wl%bjg@network-theory.co.uk> References: <87zm26a7dj.wl%bjg@network-theory.co.uk> <20070708233930.GL24297@jolt.modeemi.cs.tut.fi> <87hcoc9uvj.wl%bjg@network-theory.co.uk> Date: Wed, 31 Aug 2011 10:51:00 -0000 Message-ID: Subject: Re: request for Poisson CDF inverse From: Jason Detwiler To: Brian Gough Cc: Heikki Orsila , gsl-discuss@sources.redhat.com Content-Type: text/plain; charset=ISO-8859-1 Content-Transfer-Encoding: quoted-printable Mailing-List: contact gsl-discuss-help@sourceware.org; run by ezmlm Precedence: bulk List-Id: List-Subscribe: List-Archive: List-Post: List-Help: , Sender: gsl-discuss-owner@sourceware.org X-SW-Source: 2011-q3/txt/msg00000.txt.bz2 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) =3D 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) =3D 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 =3D mu for Poisson dist: erf((mu-k)/sqrt(2*mu))/2+0.5 ~=3D P double nSigma =3D erfInverse(2.*P-1.)*sqrt(2.); // nSigma*sqrt(mu) =3D k - mu: square both sides and solve quadratic function for mu // mu =3D [2k + nSigma^2 +/- sqrt([2k + nSigma^2]^2 - 4k^2)]/2: choose + for double b =3D k + nSigma*nSigma/2.; double mu =3D b + sqrt(b*b - k*k); // Now root-solve f(mu) using 2nd order Newton's method (converges cubicly): // mu_next =3D mu - f(mu)/f'(mu) + f''(mu)*f(mu)^2 / f'(mu)^3 / 2 double delta =3D mu; double eps =3D 1.e-8; while(fabs(delta/mu) > eps) { double f1 =3D poisson(k, mu); double f =3D (poissonCDF(k, mu) - P)/f1; double f2 =3D 1.-k/mu; delta =3D f*(1. + 0.5*f2*f); mu +=3D delta; } return mu; Thanks, Jason On Tue, Jul 10, 2007 at 7:49 AM, Brian Gough wro= te: > 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. =A0It's > mainly a question of testing it in cdf/test.c. > > -- > Brian Gough >