[Rd] qpois minor bug (PR#8058)

[mikael@amazon.com]

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