IsoSpec 2.2.1
Loading...
Searching...
No Matches
btrd.h
1/* This file was taken from Boost, as permitted by Boost licence,
2 * with slight modifications. The reason is: we don't want to introduce
3 * dependency on the whole Boost just for this one thing.
4 *
5 * Source: boost random/binomial_distribution.hpp header file, at version 1.71
6 *
7 * Copyright Steven Watanabe 2010
8 * Distributed under the Boost Software License, Version 1.0. (See
9 * accompanying file LICENSE_1_0.txt or copy at
10 * http://www.boost.org/LICENSE_1_0.txt)
11 *
12 * See http://www.boost.org for most recent version including documentation.
13 *
14 */
15
16#pragma once
17
18#include "isoMath.h"
19#include <cstdlib>
20#include <cstdint>
21#include <cmath>
22#include <limits>
23
24
25namespace IsoSpec {
26
27typedef double RealType;
28typedef int64_t IntType;
29
30
31static const RealType btrd_binomial_table[10] = {
32 0.08106146679532726,
33 0.04134069595540929,
34 0.02767792568499834,
35 0.02079067210376509,
36 0.01664469118982119,
37 0.01387612882307075,
38 0.01189670994589177,
39 0.01041126526197209,
40 0.009255462182712733,
41 0.008330563433362871
42};
43
44
63// computes the correction factor for the Stirling approximation
64// for log(k!)
65static RealType fc(IntType k)
66{
67 if(k < 10) { return btrd_binomial_table[k]; }
68 else
69 {
70 RealType ikp1 = RealType(1) / (k + 1);
71 return (RealType(1)/12
72 - (RealType(1)/360
73 - (RealType(1)/1260)*(ikp1*ikp1))*(ikp1*ikp1))*ikp1;
74 }
75}
76
77IntType btrd(IntType _t, RealType p, IntType m, std::mt19937& urng = random_gen)
78{
79 using std::floor;
80 using std::abs;
81 using std::log;
82
83 RealType btrd_r = p/(1-p);
84 RealType btrd_nr = (_t+1)*btrd_r;
85 RealType btrd_npq = _t*p*(1-p);
86 RealType sqrt_npq = sqrt(btrd_npq);
87 RealType btrd_b = 1.15 + 2.53 * sqrt_npq;
88 RealType btrd_a = -0.0873 + 0.0248*btrd_b + 0.01*p;
89 RealType btrd_c = _t*p + 0.5;
90 RealType btrd_alpha = (2.83 + 5.1/btrd_b) * sqrt_npq;
91 RealType btrd_v_r = 0.92 - 4.2/btrd_b;
92 RealType btrd_u_rv_r = 0.86*btrd_v_r;
93
94 while(true) {
95 RealType u;
96 RealType v = stdunif(urng);
97 if(v <= btrd_u_rv_r) {
98 u = v/btrd_v_r - 0.43;
99 return static_cast<IntType>(floor(
100 (2*btrd_a/(0.5 - abs(u)) + btrd_b)*u + btrd_c));
101 }
102
103 if(v >= btrd_v_r) {
104 u = stdunif(urng) - 0.5;
105 } else {
106 u = v/btrd_v_r - 0.93;
107 u = ((u < 0)? -0.5 : 0.5) - u;
108 v = stdunif(urng) * btrd_v_r;
109 }
110
111 RealType us = 0.5 - abs(u);
112 IntType k = static_cast<IntType>(floor((2*btrd_a/us + btrd_b)*u + btrd_c));
113 if(k < 0 || k > _t) continue;
114 v = v*btrd_alpha/(btrd_a/(us*us) + btrd_b);
115 RealType km = abs(k - m);
116 if(km <= 15) {
117 RealType f = 1;
118 if(m < k) {
119 IntType i = m;
120 do {
121 ++i;
122 f = f*(btrd_nr/i - btrd_r);
123 } while(i != k);
124 } else if(m > k) {
125 IntType i = k;
126 do {
127 ++i;
128 v = v*(btrd_nr/i - btrd_r);
129 } while(i != m);
130 }
131 if(v <= f) return k;
132 else continue;
133 } else {
134 // final acceptance/rejection
135 v = log(v);
136 RealType rho =
137 (km/btrd_npq)*(((km/3. + 0.625)*km + 1./6)/btrd_npq + 0.5);
138 RealType t = -km*km/(2*btrd_npq);
139 if(v < t - rho) return k;
140 if(v > t + rho) continue;
141
142 IntType nm = _t - m + 1;
143 RealType h = (m + 0.5)*log((m + 1)/(btrd_r*nm))
144 + fc(m) + fc(_t - m);
145
146 IntType nk = _t - k + 1;
147 if(v <= h + (_t+1)*log(static_cast<RealType>(nm)/nk)
148 + (k + 0.5)*log(nk*btrd_r/(k+1))
149 - fc(k)
150 - fc(_t - k))
151 {
152 return k;
153 } else {
154 continue;
155 }
156 }
157 }
158}
159
160IntType invert(IntType t, RealType p, std::mt19937& urng = random_gen)
161{
162 RealType q = 1 - p;
163 RealType s = p / q;
164 RealType a = (t + 1) * s;
165 RealType r = pow((1 - p), static_cast<RealType>(t));
166 RealType u = stdunif(urng);
167 IntType x = 0;
168 while(u > r) {
169 u = u - r;
170 ++x;
171 RealType r1 = ((a/x) - s) * r;
172 // If r gets too small then the round-off error
173 // becomes a problem. At this point, p(i) is
174 // decreasing exponentially, so if we just call
175 // it 0, it's close enough. Note that the
176 // minimum value of q_n is about 1e-7, so we
177 // may need to be a little careful to make sure that
178 // we don't terminate the first time through the loop
179 // for float. (Hence the test that r is decreasing)
180 if(r1 < std::numeric_limits<RealType>::epsilon() && r1 < r) {
181 break;
182 }
183 r = r1;
184 }
185 return x;
186}
187
188
189IntType boost_binomial_distribution_variate(IntType t_arg, RealType p_arg, std::mt19937& urng = random_gen)
190{
191 bool other_side = p_arg > 0.5;
192 RealType fake_p = other_side ? 1.0 - p_arg : p_arg;
193 IntType m = static_cast<IntType>((t_arg+1)*fake_p);
194 IntType result;
195 if(m < 11)
196 result = invert(t_arg, fake_p, urng);
197 else
198 result = btrd(t_arg, fake_p, m, urng);
199
200 if(other_side)
201 return t_arg - result;
202 else
203 return result;
204}
205
206} // namespace IsoSpec