The R Project SVN R

Rev

Rev 80266 | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 80266 Rev 80271
Line 24... Line 24...
24 
 *	double qnbinom(double p, double size, double prob,
 24 
 *	double qnbinom(double p, double size, double prob,
25 
 *                     int lower_tail, int log_p)
 25 
 *                     int lower_tail, int log_p)
26 
 *
 26 
 *
27 
 *  DESCRIPTION
 27 
 *  DESCRIPTION
28 
 *
 28 
 *
29 
 *	The quantile function of the negative binomial distribution.
 29 
 *	The quantile function of the negative binomial distribution,
- 
 
 30 
 *      for the (size, prob) parametrizations
30 
 *
 31 
 *
31 
 *  NOTES
 32 
 *  NOTES
32 
 *
 33 
 *
33 
 *	x = the number of failures before the n-th success
 34 
 *	x = the number of failures before the n-th success
34 
 *
 35 
 *
Line 42... Line 43...
42 
 */
 43 
 */
43 
 
 44 
 
44 
#include "nmath.h"
 45 
#include "nmath.h"
45 
#include "dpq.h"
 46 
#include "dpq.h"
46 
 
 47 
 
47 
/**
 - 
 
48 
   For P(x) := pnbinom(x, n, pr),  find p-quantile  y(p)   :<==>   P(y) < p <= P(y)
 - 
 
49 
 
 - 
 
50 
   @param y    current guess
 - 
 
51 
   @param *z   == pnbinom(y, ..)
 - 
 
52 
   @param p    target probability
 - 
 
53 
   @param n, pr  parameters of the negative binomial
 - 
 
54 
   @param incr increment, integer valued >= 1.
 - 
 
55 
 
 - 
 
56 
   @return  root 'y'  and  z[0] = pnbinom(y, ..)
 48 
/*-------- DEBUGGING ------------- 	make CFLAGS='-DDEBUG_qnbinom  ...'
57 
 */
 49 
 */
58 
static double
 50 
#ifdef DEBUG_qnbinom
59 
do_search(double y, double *z, double p, double n, double pr, double incr)
 - 
 
60 
{
 - 
 
61 
    if(*z >= p) {	/* search to the left */
 51 
# define R_DBG_printf(...) REprintf(__VA_ARGS__)
62 
	for(;;) {
 52 
#else
63 
	    if(y == 0 ||
 53 
# define R_DBG_printf(...)
64 
	       (*z = pnbinom(y - incr, n, pr, /*l._t.*/TRUE, /*log_p*/FALSE)) < p)
 - 
 
65 
		return y;
 54 
#endif
66 
	    y = fmax2(0, y - incr);
 - 
 
67 
	}
 55 
 
68 
    }
 - 
 
69 
    else {		/* search to the right */
 - 
 
70 
	for(;;) {
 - 
 
71 
	    y = y + incr;
 56 
#define _thisDIST_ nbinom
72 
	    if((*z = pnbinom(y, n, pr, /*l._t.*/TRUE, /*log_p*/FALSE)) >= p)
 57 
#define _dist_PARS_DECL_ double size, double prob
73 
		return y;
 58 
#define _dist_PARS_      size, prob
74 
	}
 - 
 
75 
    }
 - 
 
76 
}
 - 
 
77 
 
 59 
 
- 
 
 60 
#include "qDiscrete_search.h"
- 
 
 61 
//        ------------------>  do_search() and all called by q_DISCRETE_*() below
78 
 
 62 
 
79 
double qnbinom(double p, double size, double prob, int lower_tail, int log_p)
 63 
double qnbinom(double p, double size, double prob, int lower_tail, int log_p)
80 
{
 64 
{
81 
    double P, Q, mu, sigma, gamma, z, y;
 - 
 
82 
 
 - 
 
83 
#ifdef IEEE_754
 65 
#ifdef IEEE_754
84 
    if (ISNAN(p) || ISNAN(size) || ISNAN(prob))
 66 
    if (ISNAN(p) || ISNAN(size) || ISNAN(prob))
85 
	return p + size + prob;
 67 
	return p + size + prob;
86 
#endif
 68 
#endif
87 
 
 69 
 
88 
    /* this happens if specified via mu, size, since
 70 
    /* this happens if specified via mu, size, since
89 
       prob == size/(size+mu)
 71 
       prob == size/(size+mu)
90 
    */
 72 
    */
91 
    if (prob == 0 && size == 0) return 0;
 73 
    if (prob == 0 && size == 0) return 0;
92 
 
 - 
 
93 
    if (prob <= 0 || prob > 1 || size < 0) ML_WARN_return_NAN;
 74 
    if (prob <= 0 || prob > 1 || size < 0) ML_WARN_return_NAN;
94 
 
 - 
 
95 
    if (prob == 1 || size == 0) return 0;
 75 
    if (prob == 1 || size == 0) return 0;
96 
 
 76 
 
97 
    R_Q_P01_boundaries(p, 0, ML_POSINF);
 77 
    R_Q_P01_boundaries(p, 0, ML_POSINF);
98 
 
 78 
 
- 
 
 79 
    double
99 
    Q = 1.0 / prob;
 80 
	Q = 1.0 / prob,
100 
    P = (1.0 - prob) * Q;
 81 
	P = (1.0 - prob) * Q, // = (1 - prob) / prob  =  Q - 1
101 
    mu = size * P;
 82 
	mu = size * P,
102 
    sigma = sqrt(size * P * Q);
 83 
	sigma = sqrt(size * P * Q),
103 
    gamma = (Q + P)/sigma;
 84 
	gamma = (Q + P)/sigma;
104 
 
 85 
 
105 
    /* Note : "same" code in qpois.c, qbinom.c, qnbinom.c --
 - 
 
106 
     * FIXME: This is far from optimal [cancellation for p ~= 1, etc]: */
 86 
    R_DBG_printf("qnbinom(p=%.12g, size=%.15g, prob=%g, l.t.=%d, log=%d):"
107 
    if(!lower_tail || log_p) {
 87 
		 " mu=%g, sigma=%g, gamma=%g;\n",
108 
	p = R_DT_qIv(p); /* need check again (cancellation!): */
 - 
 
109 
	if (p == R_DT_0) return 0;
 - 
 
110 
	if (p == R_DT_1) return ML_POSINF;
 - 
 
111 
    }
 - 
 
112 
    /* temporary hack --- FIXME --- */
 - 
 
113 
    if (p + 1.01*DBL_EPSILON >= 1.) return ML_POSINF;
 - 
 
114 
 
 - 
 
115 
    /* y := approx.value (Cornish-Fisher expansion) :  */
 - 
 
116 
    z = qnorm(p, 0., 1., /*lower_tail*/TRUE, /*log_p*/FALSE);
 - 
 
117 
    y = R_forceint(mu + sigma * (z + gamma * (z*z - 1) / 6));
 - 
 
118 
 
 - 
 
119 
    z = pnbinom(y, size, prob, /*lower_tail*/TRUE, /*log_p*/FALSE);
 88 
		 p, size, prob, lower_tail, log_p, mu, sigma, gamma);
120 
 
 - 
 
121 
    /* fuzz to ensure left continuity: */
 - 
 
122 
    p *= 1 - 64*DBL_EPSILON;
 - 
 
123 
 
 - 
 
124 
    /* If the C-F value is not too large a simple search is OK */
 - 
 
125 
    if(y < 1e5) return do_search(y, &z, p, size, prob, 1);
 - 
 
126 
    /* Otherwise be a bit cleverer in the search */
 - 
 
127 
    {
 - 
 
128 
	double incr = floor(y * 0.001), oldincr;
 - 
 
129 
	do {
 - 
 
130 
	    oldincr = incr;
 - 
 
131 
	    y = do_search(y, &z, p, size, prob, incr);
 - 
 
132 
	    incr = fmax2(1, floor(incr/100));
 - 
 
133 
	} while(oldincr > 1 && incr > y*1e-15);
 - 
 
134 
	return y;
 - 
 
135 
    }
 - 
 
136 
}
 - 
 
137 
 
 89 
 
138 
double qnbinom_mu(double p, double size, double mu, int lower_tail, int log_p)
 - 
 
139 
{
 - 
 
140 
    if (size == ML_POSINF) // limit case: Poisson
 - 
 
141 
	return(qpois(p, mu, lower_tail, log_p));
 - 
 
142 
    /* FIXME!  Implement properly!! (not losing accuracy for very large size (prob ~= 1)
 - 
 
143 
       ------  --> rather invert pnbinom_mu()  */
 - 
 
144 
    double pr = size/(size+mu);
 90 
    q_DISCRETE_01_CHECKS();
145 
    if(pr == 1. && mu > 0.)
 91 
    q_DISCRETE_BODY();
146 
	MATHLIB_WARNING3(_("qnbinom_mu(%g, size=%g, mu=%g) --> prob=1: precision lost in result\n"),
 - 
 
147 
			 p, size, mu);
 - 
 
148 
    return qnbinom(p, size, /* prob = */ pr, lower_tail, log_p);
 - 
 
149 
}
 92 
}