1LoadFunctionLibrary("libv3/UtilityFunctions.bf"); 2 3/** @module random */ 4 5random.PI = 4*Arctan (1); 6 7random.LogFactorial.precomp = {{ 0.000000000000000, 8 0.000000000000000, 9 0.693147180559945, 10 1.791759469228055, 11 3.178053830347946, 12 4.787491742782046, 13 6.579251212010101, 14 8.525161361065415, 15 10.604602902745251, 16 12.801827480081469, 17 15.104412573075516, 18 17.502307845873887, 19 19.987214495661885, 20 22.552163853123421, 21 25.191221182738683, 22 27.899271383840894, 23 30.671860106080675, 24 33.505073450136891, 25 36.395445208033053, 26 39.339884187199495, 27 42.335616460753485, 28 45.380138898476908, 29 48.471181351835227, 30 51.606675567764377, 31 54.784729398112319, 32 58.003605222980518, 33 61.261701761002001, 34 64.557538627006323, 35 67.889743137181526, 36 71.257038967168000, 37 74.658236348830158, 38 78.092223553315307, 39 81.557959456115029, 40 85.054467017581516, 41 88.580827542197682, 42 92.136175603687079, 43 95.719694542143202, 44 99.330612454787428, 45 102.968198614513810, 46 106.631760260643450, 47 110.320639714757390, 48 114.034211781461690, 49 117.771881399745060, 50 121.533081515438640, 51 125.317271149356880, 52 129.123933639127240, 53 132.952575035616290, 54 136.802722637326350, 55 140.673923648234250, 56 144.565743946344900, 57 148.477766951773020, 58 152.409592584497350, 59 156.360836303078800, 60 160.331128216630930, 61 164.320112263195170, 62 168.327445448427650, 63 172.352797139162820, 64 176.395848406997370, 65 180.456291417543780, 66 184.533828861449510, 67 188.628173423671600, 68 192.739047287844900, 69 196.866181672889980, 70 201.009316399281570, 71 205.168199482641200, 72 209.342586752536820, 73 213.532241494563270, 74 217.736934113954250, 75 221.956441819130360, 76 226.190548323727570, 77 230.439043565776930, 78 234.701723442818260, 79 238.978389561834350, 80 243.268849002982730, 81 247.572914096186910, 82 251.890402209723190, 83 256.221135550009480, 84 260.564940971863220, 85 264.921649798552780, 86 269.291097651019810, 87 273.673124285693690, 88 278.067573440366120, 89 282.474292687630400, 90 286.893133295426990, 91 291.323950094270290, 92 295.766601350760600, 93 300.220948647014100, 94 304.686856765668720, 95 309.164193580146900, 96 313.652829949878990, 97 318.152639620209300, 98 322.663499126726210, 99 327.185287703775200, 100 331.717887196928470, 101 336.261181979198450, 102 340.815058870798960, 103 345.379407062266860, 104 349.954118040770250, 105 354.539085519440790, 106 359.134205369575340, 107 363.739375555563470, 108 368.354496072404690, 109 372.979468885689020, 110 377.614197873918670, 111 382.258588773060010, 112 386.912549123217560, 113 391.575988217329610, 114 396.248817051791490, 115 400.930948278915760, 116 405.622296161144900, 117 410.322776526937280, 118 415.032306728249580, 119 419.750805599544780, 120 424.478193418257090, 121 429.214391866651570, 122 433.959323995014870, 123 438.712914186121170, 124 443.475088120918940, 125 448.245772745384610, 126 453.024896238496130, 127 457.812387981278110, 128 462.608178526874890, 129 467.412199571608080, 130 472.224383926980520, 131 477.044665492585580, 132 481.872979229887900, 133 486.709261136839360, 134 491.553448223298010, 135 496.405478487217580, 136 501.265290891579240, 137 506.132825342034830, 138 511.008022665236070, 139 515.890824587822520, 140 520.781173716044240, 141 525.679013515995050, 142 530.584288294433580, 143 535.496943180169520, 144 540.416924105997740, 145 545.344177791154950, 146 550.278651724285620, 147 555.220294146894960, 148 560.169054037273100, 149 565.124881094874350, 150 570.087725725134190, 151 575.057539024710200, 152 580.034272767130800, 153 585.017879388839220, 154 590.008311975617860, 155 595.005524249382010, 156 600.009470555327430, 157 605.020105849423770, 158 610.037385686238740, 159 615.061266207084940, 160 620.091704128477430, 161 625.128656730891070, 162 630.172081847810200, 163 635.221937855059760, 164 640.278183660408100, 165 645.340778693435030, 166 650.409682895655240, 167 655.484856710889060, 168 660.566261075873510, 169 665.653857411105950, 170 670.747607611912710, 171 675.847474039736880, 172 680.953419513637530, 173 686.065407301994010, 174 691.183401114410800, 175 696.307365093814040, 176 701.437263808737160, 177 706.573062245787470, 178 711.714725802289990, 179 716.862220279103440, 180 722.015511873601330, 181 727.174567172815840, 182 732.339353146739310, 183 737.509837141777440, 184 742.685986874351220, 185 747.867770424643370, 186 753.055156230484160, 187 758.248113081374300, 188 763.446610112640200, 189 768.650616799717000, 190 773.860102952558460, 191 779.075038710167410, 192 784.295394535245690, 193 789.521141208958970, 194 794.752249825813460, 195 799.988691788643450, 196 805.230438803703120, 197 810.477462875863580, 198 815.729736303910160, 199 820.987231675937890, 200 826.249921864842800, 201 831.517780023906310, 202 836.790779582469900, 203 842.068894241700490, 204 847.352097970438420, 205 852.640365001133090, 206 857.933669825857460, 207 863.231987192405430, 208 868.535292100464630, 209 873.843559797865740, 210 879.156765776907600, 211 884.474885770751830, 212 889.797895749890240, 213 895.125771918679900, 214 900.458490711945270, 215 905.796028791646340, 216 911.138363043611210, 217 916.485470574328820, 218 921.837328707804890, 219 927.193914982476710, 220 932.555207148186240, 221 937.921183163208070, 222 943.291821191335660, 223 948.667099599019820, 224 954.046996952560450, 225 959.431492015349480, 226 964.820563745165940, 227 970.214191291518320, 228 975.612353993036210, 229 981.015031374908400, 230 986.422203146368590, 231 991.833849198223450, 232 997.249949600427840, 233 1002.670484599700300, 234 1008.095434617181700, 235 1013.524780246136200, 236 1018.958502249690200, 237 1024.396581558613400, 238 1029.838999269135500, 239 1035.285736640801600, 240 1040.736775094367400, 241 1046.192096209724900, 242 1051.651681723869200, 243 1057.115513528895000, 244 1062.583573670030100, 245 1068.055844343701400, 246 1073.532307895632800, 247 1079.012946818975000, 248 1084.497743752465600, 249 1089.986681478622400, 250 1095.479742921962700, 251 1100.976911147256000, 252 1106.478169357800900, 253 1111.983500893733000, 254 1117.492889230361000, 255 1123.006317976526100, 256 1128.523770872990800, 257 1134.045231790853000, 258 1139.570684729984800, 259 1145.100113817496100, 260 1150.633503306223700, 261 1156.170837573242400}}; 262 263lfunction random.LogFactorial(n) { 264 if (n > 254) { 265 x = n + 1; 266 return (x - 0.5) * Log (x) - x + 0.5 * Log (2 * ^"random.PI") + 1.0 / (12.0 * x); 267 } else { 268 return (^"random.LogFactorial.precomp")[n]; 269 } 270 } 271 272lfunction random.number (from, to) { 273 return Random (from,to); 274} 275 276lfunction random.exponential (lambda) { 277 return -Log (Random (0,1)) / lambda; 278} 279 280lfunction random.poisson (lambda) { 281 if (lambda < 30) { 282 L = Exp(-lambda); 283 k = 0; 284 p = 1.; 285 286 do { 287 k += 1; 288 p = p * Random(0,1); 289 290 } while (p > L); 291 292 return k - 1; 293 } else { 294 c = 0.767 - 3.36 / lambda; 295 beta = ^"random.PI" / Sqrt (3.0 * lambda); 296 alpha = beta * lambda; 297 k = Log(c) - lambda - Log(beta); 298 299 while (1) { 300 301 u = Random (0,1); 302 if (u > 0.0) { 303 x = (alpha - Log ((1.0 - u) / u)) / beta; 304 n = (x + 0.5)$1; 305 if (n < 0) { 306 continue; 307 } 308 v = Random (0,1); 309 if (v > 0.0) { 310 y = alpha - beta * x; 311 t = (1.0 + Exp(y)); 312 lhs = y + Log (v / (t*t)); 313 rhs = k + n * Log (lambda) - LogFactorial (n); 314 if (lhs <= rhs) 315 return n; 316 } 317 } 318 } 319 } 320} 321 322/*-------------------------------------------------*/ 323 324function random.gamma (alpha,beta) { 325 _sampledValue = Random (0,1); 326 FindRoot (_k,CGammaDist(_x_,alpha,beta)-_sampledValue,_x_,0,2500); 327 return _k; 328} 329 330/*-------------------------------------------------*/ 331 332function random.beta (p,q) { 333 _sampledValue = Random (0,1); 334 FindRoot (_k,IBeta(_x_,p,q)-_sampledValue,_x_,0,1); 335 return _k; 336} 337 338/*-------------------------------------------------*/ 339 340lfunction random.normal.standard () { 341 /* Box-Muller */ 342 return Sqrt (-2*Log (Random (1e-26,1)))*Cos (2 * ^"random.PI" * Random (1e-26,1)); 343} 344 345/*-------------------------------------------------*/ 346 347lfunction random.binomial (p, N) { 348 349 if (N == 0) { 350 return 0; 351 } 352 353 if (p == 1) { 354 return N; 355 } 356 q = -Log (1-p); 357 x = 0; 358 sum = 0; 359 360 do { 361 U = Random (0,1); 362 if (U > 0) { 363 if (N == x) { 364 x += 1; 365 break; 366 } 367 sum += - Log (U) / (N-x); 368 x += 1; 369 } 370 } while (sum < q); 371 372 return x - 1; 373 } 374 375 /*-------------------------------------------------*/ 376 377lfunction random.gamma_fast (alpha) { 378 // the rejection sampling algorithm with squeeze from 379 // "A Simple Method for Generating Gamma Variables" 380 // ACM Transactions on Mathematical Software, Vol. 26, No. 3, September 2000, Pages 363–372. 381 382 383 if (alpha < 1) { 384 // use the trick that Gamma (alpha) ~ Gamma (alpha+1) * Uniform (1/alpha) 385 return random.gamma_fast (alpha+1) * Random (0,1) ^ (1/alpha); 386 } else { 387 388 d = alpha - 1/3; 389 c = 1/Sqrt (9*d); 390 391 while (TRUE) { 392 x = random.normal.standard (); 393 v = (1+c*x); 394 while (v <= 0) { 395 x = random.normal.standard (); 396 v = (1+c*x); 397 } 398 v = v^3; 399 u = Random (0,1); 400 if (u < (1-0.0331 * x^4)) { 401 return d * v; 402 } 403 if (Log (u) < 0.5 * x * x + d * (1-v+Log(v))) { 404 return d * v; 405 } 406 } 407 } 408} 409 410 /*-------------------------------------------------*/ 411 412lfunction random.dirichlet (alpha) { 413 // alpha is the vector of > 0 concentration parameters 414 415 return Transpose (Random (alpha, {"PDF" : "Dirichlet"})); 416 417} 418 419lfunction random.TRUE_or_FALSE () { 420 return Random(0,1) < 0.5; 421} 422