27 typedef double RealType;
28 typedef int64_t IntType;
31 static const RealType btrd_binomial_table[10] = {
65 static RealType fc(IntType k)
67 if(k < 10) {
return btrd_binomial_table[k]; }
70 RealType ikp1 = RealType(1) / (k + 1);
71 return (RealType(1)/12
73 - (RealType(1)/1260)*(ikp1*ikp1))*(ikp1*ikp1))*ikp1;
77 IntType btrd(IntType _t, RealType p, IntType m, std::mt19937& urng = random_gen)
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;
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));
104 u = stdunif(urng) - 0.5;
106 u = v/btrd_v_r - 0.93;
107 u = ((u < 0)? -0.5 : 0.5) - u;
108 v = stdunif(urng) * btrd_v_r;
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);
122 f = f*(btrd_nr/i - btrd_r);
128 v = v*(btrd_nr/i - btrd_r);
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;
142 IntType nm = _t - m + 1;
143 RealType h = (m + 0.5)*log((m + 1)/(btrd_r*nm))
144 + fc(m) + fc(_t - m);
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))
160 IntType invert(IntType t, RealType p, std::mt19937& urng = random_gen)
164 RealType a = (t + 1) * s;
165 RealType r = pow((1 - p),
static_cast<RealType
>(t));
166 RealType u = stdunif(urng);
171 RealType r1 = ((a/x) - s) * r;
180 if(r1 < std::numeric_limits<RealType>::epsilon() && r1 < r) {
189 IntType boost_binomial_distribution_variate(IntType t_arg, RealType p_arg, std::mt19937& urng = random_gen)
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);
196 result = invert(t_arg, fake_p, urng);
198 result = btrd(t_arg, fake_p, m, urng);
201 return t_arg - result;