1 ```#include "include/legacy-distributions.h" ``` 2 3 4 ```static NPY_INLINE double legacy_double(aug_bitgen_t *aug_state) { ``` 5 1 ``` return aug_state->bit_generator->next_double(aug_state->bit_generator->state); ``` 6 ```} ``` 7 8 1 ```double legacy_gauss(aug_bitgen_t *aug_state) { ``` 9 1 ``` if (aug_state->has_gauss) { ``` 10 1 ``` const double temp = aug_state->gauss; ``` 11 1 ``` aug_state->has_gauss = false; ``` 12 1 ``` aug_state->gauss = 0.0; ``` 13 1 ``` return temp; ``` 14 ``` } else { ``` 15 ``` double f, x1, x2, r2; ``` 16 17 ``` do { ``` 18 1 ``` x1 = 2.0 * legacy_double(aug_state) - 1.0; ``` 19 1 ``` x2 = 2.0 * legacy_double(aug_state) - 1.0; ``` 20 1 ``` r2 = x1 * x1 + x2 * x2; ``` 21 1 ``` } while (r2 >= 1.0 || r2 == 0.0); ``` 22 23 ``` /* Polar method, a more efficient version of the Box-Muller approach. */ ``` 24 1 ``` f = sqrt(-2.0 * log(r2) / r2); ``` 25 ``` /* Keep for next call */ ``` 26 1 ``` aug_state->gauss = f * x1; ``` 27 1 ``` aug_state->has_gauss = true; ``` 28 1 ``` return f * x2; ``` 29 ``` } ``` 30 ```} ``` 31 32 1 ```double legacy_standard_exponential(aug_bitgen_t *aug_state) { ``` 33 ``` /* We use -log(1-U) since U is [0, 1) */ ``` 34 1 ``` return -log(1.0 - legacy_double(aug_state)); ``` 35 ```} ``` 36 37 1 ```double legacy_standard_gamma(aug_bitgen_t *aug_state, double shape) { ``` 38 ``` double b, c; ``` 39 ``` double U, V, X, Y; ``` 40 41 1 ``` if (shape == 1.0) { ``` 42 1 ``` return legacy_standard_exponential(aug_state); ``` 43 ``` } ``` 44 1 ``` else if (shape == 0.0) { ``` 45 ``` return 0.0; ``` 46 1 ``` } else if (shape < 1.0) { ``` 47 ``` for (;;) { ``` 48 1 ``` U = legacy_double(aug_state); ``` 49 1 ``` V = legacy_standard_exponential(aug_state); ``` 50 1 ``` if (U <= 1.0 - shape) { ``` 51 1 ``` X = pow(U, 1. / shape); ``` 52 1 ``` if (X <= V) { ``` 53 ``` return X; ``` 54 ``` } ``` 55 ``` } else { ``` 56 1 ``` Y = -log((1 - U) / shape); ``` 57 1 ``` X = pow(1.0 - shape + shape * Y, 1. / shape); ``` 58 1 ``` if (X <= (V + Y)) { ``` 59 ``` return X; ``` 60 ``` } ``` 61 ``` } ``` 62 ``` } ``` 63 ``` } else { ``` 64 1 ``` b = shape - 1. / 3.; ``` 65 1 ``` c = 1. / sqrt(9 * b); ``` 66 ``` for (;;) { ``` 67 ``` do { ``` 68 1 ``` X = legacy_gauss(aug_state); ``` 69 1 ``` V = 1.0 + c * X; ``` 70 1 ``` } while (V <= 0.0); ``` 71 72 1 ``` V = V * V * V; ``` 73 1 ``` U = legacy_double(aug_state); ``` 74 1 ``` if (U < 1.0 - 0.0331 * (X * X) * (X * X)) ``` 75 1 ``` return (b * V); ``` 76 1 ``` if (log(U) < 0.5 * X * X + b * (1. - V + log(V))) ``` 77 1 ``` return (b * V); ``` 78 ``` } ``` 79 ``` } ``` 80 ```} ``` 81 82 1 ```double legacy_gamma(aug_bitgen_t *aug_state, double shape, double scale) { ``` 83 1 ``` return scale * legacy_standard_gamma(aug_state, shape); ``` 84 ```} ``` 85 86 1 ```double legacy_pareto(aug_bitgen_t *aug_state, double a) { ``` 87 1 ``` return exp(legacy_standard_exponential(aug_state) / a) - 1; ``` 88 ```} ``` 89 90 1 ```double legacy_weibull(aug_bitgen_t *aug_state, double a) { ``` 91 1 ``` if (a == 0.0) { ``` 92 ``` return 0.0; ``` 93 ``` } ``` 94 1 ``` return pow(legacy_standard_exponential(aug_state), 1. / a); ``` 95 ```} ``` 96 97 1 ```double legacy_power(aug_bitgen_t *aug_state, double a) { ``` 98 1 ``` return pow(1 - exp(-legacy_standard_exponential(aug_state)), 1. / a); ``` 99 ```} ``` 100 101 1 ```double legacy_chisquare(aug_bitgen_t *aug_state, double df) { ``` 102 1 ``` return 2.0 * legacy_standard_gamma(aug_state, df / 2.0); ``` 103 ```} ``` 104 105 1 ```double legacy_noncentral_chisquare(aug_bitgen_t *aug_state, double df, ``` 106 ``` double nonc) { ``` 107 ``` double out; ``` 108 1 ``` if (nonc == 0) { ``` 109 1 ``` return legacy_chisquare(aug_state, df); ``` 110 ``` } ``` 111 1 ``` if (1 < df) { ``` 112 1 ``` const double Chi2 = legacy_chisquare(aug_state, df - 1); ``` 113 1 ``` const double n = legacy_gauss(aug_state) + sqrt(nonc); ``` 114 1 ``` return Chi2 + n * n; ``` 115 ``` } else { ``` 116 1 ``` const long i = random_poisson(aug_state->bit_generator, nonc / 2.0); ``` 117 1 ``` out = legacy_chisquare(aug_state, df + 2 * i); ``` 118 ``` /* Insert nan guard here to avoid changing the stream */ ``` 119 1 ``` if (npy_isnan(nonc)){ ``` 120 ``` return NPY_NAN; ``` 121 ``` } else { ``` 122 1 ``` return out; ``` 123 ``` } ``` 124 ``` } ``` 125 ```} ``` 126 127 1 ```double legacy_noncentral_f(aug_bitgen_t *aug_state, double dfnum, double dfden, ``` 128 ``` double nonc) { ``` 129 1 ``` double t = legacy_noncentral_chisquare(aug_state, dfnum, nonc) * dfden; ``` 130 1 ``` return t / (legacy_chisquare(aug_state, dfden) * dfnum); ``` 131 ```} ``` 132 133 1 ```double legacy_wald(aug_bitgen_t *aug_state, double mean, double scale) { ``` 134 ``` double U, X, Y; ``` 135 ``` double mu_2l; ``` 136 137 1 ``` mu_2l = mean / (2 * scale); ``` 138 1 ``` Y = legacy_gauss(aug_state); ``` 139 1 ``` Y = mean * Y * Y; ``` 140 1 ``` X = mean + mu_2l * (Y - sqrt(4 * scale * Y + Y * Y)); ``` 141 1 ``` U = legacy_double(aug_state); ``` 142 1 ``` if (U <= mean / (mean + X)) { ``` 143 ``` return X; ``` 144 ``` } else { ``` 145 1 ``` return mean * mean / X; ``` 146 ``` } ``` 147 ```} ``` 148 149 1 ```double legacy_normal(aug_bitgen_t *aug_state, double loc, double scale) { ``` 150 1 ``` return loc + scale * legacy_gauss(aug_state); ``` 151 ```} ``` 152 153 1 ```double legacy_lognormal(aug_bitgen_t *aug_state, double mean, double sigma) { ``` 154 1 ``` return exp(legacy_normal(aug_state, mean, sigma)); ``` 155 ```} ``` 156 157 1 ```double legacy_standard_t(aug_bitgen_t *aug_state, double df) { ``` 158 ``` double num, denom; ``` 159 160 1 ``` num = legacy_gauss(aug_state); ``` 161 1 ``` denom = legacy_standard_gamma(aug_state, df / 2); ``` 162 1 ``` return sqrt(df / 2) * num / sqrt(denom); ``` 163 ```} ``` 164 165 1 ```int64_t legacy_negative_binomial(aug_bitgen_t *aug_state, double n, double p) { ``` 166 1 ``` double Y = legacy_gamma(aug_state, n, (1 - p) / p); ``` 167 1 ``` return (int64_t)random_poisson(aug_state->bit_generator, Y); ``` 168 ```} ``` 169 170 1 ```double legacy_standard_cauchy(aug_bitgen_t *aug_state) { ``` 171 1 ``` return legacy_gauss(aug_state) / legacy_gauss(aug_state); ``` 172 ```} ``` 173 174 1 ```double legacy_beta(aug_bitgen_t *aug_state, double a, double b) { ``` 175 ``` double Ga, Gb; ``` 176 177 1 ``` if ((a <= 1.0) && (b <= 1.0)) { ``` 178 ``` double U, V, X, Y; ``` 179 ``` /* Use Johnk's algorithm */ ``` 180 181 ``` while (1) { ``` 182 1 ``` U = legacy_double(aug_state); ``` 183 1 ``` V = legacy_double(aug_state); ``` 184 1 ``` X = pow(U, 1.0 / a); ``` 185 1 ``` Y = pow(V, 1.0 / b); ``` 186 187 1 ``` if ((X + Y) <= 1.0) { ``` 188 1 ``` if (X + Y > 0) { ``` 189 1 ``` return X / (X + Y); ``` 190 ``` } else { ``` 191 1 ``` double logX = log(U) / a; ``` 192 1 ``` double logY = log(V) / b; ``` 193 1 ``` double logM = logX > logY ? logX : logY; ``` 194 1 ``` logX -= logM; ``` 195 1 ``` logY -= logM; ``` 196 197 1 ``` return exp(logX - log(exp(logX) + exp(logY))); ``` 198 ``` } ``` 199 ``` } ``` 200 ``` } ``` 201 ``` } else { ``` 202 1 ``` Ga = legacy_standard_gamma(aug_state, a); ``` 203 1 ``` Gb = legacy_standard_gamma(aug_state, b); ``` 204 1 ``` return Ga / (Ga + Gb); ``` 205 ``` } ``` 206 ```} ``` 207 208 1 ```double legacy_f(aug_bitgen_t *aug_state, double dfnum, double dfden) { ``` 209 1 ``` return ((legacy_chisquare(aug_state, dfnum) * dfden) / ``` 210 1 ``` (legacy_chisquare(aug_state, dfden) * dfnum)); ``` 211 ```} ``` 212 213 1 ```double legacy_exponential(aug_bitgen_t *aug_state, double scale) { ``` 214 1 ``` return scale * legacy_standard_exponential(aug_state); ``` 215 ```} ``` 216 217 218 1 ```static RAND_INT_TYPE legacy_random_binomial_original(bitgen_t *bitgen_state, ``` 219 ``` double p, ``` 220 ``` RAND_INT_TYPE n, ``` 221 ``` binomial_t *binomial) { ``` 222 ``` double q; ``` 223 224 1 ``` if (p <= 0.5) { ``` 225 1 ``` if (p * n <= 30.0) { ``` 226 1 ``` return random_binomial_inversion(bitgen_state, n, p, binomial); ``` 227 ``` } else { ``` 228 1 ``` return random_binomial_btpe(bitgen_state, n, p, binomial); ``` 229 ``` } ``` 230 ``` } else { ``` 231 1 ``` q = 1.0 - p; ``` 232 1 ``` if (q * n <= 30.0) { ``` 233 1 ``` return n - random_binomial_inversion(bitgen_state, n, q, binomial); ``` 234 ``` } else { ``` 235 1 ``` return n - random_binomial_btpe(bitgen_state, n, q, binomial); ``` 236 ``` } ``` 237 ``` } ``` 238 ```} ``` 239 240 241 1 ```int64_t legacy_random_binomial(bitgen_t *bitgen_state, double p, ``` 242 ``` int64_t n, binomial_t *binomial) { ``` 243 1 ``` return (int64_t) legacy_random_binomial_original(bitgen_state, p, ``` 244 ``` (RAND_INT_TYPE) n, ``` 245 ``` binomial); ``` 246 ```} ``` 247 248 249 1 ```static RAND_INT_TYPE random_hypergeometric_hyp(bitgen_t *bitgen_state, ``` 250 ``` RAND_INT_TYPE good, ``` 251 ``` RAND_INT_TYPE bad, ``` 252 ``` RAND_INT_TYPE sample) { ``` 253 ``` RAND_INT_TYPE d1, k, z; ``` 254 ``` double d2, u, y; ``` 255 256 1 ``` d1 = bad + good - sample; ``` 257 1 ``` d2 = (double)MIN(bad, good); ``` 258 259 1 ``` y = d2; ``` 260 1 ``` k = sample; ``` 261 1 ``` while (y > 0.0) { ``` 262 1 ``` u = next_double(bitgen_state); ``` 263 1 ``` y -= (RAND_INT_TYPE)floor(u + y / (d1 + k)); ``` 264 1 ``` k--; ``` 265 1 ``` if (k == 0) ``` 266 ``` break; ``` 267 ``` } ``` 268 1 ``` z = (RAND_INT_TYPE)(d2 - y); ``` 269 1 ``` if (good > bad) ``` 270 1 ``` z = sample - z; ``` 271 1 ``` return z; ``` 272 ```} ``` 273 274 ```/* D1 = 2*sqrt(2/e) */ ``` 275 ```/* D2 = 3 - 2*sqrt(3/e) */ ``` 276 ```#define D1 1.7155277699214135 ``` 277 ```#define D2 0.8989161620588988 ``` 278 1 ```static RAND_INT_TYPE random_hypergeometric_hrua(bitgen_t *bitgen_state, ``` 279 ``` RAND_INT_TYPE good, ``` 280 ``` RAND_INT_TYPE bad, ``` 281 ``` RAND_INT_TYPE sample) { ``` 282 ``` RAND_INT_TYPE mingoodbad, maxgoodbad, popsize, m, d9; ``` 283 ``` double d4, d5, d6, d7, d8, d10, d11; ``` 284 ``` RAND_INT_TYPE Z; ``` 285 ``` double T, W, X, Y; ``` 286 287 1 ``` mingoodbad = MIN(good, bad); ``` 288 1 ``` popsize = good + bad; ``` 289 1 ``` maxgoodbad = MAX(good, bad); ``` 290 1 ``` m = MIN(sample, popsize - sample); ``` 291 1 ``` d4 = ((double)mingoodbad) / popsize; ``` 292 1 ``` d5 = 1.0 - d4; ``` 293 1 ``` d6 = m * d4 + 0.5; ``` 294 1 ``` d7 = sqrt((double)(popsize - m) * sample * d4 * d5 / (popsize - 1) + 0.5); ``` 295 1 ``` d8 = D1 * d7 + D2; ``` 296 1 ``` d9 = (RAND_INT_TYPE)floor((double)(m + 1) * (mingoodbad + 1) / (popsize + 2)); ``` 297 1 ``` d10 = (random_loggam(d9 + 1) + random_loggam(mingoodbad - d9 + 1) + ``` 298 1 ``` random_loggam(m - d9 + 1) + random_loggam(maxgoodbad - m + d9 + 1)); ``` 299 1 ``` d11 = MIN(MIN(m, mingoodbad) + 1.0, floor(d6 + 16 * d7)); ``` 300 ``` /* 16 for 16-decimal-digit precision in D1 and D2 */ ``` 301 302 ``` while (1) { ``` 303 1 ``` X = next_double(bitgen_state); ``` 304 1 ``` Y = next_double(bitgen_state); ``` 305 1 ``` W = d6 + d8 * (Y - 0.5) / X; ``` 306 307 ``` /* fast rejection: */ ``` 308 1 ``` if ((W < 0.0) || (W >= d11)) ``` 309 1 ``` continue; ``` 310 311 1 ``` Z = (RAND_INT_TYPE)floor(W); ``` 312 1 ``` T = d10 - (random_loggam(Z + 1) + random_loggam(mingoodbad - Z + 1) + ``` 313 1 ``` random_loggam(m - Z + 1) + random_loggam(maxgoodbad - m + Z + 1)); ``` 314 315 ``` /* fast acceptance: */ ``` 316 1 ``` if ((X * (4.0 - X) - 3.0) <= T) ``` 317 ``` break; ``` 318 319 ``` /* fast rejection: */ ``` 320 1 ``` if (X * (X - T) >= 1) ``` 321 1 ``` continue; ``` 322 ``` /* log(0.0) is ok here, since always accept */ ``` 323 1 ``` if (2.0 * log(X) <= T) ``` 324 ``` break; /* acceptance */ ``` 325 ``` } ``` 326 327 ``` /* this is a correction to HRUA* by Ivan Frohne in rv.py */ ``` 328 1 ``` if (good > bad) ``` 329 1 ``` Z = m - Z; ``` 330 331 ``` /* another fix from rv.py to allow sample to exceed popsize/2 */ ``` 332 1 ``` if (m < sample) ``` 333 1 ``` Z = good - Z; ``` 334 335 1 ``` return Z; ``` 336 ```} ``` 337 ```#undef D1 ``` 338 ```#undef D2 ``` 339 340 1 ```static RAND_INT_TYPE random_hypergeometric_original(bitgen_t *bitgen_state, ``` 341 ``` RAND_INT_TYPE good, ``` 342 ``` RAND_INT_TYPE bad, ``` 343 ``` RAND_INT_TYPE sample) ``` 344 ```{ ``` 345 1 ``` if (sample > 10) { ``` 346 1 ``` return random_hypergeometric_hrua(bitgen_state, good, bad, sample); ``` 347 1 ``` } else if (sample > 0) { ``` 348 1 ``` return random_hypergeometric_hyp(bitgen_state, good, bad, sample); ``` 349 ``` } else { ``` 350 ``` return 0; ``` 351 ``` } ``` 352 ```} ``` 353 354 355 ```/* ``` 356 ``` * This is a wrapper function that matches the expected template. In the legacy ``` 357 ``` * generator, all int types are long, so this accepts int64 and then converts ``` 358 ``` * them to longs. These values must be in bounds for long and this is checked ``` 359 ``` * outside this function ``` 360 ``` * ``` 361 ``` * The remaining are included for the return type only ``` 362 ``` */ ``` 363 1 ```int64_t legacy_random_hypergeometric(bitgen_t *bitgen_state, int64_t good, ``` 364 ``` int64_t bad, int64_t sample) { ``` 365 1 ``` return (int64_t)random_hypergeometric_original(bitgen_state, ``` 366 ``` (RAND_INT_TYPE)good, ``` 367 ``` (RAND_INT_TYPE)bad, ``` 368 ``` (RAND_INT_TYPE)sample); ``` 369 ```} ``` 370 371 372 1 ```int64_t legacy_random_logseries(bitgen_t *bitgen_state, double p) { ``` 373 1 ``` return (int64_t)random_logseries(bitgen_state, p); ``` 374 ```} ``` 375 376 1 ``` int64_t legacy_random_poisson(bitgen_t *bitgen_state, double lam) { ``` 377 1 ``` return (int64_t)random_poisson(bitgen_state, lam); ``` 378 ```} ``` 379 380 1 ``` int64_t legacy_random_zipf(bitgen_t *bitgen_state, double a) { ``` 381 1 ``` return (int64_t)random_zipf(bitgen_state, a); ``` 382 ```} ``` 383 384 1 ``` int64_t legacy_random_geometric(bitgen_t *bitgen_state, double p) { ``` 385 1 ``` return (int64_t)random_geometric(bitgen_state, p); ``` 386 ```} ``` 387 388 1 ``` void legacy_random_multinomial(bitgen_t *bitgen_state, RAND_INT_TYPE n, ``` 389 ``` RAND_INT_TYPE *mnix, double *pix, npy_intp d, ``` 390 ``` binomial_t *binomial) { ``` 391 1 ``` return random_multinomial(bitgen_state, n, mnix, pix, d, binomial); ``` 392 ```} ```

Read our documentation on viewing source code .