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