#include #include #include /* for FLT_MAX and FLT_MIN*/ #include "gen_beta.h" #include "Utilities/Random.h" /* * gen_beta.cc generating beta-distributed random numbers * Copyright (C) 2000,2005 Kevin Karplus * * * This library is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public * License as published by the Free Software Foundation; * version 2.1 of the License. * * This library is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * Lesser General Public License for more details. * * See http://www.gnu.org/copyleft/lesser.html for the license details, * or write to the Free Software Foundation, Inc., 59 Temple Place, Suite * 330, Boston, MA 02111-1307 USA */ /* This code was initially written by Kevin Karplus 2-6 Nov 2000 * * Some inspiration was taken from genbet in ranlib.c, which was * written by Barry Brown and James Lovato. * Their code was full of completely unnecessary goto statements, * since it was generated by a Fortran-to-c conversion. * They also did not provide for generating interleaved streams from * multiple beta distributions, and one of their cases did not seem to * work correctly after cleaning up the code. * After failing to fix their code, I started over from scratch. * * The current version uses 3 different beta generators, based on the * parameter values. * This choice of generators was based in part on the information in * Dagpunar's book "Principles of Random Variate Generation", * and partly on Chapter 21 of the documentation for the Numerical * Analysis Group's software: routine nag_rand_beta. * www.nag.com/numeric/fbfn/fn/Ctwentyone/Ctwentyone_txt.html * * I could not get any of the implementations of Cheng's Algorithm BC to * produce the right moments, so I stopped using it, and used Atkinson's * switching method for 0.5<= max_ab < 1 instead. The bugs were * probably in my translation of the method, not in the method itself, * as I was working exclusively from secondary sources. * * The implementations finally chosen were selected more for * robustness than absolute efficiency: * 1 < min_ab , using Cheng's algortihm BB * min_ab <= 1 <= max_ab, using Atkinson's switching algorithm * max_ab < 0.5, using Joehnk's algorithm * 0.5<= max_ab < 1, using Atkinson's switching algorithm */ // Wed Jul 20 15:20:11 PDT 2005 Kevin Karplus // Changed from C to C++ for inclusion in ultimate library /* Initialize the variables of the already allocated generator in gen . * This must be called before gen can be used. */ void gen_beta::initialize(double aa, double bb) { double t; assert (aa>0 && bb >0); a = aa; b=bb; if (a< b) {min_ab=a; max_ab=b;} else {min_ab=b; max_ab=a;} sum_ab = a+b; if (max_ab < 0.5) { /* use Joehnk's algorithm, no precomputation */ /* BUG fix by Vladimir Valenta, vladimir.valenta@bankofamerica.com * 28 April 2003 (the test had been incorrectly <=0.5) */ } else if (min_ab>1.0) { /* use Cheng's Algorithm BB */ param[0] = sqrt((sum_ab-2.0)/(2.0*a*b-sum_ab)); param[1] = min_ab+1.0/param[0]; } else if (max_ab > 1.0) { /* use Atkinson's switching algorithm * p=min_ab, q=max_ab * t stored as param[0], r as param[1] */ t = (1-min_ab)/(1+max_ab - min_ab); param[0] = t; param[1] = max_ab * t / ( max_ab * t + min_ab*pow(1-param[0], max_ab)); } else { /* use Atkinson's switching algorithm 0.5 <= max_ab <= 1.0 * p=min_ab, q=max_ab * */ if (min_ab == 1.0) { param[0] = param[1] = 0.5; } else { t = 1/(1+sqrt(max_ab*(1-max_ab)/ (min_ab*(1-min_ab)))); param[0] = t ; param[1] = max_ab * t / (max_ab * t + min_ab * (1-t)); } } } /* define a safe version of a*exp(v), with "infinity" as FLT_MAX * (approx 3.40282347e+38f) */ #define aexp(a,v) (v> 88.7? FLT_MAX: a*exp(v)) /* define a return function, taking into account the pseudosymmetry of the * Beta distribution */ #define ret(a,mina, b,w) ((a==mina)? w/(b+w) : b/(b+w)) /* Generate one number from the distribution specified by gen. Returns a single random deviate from the beta distribution with parameters A and B. The density of the beta is x^(a-1) * (1-x)^(b-1) / B(a,b) for 0 < x < 1 Method R. C. H. Cheng Generating Beta Variates with Nonintegral Shape Parameters Communications of the ACM, 21:317-322 (1978) (Algorithms BB and BC) Atkinson, A. C. A family of switching algorithms for the computer generation of beta random variables. Biometrika 66:141-145. (1979) Joehnk, M.D. Erzeugung von Betaverteilten und Param[1]verteilten Zufallszahlen. Metrika, 8:5-15. (1964) Cited in Dagpunar, John. Principles of Random Variate Generation. Oxford University Press, 1988. */ double gen_beta::generate(void) const { /* The following temporaries are recomputed on each iteration, * or restored from param array. */ double c,r,s,t,u1,u2,v,w,z,lambda; double logv, logw, log_sum; if (max_ab<0.5) { /* Use Joehnk's algorithm. * Use logv and logw, rather than v and w, to avoid * floating-point underflow with very small a or b values. */ do { u1 = drandom(); u2 = drandom(); logv = log(u1)/a; logw = log(u2)/b; log_sum = logv>logw? logv + log(1+ exp(logw-logv)) : logw + log(1+ exp(logv-logw)); } while (log_sum>0.0); assert(logv<=log_sum); return exp(logv - log_sum); } if (min_ab > 1.0) { /* use Algorithm BB */ lambda = param[0]; c = param[1]; do { u1 = drandom(); u2 = drandom(); v = lambda*log(u1/(1.0-u1)); w = aexp(min_ab,v); z = u1*u1*u2; r = c*v-1.38629436112; s = min_ab+r-w; if(s+2.609438 >= 5.0*z) break; t = log(z); } while ( /* s<=t && */ r+sum_ab*log(sum_ab/(max_ab+w)) < t); return ret(a,min_ab, max_ab, w); } if (max_ab>= 1.0) { /* use Atkinson's switching method, as * described in Dagpunar's book * p=min_ab, q=max_ab * t stored as param[0], r as param[1] */ t = param[0]; r = param[1]; for(;;) { u1 = drandom(); u2 = drandom(); if (u1 < r) { w = t*pow(u1/r, 1/min_ab); if (log(u2) < (max_ab -1)*log(1-w)) break; } else { w = 1- (1-t)*pow((1-u1)/(1-r), 1/max_ab); if (log(u2) < (min_ab -1) * log(w/t)) break; } } return (a==min_ab)? w : 1-w; } else { /* use Atkinson's Algorithm */ t = param[0]; r = param[1]; for(;;) { u1 = drandom(); u2 = drandom(); if (u1 < r) { w = t*pow(u1/r, 1/min_ab); if (log(u2) < (max_ab -1)*log((1-w)/(1-t))) break; } else { w = 1- (1-t)*pow((1-u1)/(1-r), 1/max_ab); if (log(u2) < (min_ab -1) * log(w/t)) break; } } return (a==min_ab)? w : 1-w; } }