rmathd.binomial source code

1 module rmathd.binomial;
2 
3 public import rmathd.common;
4 public import rmathd.beta;
5 //import toms708;
6 //import normal;
7 
8 /*
9 ** dmd binomial.d common.d toms708.d beta.d normal.d && ./binomial
10 */
11 
12 T dbinom(T: double)(T x, T n, T p, int give_log)
13 {
14     /* NaNs propagated correctly */
15     if (isNaN(x) || isNaN(n) || isNaN(p))
16     	return x + n + p;
17 
18     mixin R_D__0!give_log;
19 
20     if (p < 0 || p > 1 || R_D_negInonint!T(n))
21 	    return T.nan;
22     mixin R_D_nonint_check!(x);
23     if (x < 0 || !isFinite(x))
24         return R_D__0;
25 
26     n = nearbyint(n);
27     x = nearbyint(x);
28 
29     return dbinom_raw(x, n, p, 1 - p, give_log);
30 }
31 
32 T pbinom(T: double)(T x, T n, T p, int lower_tail, int log_p)
33 {
34     if (isNaN(x) || isNaN(n) || isNaN(p))
35         return x + n + p;
36     if (!isFinite(n) || !isFinite(p))
37         return T.nan;
38 
39     if(R_nonint!T(n)) {
40         return T.nan;
41     }
42     n = nearbyint(n);
43     if(n < 0 || p < 0 || p > 1)
44         return T.nan;
45 
46     if (x < 0)
47         return R_DT_0!T(lower_tail, log_p);
48     x = floor(x + 1e-7);
49     if (n <= x)
50         return R_DT_1!T(lower_tail, log_p);
51     return pbeta!T(p, x + 1, n - x, !lower_tail, log_p);
52 }
53 
54 /* toms708.c */
55 static T do_search(T)(T y, T *z, T p, T n, T pr, T incr)
56 {
57     if(*z >= p) {
58             /* search to the left */
59         //#ifdef DEBUG_qbinom
60         //    REprintf("\tnew z=%7g >= p = %7g  --> search to left (y--) ..\n", z,p);
61         //#endif
62         for(;;) {
63             T newz;
64             if(y == 0 || (newz = pbinom!T(y - incr, n, pr, /*l._t.*/1, /*log_p*/0)) < p)
65                 return y;
66             y = fmax2!T(0, y - incr);
67             *z = newz;
68         }
69     } else {      /* search to the right */
70         //#ifdef DEBUG_qbinom
71         //    REprintf("\tnew z=%7g < p = %7g  --> search to right (y++) ..\n", z,p);
72         //#endif
73         for(;;) {
74             y = fmin2!T(y + incr, n);
75             if(y == n || (*z = pbinom!T(y, n, pr, /*l._t.*/1, /*log_p*/0)) >= p)
76                 return y;
77         }
78     }
79 }
80 
81 
82 T qbinom(T: double)(T p, T n, T pr, int lower_tail, int log_p)
83 {
84     T q, mu, sigma, gamma, z, y;
85     
86     if (isNaN(p) || isNaN(n) || isNaN(pr))
87         return p + n + pr;
88     
89     if(!isFinite(n) || !isFinite(pr))
90         return T.nan;
91     /* if log_p is true, p = -Inf is a legitimate value */
92     if(!isFinite(p) && !log_p)
93         return T.nan;
94 
95     if(n != floor(n + 0.5))
96         return T.nan;
97     if (pr < 0 || pr > 1 || n < 0)
98         return T.nan;
99 
100     mixin (R_Q_P01_boundaries!(p, 0, n)());
101 
102     if (pr == 0. || n == 0)
103         return 0.;
104 
105     q = 1 - pr;
106     if(q == 0.) return n; /* covers the full range of the distribution */
107     mu = n * pr;
108     sigma = sqrt(n * pr * q);
109     gamma = (q - pr) / sigma;
110 
111     //#ifdef DEBUG_qbinom
112     //    REprintf("qbinom(p=%7g, n=%g, pr=%7g, l.t.=%d, log=%d): sigm=%g, gam=%g\n",
113     //         p,n,pr, lower_tail, log_p, sigma, gamma);
114     //#endif
115 
116     /* Note : "same" code in qpois.c, qbinom.c, qnbinom.c --
117      * FIXME: This is far from optimal [cancellation for p ~= 1, etc]: */
118      mixin R_DT_qIv!(p);
119     if(!lower_tail || log_p) {
120         p = R_DT_qIv; /* need check again (cancellation!): */
121         if (p == 0.)
122             return 0.;
123         if (p == 1.)
124             return n;
125     }
126     /* temporary hack --- FIXME --- */
127     if (p + 1.01*DBL_EPSILON >= 1.) return n;
128 
129     /* y := approx.value (Cornish-Fisher expansion) :  */
130     z = qnorm!T(p, 0., 1., /*lower_tail*/1, /*log_p*/0);
131     y = floor(mu + sigma * (z + gamma * (z*z - 1) / 6) + 0.5);
132 
133     if(y > n) /* way off */ y = n;
134 
135     //#ifdef DEBUG_qbinom
136     //    REprintf("  new (p,1-p)=(%7g,%7g), z=qnorm(..)=%7g, y=%5g\n", p, 1-p, z, y);
137     //#endif
138     z = pbinom!T(y, n, pr, /*lower_tail*/1, /*log_p*/0);
139 
140     /* fuzz to ensure left continuity: */
141     p *= 1 - 64*DBL_EPSILON;
142 
143     if(n < 1e5) return do_search!T(y, &z, p, n, pr, 1);
144     /* Otherwise be a bit cleverer in the search */
145     {
146     T incr = floor(n * 0.001), oldincr;
147     do {
148         oldincr = incr;
149         y = do_search!T(y, &z, p, n, pr, incr);
150         incr = fmax2!T(1, floor(incr/100));
151     } while(oldincr > 1 && incr > n*1e-15);
152     return y;
153     }
154 }
155 
156 T rbinom(T: double)(T nin, T pp)
157 {
158     /* FIXME: These should become THREAD_specific globals : */
159 
160     static T c, fm, npq, p1, p2, p3, p4, qn;
161     static T xl, xll, xlr, xm, xr;
162 
163     static T psave = -1.0;
164     static int nsave = -1;
165     static int m;
166 
167     T f, f1, f2, u, v, w, w2, x, x1, x2, z, z2;
168     T p, q, np, g, r, al, alv, amaxp, ffm, ynorm;
169     int i, ix, k, n;
170 
171     if (!isFinite(nin))
172         return T.nan;
173     r = nearbyint(nin);
174     if (r != nin)
175         return T.nan;
176     if (!isFinite(pp) ||
177     /* n=0, p=0, p=1 are not errors <TSL>*/
178     r < 0 || pp < 0. || pp > 1.)
179         return T.nan;
180 
181     if (r == 0 || pp == 0.) return 0;
182     if (pp == 1.) return r;
183 
184     if (r >= INT_MAX)/* evade integer overflow,
185             and r == INT_MAX gave only even values */
186     return qbinom!T(unif_rand!T(), r, pp, /*lower_tail*/ 0, /*log_p*/ 0);
187     /* else */
188     n = cast(int) r;
189 
190     p = fmin2(pp, 1. - pp);
191     q = 1. - p;
192     np = n * p;
193     r = p / q;
194     g = r * (n + 1);
195 
196     /* Setup, perform only when parameters change [using static (globals): */
197 
198     /* FIXING: Want this thread safe
199        -- use as little (thread globals) as possible
200     */
201     if (pp != psave || n != nsave) {
202     psave = pp;
203     nsave = n;
204     if (np < 30.0) {
205         /* inverse cdf logic for mean less than 30 */
206         qn = R_pow_di!T(q, n);
207         goto L_np_small;
208     } else {
209         ffm = np + p;
210         m = cast(int) ffm;
211         fm = m;
212         npq = np*q;
213         p1 = cast(int)(2.195 * sqrt(npq) - 4.6 * q) + 0.5;
214         xm = fm + 0.5;
215         xl = xm - p1;
216         xr = xm + p1;
217         c = 0.134 + 20.5 / (15.3 + fm);
218         al = (ffm - xl) / (ffm - xl * p);
219         xll = al * (1.0 + 0.5 * al);
220         al = (xr - ffm) / (xr * q);
221         xlr = al * (1.0 + 0.5 * al);
222         p2 = p1 * (1.0 + c + c);
223         p3 = p2 + c / xll;
224         p4 = p3 + c / xlr;
225     }
226     } else if (n == nsave) {
227     if (np < 30.0)
228         goto L_np_small;
229     }
230 
231     /*-------------------------- np = n*p >= 30 : ------------------- */
232     for(;;) {
233       u = unif_rand!T() * p4;
234       v = unif_rand!T();
235       /* triangular region */
236       if (u <= p1) {
237       ix = cast(int)(xm - p1 * v + u);
238       goto finis;
239       }
240       /* parallelogram region */
241       if (u <= p2) {
242       x = xl + (u - p1) / c;
243       v = v * c + 1.0 - fabs(xm - x) / p1;
244       if (v > 1.0 || v <= 0.)
245           continue;
246       ix = cast(int) x;
247       } else {
248       if (u > p3) { /* right tail */
249           ix = cast(int)(xr - log(v) / xlr);
250           if (ix > n)
251           continue;
252           v = v * (u - p3) * xlr;
253       } else {/* left tail */
254           ix = cast(int)(xl + log(v) / xll);
255           if (ix < 0)
256           continue;
257           v = v * (u - p2) * xll;
258       }
259       }
260       /* determine appropriate way to perform accept/reject test */
261       k = abs(ix - m);
262       if (k <= 20 || k >= npq / 2 - 1) {
263       /* explicit evaluation */
264       f = 1.0;
265       if (m < ix) {
266           for (i = m + 1; i <= ix; i++)
267           f *= (g / i - r);
268       } else if (m != ix) {
269           for (i = ix + 1; i <= m; i++)
270           f /= (g / i - r);
271       }
272       if (v <= f)
273           goto finis;
274       } else {
275       /* squeezing using upper and lower bounds on log(f(x)) */
276       amaxp = (k / npq) * ((k * (k / 3. + 0.625) + 0.1666666666666) / npq + 0.5);
277       ynorm = -k * k / (2.0 * npq);
278       alv = log(v);
279       if (alv < ynorm - amaxp)
280           goto finis;
281       if (alv <= ynorm + amaxp) {
282           /* stirling's formula to machine accuracy */
283           /* for the final acceptance/rejection test */
284           x1 = ix + 1;
285           f1 = fm + 1.0;
286           z = n + 1 - fm;
287           w = n - ix + 1.0;
288           z2 = z * z;
289           x2 = x1 * x1;
290           f2 = f1 * f1;
291           w2 = w * w;
292           if (alv <= xm * log(f1 / x1) + (n - m + 0.5) * log(z / w) + (ix - m) * log(w * p / (x1 * q)) + (13860.0 - (462.0 - (132.0 - (99.0 - 140.0 / f2) / f2) / f2) / f2) / f1 / 166320.0 + (13860.0 - (462.0 - (132.0 - (99.0 - 140.0 / z2) / z2) / z2) / z2) / z / 166320.0 + (13860.0 - (462.0 - (132.0 - (99.0 - 140.0 / x2) / x2) / x2) / x2) / x1 / 166320.0 + (13860.0 - (462.0 - (132.0 - (99.0 - 140.0 / w2) / w2) / w2) / w2) / w / 166320.)
293           goto finis;
294       }
295       }
296   }
297 
298  L_np_small:
299     /*---------------------- np = n*p < 30 : ------------------------- */
300 
301   for(;;) {
302      ix = 0;
303      f = qn;
304      u = unif_rand!T();
305      for(;;) {
306      if (u < f)
307          goto finis;
308      if (ix > 110)
309          break;
310      u -= f;
311      ix++;
312      f *= (g / ix - r);
313      }
314   }
315  finis:
316     if (psave > 0.5)
317      ix = n - ix;
318   return cast(T)ix;
319 }
320 
321 
322 
323 void test_binomial()
324 {
325     import std.stdio: writeln;
326     writeln("dbinom: ", dbinom(2., 4., .5, 0));
327     writeln("pbinom: ", pbinom(2., 4., .5, 1, 0));
328     writeln("qbinom: ", qbinom(.8, 4., .5, 1, 0));
329     writeln("rbinom: ", rbinom(20, .3), ", rbinom: ", rbinom(20, .3), ", rbinom: ", rbinom(20, .3));
330 }