ADMB Documentation  11.1x.2711
 All Classes Files Functions Variables Typedefs Friends Defines
drangam.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: drangam.cpp 2259 2014-08-27 23:16:51Z 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 static double q1 = 4.166669E-2;
00099 static double q2 = 2.083148E-2;
00100 static double q3 = 8.01191E-3;
00101 static double q4 = 1.44121E-3;
00102 static double q5 = -7.388E-5;
00103 static double q6 = 2.4511E-4;
00104 static double q7 = 2.424E-4;
00105 static double a1 = 0.3333333;
00106 static double a2 = -0.250003;
00107 static double a3 = 0.2000062;
00108 static double a4 = -0.1662921;
00109 static double a5 = 0.1423657;
00110 static double a6 = -0.1367177;
00111 static double a7 = 0.1233795;
00112 static double e1 = 1.0;
00113 static double e2 = 0.4999897;
00114 static double e3 = 0.166829;
00115 static double e4 = 4.07753E-2;
00116 static double e5 = 1.0293E-2;
00117 static double aa = 0.0;
00118 static double aaa = 0.0;
00119 static double sqrt32 = 5.656854;
00120 /* JJV added b0 to fix rare and subtle bug */
00121 static double sgamma,s2,s,d,t,x,u,r,q0,b,b0,si,c,v,q,e,w,p;
00122     if(a == aa) goto S10;
00123     if(a < 1.0) goto S120;
00124 /*
00125      STEP  1:  RECALCULATIONS OF S2,S,D IF A HAS CHANGED
00126 */
00127     aa = a;
00128     s2 = a-0.5;
00129     s = sqrt(s2);
00130     d = sqrt32-12.0*s;
00131 S10:
00132 /*
00133      STEP  2:  T=STANDARD NORMAL DEVIATE,
00134                X=(S,1/2)-NORMAL DEVIATE.
00135                IMMEDIATE ACCEPTANCE (I)
00136 */
00137     t = gasdev(rng);
00138     x = s+0.5*t;
00139     sgamma = x*x;
00140     if(t >= 0.0)
00141       return sgamma;
00142 /*
00143      STEP  3:  U= 0,1 -UNIFORM SAMPLE. SQUEEZE ACCEPTANCE (S)
00144 */
00145     u = rng.better_rand();
00146     //u=ranf();
00147     if(d*u <= t*t*t) return sgamma;
00148 /*
00149      STEP  4:  RECALCULATIONS OF Q0,B,SI,C IF NECESSARY
00150 */
00151     if(a == aaa) goto S40;
00152     aaa = a;
00153     r = 1.0/ a;
00154     q0 = ((((((q7*r+q6)*r+q5)*r+q4)*r+q3)*r+q2)*r+q1)*r;
00155 /*
00156                APPROXIMATION DEPENDING ON SIZE OF PARAMETER A
00157                THE CONSTANTS IN THE EXPRESSIONS FOR B, SI AND
00158                C WERE ESTABLISHED BY NUMERICAL EXPERIMENTS
00159 */
00160     if(a <= 3.686) goto S30;
00161     if(a <= 13.022) goto S20;
00162 /*
00163                CASE 3:  A .GT. 13.022
00164 */
00165     b = 1.77;
00166     si = 0.75;
00167     c = 0.1515/s;
00168     goto S40;
00169 S20:
00170 /*
00171                CASE 2:  3.686 .LT. A .LE. 13.022
00172 */
00173     b = 1.654+7.6E-3*s2;
00174     si = 1.68/s+0.275;
00175     c = 6.2E-2/s+2.4E-2;
00176     goto S40;
00177 S30:
00178 /*
00179                CASE 1:  A .LE. 3.686
00180 */
00181     b = 0.463+s+0.178*s2;
00182     si = 1.235;
00183     c = 0.195/s-7.9E-2+1.6E-1*s;
00184 S40:
00185 /*
00186      STEP  5:  NO QUOTIENT TEST IF X NOT POSITIVE
00187 */
00188     if(x <= 0.0) goto S70;
00189 /*
00190      STEP  6:  CALCULATION OF V AND QUOTIENT Q
00191 */
00192     v = t/(s+s);
00193     if(fabs(v) <= 0.25) goto S50;
00194     q = q0-s*t+0.25*t*t+(s2+s2)*log(1.0+v);
00195     goto S60;
00196 S50:
00197     q = q0+0.5*t*t*((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v;
00198 S60:
00199 /*
00200      STEP  7:  QUOTIENT ACCEPTANCE (Q)
00201 */
00202     if(log(1.0-u) <= q) return sgamma;
00203 S70:
00204 /*
00205      STEP  8:  E=STANDARD EXPONENTIAL DEVIATE
00206                U= 0,1 -UNIFORM DEVIATE
00207                T=(B,SI)-DOUBLE EXPONENTIAL (LAPLACE) SAMPLE
00208 */
00209     e = expdev(rng);
00210     u = rng.better_rand();
00211     // u = ranf();
00212     u += (u-1.0);
00213     t = b+fsign(si*e,u);
00214 /*
00215      STEP  9:  REJECTION IF T .LT. TAU(1) = -.71874483771719
00216 */
00217     if(t < -0.7187449) goto S70;
00218 /*
00219      STEP 10:  CALCULATION OF V AND QUOTIENT Q
00220 */
00221     v = t/(s+s);
00222     if(fabs(v) <= 0.25) goto S80;
00223     q = q0-s*t+0.25*t*t+(s2+s2)*log(1.0+v);
00224     goto S90;
00225 S80:
00226     q = q0+0.5*t*t*((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v;
00227 S90:
00228 /*
00229      STEP 11:  HAT ACCEPTANCE (H) (IF Q NOT POSITIVE GO TO STEP 8)
00230 */
00231     if(q <= 0.0) goto S70;
00232     if(q <= 0.5) goto S100;
00233 /*
00234  * JJV modified the code through line 115 to handle large Q case
00235  */
00236     if(q < 15.0) goto S95;
00237 /*
00238  * JJV Here Q is large enough that Q = log(exp(Q) - 1.0) (for real Q)
00239  * JJV so reformulate test at 110 in terms of one EXP, if not too big
00240  * JJV 87.49823 is close to the largest real which can be
00241  * JJV exponentiated (87.49823 = log(1.0E38))
00242  */
00243     if((q+e-0.5*t*t) > 87.49823) goto S115;
00244     if(c*fabs(u) > exp(q+e-0.5*t*t)) goto S70;
00245     goto S115;
00246 S95:
00247     w = exp(q)-1.0;
00248     goto S110;
00249 S100:
00250     w = ((((e5*q+e4)*q+e3)*q+e2)*q+e1)*q;
00251 S110:
00252 /*
00253                IF T IS REJECTED, SAMPLE AGAIN AT STEP 8
00254 */
00255     if(c*fabs(u) > w*exp(e-0.5*t*t)) goto S70;
00256 S115:
00257     x = s+0.5*t;
00258     sgamma = x*x;
00259     return sgamma;
00260 S120:
00261 /*
00262      ALTERNATE METHOD FOR PARAMETERS A BELOW 1  (.3678794=EXP(-1.))
00263 
00264      JJV changed B to B0 (which was added to declarations for this)
00265      JJV in 120 to END to fix rare and subtle bug.
00266      JJV Line: 'aa = 0.0' was removed (unnecessary, wasteful).
00267      JJV Reasons: the state of AA only serves to tell the A >= 1.0
00268      JJV case if certain A-dependent constants need to be recalculated.
00269      JJV The A < 1.0 case (here) no longer changes any of these, and
00270      JJV the recalculation of B (which used to change with an
00271      JJV A < 1.0 call) is governed by the state of AAA anyway.
00272     aa = 0.0;
00273 */
00274     b0 = 1.0+0.3678794*a;
00275 S130:
00276     p = b0*rng.better_rand();
00277     // p = b0*ranf();
00278     if(p >= 1.0) goto S140;
00279     sgamma = exp(log(p)/ a);
00280     if(expdev(rng) < sgamma) goto S130;
00281     return sgamma;
00282 S140:
00283     sgamma = -log((b0-p)/ a);
00284     if(expdev(rng) < (1.0-a)*log(sgamma)) goto S130;
00285     return sgamma;
00286 }
00287 
00292 double gasdev(const random_number_generator& _rng)
00293 {
00294   random_number_generator& rng=(random_number_generator&) _rng;
00295   double fac,rsq,v1,v2;
00296 
00297   do
00298   {
00299     v1=2.0* rng.better_rand() -1.0;
00300     v2=2.0*rng.better_rand()-1.0;
00301     rsq=v1*v1+v2*v2;
00302   } while (rsq >= 1.0 || rsq == 0.0);
00303   fac=sqrt(-2.0*log(rsq)/rsq);
00304   return v1*fac;
00305 }
00306 
00311 double expdev(const random_number_generator& _rng)
00312 {
00313   random_number_generator& rng=(random_number_generator&) _rng;
00314   double dum;
00315 
00316   do
00317     dum=rng.better_rand();
00318   while (dum == 0.0);
00319   return -log(dum);
00320 }