ADMB Documentation  11.1.2247
 All Classes Files Functions Variables Typedefs Friends Defines
drangam.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: drangam.cpp 1919 2014-04-22 22:02:01Z johnoel $
00003  *
00004  * Authors: H. H. Ahrens, U. Dieter and B. W. Brown.
00005  * License: Public domain.
00006  */
00011 #include <fvar.hpp>
00012 
00013 /*
00014 void main()
00015 {
00016   int n=1000;
00017   dvector v(1,n);
00018   random_number_generator rng(1011);
00019 
00020   v(1)=sgamma(18.0,rng);
00021   for (int i=1;i<=n;i++)
00022   {
00023     //v(i)=sgamma(double(i),rng);
00024     v(i)=sgamma(.5,rng);
00025   }
00026   cout << mean(v) << endl;
00027 }
00028 */
00029 
00033 static double fsign( double num, double sign )
00034 {
00035 if ( ( sign>0.0f && num<0.0f ) || ( sign<0.0f && num>0.0f ) )
00036     return -num;
00037 else return num;
00038 }
00039 
00044 double sgamma(double a,const random_number_generator& _rng)
00045 /*
00046 **********************************************************************
00047 
00048 
00049      (STANDARD-)  G A M M A  DISTRIBUTION
00050 
00051 
00052 **********************************************************************
00053 **********************************************************************
00054 
00055                PARAMETER  A >= 1.0  !
00056 
00057 **********************************************************************
00058 
00059      FOR DETAILS SEE:
00060 
00061                AHRENS, J.H. AND DIETER, U.
00062                GENERATING GAMMA VARIATES BY A
00063                MODIFIED REJECTION TECHNIQUE.
00064                COMM. ACM, 25,1 (JAN. 1982), 47 - 54.
00065 
00066      STEP NUMBERS CORRESPOND TO ALGORITHM 'GD' IN THE ABOVE PAPER
00067                                  (STRAIGHTFORWARD IMPLEMENTATION)
00068 
00069      Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of
00070      SUNIF.  The argument IR thus goes away.
00071 
00072 **********************************************************************
00073 
00074                PARAMETER  0.0 < A < 1.0  !
00075 
00076 **********************************************************************
00077 
00078      FOR DETAILS SEE:
00079 
00080                AHRENS, J.H. AND DIETER, U.
00081                COMPUTER METHODS FOR SAMPLING FROM GAMMA,
00082                BETA, POISSON AND BINOMIAL DISTRIBUTIONS.
00083                COMPUTING, 12 (1974), 223 - 246.
00084 
00085      (ADAPTED IMPLEMENTATION OF ALGORITHM 'GS' IN THE ABOVE PAPER)
00086 
00087 **********************************************************************
00088      INPUT: A =PARAMETER (MEAN) OF THE STANDARD GAMMA DISTRIBUTION
00089      OUTPUT: SGAMMA = SAMPLE FROM THE GAMMA-(A)-DISTRIBUTION
00090      COEFFICIENTS Q(K) - FOR Q0 = SUM(Q(K)*A**(-K))
00091      COEFFICIENTS A(K) - FOR Q = Q0+(T*T/2)*SUM(A(K)*V**K)
00092      COEFFICIENTS E(K) - FOR EXP(Q)-1 = SUM(E(K)*Q**K)
00093      PREVIOUS A PRE-SET TO ZERO - AA IS A', AAA IS A"
00094      SQRT32 IS THE SQUAREROOT OF 32 = 5.656854249492380
00095 */
00096 {
00097   random_number_generator& rng=(random_number_generator&) _rng;
00098 #ifndef __OPENCC__
00099 extern double fsign( double num, double sign );
00100 #endif
00101 static double q1 = 4.166669E-2;
00102 static double q2 = 2.083148E-2;
00103 static double q3 = 8.01191E-3;
00104 static double q4 = 1.44121E-3;
00105 static double q5 = -7.388E-5;
00106 static double q6 = 2.4511E-4;
00107 static double q7 = 2.424E-4;
00108 static double a1 = 0.3333333;
00109 static double a2 = -0.250003;
00110 static double a3 = 0.2000062;
00111 static double a4 = -0.1662921;
00112 static double a5 = 0.1423657;
00113 static double a6 = -0.1367177;
00114 static double a7 = 0.1233795;
00115 static double e1 = 1.0;
00116 static double e2 = 0.4999897;
00117 static double e3 = 0.166829;
00118 static double e4 = 4.07753E-2;
00119 static double e5 = 1.0293E-2;
00120 static double aa = 0.0;
00121 static double aaa = 0.0;
00122 static double sqrt32 = 5.656854;
00123 /* JJV added b0 to fix rare and subtle bug */
00124 static double sgamma,s2,s,d,t,x,u,r,q0,b,b0,si,c,v,q,e,w,p;
00125     if(a == aa) goto S10;
00126     if(a < 1.0) goto S120;
00127 /*
00128      STEP  1:  RECALCULATIONS OF S2,S,D IF A HAS CHANGED
00129 */
00130     aa = a;
00131     s2 = a-0.5;
00132     s = sqrt(s2);
00133     d = sqrt32-12.0*s;
00134 S10:
00135 /*
00136      STEP  2:  T=STANDARD NORMAL DEVIATE,
00137                X=(S,1/2)-NORMAL DEVIATE.
00138                IMMEDIATE ACCEPTANCE (I)
00139 */
00140     t = gasdev(rng);
00141     x = s+0.5*t;
00142     sgamma = x*x;
00143     if(t >= 0.0)
00144       return sgamma;
00145 /*
00146      STEP  3:  U= 0,1 -UNIFORM SAMPLE. SQUEEZE ACCEPTANCE (S)
00147 */
00148     u = rng.better_rand();
00149     //u=ranf();
00150     if(d*u <= t*t*t) return sgamma;
00151 /*
00152      STEP  4:  RECALCULATIONS OF Q0,B,SI,C IF NECESSARY
00153 */
00154     if(a == aaa) goto S40;
00155     aaa = a;
00156     r = 1.0/ a;
00157     q0 = ((((((q7*r+q6)*r+q5)*r+q4)*r+q3)*r+q2)*r+q1)*r;
00158 /*
00159                APPROXIMATION DEPENDING ON SIZE OF PARAMETER A
00160                THE CONSTANTS IN THE EXPRESSIONS FOR B, SI AND
00161                C WERE ESTABLISHED BY NUMERICAL EXPERIMENTS
00162 */
00163     if(a <= 3.686) goto S30;
00164     if(a <= 13.022) goto S20;
00165 /*
00166                CASE 3:  A .GT. 13.022
00167 */
00168     b = 1.77;
00169     si = 0.75;
00170     c = 0.1515/s;
00171     goto S40;
00172 S20:
00173 /*
00174                CASE 2:  3.686 .LT. A .LE. 13.022
00175 */
00176     b = 1.654+7.6E-3*s2;
00177     si = 1.68/s+0.275;
00178     c = 6.2E-2/s+2.4E-2;
00179     goto S40;
00180 S30:
00181 /*
00182                CASE 1:  A .LE. 3.686
00183 */
00184     b = 0.463+s+0.178*s2;
00185     si = 1.235;
00186     c = 0.195/s-7.9E-2+1.6E-1*s;
00187 S40:
00188 /*
00189      STEP  5:  NO QUOTIENT TEST IF X NOT POSITIVE
00190 */
00191     if(x <= 0.0) goto S70;
00192 /*
00193      STEP  6:  CALCULATION OF V AND QUOTIENT Q
00194 */
00195     v = t/(s+s);
00196     if(fabs(v) <= 0.25) goto S50;
00197     q = q0-s*t+0.25*t*t+(s2+s2)*log(1.0+v);
00198     goto S60;
00199 S50:
00200     q = q0+0.5*t*t*((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v;
00201 S60:
00202 /*
00203      STEP  7:  QUOTIENT ACCEPTANCE (Q)
00204 */
00205     if(log(1.0-u) <= q) return sgamma;
00206 S70:
00207 /*
00208      STEP  8:  E=STANDARD EXPONENTIAL DEVIATE
00209                U= 0,1 -UNIFORM DEVIATE
00210                T=(B,SI)-DOUBLE EXPONENTIAL (LAPLACE) SAMPLE
00211 */
00212     e = expdev(rng);
00213     u = rng.better_rand();
00214     // u = ranf();
00215     u += (u-1.0);
00216     t = b+fsign(si*e,u);
00217 /*
00218      STEP  9:  REJECTION IF T .LT. TAU(1) = -.71874483771719
00219 */
00220     if(t < -0.7187449) goto S70;
00221 /*
00222      STEP 10:  CALCULATION OF V AND QUOTIENT Q
00223 */
00224     v = t/(s+s);
00225     if(fabs(v) <= 0.25) goto S80;
00226     q = q0-s*t+0.25*t*t+(s2+s2)*log(1.0+v);
00227     goto S90;
00228 S80:
00229     q = q0+0.5*t*t*((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v;
00230 S90:
00231 /*
00232      STEP 11:  HAT ACCEPTANCE (H) (IF Q NOT POSITIVE GO TO STEP 8)
00233 */
00234     if(q <= 0.0) goto S70;
00235     if(q <= 0.5) goto S100;
00236 /*
00237  * JJV modified the code through line 115 to handle large Q case
00238  */
00239     if(q < 15.0) goto S95;
00240 /*
00241  * JJV Here Q is large enough that Q = log(exp(Q) - 1.0) (for real Q)
00242  * JJV so reformulate test at 110 in terms of one EXP, if not too big
00243  * JJV 87.49823 is close to the largest real which can be
00244  * JJV exponentiated (87.49823 = log(1.0E38))
00245  */
00246     if((q+e-0.5*t*t) > 87.49823) goto S115;
00247     if(c*fabs(u) > exp(q+e-0.5*t*t)) goto S70;
00248     goto S115;
00249 S95:
00250     w = exp(q)-1.0;
00251     goto S110;
00252 S100:
00253     w = ((((e5*q+e4)*q+e3)*q+e2)*q+e1)*q;
00254 S110:
00255 /*
00256                IF T IS REJECTED, SAMPLE AGAIN AT STEP 8
00257 */
00258     if(c*fabs(u) > w*exp(e-0.5*t*t)) goto S70;
00259 S115:
00260     x = s+0.5*t;
00261     sgamma = x*x;
00262     return sgamma;
00263 S120:
00264 /*
00265      ALTERNATE METHOD FOR PARAMETERS A BELOW 1  (.3678794=EXP(-1.))
00266 
00267      JJV changed B to B0 (which was added to declarations for this)
00268      JJV in 120 to END to fix rare and subtle bug.
00269      JJV Line: 'aa = 0.0' was removed (unnecessary, wasteful).
00270      JJV Reasons: the state of AA only serves to tell the A >= 1.0
00271      JJV case if certain A-dependent constants need to be recalculated.
00272      JJV The A < 1.0 case (here) no longer changes any of these, and
00273      JJV the recalculation of B (which used to change with an
00274      JJV A < 1.0 call) is governed by the state of AAA anyway.
00275     aa = 0.0;
00276 */
00277     b0 = 1.0+0.3678794*a;
00278 S130:
00279     p = b0*rng.better_rand();
00280     // p = b0*ranf();
00281     if(p >= 1.0) goto S140;
00282     sgamma = exp(log(p)/ a);
00283     if(expdev(rng) < sgamma) goto S130;
00284     return sgamma;
00285 S140:
00286     sgamma = -log((b0-p)/ a);
00287     if(expdev(rng) < (1.0-a)*log(sgamma)) goto S130;
00288     return sgamma;
00289 }
00290 
00295 double gasdev(const random_number_generator& _rng)
00296 {
00297   random_number_generator& rng=(random_number_generator&) _rng;
00298   double fac,rsq,v1,v2;
00299 
00300   do
00301   {
00302     v1=2.0* rng.better_rand() -1.0;
00303     v2=2.0*rng.better_rand()-1.0;
00304     rsq=v1*v1+v2*v2;
00305   } while (rsq >= 1.0 || rsq == 0.0);
00306   fac=sqrt(-2.0*log(rsq)/rsq);
00307   return v1*fac;
00308 }
00309 
00314 double expdev(const random_number_generator& _rng)
00315 {
00316   random_number_generator& rng=(random_number_generator&) _rng;
00317   double dum;
00318 
00319   do
00320     dum=rng.better_rand();
00321   while (dum == 0.0);
00322   return -log(dum);
00323 }