tao updates

CFiles/rand.c

@@ -143,12 +143,15 @@ double gaussrnd(void)
    When a<1.0, Johnk's generator (Devroye, page 418).
    When a>1.0, a rejection method or Best's algorithm (Devroye, page 410).
 
-   A general gamma variate can be obtained as follows.  Let z=b*x.  Then,
+   A general gamma variate can be obtained as follows.  Let z=x/b.  Then,
    z is drawn from a general gamma distribution whose density is
 
-                                      1
-                         p(z) = -------------- z^(a-1) exp(-z/b).
-                                 gamma(a) b^a
+                                   b^a
+                         p(z) = ---------- z^(a-1) exp(-b*z).
+                                 gamma(a)
+
+   This density is the same as Matlab where b is an inverse scale parameter (same
+   as Gelman's book but different from Matlab).
 
    Uses algorithm translated by Iskander Karibzhanov from the Matlab function
    gamrnd.m, which follows Johnk's generator in Devroye ("Non-Uniform Random
@@ -158,16 +161,14 @@ double gammrnd(double a) {
    double b = a-1.0,
           u, v, w, x, y, z;
 
-   if (a <= 0.0) {
+   if (a <= 0.0)
+   {
       //** When used with a C MEX-file for MATLAB.
       printf("\nThe input argument x for gammrnd(x) in rand.c must be a positive double.");
-      #ifdef WIN_MATLABAPI
-         mexErrMsgTxt("Error! ???.");       // This terminates the entire Matlab program.
-      #else
+
       //@ When used entirely with C.
       getchar();
       exit(1);                          // This exits the entire C program.
-      #endif
    }
    if (a==1.0)
       return -log(unirnd());
@@ -198,6 +199,24 @@ double gammrnd(double a) {
 }
 
 
+/*
+   Returns a random draw from beta(a,b).  The density function is
+
+                                     gamma(a+b)
+                         p(x) =  ------------------- x^(a-1) (1-x)^(b-1).
+                                  gamma(a)*gamma(b)
+
+   where a>0, b>0 and gamma denotes the gamma function, which is defined as
+   the integral with respect to t from 0 to infinity of exp(-t)*t^(z-1).
+*/
+double betarnd(double a, double b)
+{
+   double xa=2.0*gammrnd(a),
+          xb=2.0*gammrnd(b);
+
+   return (xa/(xa+xb));
+}
+
 
 /*
    Returns the integral from -infinity to x of 1/sqrt(2*PI)*exp(-y^2/2).
@@ -276,12 +295,19 @@ double StandardNormalDouble(void) {
       double x;
       imsls_d_random_normal(1, IMSLS_RETURN_USER, &x, 0);
       return x;
-   #else 
-      // CASE2_RANDOMNUMBERGENERATOR Imported from the C recipe book -- Case 2 and my own (Iskander) code for generating a gamma distribution.
+   #else CASE2_RANDOMNUMBERGENERATOR   //Imported from the C recipe book -- Case 2 and my own (Iskander) code for generating a gamma distribution.
       return gaussrnd();
    #endif
 }
 
+double NormalDouble(double m, double s)
+{
+   //m: mean;
+   //s: standard deviation.
+
+   return (m+s*StandardNormalDouble());
+}
+
 double LogNormalDouble(double mean, double std)
 {
    //Return x where log(x) is normal(mean, std^2).
@@ -290,8 +316,7 @@ double LogNormalDouble(double mean, double std)
       double x;
       imsls_d_random_lognormal(1, mean, std, IMSLS_RETURN_USER, &x, 0);
       return x;
-   #else 
-      // CASE2_RANDOMNUMBERGENERATOR Imported from the C recipe book -- Case 2 and my own (Iskander) code for generating a gamma distribution.
+   #else CASE2_RANDOMNUMBERGENERATOR   //Imported from the C recipe book -- Case 2 and my own (Iskander) code for generating a gamma distribution.
       fn_DisplayError(".../rand.c/LogNormalDouble(): NO defaul available for log normal draws");
       return (0);
    #endif
@@ -299,20 +324,30 @@ double LogNormalDouble(double mean, double std)
 
 
 double GammaDouble(const double a, const double b) {
-   //The probability density of x is of the form: ( x^(a-1) exp(-x/b) )/( Gamma(a) * b^a ).
-   //The same form as in the MATLAB (not Gelman et al's) notation.
+   //The probability density of x is of the form: ( x^(a-1) exp(-bx) )/( b^a /Gamma(a) ).
+   //The same form as in the Gelman et al's (not MATLAB) notation.
    #ifdef IMSL_RANDOMNUMBERGENERATOR             // IMSL optimization routines.
       //The GFSR algorithm to generate random numbers.
       double x;
       imsls_d_random_gamma(1, a, IMSLS_RETURN_USER, &x, 0);
-      if (b != 1.0)  x *= b;
+      if (b != 1.0)  x /= b;
       return x;
-   #else 
-      // CASE2_RANDOMNUMBERGENERATOR Imported from the C recipe book -- Case 2 and my own (Iskander) code for generating a gamma distribution.
-      return ( (b==1.0) ? gammrnd(a) : b*gammrnd(a) );
+   #else CASE2_RANDOMNUMBERGENERATOR   //Imported from the C recipe book -- Case 2 and my own (Iskander) code for generating a gamma distribution.
+      return ( (b==1.0) ? gammrnd(a) : gammrnd(a)/b );
    #endif
 }
 
+double InverseGammaDouble(const double a, const double b)
+{
+   //p(x) = ( b^a/Gamma(a) ) x^(-a-1) exp(-b/x) for a>0 and b>0.
+   //   where a is shape and b is scale parameter.
+   //E(x) = b/(a-1) for a>1;  var(x) = b^2/( (a-1)^2*(a-2) ) for a>2;
+   //Noninformative distribution: a,b -> 0.
+   //How to draw: (1) draw z from Gamma(a,b); (2) let x=1/z.
+
+   return (1.0/GammaDouble(a,b));
+}
+
 
 void UniformVector(TSdvector *x_dv)
 {
@@ -419,8 +454,7 @@ void GammaMatrix(TSdmatrix *X_dm, TSdmatrix *A_dm, TSdmatrix *B_dm) {
       }
 
       X_dm->flag = M_GE;
-   #else 
-      //CASE2_RANDOMNUMBERGENERATOR Imported from the C recipe book -- Case 2 and my own (Iskander) code for generating a gamma distribution.
+   #else CASE2_RANDOMNUMBERGENERATOR   //Imported from the C recipe book -- Case 2 and my own (Iskander) code for generating a gamma distribution.
       int _i, nels, nrows, ncols;
       double *X, *A, *B;
 
@@ -456,8 +490,7 @@ double ChisquareDouble(const double v) {
       double x;
       imsls_d_random_chi_squared(1, v, IMSLS_RETURN_USER, &x, 0);
       return x;
-   #else 
-      // CASE2_RANDOMNUMBERGENERATOR Imported from the C recipe book -- Case 2 and my own (Iskander) code for generating a gamma distribution.
+   #else CASE2_RANDOMNUMBERGENERATOR   //Imported from the C recipe book -- Case 2 and my own (Iskander) code for generating a gamma distribution.
       return ( 2.0*gammrnd(0.5*v) );
    #endif
 }