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