# [Rd] qpois minor bug (PR#8058)

From: <mikael_at_amazon.com>
Date: Tue 09 Aug 2005 - 17:59:01 GMT

Full_Name: Mikael Weigelt
Version: 2.0
OS: windows
Submission from: (NULL) (207.171.180.101)

The calculation of the qpois attempts to use the Cornish-Fisher expansion as a starting approximation. The definition of the expansion is incorrect. However, since this approximation just gives an initial solution, the end result of the function is still correct.

To fix the approximation, in the snippet below the line

gamma = sigma;
should be replaced by

gamma = 1.0/sigma; /* the skewness */

The reference is Abramowitz and Stegun 'Handbook of Mathmatical Functions' pages 935 and 928

Mikael

double qpois(double p, double lambda, int lower_tail, int log_p) {

double mu, sigma, gamma, z, y;
#ifdef IEEE_754

if (ISNAN(p) || ISNAN(lambda))

return p + lambda;
#endif

if(!R_FINITE(lambda))

ML_ERR_return_NAN;

R_Q_P01_boundaries(p, 0, ML_POSINF);

if(lambda < 0) ML_ERR_return_NAN;
if(lambda == 0) return 0;

mu = lambda;
sigma = sqrt(lambda);
gamma = sigma;

/* Note : "same" code in qpois.c, qbinom.c, qnbinom.c --

• FIXME: This is far from optimal [cancellation for p ~= 1, etc]: */ if(!lower_tail || log_p) { p = R_DT_qIv(p); /* need check again (cancellation!): */ if (p == 0.) return 0; if (p == 1.) return ML_POSINF; }
/* temporary hack --- FIXME --- */
if (p + 1.01*DBL_EPSILON >= 1.) return ML_POSINF;

/* y := approx.value (Cornish-Fisher expansion) : */
z = qnorm(p, 0., 1., /*lower_tail*/TRUE, /*log_p*/FALSE);     y = floor(mu + sigma * (z + gamma * (z*z - 1) / 6) + 0.5);

R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel Received on Wed Aug 10 04:05:21 2005

This archive was generated by hypermail 2.1.8 : Mon 24 Oct 2005 - 22:27:37 GMT