43 #include "../NumConstants.h"
44 #include "../VectorTools.h"
58 DEFAULT_GENERATOR->setSeed(seed);
83 assert(dblAlpha > 0.0);
86 else if (dblAlpha > 1.0)
107 vector<size_t> sample(n);
108 for (
unsigned int i = 0; i < n; i++)
113 for (
unsigned int j = 0; test& (j < probs.size()); j++)
115 cumprob += probs[j] / s;
124 sample[i] = probs.size();
145 rgdbl[1] = dblAlpha - 1.0;
146 rgdbl[2] = (dblAlpha - (1.0 / (6.0 * dblAlpha))) / rgdbl[1];
147 rgdbl[3] = 2.0 / rgdbl[1];
148 rgdbl[4] = rgdbl[3] + 2.0;
149 rgdbl[5] = 1.0 / sqrt(dblAlpha);
160 dblRand1 = dblRand2 + rgdbl[5] * (1.0 - 1.86 * dblRand1);
162 while (!(0.0 < dblRand1 && dblRand1 < 1.0));
164 double dblTemp = rgdbl[2] * dblRand2 / dblRand1;
166 if (rgdbl[3] * dblRand1 + dblTemp + 1.0 / dblTemp <= rgdbl[4] ||
167 rgdbl[3] * log(dblRand1) + dblTemp - log(dblTemp) < 1.0)
169 return dblTemp * rgdbl[1];
182 const double dblexp = exp(1.0);
185 double dblRand0 = giveRandomNumberBetweenZeroAndEntry(1.0, generator);
186 double dblRand1 = giveRandomNumberBetweenZeroAndEntry(1.0, generator);
187 if (dblRand0 <= (dblexp / (dblAlpha + dblexp)))
189 dblTemp = pow(((dblAlpha + dblexp) * dblRand0) /
190 dblexp, 1.0 / dblAlpha);
191 if (dblRand1 <= exp(-1.0 * dblTemp))
196 dblTemp = -1.0 * log((dblAlpha + dblexp) * (1.0 - dblRand0) / (dblAlpha * dblexp));
197 if (dblRand1 <= pow(dblTemp, dblAlpha - 1.0))
213 double a0 = -.322232431088, a1 = -1, a2 = -.342242088547, a3 = -.0204231210245;
214 double a4 = -.453642210148e-4, b0 = .0993484626060, b1 = .588581570495;
215 double b2 = .531103462366, b3 = .103537752850, b4 = .0038560700634;
216 double y, z = 0, p = prob, p1;
218 p1 = (p < 0.5 ? p : 1 - p);
222 y = sqrt (log(1 / (p1 * p1)));
223 z = y + ((((y * a4 + a3) * y + a2) * y + a1) * y + a0) / ((((y * b4 + b3) * y + b2) * y + b1) * y + b0);
224 return p < 0.5 ? -z : z;
235 double p = alpha, g = ln_gamma_alpha;
236 double accurate = 1e-8, overflow = 1e30;
237 double factor, gin = 0, rn = 0, a = 0, b = 0, an = 0, dif = 0, term = 0;
238 vector<double> pn(6);
245 factor = exp(p * log(x) - x - g);
249 gin = 1; term = 1; rn = p;
252 term *= x / rn; gin += term;
260 a = 1 - p; b = a + x + 1; term = 0;
261 pn[0] = 1; pn[1] = x; pn[2] = x + 1; pn[3] = x * b;
264 a++; b += 2; term++; an = a * term;
265 for (i = 0; i < 2; i++)
267 pn[i + 4] = b * pn[i + 2] - an * pn[i];
271 rn = pn[4] / pn[5]; dif = fabs(gin - rn);
274 if (dif <= accurate * rn)
279 for (i = 0; i < 4; i++)
283 if (fabs(pn[4]) < overflow)
285 for (i = 0; i < 4; i++)
291 gin = 1 - factor * gin;
301 double e = .5e-6, aa = .6931471805, p = prob, g;
302 double xx, c, ch, a = 0, q = 0, p1 = 0, p2 = 0, t = 0, x = 0, b = 0, s1, s2, s3, s4, s5, s6;
304 if (p < .000002 || p > .999998 || v <= 0)
308 xx = v / 2; c = xx - 1;
309 if (v >= -1.24 * log(p))
312 ch = pow((p * xx * exp(g + xx * aa)), 1 / xx);
319 ch = 0.4; a = log(1 - p);
321 q = ch; p1 = 1 + ch * (4.67 + ch); p2 = ch * (6.73 + ch * (6.66 + ch));
322 t = -0.5 + (4.67 + 2 * ch) / p1 - (6.73 + ch * (13.32 + 3 * ch)) / p2;
323 ch -= (1 - exp(a + g + .5 * ch + c * aa) * p2 / p1) / t;
324 if (fabs(q / ch - 1) - .01 <= 0)
331 p1 = 0.222222 / v; ch = v * pow((x * sqrt(p1) + 1 - p1), 3.0);
332 if (ch > 2.2 * v + 6)
333 ch = -2 * (log(1 - p) - c * log(.5 * ch) + g);
335 q = ch; p1 = .5 * ch;
336 if ((t = incompleteGamma (p1, xx, g)) < 0)
338 std::cerr <<
"err IncompleteGamma" << std::endl;
342 t = p2 * exp(xx * aa + g + p1 - c * log(ch));
343 b = t / ch; a = 0.5 * t - b * c;
345 s1 = (210 + a * (140 + a * (105 + a * (84 + a * (70 + 60 * a))))) / 420;
346 s2 = (420 + a * (735 + a * (966 + a * (1141 + 1278 * a)))) / 2520;
347 s3 = (210 + a * (462 + a * (707 + 932 * a))) / 2520;
348 s4 = (252 + a * (672 + 1182 * a) + c * (294 + a * (889 + 1740 * a))) / 5040;
349 s5 = (84 + 264 * a + c * (175 + 606 * a)) / 2520;
350 s6 = (120 + c * (346 + 127 * c)) / 5040;
351 ch += t * (1 + 0.5 * t * s1 - b * c * (s1 - b * (s2 - b * (s3 - b * (s4 - b * (s5 - b * s6))))));
352 if (fabs(q / ch - 1) > e)
367 const static double a[5] = {
368 2.2352520354606839287,
369 161.02823106855587881,
370 1067.6894854603709582,
371 18154.981253343561249,
372 0.065682337918207449113
374 const static double b[4] = {
375 47.20258190468824187,
376 976.09855173777669322,
377 10260.932208618978205,
378 45507.789335026729956
380 const static double c[9] = {
381 0.39894151208813466764,
382 8.8831497943883759412,
383 93.506656132177855979,
384 597.27027639480026226,
385 2494.5375852903726711,
386 6848.1904505362823326,
387 11602.651437647350124,
388 9842.7148383839780218,
389 1.0765576773720192317e-8
391 const static double d[8] = {
392 22.266688044328115691,
393 235.38790178262499861,
394 1519.377599407554805,
395 6485.558298266760755,
396 18615.571640885098091,
397 34900.952721145977266,
398 38912.003286093271411,
399 19685.429676859990727
401 const static double p[6] = {
403 0.1274011611602473639,
404 0.022235277870649807,
405 0.001421619193227893466,
406 2.9112874951168792e-5,
407 0.02307344176494017303
409 const static double q[5] = {
411 0.468238212480865118,
412 0.0659881378689285515,
413 0.00378239633202758244,
414 7.29751555083966205e-5
417 double xden, xnum, temp, del, eps, xsq, y, cum;
430 for (i = 0; i < 3; ++i)
432 xnum = (xnum + a[i]) * xsq;
433 xden = (xden + b[i]) * xsq;
439 temp = x * (xnum + a[3]) / (xden + b[3]);
442 else if (y <= sqrt(32))
448 for (i = 0; i < 7; ++i)
450 xnum = (xnum + c[i]) * y;
451 xden = (xden + d[i]) * y;
453 temp = (xnum + c[7]) / (xden + d[7]);
455 xsq = trunc(y * 16) / 16;
456 del = (y - xsq) * (y + xsq);
457 cum = exp(-xsq * xsq * 0.5) * exp(-del * 0.5) * temp;
462 else if (-37.5193 < x && x < 8.2924)
467 for (i = 0; i < 4; ++i)
469 xnum = (xnum + p[i]) * xsq;
470 xden = (xden + q[i]) * xsq;
472 temp = xsq * (xnum + p[4]) / (xden + q[4]);
473 temp = (1 / sqrt(2 * M_PI) - temp) / y;
475 xsq = trunc(x * 16) / 16;
476 del = (x - xsq) * (x + xsq);
478 cum = exp(-xsq * xsq * 0.5) * exp(-del * 0.5) * temp;
496 return lnGamma(alpha) + lnGamma(beta) - lnGamma(alpha + beta);
509 double const1 = 2.30753;
510 double const2 = 0.27061;
511 double const3 = 0.99229;
512 double const4 = 0.04481;
515 int swap_tail, i_pb, i_inn;
516 double a, adj, logbeta, g, h, pp, prev, qq, r, s, t, tx, w, y, yprev;
518 volatile double xinbta;
520 if (alpha <= 0. || beta < 0.)
521 throw (
"RandomTools::qBeta with non positive parameters");
523 if (prob < 0. || prob > 1.)
524 throw (
"RandomTools::qBeta with bad probability");
527 logbeta = lnBeta(alpha, beta);
532 a = prob; pp = alpha; qq = beta; swap_tail = 0;
537 pp = beta; qq = alpha; swap_tail = 1;
543 r = sqrt(-2 * log(a));
544 y = r - (const1 + const2 * r) / (1. + (const3 + const4 * r) * r);
545 if (pp > 1 && qq > 1)
547 r = (y * y - 3.) / 6.;
548 s = 1. / (pp + pp - 1.);
549 t = 1. / (qq + qq - 1.);
551 w = y * sqrt(h + r) / h - (t - s) * (r + 5. / 6. - 2. / (3. * h));
552 xinbta = pp / (pp + qq * exp(w + w));
558 t = r * pow(1. - t + y * sqrt(t), 3.0);
560 xinbta = 1. - exp((log1p(-a) + log(qq) + logbeta) / qq);
563 t = (4. * pp + r - 2.) / t;
565 xinbta = exp((log(a * pp) + logbeta) / pp);
567 xinbta = 1. - 2. / (t + 1.);
581 else if (xinbta > upper)
594 double po = pow(10., -13 - 2.5 / (pp * pp) - 0.5 / (a * a));
595 acu = (lower > po) ? lower : po;
599 for (i_pb = 0; i_pb < 1000; i_pb++)
601 y = incompleteBeta(xinbta, pp, qq);
610 exp(logbeta + r * log(xinbta) + t * log1p(-xinbta));
612 prev = (fabs(adj) > lower) ? fabs(adj) : lower;
614 for (i_inn = 0; i_inn < 1000; i_inn++)
617 if (fabs(adj) < prev)
620 if (tx >= 0. && tx <= 1)
622 if ((prev <= acu) || (fabs(y) <= acu))
623 return swap_tail ? 1 - xinbta : xinbta;
624 if (tx != 0. && tx != 1)
630 if (fabs(tx - xinbta) < 1e-15 * xinbta)
631 return swap_tail ? 1 - xinbta : xinbta;
638 return swap_tail ? 1 - xinbta : xinbta;
654 big = 4.503599627370496e15;
655 biginv = 2.22044604925031308085e-16;
656 maxgam = 171.624376956302725;
660 if ((alpha <= 0) || (beta <= 0))
661 throw Exception(
"RandomTools::incompleteBeta not valid with non-positive parameters");
663 if ((x < 0) || (x > 1))
664 throw Exception(
"RandomTools::incompleteBeta out of bounds limit");
673 if ((beta * x <= 1.0) && (x <= 0.95))
675 return incompletebetaps(alpha, beta, x, maxgam);
679 if (x > alpha / (alpha + beta))
692 if (flag == 1 && (beta * x <= 1.0) && (x <= 0.95) )
694 t = incompletebetaps(alpha, beta, x, maxgam);
701 y = x * (alpha + beta - 2.0) - (alpha - 1.0);
704 w = incompletebetafe(alpha, beta, x, big, biginv);
708 w = incompletebetafe2(alpha, beta, x, big, biginv) / xc;
712 if ( (alpha + beta < maxgam) && (fabs(y) < maxlog) && (fabs(t) < maxlog) )
715 t = t * pow(x, alpha);
718 t = t * exp(lnGamma(alpha + beta) - (lnGamma(alpha) + lnGamma(beta)));
729 y = y + t + lnGamma(alpha + beta) - lnGamma(alpha) - lnGamma(beta);
730 y = y + log(w / alpha);
798 xk = -x * k1 * k2 / (k3 * k4);
799 pk = pkm1 + pkm2 * xk;
800 qk = qkm1 + qkm2 * xk;
805 xk = x * k5 * k6 / (k7 * k8);
806 pk = pkm1 + pkm2 * xk;
807 qk = qkm1 + qkm2 * xk;
818 t = fabs((ans - r) / r);
837 if (fabs(qk) + fabs(pk) > big)
839 pkm2 = pkm2 * biginv;
840 pkm1 = pkm1 * biginv;
841 qkm2 = qkm2 * biginv;
842 qkm1 = qkm1 * biginv;
844 if ((fabs(qk) < biginv) || (fabs(pk) < biginv))
914 xk = -z * k1 * k2 / (k3 * k4);
915 pk = pkm1 + pkm2 * xk;
916 qk = qkm1 + qkm2 * xk;
921 xk = z * k5 * k6 / (k7 * k8);
922 pk = pkm1 + pkm2 * xk;
923 qk = qkm1 + qkm2 * xk;
934 t = fabs((ans - r) / r);
953 if (fabs(qk) + fabs(pk) > big)
955 pkm2 = pkm2 * biginv;
956 pkm1 = pkm1 * biginv;
957 qkm2 = qkm2 * biginv;
958 qkm1 = qkm1 * biginv;
960 if ((fabs(qk) < biginv) || (fabs(pk) < biginv))
1004 u = (n - b) * x / n;
1015 t = exp(lnGamma(a + b) - (lnGamma(a) + lnGamma(b)));
1016 s = s * t * pow(x, a);
1020 t = lnGamma(a + b) - lnGamma(a) - lnGamma(b) + u + log(s);
Exception base class. Overload exception constructor (to control the exceptions mechanism)....
static double VERY_TINY()
This is the interface for the Random Number Generators.
virtual double drawNumber() const =0
Return a random number.