rmathd.geom source code

1 module rmathd.geom;
2 
3 public import rmathd.common;
4 public import rmathd.poisson;
5 public import rmathd.exponential;
6 
7 
8 /*
9 ** dmd geom.d common.d normal.d poisson.d exponential.d && ./geom
10 */
11 
12 
13 T dgeom(T: double)(T x, T p, int give_log)
14 { 
15     T prob;
16     mixin R_D__0!give_log;
17     
18     if (isNaN(x) || isNaN(p))
19     	return x + p;
20 
21     if (p <= 0 || p > 1)
22     	return T.nan;
23 
24     mixin (R_D_nonint_check!(x));
25     if (x < 0 || !isFinite(x) || p == 0)
26     	return R_D__0;
27     x = nearbyint(x);
28 
29     /* prob = (1-p)^x, stable for small p */
30     prob = dbinom_raw!T(0., x, p, 1 - p, give_log);
31 
32     return((give_log) ? log(p) + prob : p*prob);
33 }
34 
35 
36 T pgeom(T: double)(T x, T p, int lower_tail, int log_p)
37 {
38     if (isNaN(x) || isNaN(p))
39 	    return x + p;
40     
41     if(p <= 0 || p > 1)
42     	return T.nan;
43 
44     if (x < 0.)
45     	return R_DT_0!T(lower_tail, log_p);
46     if (!isFinite(x))
47     	return R_DT_1!T(lower_tail, log_p);
48     x = floor(x + 1e-7);
49 
50     if(p == 1.) { /* we cannot assume IEEE */
51 	    x = lower_tail ? 1: 0;
52 	    return log_p ? log(x) : x;
53     }
54     x = log1p(-p) * (x + 1);
55     if (log_p)
56 	    return R_DT_Clog!T(x, lower_tail, log_p);
57     else
58 	    return lower_tail ? -expm1(x) : exp(x);
59 }
60 
61 
62 T qgeom(T: double)(T p, T prob, int lower_tail, int log_p)
63 {
64     if (isNaN(p) || isNaN(prob))
65 	    return p + prob;
66     
67     if (prob <= 0 || prob > 1)
68     	return T.nan;
69 
70     mixin (R_Q_P01_check!(p)());
71     if (prob == 1)
72     	return 0;
73     immutable(T) INF = T.infinity;
74     mixin (R_Q_P01_boundaries!(p, 0, INF)());
75 
76 /* add a fuzz to ensure left continuity, but value must be >= 0 */
77     return fmax2!T(0, ceil(R_DT_Clog!T(p, lower_tail, log_p) / log1p(- prob) - 1 - 1e-12));
78 }
79 
80 
81 T rgeom(T: double)(T p)
82 {
83     if (!isFinite(p) || p <= 0 || p > 1)
84     	return T.nan;
85 
86     return rpois!T(exp_rand!T() * ((1 - p) / p));
87 }
88 
89 
90 
91 void test_geom()
92 {
93 	import std.stdio: writeln;
94 	writeln("dgeom: ", dgeom(1., .5, 0));
95 	writeln("pgeom: ", pgeom(1., .5, 1, 0));
96 	writeln("qgeom: ", qgeom(0.8, 0.2, 1, 0));
97 	writeln("rgeom: ", rgeom(0.2), ", rgeom: ", rgeom(0.2), ", rgeom: ", rgeom(0.2), ", rgeom: ", rgeom(0.2));
98 }
99 
100 
101