The R Project SVN R

Rev

Rev 8872 | Rev 35246 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 8872 Rev 10470
Line 63... Line 63...
63 
    }
 63 
    }
64 
    /* temporary hack --- FIXME --- */
 64 
    /* temporary hack --- FIXME --- */
65 
    if (p + 1.01*DBL_EPSILON >= 1.) return ML_POSINF;
 65 
    if (p + 1.01*DBL_EPSILON >= 1.) return ML_POSINF;
66 
 
 66 
 
67 
    /* y := approx.value (Cornish-Fisher expansion) :  */
 67 
    /* y := approx.value (Cornish-Fisher expansion) :  */
68 
    z = qnorm(p, 0., 1., /*lower_tail*/LTRUE, /*log_p*/LFALSE);
 68 
    z = qnorm(p, 0., 1., /*lower_tail*/TRUE, /*log_p*/FALSE);
69 
    y = floor(mu + sigma * (z + gamma * (z*z - 1) / 6) + 0.5);
 69 
    y = floor(mu + sigma * (z + gamma * (z*z - 1) / 6) + 0.5);
70 
 
 70 
 
71 
    z = ppois(y, lambda, /*lower_tail*/LTRUE, /*log_p*/LFALSE);
 71 
    z = ppois(y, lambda, /*lower_tail*/TRUE, /*log_p*/FALSE);
72 
 
 72 
 
73 
    /* fuzz to ensure left continuity; 1 - 1e-7 may lose too much : */
 73 
    /* fuzz to ensure left continuity; 1 - 1e-7 may lose too much : */
74 
    p *= 1 - 64*DBL_EPSILON;
 74 
    p *= 1 - 64*DBL_EPSILON;
75 
 
 75 
 
76 
/*-- Fixme, here y can be way off --
 76 
/*-- Fixme, here y can be way off --
Line 82... Line 82...
82 
    if(z >= p) {
 82 
    if(z >= p) {
83 
#endif
 83 
#endif
84 
			/* search to the left */
 84 
			/* search to the left */
85 
	for(;;) {
 85 
	for(;;) {
86 
	    if(y == 0 ||
 86 
	    if(y == 0 ||
87 
	       (z = ppois(y - 1, lambda, /*l._t.*/LTRUE, /*log_p*/LFALSE)) < p)
 87 
	       (z = ppois(y - 1, lambda, /*l._t.*/TRUE, /*log_p*/FALSE)) < p)
88 
		return y;
 88 
		return y;
89 
	    y = y - 1;
 89 
	    y = y - 1;
90 
	}
 90 
	}
91 
    }
 91 
    }
92 
    else {		/* search to the right */
 92 
    else {		/* search to the right */
93 
	for(;;) {
 93 
	for(;;) {
94 
	    y = y + 1;
 94 
	    y = y + 1;
95 
	    if((z = ppois(y, lambda, /*l._t.*/LTRUE, /*log_p*/LFALSE)) >= p)
 95 
	    if((z = ppois(y, lambda, /*l._t.*/TRUE, /*log_p*/FALSE)) >= p)
96 
		return y;
 96 
		return y;
97 
	}
 97 
	}
98 
    }
 98 
    }
99 
}
 99 
}