1 /* random.c
2
3 Copyright (c) 1992-2018 Free Software Foundation, Inc.
4
5 This file is part of GNU MCSim.
6
7 GNU MCSim is free software; you can redistribute it and/or
8 modify it under the terms of the GNU General Public License
9 as published by the Free Software Foundation; either version 3
10 of the License, or (at your option) any later version.
11
12 GNU MCSim is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16
17 You should have received a copy of the GNU General Public License
18 along with GNU MCSim; if not, see <http://www.gnu.org/licenses/>
19
20
21 * Random number generator module: Two versions are provided. One uses the
22 GNU Scientific Library (GSL). The other is independent of GSL.
23
24 The GSL version provides:
25 Randoms(), which in fact call gsl_rng_uniform()
26 InitRandom(rank, seed, f), must be called to initializes the generator
27 SetSeed(seed), which sets the seed of Randoms()
28 GetSeed(seed), not implemented in fact and exits with an error message
29
30 The independent (non-GSL) code provides two basic uniform random number
31 generators and associated functions:
32 Randoms(), yielding uniform random variates between 0 and 1
33 InitRandom(rank, seed, f), which initializes the seed
34 SetSeed(seed), which sets the seed of Randoms()
35 GetSeed(seed), which gets the current seed
36
37 * Other types of random variates generators are available. They call
38 Randoms() if they require a uniform variate:
39
40 BetaRandom(alpha, beta, a, b) -- Beta(alpha, beta) over [a, b]
41 BinomialBetaRandom(E, a, b) -- BetaBinomial(n = E + E * a / b)
42 BinomialRandom(p, n) -- Binomial of n trials, P(each) = p
43 CauchyRandom(s) -- Cauchy, with scale s
44 Chi2Random(dof) -- Chi-squared w/dof degrees of freedom
45 ExpRandom(beta) -- Exponential of inverse scale beta
46 GammaRandom(alpha) -- Gamma variate
47 GenLogNormalRandom(m, sn, sln) -- Generalized log transformation of the
48 two-component error model
49 GGammaRandom(alpha, beta) -- General gamma variate
50 InvGGammaRandom(alpha, beta) -- General inverse gamma variate
51 LogNormalRandom(m, s) -- exp(Normal)
52 LogUniformRandom(a, b) -- LogUniform over [a, b]
53 Multinomial(...) -- Multinomial variates
54 NegativeBinomialRandom(r, p) -- Negative Binomial variates
55 NormalRandom(m, s) -- General Normal
56 PiecewiseRandom(min, a, b, max) -- Draws from a Mayan pyramid !
57 PiecewiseVariate(n, x[], p[], Cdf[], o, p) -- Draws from a tabulated PDF
58 PoissonRandom(mu) -- Poisson with rate mu
59 StudentTRandom(dof, m, s) -- Student-t with dof degrees of freedom
60 location m and scale s
61 TruncInvGGammaRandom(alpha, beta, a, b) -- Truncated general inv-gamma
62 TruncLogNormalRandom(m, s, a, b) -- Truncated log normal
63 TruncNormalRandom(m, s, a, b) -- Truncated general normal
64 UniformRandom(a, b) -- Uniform over [a, b]
65 WishartRandom(...) -- Wishart(matrix) variates
66
67 * And utility functions:
68
69 and() -- Logical "and" function
70 CDFNormal(z) -- Integral of the normal at z
71 CalcCumulative(n, x[], p[], Cdf[], o) -- Constructs a CDF given a PDF
72 DFNormal(x, mu, sd) -- Normal density
73 erfc(x) -- Complementary error function
74 leq() -- Logical lower or equal function
75 lnDFNormal(x, mu, sd) -- Log of the normal density
76 lnGamma(x) -- Natural log of gamma function
77 lnDFBeta(x, alpha, beta, min, max) -- Log of beta density
78 piecewise(...) -- Piecewise function
79
80 */
81
82 #include <float.h>
83 #include <assert.h>
84 #include <limits.h>
85 #include <math.h>
86 #include <stdio.h>
87 #include <stdlib.h>
88
89 #include "hungtype.h"
90 #include "lex.h"
91 #include "lexerr.h"
92 #include "random.h"
93
94
95 /* ----------------------------------------------------------------------------
96 Globals/Externals
97 */
98
99 static BOOL vbSwitchGauss = FALSE; /* Flag to reset switchG in NormalRandom */
100
101
102 /* ----------------------------------------------------------------------------
103 We have two versions of Randoms() depending whether we use gsl or not
104 */
105
106 #ifdef HAVE_LIBGSL
107
108 #include <gsl/gsl_rng.h>
109
110 /* ----------------------------------------------------------------------------
111 GSL version, global definitions, private
112 */
113
114 static const gsl_rng_type * genType;
115 static gsl_rng * rGenerator;
116
117 /* ----------------------------------------------------------------------------
118 GetSeed
119
120 Not available if GSL is used. Issue an error message.
121 */
GetSeed(void)122 double GetSeed (void)
123 {
124 printf ("Error: the GetSeed() function is not available if GSL is used "
125 "- Exiting\n\n");
126 exit (0);
127 }
128
129
130 /* ----------------------------------------------------------------------------
131 SetSeed
132
133 Sets vRandRec.seed to given dSeed. Note that the seed is a long int.
134 */
SetSeed(double dSeed)135 void SetSeed (double dSeed)
136 {
137
138 gsl_rng_set(rGenerator, (unsigned long int) dSeed);
139 vbSwitchGauss = FALSE; /* Force a reset of the normal variate sampler */
140
141 } /* SetSeed */
142
143
144 /* ----------------------------------------------------------------------------
145 Randoms, gsl version
146 */
Randoms(void)147 double Randoms (void)
148 {
149
150 return(gsl_rng_uniform(rGenerator));
151
152 } /* Randoms, gsl version */
153
154
155 /* ----------------------------------------------------------------------------
156 InitRandom, gsl version
157
158 Selects a generator and sets the gsl seed to given dSeed, silently
159 corrects an invalid dSeed.
160 */
InitRandom(int rank,double dSeed,int rdm_gen_name)161 void InitRandom (int rank, double dSeed, int rdm_gen_name)
162 {
163 static int bInit = 0;
164
165 if (!bInit) { /* initialize */
166
167 /* set the type of random generator, in fact we should use an enum */
168 switch (rdm_gen_name) {
169 case TRUE:
170 genType = gsl_rng_default; /* mt19937 */
171 break;
172 case FALSE:
173 genType = gsl_rng_taus2;
174 break;
175 default:
176 genType = gsl_rng_default;
177 }
178
179 /* create an instance of a generator of type genType */
180 rGenerator = gsl_rng_alloc(genType);
181
182 /* set the seed */
183 gsl_rng_set(rGenerator, (unsigned long int) dSeed);
184
185 vbSwitchGauss = FALSE; /* Force a reset of the normal variate sampler */
186
187 /* next calls to InitRandoms will do nothing */
188 bInit = 1;
189 }
190
191 } /* InitRandom, gsl version */
192
193
194 #else /* non-gsl version */
195
196 /* ----------------------------------------------------------------------------
197 Non-GSL version, global definitions, private
198 */
199
200 typedef struct tagRANDREC {
201 double seed;
202 } RANDREC, *PRANDREC;
203
204
205 static BOOL vbNoSeed = TRUE; /* Flag to prevent use without seed */
206 static BOOL vbNotInitd = TRUE; /* Flag to prevent use without initing */
207 static RANDREC vRandRec; /* Global random information shared by */
208 /* all random number functions */
209
210 /* ----------------------------------------------------------------------------
211 GetSeed
212
213 Returns the current value of vRandRec.seed.
214 */
GetSeed(void)215 double GetSeed (void)
216 {
217 return (vRandRec.seed);
218 }
219
220
221 /* ----------------------------------------------------------------------------
222 SetSeed
223
224 Sets vRandRec.seed to given dSeed, silently corrects an invalid dSeed.
225 */
SetSeed(double dSeed)226 void SetSeed (double dSeed)
227 {
228 int bCorrected = 0;
229
230 if (dSeed == 0.0) {
231 dSeed = SEED_DEFAULT;
232 bCorrected++;
233 }
234
235 if (dSeed < 0)
236 dSeed = -dSeed; /* Don't announce this correction */
237
238 if (dSeed < SEED_MIN) {
239 dSeed = SEED_MIN + (dSeed/SEED_MIN) / (SEED_MAX-SEED_MIN);
240 bCorrected++;
241 }
242
243 if (dSeed > SEED_MAX) {
244 dSeed = SEED_MIN + (SEED_MAX/dSeed) / (SEED_MAX-SEED_MIN);
245 bCorrected++;
246 }
247
248 assert ((/* Invalid Seed */ dSeed >= SEED_MIN && dSeed <= SEED_MAX));
249
250 /* Assign valid seed */
251
252 if (bCorrected)
253 printf ("SetSeed(): corrected out of range random number seed\n"
254 "Seed must lie in the range [%g, %g]\n"
255 "New seed --> %g\n", SEED_MIN, SEED_MAX, dSeed);
256
257 vRandRec.seed = dSeed;
258 vbNoSeed = FALSE; /* Flag that seed has been set */
259 vbSwitchGauss = FALSE; /* Force a reset of the normal variate sampler */
260
261 } /* SetSeed */
262
263
264 /* ----------------------------------------------------------------------------
265 InitRandom
266
267 initializes the random generator with the given seed.
268 If an invalid seed is given, SetSeed() silently corrects it.
269
270 If the boolean bWarmUp is non-zero, the random number generator is
271 "warmed up" by running it a number of times.
272 Also, a flag used by the Normal() routine is initialized.
273 */
InitRandom(int rank,double dSeed,int bWarmUp)274 void InitRandom (int rank, double dSeed, int bWarmUp)
275 {
276 long i;
277
278 /* Prevent nuking user's seed if not initd */
279 if (vbNoSeed || dSeed != SEED_DEFAULT)
280 SetSeed(dSeed);
281
282 if (bWarmUp) {
283 /* Warm up generator */
284 for (i = 0; i < 50; i++) (void) Randoms();
285
286 vbNotInitd = FALSE; /* Flag as initialized */
287 }
288
289 } /* InitRandom, non-gsl version */
290
291
292 /* ----------------------------------------------------------------------------
293 Randoms
294
295 An alternative random number generator, so you don't have to use
296 the (probably not so good) system supplied standard C version.
297
298 Randoms() returns random numbers between 0 and 1. The minimum
299 returned value is 1/m and the maximum 1 - 1/m. The generator can
300 be initialized with InitRandom(). If not, it auto-initializes its seed.
301
302 This generator should be correct on any system for which the
303 representattion of reals uses at least a 32-bit mantissa, including
304 the sign bit.
305
306 From PARK SK, MILLER KW: Random number generators: good ones are
307 hard to find. Commun. ACM 1988; 31: 1192. (Version Real 2).
308 */
Randoms(void)309 double Randoms (void)
310 {
311 #define a 16807.0
312 #define m 2147483647.0
313 #define q 127773.0 /* m Div a */
314 #define r 2836.0 /* m Mod a */
315
316 double hi, test;
317
318 if (vbNoSeed)
319 SetSeed (SEED_DEFAULT);
320
321 hi = (long)(vRandRec.seed / q);
322 test = a * (vRandRec.seed - q * hi) - r * hi;
323
324 if (test > 0.0)
325 vRandRec.seed = test;
326 else
327 vRandRec.seed = test + m;
328
329 return (vRandRec.seed / m);
330
331 #undef a
332 #undef m
333 #undef q
334 #undef r
335
336 } /* Randoms, non-gsl version */
337
338 #endif
339
340 /* End of the non-gsl random uniform generators */
341
342
343 /* ----------------------------------------------------------------------------
344 BetaRandom
345
346 returns a variate that is Beta distributed on the interval [a,b]
347 with shape parameters alpha and beta.
348
349 The Beta function has two shaping parameters, alpha and beta.
350 Setting these parameters to 1.5 and 1.5 yields a normal-like
351 distribution, but without tails. If alpha and beta are equal to
352 1 it is a uniform distribution.
353
354 If alpha and beta are less than 1, use a rejection algorithm;
355 Otherwise use the fact that if x is distributed Gamma(alpha) and y
356 Gamma(beta) then x/(x+y) is Beta(alpha, beta).
357
358 For the rejection algorithm first a Beta variate is found over the
359 interval [0, 1] with not the most efficient algorithm. This is then
360 scaled at the end to desired range.
361
362 It may be tempting to re-use the second number drawn as the first
363 random number of the next iteration, and simply draw one more.
364 *** Don't do it. You will produce an incorrect distribution. You
365 must draw two new numbers for the rejection sampling to be correct.
366
367 References:
368 - Ripley, Stochastic Simulations, John Wiley and Sons, 1987, p 90.
369 - J.H.Maindonald, Statistical Computation, John Wiley and Sons,
370 1984, p 370.
371 */
BetaRandom(double alpha,double beta,double a,double b)372 double BetaRandom (double alpha, double beta, double a, double b)
373 {
374 double u1, u2, w;
375
376 if (b <= a || alpha <= 0 || beta <= 0) {
377 printf ("Error: bad shape or range for a beta variate - Exiting\n\n");
378 exit (0);
379 }
380
381 if ((alpha < 1) && (beta < 1))
382 /* use rejection */
383 do {
384 u1 = Randoms(); /* Draw two numbers */
385 u2 = Randoms();
386
387 u1 = pow(u1, 1/alpha); /* alpha and beta are > 0 */
388 u2 = pow(u2, 1/beta);
389
390 w = u1 + u2;
391
392 } while (w > 1.0);
393
394 else {
395 /* use relation to Gamma */
396 u1 = GammaRandom(alpha);
397 u2 = GammaRandom(beta);
398 w = u1 + u2;
399 }
400
401 return (a + (u1/w) * (b - a)); /* Scale to interval [a, b] */
402
403 } /* BetaRandom */
404
405
406 /* ----------------------------------------------------------------------------
407 BinomialBetaRandom
408
409 Return as a double floating-point number an integer value that is a random
410 deviate drawn from a binomial distribution of n trials each of
411 probability p, p being beta distributed with parameters alpha and beta.
412 I use the expectation in input. The classical N is equal to
413 E + E * beta / alpha.
414 See Bernardo & Smith "Bayesian Theory"
415 */
BinomialBetaRandom(double Expectation,double alpha,double beta)416 double BinomialBetaRandom (double Expectation, double alpha, double beta)
417 {
418 double dTmp = Expectation + Expectation * beta / alpha;
419
420 if (dTmp < LONG_MAX)
421 /* parameters will be checked in BinomialRandom */
422 return BinomialRandom (BetaRandom (alpha, beta, 0, 1), (long) dTmp);
423 else {
424 printf("BinomialBetaRandom: N (= %g) too large - Exiting...", dTmp);
425 exit (0);
426 }
427
428 } /* BinomialBetaRandom */
429
430
431 /* ----------------------------------------------------------------------------
432 BinomialRandom
433
434 Return as a double floating-point number an integer value that is a random
435 deviate drawn from a binomial distribution of n trials each of
436 probability p, using Randoms() as a source of uniform random deviates.
437 Adapted from the algorithm described in the book Numerical Recipes by
438 Press et al.
439 */
BinomialRandom(double p,long N)440 double BinomialRandom (double p, long N)
441 {
442 long j;
443 static long iOldN = -1;
444 double dAngle, dDeviate, dMean, dPtemp, dSqrt, dTangent, dTemp1, dTemp2;
445 static double dLnFactN, dPold = -1, dLnP, dQ, dLnQ;
446
447 if (p < 0 || p > 1 || N < 0) {
448 printf ("Error: parameters out of bounds for a binomial variate "
449 "- Exiting\n\n");
450 exit (0);
451 }
452
453 dPtemp = ( p <= 0.5 ? p : 1 - p);
454 dMean = N * dPtemp; /* mean of the deviate to be produced. */
455
456 /* Use the direct method if N is not too large.
457 This can require up to 25 calls to random */
458
459 if (N < 25) {
460 dDeviate = 0;
461 for (j = 0; j < N; j++)
462 if (Randoms() < dPtemp)
463 dDeviate = dDeviate + 1;
464 }
465 else
466 if (dMean < 1) {
467 /* if less than one event is expected out of 25 or more trials,then the
468 distribution is quite accurately Poisson. Use direct method. */
469 dTemp1 = exp(-dMean);
470 dTemp2 = 1.0;
471 for (j = 0; j <= N; j++) {
472 dTemp2 = dTemp2 * Randoms();
473 if (dTemp2 < dTemp1) break;
474 }
475
476 dDeviate = (j <= N ? j : N);
477 }
478 else { /* Use rejection */
479
480 if (N != iOldN) {
481 /* if N has changed or it's the first call, initialize */
482 dLnFactN = lnGamma((double) N + 1);
483 iOldN = N;
484 }
485
486 if (dPtemp != dPold) {
487 /* if dPtemp has changed or it's the first call, initialize. */
488 dPold = dPtemp;
489 dQ = 1 - dPtemp;
490 dLnP = log(dPtemp);
491 dLnQ = log(dQ);
492 } /* if */
493
494 dSqrt = sqrt(2 * dMean * dQ);
495
496 /* Rejection method with a Lorentzian comparison function. */
497
498 do {
499 do {
500 dAngle = PI * Randoms();
501 dTangent = tan(dAngle);
502 dTemp1 = dSqrt * dTangent + dMean;
503 } while (dTemp1 < 0 || dTemp1 >= (N + 1)); /* Reject */
504
505 dTemp1 = floor(dTemp1); /* discrete distribution */
506
507 dTemp2 = 1.2 * dSqrt * (1 + dTangent * dTangent) *
508 exp(dLnFactN - lnGamma(dTemp1 + 1) -
509 lnGamma(N - dTemp1 + 1) +
510 dTemp1 * dLnP + (N - dTemp1) * dLnQ);
511
512 } while (Randoms() > dTemp2);
513
514 /* Reject on average about 1.5 time per deviate */
515
516 dDeviate = dTemp1;
517
518 } /* else */ /* end of rejection */
519
520 if (dPtemp != p)
521 dDeviate = N - dDeviate; /* undo the symmetry tranformation */
522
523 return (dDeviate);
524
525 } /* BinomialRandom */
526
527
528 /* ----------------------------------------------------------------------------
529 CauchyRandom
530
531 Returns a random variate from Cauchy's distribution (generated by dividing a
532 standard Normal by a Chi-square with 1 degree of freedom.
533 */
CauchyRandom(double dScale)534 double CauchyRandom (double dScale)
535 {
536 double z, x;
537
538 z = NormalRandom(0, dScale);
539 x = GGammaRandom (0.5, 0.5); /* Chi-2, see Chi2Random */
540
541 return (z / sqrt(x));
542
543 } /* CauchyRandom */
544
545
546 /* ----------------------------------------------------------------------------
547 Chi2Random
548
549 returns a chi-squared random variate, which is a gamma(dof/2, 1/2).
550 */
Chi2Random(double dof)551 double Chi2Random (double dof)
552 {
553
554 return (GGammaRandom (dof / 2.0, 0.5));
555
556 } /* Chi2Random */
557
558
559 /* ----------------------------------------------------------------------------
560 ExpRandom
561
562 returns an exponential variate with inverse scale beta
563
564 Algorithm 3.2 from Ripley "Stochastic Simulations" Wiley 1987, p. 55.
565 */
ExpRandom(double beta)566 double ExpRandom (double beta)
567 {
568
569 if (beta <= 0) {
570 printf ("Error: negative or null inverse scale for an exponential variate "
571 "- Exiting\n\n");
572 exit (0);
573 }
574
575 return -log(Randoms()) / beta;
576
577 } /* ExpRandom */
578
579
580 /* ----------------------------------------------------------------------------
581 GammaRandom
582
583 returns a gamma distributed random variate with shape parameter
584 alpha.
585
586 If alpha < 1 uses algorithm 3.19 of Ripley; if alpha > 1 uses
587 algorithm 3.20 of Ripley; if alpha is 1 uses returns an
588 ExpRandom variate.
589
590 Reference:
591 - Ripley, Stochastic Simulations, John Wiley and Sons, 1987, pp 88-90.
592 */
GammaRandom(double alpha)593 double GammaRandom (double alpha)
594 {
595 #define E 2.718281828459
596 static double aprev = 0.0, c1, c2, c3, c4, c5;
597 double b, u1, u2, w, x;
598
599 if (alpha <= 0) {
600 printf ("Error: negative or null shape parameter for a gamma variate "
601 "- Exiting\n\n");
602 exit (0);
603 }
604 else if (alpha < 1) {
605
606 b = (alpha + E) / E;
607
608 do {
609 u1 = b * Randoms();
610 if (u1 <= 1.0) {
611 x = pow(u1, 1./alpha);
612 /* problem: if alpha is too small, x will be about zero and that's
613 bad, particularly for inverse-gamma variates.
614 Fixed by blocking zeros - FB 5/11/1997 */
615 if ((x > DBL_MIN) && (x <= -log(Randoms())))
616 return(x);
617 }
618 else {
619 x = -log((b - u1) / alpha);
620 if (pow(x, alpha - 1) >= Randoms()) return(x);
621 }
622 } while (1 == 1);
623
624 } /* end if alpha < 1 */
625
626 else {
627 if (alpha > 1) {
628
629 if (alpha != aprev) {
630 /* initialize */
631 aprev = alpha;
632 c1 = alpha - 1;
633 b = 1.0 / c1;
634 c2 = b * (alpha - (1 / (6.0 * alpha)));
635 c3 = 2 * b;
636 c4 = c3 + 2.0;
637 if (alpha > 2.5)
638 c5 = 1.0 / sqrt(alpha);
639 }
640
641 do {
642 do {
643 u1 = Randoms();
644 u2 = Randoms();
645 if (alpha > 2.5)
646 u1 = u2 + c5 * (1 - 1.86 * u1);
647 } while ((u1 >= 1) || ( u1 <= 0));
648
649 w = c2 * u2 / u1;
650 if (((c3 * u1 + w + 1 / w) <= c4) ||
651 ((c3 * log(u1) - log(w) + w) < 1))
652 return(c1 * w);
653 } while (1 == 1);
654
655 }
656 else
657 return ExpRandom(1.0);
658 }
659
660 #undef E
661
662 } /* GammaRandom */
663
664
665 /* ----------------------------------------------------------------------------
666 GenLogNormalRandom
667
668 Returns a variate drawn from the generalized lognormal distribution,
669 which is an approximation of the two-component error model of Rocke
670 and Lorenzato 1995. At low values of the mean, it is approximately
671 normal and at high values of the mean, it is approximately lognormal
672
673 dMean and dSDNorn are in natural space.
674 Note that here SDLogNorm is sigma in log space; this convention
675 is different from that for the LogNormal.
676 */
GenLogNormalRandom(double dMean,double dSDNorm,double dSDLogNorm)677 double GenLogNormalRandom (double dMean, double dSDNorm,
678 double dSDLogNorm)
679 {
680 double dmuz, dSLogNorm, dLambda, dz;
681
682 if (dMean < 0) { /* "True value must be >= 0 */
683 char str[10];
684 sprintf(str, "%5.2e", dMean);
685 ReportRunTimeError(NULL, RE_BADLOGNORMALMEAN | RE_FATAL,
686 "", str, "GenLogNormalRandom");
687 }
688 else
689 if (dSDLogNorm <= 0) {
690 char str[10];
691 sprintf(str, "%5.2e", dSDLogNorm);
692 ReportRunTimeError(NULL, RE_BADLOGNORMALSD | RE_FATAL,
693 "", str, "GenLogNormalRandom");
694 }
695
696 /* Relative Standard Deviation of Lognormal */
697 dSLogNorm = sqrt(exp(pow(dSDLogNorm,2)) * (exp(pow(dSDLogNorm,2)) - 1));
698 /* Transformation parameter */
699 dLambda = pow(dSDNorm/dSLogNorm,2);
700 /* Transformation of mean */
701 dmuz = log(dMean + sqrt(pow(dMean,2) + dLambda));
702 /* Random variate of transformed variable */
703 dz = NormalRandom(dmuz, dSLogNorm);
704
705 return (exp(dz) - dLambda * exp(-dz))/2; /* return untransformed variable */
706
707 } /* GenLogNormalRandom */
708
709
710 /* ----------------------------------------------------------------------------
711 GGammaRandom
712
713 Returns a gamma distributed random variate with shaping parameter
714 alpha and inverse scale (rate) parameter beta (> 0).
715 */
GGammaRandom(double alpha,double beta)716 double GGammaRandom (double alpha, double beta)
717 {
718
719 if (beta <= 0) {
720 printf ("Error: negative or null inverse scale for a gamma variate "
721 "- Exiting\n\n");
722 exit (0);
723 }
724 return GammaRandom(alpha) / beta;
725
726 } /* GGammaRandom */
727
728
729 /* ----------------------------------------------------------------------------
730 InvGGammaRandom
731
732 Returns an inverse gamma distributed random variate with shape parameter
733 alpha and scale parameter beta.
734 This just gets a general gamma variate and returns its inverse
735 See Gelman et al. "Bayesian Data Analysis"
736 */
InvGGammaRandom(double alpha,double beta)737 double InvGGammaRandom (double alpha, double beta)
738 {
739
740 /* parameters will be checked in GGammaRandom */
741 if (beta <= 0) {
742 printf ("Error: negative or null scale for an inverse gamma variate "
743 "- Exiting\n\n");
744 exit (0);
745 }
746
747 return beta / GammaRandom(alpha);
748
749 } /* InvGGammaRandom */
750
751
752 /* ----------------------------------------------------------------------------
753 LogNormalRandom
754
755 returns a variate such that the log of the variate is normally
756 distributed.
757 Uses the geometric mean and geometric standard deviation.
758 */
LogNormalRandom(double dGeoMean,double dGeoSD)759 double LogNormalRandom (double dGeoMean, double dGeoSD)
760 {
761
762 if (dGeoMean <= 0) {
763 char str[10];
764 sprintf(str, "%5.2e", dGeoMean);
765 ReportRunTimeError(NULL, RE_BADLOGNORMALMEAN | RE_FATAL,
766 "", str, "LogNormalRandom");
767 }
768 else {
769 if (dGeoSD < 1) {
770 char str[10];
771 sprintf(str, "%5.2e", dGeoSD);
772 ReportRunTimeError(NULL, RE_BADLOGNORMALSD | RE_FATAL,
773 "", str, "LogNormalRandom");
774 }
775 }
776
777 return exp(NormalRandom (log(dGeoMean), log(dGeoSD)));
778
779 } /* LogNormalRandom */
780
781
782 /* ----------------------------------------------------------------------------
783 LogUniformRandom
784
785 returns a variate that is log-uniformly distributed on the interval
786 [a,b].
787 */
LogUniformRandom(double a,double b)788 double LogUniformRandom (double a, double b)
789 {
790
791 if (b < a) {
792 printf ("Error: bad range a for uniform variate - Exiting\n\n");
793 exit (0);
794 }
795
796 return ( a * pow(b/a, Randoms()) );
797
798 } /* LogUniformRandom */
799
800
801 /* ----------------------------------------------------------------------------
802 Multinomial
803
804 Generates multinomial deviates.
805 N is the number of trials,
806 p the array of probabilities,
807 dim the dimension of the array (number of possible events),
808 x the array of event occurences which is returned.
809
810 From Devroye "Non-Uniform Random Numbers...".
811 */
Multinomial(long n,int dim,double * p,double * x)812 void Multinomial (long n, int dim, double *p, double *x)
813 {
814 int i;
815 double sum, ptemp;
816
817 sum = 1;
818
819 for (i = 1; i <= dim; i++) {
820 if (p[i]) {
821 ptemp = p[i] / sum;
822 x[i] = BinomialRandom (ptemp, n);
823 n = n - (long)x[i];
824 sum = sum - p[i];
825 }
826 else x[i] = 0.0;
827
828 } /* for */
829
830 } /* Multinomial */
831
832
833 /* ----------------------------------------------------------------------------
834 NegativeBinomialRandom
835
836 Generates negative binomial deviates.
837 r (> 0) is number of failures until experiments are stopped
838 p (in [0,1]) is the success probability in each experiment (real)
839
840 Adapted from Devroye "Non-Uniform Random Numbers..." page 543.
841 */
NegativeBinomialRandom(double r,double p)842 long NegativeBinomialRandom (double r, double p)
843 {
844 if (p < 0) {
845 printf ("Error: parameter p negative for a negative binomial variate "
846 "- Exiting\n\n");
847 exit (0);
848 }
849
850 if (p >= 1) {
851 printf ("Error: parameter p >= 1 for a negative binomial variate "
852 "- Exiting\n\n");
853 exit (0);
854 }
855
856 if (r < 0) {
857 printf ("Error: parameter r negative for a negative binomial variate "
858 "- Exiting\n\n");
859 exit (0);
860 }
861
862 if (r == 0 || p == 0) return(0);
863
864 return(PoissonRandom(GGammaRandom(r, (1 - p) / p)));
865
866 } /* NegativeBinomialRandom */
867
868
869 /* ----------------------------------------------------------------------------
870 NormalRandom
871
872 Returns a Normal random variate based on a unit variate,
873 using a random generator as a source of uniform deviates.
874 Adapted from the algorithm described in the book Numerical Recipes by
875 Press et al.
876 */
NormalRandom(double dMean,double dSD)877 double NormalRandom (double dMean, double dSD)
878 {
879 double dRacine, dTemp1, dTemp2, dTemp3;
880 static double memGauss;
881
882 if (vbSwitchGauss) { /* used stored variate */
883 vbSwitchGauss = FALSE;
884 return (dMean + dSD * memGauss);
885 }
886
887 do {
888 dTemp1 = 2 * Randoms() - 1;
889 dTemp2 = 2 * Randoms() - 1;
890 dRacine = dTemp1 * dTemp1 + dTemp2 * dTemp2;
891 } while ((dRacine >= 1) || (dRacine == 0));
892
893 dTemp3 = sqrt(-2 * log(dRacine) / dRacine);
894 vbSwitchGauss = TRUE;
895 memGauss = dTemp1 * dTemp3;
896 return (dMean + dSD * (dTemp2 * dTemp3));
897
898 } /* NormalRandom */
899
900
901 /* ----------------------------------------------------------------------------
902 PiecewiseRandom
903 __
904 returns a variate that is distributed along a / \ shaped distribution.
905 */
PiecewiseRandom(double min,double a,double b,double max)906 double PiecewiseRandom (double min, double a, double b, double max)
907 {
908 double dTemp;
909 static double Grille[4];
910 static double densite[4];
911 static double densiteCum[4];
912 double nvlle_densite;
913
914 Grille[0] = min;
915 Grille[1] = a;
916 Grille[2] = b;
917 Grille[3] = max;
918 densite[0] = 0;
919 densite[1] = 1/(max/2+b/2-a/2-min/2);
920 densite[2] = 1/(max/2+b/2-a/2-min/2);
921 densite[3] = 0;
922
923 CalcCumulative (4, Grille, densite, densiteCum, 1);
924
925 dTemp = PiecewiseVariate (4, Grille, densite, densiteCum, 1,
926 &nvlle_densite);
927 return (dTemp);
928
929 } /* PiecewiseRandom */
930
931
932 /* ----------------------------------------------------------------------------
933 PiecewiseVariate
934
935 Returns a variate drawn by tabulated inversion from the
936 cumulative, Cdf[], calculated to order iOrder
937 (0 = piecewise-constant, etc.)
938
939 Inputs: dimension of the table, x values, cdf values at x,
940 pdf values at x, order of the interpolation.
941
942 Returns the sampled variate as its value. If a pointer is given,
943 the value of pdf[] at the sampled variate is returned in *pVal_pdf.
944
945 Note: For piecewise-constant variates, the grid is corrected with
946 CorrectPWConstantGrid() to center the intervals around
947 sampled points.
948 */
PiecewiseVariate(long cDim,double rg_x[],double rg_pdf[],double rg_Cdf[],int iOrder,double * pVal_pdf)949 double PiecewiseVariate (long cDim, double rg_x[], double rg_pdf[],
950 double rg_Cdf[], int iOrder, double *pVal_pdf)
951 {
952 double dPWVariate; /* the variate chosen */
953 double dValPdf; /* the value of the pdf at variate */
954 double dUniform = UniformRandom(0, rg_Cdf[cDim - 1]);
955 long lUpper, lLower, lIndex;
956
957 if (iOrder > 1) {
958 printf ("CalcCumulative: Order %d not supported"
959 "-> using piecewise-linear\n", iOrder);
960 iOrder = 1;
961 }
962
963 /* Find bounding Xs by a binary search of the ordered rg_x's
964 */
965 lUpper = cDim;
966 lLower = 0;
967 lIndex = 0;
968
969 while (lUpper - lLower > 1) {
970 lIndex = (lUpper + lLower)/2;
971
972 if (dUniform > rg_Cdf[lIndex]) lLower = lIndex; /* Move to right half */
973 else
974 if (dUniform < rg_Cdf[lIndex])
975 lUpper = lIndex; /* Move to left half */
976 else lUpper = lLower = lIndex;
977 }
978
979 /* If we are exactly on a cumulative data point (unlikely),
980 the value of the pdf is known and the variate is a grid point.
981 */
982 if (lUpper == lLower) {
983 dValPdf = rg_pdf[lLower];
984 dPWVariate = rg_x[lLower];
985 }
986
987 /* Otherwise we do the appropriate interpolation
988 */
989 else
990 switch (iOrder) {
991
992 /* Piecewise-constant pdf: the Cdf is piecewise-linear
993 */
994 case 0:
995 dValPdf = rg_pdf[lLower];
996 dPWVariate = InterpolateX (rg_x, rg_Cdf, lLower, dUniform);
997 break;
998
999 /* Piecewise-linear pdf: the Cdf is piecewise-quadratic
1000 */
1001 case 1: {
1002
1003 if (rg_pdf[lLower] == rg_pdf[lUpper]) { /* A linear segment */
1004 dValPdf = rg_pdf[lLower];
1005 dPWVariate = InterpolateX (rg_x, rg_Cdf, lLower, dUniform);
1006 }
1007
1008 else { /* Interpolate a quadratic */
1009
1010 double a, b, c, dRadical;
1011
1012 /* Find a, b, and c from the quadratic equation.
1013 a is guaranteed not zero from if() above. */
1014
1015 a = (rg_pdf[lUpper] - rg_pdf[lLower]) /(rg_x[lUpper] - rg_x[lLower]);
1016
1017 b = rg_pdf[lLower] - a*rg_x[lLower];
1018
1019 c = rg_Cdf[lLower] - (a*rg_x[lLower]/2.0 + b) * rg_x[lLower];
1020
1021 dRadical = sqrt(b*b - 2*a*(c - dUniform));
1022
1023 dPWVariate = (-b + dRadical) / a;
1024
1025 assert (dPWVariate >= rg_x[lLower] && dPWVariate <= rg_x[lUpper]);
1026
1027 dValPdf = a * dPWVariate + b;
1028
1029 if (a > 0)
1030 assert (dValPdf >= rg_pdf[lLower] && dValPdf <= rg_pdf[lUpper]);
1031 else
1032 assert (dValPdf <= rg_pdf[lLower] && dValPdf >= rg_pdf[lUpper]);
1033
1034 } /* else */
1035
1036 } /* case block */
1037 break;
1038
1039 default:
1040 dValPdf = 0;
1041 dPWVariate = 0;
1042 assert(0);
1043 break;
1044
1045 } /* switch */
1046
1047 if (pVal_pdf) *pVal_pdf = dValPdf; /* Return the value if requested */
1048
1049 return dPWVariate;
1050
1051 } /* PiecewiseVariate */
1052
1053
1054 /* ----------------------------------------------------------------------------
1055 PoissonRandom
1056
1057 returns a Poisson random variate, with rate mu.
1058
1059 If mu is less than 60, uses inversion; otherwise uses the rejection
1060 method of Atkinson (as presented by Ripley "Stochastic Simulations",
1061 Wiley 1987, p 79).
1062 */
PoissonRandom(double mu)1063 long PoissonRandom(double mu)
1064 {
1065
1066 double u1, x, u2, lnfact, s, t;
1067 static double prev_mu = 0, c, beta, alpha, k;
1068 long n = 0;
1069
1070 if (mu <= 0) {
1071 printf("Error: negative or null rate for a Poisson variate "
1072 "- Exiting\n\n");
1073 exit(0);
1074 }
1075
1076 if (mu <= 60) {
1077 /* inversion */
1078 s = 1;
1079 t = 1;
1080 u1 = Randoms() * exp(mu);
1081 while(s < u1){
1082 n++;
1083 t = t * mu / n;
1084 s = s + t;
1085 }
1086 }
1087 else {
1088 /* rejection */
1089 if (mu != prev_mu) {
1090 c = 0.767 - 3.36 / mu;
1091 beta = PI / sqrt(3 * mu);
1092 alpha = beta * mu;
1093 k = log(c) - mu - log(beta);
1094 }
1095
1096 do {
1097 do {
1098 u1 = Randoms();
1099 x = (alpha - log((1 - u1) / u1)) / beta;
1100 } while (x <= -0.5);
1101
1102 n = (long)(x + 0.5);
1103 u2 = Randoms();
1104
1105 /* calculate log n factorial using Stirling's formula */
1106 lnfact = 0.918938533 - n + (n + 0.5) * log(n);
1107 } while (alpha - beta * x + log(u2 / pow((1 + exp(alpha - beta * x)), 2))
1108 > k + n * log(mu) - lnfact);
1109 }
1110
1111 return n;
1112
1113 } /* PoissonRandom */
1114
1115
1116 /* ----------------------------------------------------------------------------
1117 StudentTRandom
1118
1119 Returns a random variate from Student T distribution with dof
1120 degrees of freedom, location dMean, and scale dSD.
1121 */
StudentTRandom(double dof,double dMean,double dSD)1122 double StudentTRandom (double dof, double dMean, double dSD)
1123 {
1124
1125 double z, x;
1126
1127 if (dof <= 0) {
1128 printf ("Error: StudentTRandom: dof <= 0\n");
1129 exit (0);
1130 }
1131 z = NormalRandom(0,1);
1132 x = Chi2Random(dof);
1133
1134 return (dMean + dSD*z*sqrt(dof/x));
1135
1136 } /* StudentTRandom */
1137
1138
1139 /* ----------------------------------------------------------------------------
1140 TruncInvGGammaRandom
1141
1142 returns a truncated general inverse gamme variate in the range [a, b].
1143 */
TruncInvGGammaRandom(double alpha,double beta,double a,double b)1144 double TruncInvGGammaRandom (double alpha, double beta, double a, double b)
1145 {
1146 double X = 0.0;
1147 int iter = 0;
1148
1149 if (a >= b)
1150 printf ("TruncLogNormalRandom: min >= max [%g %g]\n", a, b);
1151
1152 else do {
1153 if(++iter == 25) {
1154 printf("TruncInvGGammaRandom: problem with range: ");
1155 printf("min %g, max %g, alpha %g, beta %g\n", a, b, alpha, beta);
1156 }
1157 X = InvGGammaRandom(alpha, beta);
1158 }
1159 while (X < a || X > b);
1160
1161 return X;
1162
1163 } /* TruncInvGGammaRandom */
1164
1165
1166 /* ----------------------------------------------------------------------------
1167 TruncLogNormalRandom
1168
1169 returns a truncated LogNormal variate in the range [a, b].
1170 All parameters are in natural space.
1171 */
TruncLogNormalRandom(double dGeoMean,double dGeoSD,double a,double b)1172 double TruncLogNormalRandom (double dGeoMean, double dGeoSD,
1173 double a, double b)
1174 {
1175
1176 return exp(TruncNormalRandom(log(dGeoMean), log(dGeoSD), log(a), log(b)));
1177
1178 } /* TruncLogNormalRandom */
1179
1180
1181 /* ----------------------------------------------------------------------------
1182 TruncLogNormalRandom_old
1183
1184 returns a truncated LogNormal variate in the range [a, b].
1185 All parameters are in natural space.
1186 Deprecated.
1187 */
TruncLogNormalRandom_old(double dGeoMean,double dGeoSD,double a,double b)1188 double TruncLogNormalRandom_old (double dGeoMean, double dGeoSD,
1189 double a, double b)
1190 {
1191 double X = 0.0;
1192 int iter = 0;
1193
1194 if (a >= b)
1195 printf ("TruncLogNormalRandom: min >= max [%g %g]\n", a, b);
1196
1197 else do {
1198 if(++iter == 25) {
1199 printf("TruncLogNormalRandom: problem with range: ");
1200 printf("min %g, max %g, ave %g, sd %g\n", a, b, dGeoMean, dGeoSD);
1201 }
1202 X = LogNormalRandom(dGeoMean, dGeoSD);
1203 }
1204 while (X < a || X > b);
1205
1206 return X;
1207
1208 } /* TruncLogNormalRandom_old */
1209
1210
1211 /* ----------------------------------------------------------------------------
1212 TruncNormalRandom
1213
1214 returns a truncated Normal variate in the range [a, b].
1215 Uses different rejection sampling algorithms, depending on where the
1216 upper and lower limits lie. See Robert, Stat. Comp (1995) 5:121-5.
1217 */
TruncNormalRandom(double dMean,double dSD,double a,double b)1218 double TruncNormalRandom (double dMean, double dSD, double a, double b)
1219 {
1220 /* the algorithm works on mean 0, sd 1 scale */
1221 double lower = (a - dMean) / dSD;
1222 double upper = (b - dMean) / dSD;
1223 double dTmp, y;
1224
1225 if (a >= b) {
1226 printf ("Error: TruncNormalRandom: min >= max [%g %g]\n", a, b);
1227 exit (0);
1228 }
1229
1230 if (((lower < 0) && (upper > 0) && (upper - lower > SQRT_2PI))) {
1231 /* wide range: use the standard normal proposal */
1232 do {
1233 y = NormalRandom(0, 1);
1234 } while ((y < lower) || (y > upper));
1235 }
1236 else {
1237 dTmp = lower + sqrt(lower * lower + 4);
1238 if ((lower >= 0) &&
1239 (upper > lower + TWO_SQRTEXP1 / dTmp *
1240 exp(lower * (2 - sqrt(lower*lower + 4)) * 0.25))) {
1241 /* lower >> mean: rejection sampling with exponential proposal */
1242 dTmp = dTmp * 0.5;
1243 do {
1244 y = ExpRandom(dTmp) + lower;
1245 } while ((Randoms() > exp(-(y - dTmp) * (y - dTmp) * 0.5)) ||
1246 (y > upper));
1247 }
1248 else {
1249 dTmp = -upper + sqrt(upper * upper + 4);
1250 if ((upper <= 0) &&
1251 (lower < upper - TWO_SQRTEXP1 / dTmp *
1252 exp(upper * (2 + sqrt(upper*upper + 4)) * 0.25))) {
1253 /* upper << mean: rejection sampling with exponential proposal */
1254 dTmp = dTmp * 0.5;
1255 do {
1256 y = ExpRandom(dTmp) - upper;
1257 } while ((Randoms() > exp(-(y - dTmp) * (y - dTmp) * 0.5)) ||
1258 (y > -lower));
1259 y = -y;
1260 }
1261 else {
1262 /* range narrow and central: rejection with uniform proposal */
1263 do {
1264 y = UniformRandom(lower, upper);
1265 if (lower > 0) {
1266 dTmp = exp((lower * lower - y * y) * 0.5);
1267 }
1268 else {
1269 if (upper < 0) {
1270 dTmp = exp((upper * upper - y * y) * 0.5);
1271 } else {
1272 dTmp = exp(-y * y * 0.5);
1273 }
1274 }
1275 } while (Randoms() > dTmp);
1276 }
1277 }
1278 }
1279
1280 return(y * dSD + dMean);
1281
1282 } /* TruncNormalRandom */
1283
1284
1285 /* ----------------------------------------------------------------------------
1286 TruncNormalRandom_old
1287
1288 returns a truncated Normal variate in the range [a, b].
1289 Deprecated.
1290 */
TruncNormalRandom_old(double dMean,double dSD,double a,double b)1291 double TruncNormalRandom_old (double dMean, double dSD, double a, double b)
1292 {
1293 double X = 0.0;
1294 double unif_density, normConstant, imp_ratio;
1295
1296 if (a >= b) {
1297 printf ("Error: TruncNormalRandom: min >= max [%g %g]\n", a, b);
1298 exit (0);
1299 }
1300
1301 else { /* proceed */
1302 if ((b - a) / dSD > 1.5) { /* wide truncation range relative to SD */
1303 do { /* rejection sampling from normal proposal */
1304 X = NormalRandom(dMean, dSD);
1305 } while (X < a || X > b);
1306 }
1307 else { /* narrow truncation range relative to SD */
1308 /* normalization constant of the truncated normal */
1309 normConstant = CDFNormal((b - dMean) / dSD) -
1310 CDFNormal((a - dMean) / dSD);
1311 /* density of any X under a normalized uniform */
1312 unif_density = (dMean < a ?
1313 DFNormal(a, dMean, dSD) / normConstant :
1314 (dMean < b ?
1315 DFNormal(dMean, dMean, dSD) / normConstant :
1316 DFNormal(b, dMean, dSD) / normConstant));
1317 do { /* rejection sampling from uniform proposal */
1318 X = UniformRandom (a, b);
1319 /* density of X under the normalized truncated normal */
1320 imp_ratio = DFNormal(X, dMean, dSD) /
1321 (normConstant * unif_density);
1322 } while ((imp_ratio < 1) && (Randoms() > imp_ratio));
1323 }
1324 }
1325
1326 return X;
1327
1328 } /* TruncNormalRandom_old */
1329
1330
1331 /* ----------------------------------------------------------------------------
1332 UniformRandom
1333
1334 returns a variate that is uniformly distributed on the interval [a,b].
1335 */
UniformRandom(double a,double b)1336 double UniformRandom (double a, double b)
1337 {
1338
1339 if (b < a) {
1340 printf ("Error: bad range a for uniform variate - Exiting\n\n");
1341 exit (0);
1342 }
1343
1344 return (Randoms() * (b - a) + a);
1345
1346 } /* UniformRandom */
1347
1348
1349 /* ----------------------------------------------------------------------------
1350 Wishart
1351
1352 samples a matrix according to the Wishart distribution by the method
1353 of Odell and Feiveson (1966).
1354
1355 Paramters are:
1356 n (degrees of freedom); p (dimension of Wishart matrix);
1357 t (pointer to a Cholesky decomposition of a covariance matrix);
1358 w (pointer to the sampled Wishart matrix, in
1359 triangular form; work (pointer to a work space, of length p*p).
1360
1361 Triangular matrices are stored in order
1362 0 1 3
1363 2 4
1364 5 etc.
1365 */
WishartRandom(long n,long p,double * t,double * w,double * work)1366 void WishartRandom (long n, long p, double *t, double *w, double *work)
1367 {
1368 double eta, sum;
1369 long i, j, k, m, k1, k2, k3;
1370
1371 printf ("WishartRandom not tested - Exiting...");
1372 exit(0);
1373
1374 /* generate random quantities for Bartlett's decomposition */
1375 for (j = 0, k = 0; j < p; j++) {
1376 for (i = 0; i < j; i++)
1377 w[k++] = NormalRandom(0, 1);
1378
1379 /* Chi-square with n-i degrees of freedom */
1380 w[k++] = GGammaRandom((n - i) / 2.0, 0.5);
1381 }
1382
1383 /* generate a standard Wishart */
1384 for (j = p - 1, m = k - 1, k2 = (p * (p - 1)) / 2; j >= 0; k2 = k2 - (j--)) {
1385 eta = w[m];
1386 for (i = j, k1 = (i * (i + 1)) / 2; i >= 0; k1 = k1 - (i--), m--) {
1387 for (k = 0, sum = 0.0; k < i; k++)
1388 sum = sum + w[k1+k] * w[k2+k];
1389
1390 if (i == j)
1391 w[m] = sum + eta;
1392 else
1393 w[m] = sum + sqrt(eta) * w[m];
1394 }
1395 }
1396
1397 /* form product L * W * L' */
1398 for (i = 0, k1 = 0, m = 0; i < p; k1 = k1 + (++i)) {
1399 for (j = 0, k2 = 0; j < p; k2 = k2 + (++j), m++) {
1400 for (k = 0, sum = 0.0; k < j; k++)
1401 sum = sum + t[k1+k] * w[k2+k];
1402
1403 for (k = j, k3 = j; k <= i; k3 = k3 + (++k))
1404 sum = sum + t[k1+k] * w[k2+k3];
1405
1406 work[m] = sum;
1407 }
1408 }
1409
1410 for (i = 0, m = 0, k1 = 0; i < p; i++, k1 = k1 + p) {
1411 for (j = 0, k2 = 0; j <= i; k2 = k2 + (++j), m++) {
1412 for (k = 0, sum = 0.0; k <= j; k++)
1413 sum = sum + work[k1+k] * t[k2+k];
1414
1415 w[m] = sum;
1416 }
1417 }
1418
1419 } /* WishartRandom */
1420
1421
1422 /* ----------------------------------------------------------------------------
1423 Utility functions
1424 */
1425
1426 /* ----------------------------------------------------------------------------
1427 Boolean "AND" logical function.
1428
1429 Return TRUE if both its two arguments are TRUE. This is just for
1430 compatibility with SBML.
1431 */
and(BOOL A,BOOL B)1432 BOOL and (BOOL A, BOOL B)
1433 {
1434 return (A && B);
1435 } /* piecewise */
1436
1437
1438 /* ----------------------------------------------------------------------------
1439 CalcCumulative
1440
1441 Approximates to an iOrder the cumulative distribution rg_Cdf
1442 given a sampling grid (rg_x) of dimension cDim, and the
1443 sampled pdf (rg_pdf) points.
1444
1445 Supports piecewise-constant (order 0) and piecewise-linear
1446 (order 1) pdfs.
1447 */
CalcCumulative(long cDim,double * rg_x,double * rg_pdf,double * rg_Cdf,int iOrder)1448 void CalcCumulative (long cDim, double *rg_x, double *rg_pdf,
1449 double *rg_Cdf, int iOrder)
1450 {
1451 long i; /* Index for the samples */
1452
1453 if (iOrder > 1) {
1454 printf ("CalcCumulative: Order %d not supported"
1455 "-> using piecewise-linear\n", iOrder);
1456 iOrder = 1;
1457 }
1458
1459 rg_Cdf[0] = 0.0; /* Cumulative starts at 0.0 */
1460 switch (iOrder) {
1461
1462 /* Piecewise Constant: sum of rectangles */
1463 case 0:
1464 for (i = 1; i < cDim; i++)
1465 rg_Cdf[i] = rg_Cdf[i-1] + rg_pdf[i]*(rg_x[i] - rg_x[i-1]);
1466 break;
1467
1468 /* Piecewise Linear: sum of trapezoids */
1469 case 1:
1470 for (i = 1; i < cDim; i++)
1471 rg_Cdf[i] = rg_Cdf[i-1] + ((rg_x[i] - rg_x[i - 1]) *
1472 (rg_pdf[i] + rg_pdf[i - 1]) / 2);
1473 break;
1474
1475 default:
1476 assert (0); /* This is an error condition */
1477 break;
1478
1479 } /* switch */
1480
1481 } /* CalcCumulative */
1482
1483
1484 /* ----------------------------------------------------------------------------
1485 CDFNormal
1486
1487 the probability for [-inf;Z] under the normal distribution
1488 */
CDFNormal(double z)1489 double CDFNormal (double z)
1490 {
1491 register double tmp;
1492
1493 tmp = z/SQRT_2;
1494
1495 /* avoid roundoff errors generated by subtracting a small number from
1496 a much bigger one */
1497 if (tmp >= 0)
1498 return ( 0.5 * (2 - erfc(tmp)) );
1499 else
1500 return ( 0.5 * erfc(-tmp));
1501 }
1502
1503
1504 /* ----------------------------------------------------------------------------
1505 erfc
1506
1507 the complementary error function of z (= 1 - erf(z))
1508
1509 Adapted from the algorithm described in the book Numerical Recipes by
1510 Press et al.
1511 */
erfc(double x)1512 double erfc (double x)
1513 {
1514 double dAbsX, t, dVal;
1515
1516 dAbsX = fabs(x);
1517 if (dAbsX > 20) { /* FB 16/2/99) */
1518 return ( x >= 0 ? 0 : 2 );
1519 }
1520 else {
1521 t = 1 / (1 + 0.5 * dAbsX);
1522 dVal = t * exp(-dAbsX*dAbsX - 1.26551223 + t*(1.00002368 + t*(0.37409196 +
1523 t*(0.09678418 + t*(-0.18628806 + t*(0.27886807 + t*(-1.13520398 +
1524 t*(1.48851587 + t*(-0.82215223 + t*(0.17087277))))))))));
1525 return ( x >= 0 ? dVal : 2 - dVal );
1526 }
1527 } /* erfc */
1528
1529
1530 /* ----------------------------------------------------------------------------
1531 DFNormal
1532 Normal density function
1533 */
DFNormal(double x,double mu,double sd)1534 double DFNormal (double x, double mu, double sd)
1535 {
1536 if (sd <= 0.0) {
1537 printf ("Error: negative or null SD in DFNormal\n");
1538 exit (0);
1539 }
1540
1541 double tmp = (mu - x) / sd;
1542 return ( (INV_SQRT_2PI / sd) * exp(-0.5 * tmp * tmp) );
1543 }
1544
1545
1546 /* ----------------------------------------------------------------------------
1547 InterpolateX
1548
1549 Do a linear interpolation to return x
1550 */
InterpolateX(double rgX[],double rgY[],long lLower,double dY)1551 double InterpolateX (double rgX[], double rgY[], long lLower, double dY)
1552 {
1553 return rgX[lLower] + (dY - rgY[lLower]) *
1554 (rgX[lLower + 1] - rgX[lLower]) /
1555 (rgY[lLower + 1] - rgY[lLower]);
1556
1557 } /* InterpolateX */
1558
1559
1560 /* ----------------------------------------------------------------------------
1561 lnDFBeta
1562 the log of the beta density function
1563 FB 08/07/97
1564 */
lnDFBeta(double x,double alpha,double beta,double min,double max)1565 double lnDFBeta (double x, double alpha, double beta, double min, double max)
1566 {
1567 if (max <= min) {
1568 printf ("Error: bad range for beta variate in lnDFBeta\n");
1569 exit (0);
1570 }
1571 if (alpha <= 0) {
1572 printf ("Error: bad alpha for beta variate in LnDensity\n");
1573 exit (0);
1574 }
1575 if (beta <= 0) {
1576 printf ("Error: bad beta for beta variate in LnDensity\n");
1577 exit (0);
1578 }
1579
1580 x = (x - min) / (max - min);
1581 return (alpha - 1) * log (x) + (beta - 1) * log (1 - x) +
1582 lnGamma (alpha + beta) - lnGamma (alpha) - lnGamma(beta) -
1583 log (max - min);
1584 }
1585
1586
1587 /* ----------------------------------------------------------------------------
1588 lnDFNormal
1589 log of the normal density function
1590 */
lnDFNormal(double x,double mu,double sd)1591 double lnDFNormal (double x, double mu, double sd)
1592 {
1593 if (sd <= 0.0) {
1594 printf ("Error: negative or null SD in lnDFNormal\n");
1595 exit (0);
1596 }
1597
1598 double tmp = (mu - x) / sd;
1599 return ( -LOG_SQRT_2PI - log(sd) - 0.5 * tmp * tmp );
1600 }
1601
1602
1603 /* ----------------------------------------------------------------------------
1604 lnGamma
1605
1606 A function to return the natural log of the Gamma function of x.
1607 Adapted from the algorithm described in the book Numerical Recipes by
1608 Press et al.
1609 It can be used to compute factorials since ln(n!) = lnGamma(n + 1)
1610 */
lnGamma(double x)1611 double lnGamma (double x)
1612 {
1613 double dSeries, dTemp;
1614
1615 if (x <= 0.0) {
1616 printf ("Error: negative or null parameter for lnGamma function\n");
1617 exit (0);
1618 }
1619
1620 dSeries = 1.000000000190015 +
1621 76.18009172947146 / x -
1622 86.50532032141677 / (x + 1) +
1623 24.01409824083091 / (x + 2) -
1624 1.231739572450155 / (x + 3) +
1625 1.20865097386617E-3 / (x + 4) -
1626 5.39523938495E-6 / (x + 5);
1627
1628 dTemp = x + 4.5;
1629 dTemp = -dTemp + (x - 0.5) * log (dTemp) + log (2.50662827465 * dSeries);
1630 return dTemp;
1631
1632 } /* lnGamma */
1633
1634
1635 #ifdef ndef
1636 /* ----------------------------------------------------------------------------
1637 Piecewise math function.
1638
1639 Takes an unknwn number of pairs of (double, boolean), and eventually a
1640 terminal double. Returns the first double value for which the associated
1641 condition is true. If none is true it returns FALSE or the terminal double
1642 if one is given.
1643 Example: {double am = 1., pm = 2., noon = 3., time = 55.888;
1644 return piecewise (am, (time < 12), pm, (time > 12), noon);}
1645 */
piecewise(double X1...)1646 void piecewise (double X1 ...)
1647 {
1648 va_list argptr;
1649
1650 va_start(argprt, X1);
1651
1652 if ((arg1 = va_arg (argptr, char*)) != NULL) {
1653 test = arg1;
1654 while ((arg[++nargs] = va_arg(ap, char*)) != NULL) {};
1655 }
1656 va_end (ap);
1657
1658 /* get a target - clause pair
1659 szMsg1 = va_arg(ap, PSTR);
1660 szMsg2 = va_arg(ap, PSTR); */
1661 va_end(ap);
1662
1663 } /* piecewise */
1664 #endif
1665
1666
1667 /* End of random module */
1668