rmathd.nbinomial source code

1 module rmathd.nbinomial;
2 
3 public import rmathd.common;
4 public import rmathd.beta;
5 public import rmathd.poisson;
6 public import rmathd.toms708;
7 public import rmathd.normal;
8 public import rmathd.gamma;
9 
10 /* 
11 ** 
12 ** dmd nbinomial.d common.d beta.d poisson.d toms708 normal.d gamma.d && ./nbinomial
13 */
14 
15 
16 T dnbinom(T: double)(T x, T size, T prob, int give_log)
17 {
18     T ans, p;
19     mixin R_D__0!give_log;
20     mixin R_D__1!give_log;
21     
22     if (isNaN(x) || isNaN(size) || isNaN(prob))
23         return x + size + prob;
24 
25     if (prob <= 0 || prob > 1 || size < 0)
26     	return T.nan;
27     mixin (R_D_nonint_check!(x));
28     if (x < 0 || !isFinite(x))
29     	return R_D__0;
30     /* limiting case as size approaches zero is point mass at zero */
31     if (x == 0 && size == 0)
32     	return R_D__1;
33     x = nearbyint(x);
34     if(!isFinite(size)) size = DBL_MAX;
35 
36     ans = dbinom_raw!T(size, x + size, prob, 1 - prob, give_log);
37     p = (cast(T)size)/(size + x);
38     return((give_log) ? log(p) + ans : p * ans);
39 }
40 
41 
42 T dnbinom_mu(T: double)(T x, T size, T mu, int give_log)
43 {
44     /* originally, just set  prob :=  size / (size + mu)  and called dbinom_raw(),
45      * but that suffers from cancellation when   mu << size  */
46 
47     mixin R_D__0!give_log;
48     mixin R_D__1!give_log;
49 
50     if (isNaN(x) || isNaN(size) || isNaN(mu))
51         return x + size + mu;
52 
53     if (mu < 0 || size < 0)
54         return T.nan;
55     mixin (R_D_nonint_check!(x));
56     if (x < 0 || !isFinite(x))
57         return R_D__0;
58 
59     /* limiting case as size approaches zero is point mass at zero,
60      * even if mu is kept constant. limit distribution does not
61      * have mean mu, though.
62      */
63     if (x == 0 && size == 0)
64         return R_D__1;
65     x = nearbyint(x);
66     if(!isFinite(size)) // limit case: Poisson
67         return(dpois_raw!T(x, mu, give_log));
68 
69     if(x == 0)/* be accurate, both for n << mu, and n >> mu :*/
70         return R_D_exp!T(size * (size < mu ? log(size/(size + mu)) : log1p(-mu/(size + mu))), give_log);
71     if(x < 1e-10 * size) { /* don't use dbinom_raw() but MM's formula: */
72         /* FIXME --- 1e-8 shows problem; rather use algdiv() from ./toms708.c */
73         T p = (size < mu ? log(size/(1 + size/mu)) : log(mu / (1 + mu/size)));
74         return R_D_exp!T(x * p - mu - lgamma!T(x + 1) +
75                    log1p(x*(x - 1)/(2*size)), give_log);
76     } else {
77         /* no unnecessary cancellation inside dbinom_raw, when
78          * x_ = size and n_ = x+size are so close that n_ - x_ loses accuracy */
79         T p = (cast(T)size)/(size + x),
80             ans = dbinom_raw!T(size, x + size, size/(size + mu), mu/(size + mu), give_log);
81         return((give_log) ? log(p) + ans : p * ans);
82     }
83 }
84 
85 
86 T pnbinom(T: double)(T x, T size, T prob, int lower_tail, int log_p)
87 {
88     if (isNaN(x) || isNaN(size) || isNaN(prob))
89         return x + size + prob;
90     if(!isFinite(size) || !isFinite(prob))
91         return T.nan;
92 
93     if (size < 0 || prob <= 0 || prob > 1)
94         return T.nan;
95 
96     /* limiting case: point mass at zero */
97     if (size == 0)
98         return (x >= 0) ? R_DT_1!T(lower_tail, log_p): R_DT_0!T(lower_tail, log_p);
99 
100     if (x < 0)
101         return R_DT_0!T(lower_tail, log_p);
102     if (!isFinite(x))
103         return R_DT_1!T(lower_tail, log_p);
104     x = floor(x + 1e-7);
105     return pbeta!T(prob, size, x + 1, lower_tail, log_p);
106 }
107 
108 
109 T pnbinom_mu(T: double)(T x, T size, T mu, int lower_tail, int log_p)
110 {
111     if (isNaN(x) || isNaN(size) || isNaN(mu))
112         return x + size + mu;
113     if(!isFinite(mu))
114         return T.nan;
115     
116     if (size < 0 || mu < 0)
117         return T.nan;
118 
119     /* limiting case: point mass at zero */
120     if (size == 0)
121         return (x >= 0) ? R_DT_1!T(lower_tail, log_p) : R_DT_0!T(lower_tail, log_p);
122 
123     if (x < 0)
124         return R_DT_0!T(lower_tail, log_p);
125     if (!isFinite(x))
126         return R_DT_1!T(lower_tail, log_p);
127     if (!isFinite(size)) // limit case: Poisson
128         return(ppois!T(x, mu, lower_tail, log_p));
129 
130     x = floor(x + 1e-7);
131     /* return
132      * pbeta(pr, size, x + 1, lower_tail, log_p);  pr = size/(size + mu), 1-pr = mu/(size+mu)
133      *
134      *= pbeta_raw(pr, size, x + 1, lower_tail, log_p)
135      *            x.  pin   qin
136      *=  bratio (pin,  qin, x., 1-x., &w, &wc, &ierr, log_p),  and return w or wc ..
137      *=  bratio (size, x+1, pr, 1-pr, &w, &wc, &ierr, log_p) */
138     {
139     int ierr;
140     T w, wc;
141     bratio!T(size, x + 1, size/(size + mu), mu/(size + mu), &w, &wc, &ierr, log_p);
142     if(ierr){
143         //MATHLIB_WARNING(_("pnbinom_mu() -> bratio() gave error code %d"), ierr);
144         doNothing();
145     }
146     return lower_tail ? w : wc;
147     }
148 }
149 
150 
151 static T do_search(T)(T y, T *z, T p, T n, T pr, T incr)
152 {
153     if(*z >= p) {   /* search to the left */
154         for(;;) {
155             if(y == 0 || (*z = pnbinom!T(y - incr, n, pr, /*l._t.*/1, /*log_p*/0)) < p)
156                 return y;
157             y = fmax2!T(0, y - incr);
158         }
159     } else {      /* search to the right */
160         for(;;) {
161             y = y + incr;
162             if((*z = pnbinom!T(y, n, pr, /*l._t.*/1, /*log_p*/0)) >= p)
163                 return y;
164         }
165     }
166 }
167 
168 
169 T qnbinom(T: double)(T p, T size, T prob, int lower_tail, int log_p)
170 {
171     T P, Q, mu, sigma, gamma, z, y;
172 
173     if (isNaN(p) || isNaN(size) || isNaN(prob))
174     return p + size + prob;
175 
176     /* this happens if specified via mu, size, since
177        prob == size/(size+mu)
178     */
179     if (prob == 0 && size == 0)
180         return 0;
181 
182     if (prob <= 0 || prob > 1 || size < 0)
183         return T.nan;
184 
185     if (prob == 1 || size == 0)
186         return 0;
187 
188     immutable(T) INF = T.infinity;
189     mixin (R_Q_P01_boundaries!(p, 0, INF));
190 
191     Q = 1.0 / prob;
192     P = (1.0 - prob) * Q;
193     mu = size * P;
194     sigma = sqrt(size * P * Q);
195     gamma = (Q + P)/sigma;
196 
197     mixin R_DT_qIv!p;
198     /* Note : "same" code in qpois.c, qbinom.c, qnbinom.c --
199      * FIXME: This is far from optimal [cancellation for p ~= 1, etc]: */
200     if(!lower_tail || log_p) {
201     p = R_DT_qIv; /* need check again (cancellation!): */
202     if (p == R_DT_0!T(lower_tail, log_p))
203         return 0;
204     if (p == R_DT_1!T(lower_tail, log_p))
205         return T.infinity;
206     }
207     /* temporary hack --- FIXME --- */
208     if (p + 1.01*DBL_EPSILON >= 1.)
209         return T.infinity;
210 
211     /* y := approx.value (Cornish-Fisher expansion) :  */
212     z = qnorm!T(p, 0., 1., /*lower_tail*/1, /*log_p*/0);
213     y = nearbyint(mu + sigma * (z + gamma * (z*z - 1) / 6));
214 
215     z = pnbinom!T(y, size, prob, /*lower_tail*/1, /*log_p*/0);
216 
217     /* fuzz to ensure left continuity: */
218     p *= 1 - 64*DBL_EPSILON;
219 
220     /* If the C-F value is not too large a simple search is OK */
221     if(y < 1e5)
222         return do_search!T(y, &z, p, size, prob, 1);
223     /* Otherwise be a bit cleverer in the search */
224     {
225     T incr = floor(y * 0.001), oldincr;
226     do {
227         oldincr = incr;
228         y = do_search!T(y, &z, p, size, prob, incr);
229         incr = fmax2!T(1, floor(incr/100));
230     } while(oldincr > 1 && incr > y*1e-15);
231     return y;
232     }
233 }
234 
235 T qnbinom_mu(T: double)(T p, T size, T mu, int lower_tail, int log_p)
236 {
237     if (size == T.infinity) // limit case: Poisson
238         return(qpois!T(p, mu, lower_tail, log_p));
239     /* FIXME!  Implement properly!! (not losing accuracy for very large size (prob ~= 1)*/
240     return qnbinom!T(p, size, /* prob = */ size/(size + mu), lower_tail, log_p);
241 }
242 
243 
244 T rnbinom(T: double)(T size, T prob)
245 {
246     if(!isFinite(prob) || isNaN(size) || size <= 0 || prob <= 0 || prob > 1)
247     /* prob = 1 is ok, PR#1218 */
248         return T.nan;
249     if(!isFinite(size)) size = DBL_MAX / 2.; // '/2' to prevent rgamma() returning Inf
250     return (prob == 1) ? 0 : rpois(rgamma!T(size, (1 - prob) / prob));
251 }
252 
253 T rnbinom_mu(T: double)(T size, T mu)
254 {
255     if(!isFinite(mu) || isNaN(size) || size <= 0 || mu < 0)
256         return T.nan;
257     if(!isFinite(size)) size = DBL_MAX / 2.;
258     return (mu == 0) ? 0 : rpois!T(rgamma!T(size, mu / size));
259 }
260 
261 
262 void test_nbinomial()
263 {
264     import std.stdio: writeln;
265 	writeln("dnbinom: ", dnbinom(2., 8., 0.7, 0));
266     writeln("dnbinom_mu: ", dnbinom_mu(2., 8., 0.7, 0));
267     writeln("pnbinom: ", pnbinom(2., 8., 0.7, 1, 0));
268     writeln("pnbinom_mu: ", pnbinom_mu(2., 8., 0.7, 1, 0));
269     writeln("qnbinom: ", qnbinom(.8, 8., 0.7, 1, 0));
270     writeln("qnbinom_mu: ", qnbinom_mu(.8, 8., 0.7, 1, 0));
271     writeln("rnbinom: ", rnbinom(8., 0.2), ", rnbinom: ", rnbinom(8., 0.2), ", rnbinom: ", rnbinom(8., 0.2));
272     writeln("rnbinom_mu: ", rnbinom_mu(8., 0.7), ", rnbinom_mu: ", rnbinom_mu(8., 0.7), ", rnbinom_mu: ", rnbinom_mu(8., 0.7));
273 }
274