IsoSpec  2.1.2
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 
25 namespace IsoSpec {
26 
27 typedef double RealType;
28 typedef int64_t IntType;
29 
30 
31 static 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!)
65 static 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 
77 IntType 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 
160 IntType 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 
189 IntType 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
IsoSpec
Definition: allocator.cpp:20