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);
