#include #include #include #include /* for FLT_MAX and FLT_MIN*/ #include "gen_beta.h" /* * gen_beta.c generating beta-distributed random numbers * Copyright (C) 2000 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 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 */ /* pick which underlying uniform random number generator to use. */ #define DRAND() ((random()+0.5)/ 2147483648.0) /* Initialize the variables of the already allocated generator in gen . * This must be called before gen can be used. */ void gen_beta_initialize(gen_beta_param *gen, double a, double b) { double t; assert (a>0 && b >0); gen->a = a; gen->b=b; if (a< b) {gen->min_ab=a; gen->max_ab=b;} else {gen->min_ab=b; gen->max_ab=a;} gen->sum_ab = a+b; if (gen->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 (gen->min_ab>1.0) { /* use Cheng's Algorithm BB */ gen->param[0] = sqrt((gen->sum_ab-2.0)/(2.0*a*b-gen->sum_ab)); gen->param[1] = gen->min_ab+1.0/gen->param[0]; } else if (gen->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-gen->min_ab)/(1+gen->max_ab - gen->min_ab); gen->param[0] = t; gen->param[1] = gen->max_ab * t / ( gen->max_ab * t + gen->min_ab*pow(1-gen->param[0], gen->max_ab)); } else { /* use Atkinson's switching algorithm 0.5 <= max_ab <= 1.0 * p=min_ab, q=max_ab * */ if (gen->min_ab == 1.0) { gen->param[0] = gen->param[1] = 0.5; } else { t = 1/(1+sqrt(gen->max_ab*(1-gen->max_ab)/ (gen->min_ab*(1-gen->min_ab)))); gen->param[0] = t ; gen->param[1] = gen->max_ab * t / (gen->max_ab * t + gen->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(const gen_beta_param *gen) { /* The following temporaries are recomputed on each iteration, * or restored from gen->param array. */ double c,r,s,t,u1,u2,v,w,z,lambda; double logv, logw, log_sum; double a = gen->a; double b = gen->b; double min_ab = gen->min_ab; double max_ab = gen->max_ab; double sum_ab = gen->sum_ab; 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 = DRAND(); u2 = DRAND(); 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 = gen->param[0]; c = gen->param[1]; do { u1 = DRAND(); u2 = DRAND(); 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 gen->param[0], r as gen->param[1] */ t = gen->param[0]; r = gen->param[1]; for(;;) { u1 = DRAND(); u2 = DRAND(); 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 = gen->param[0]; r = gen->param[1]; for(;;) { u1 = DRAND(); u2 = DRAND(); 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; } }