From b64f9fe088c534a4b405921fafb90d248384ae6f Mon Sep 17 00:00:00 2001
From: Jacob Smith <jake@laptop.(none)>
Date: Mon, 19 Apr 2010 04:03:25 -0400
Subject: [PATCH] tao updates
---
CFiles/rand.c | 85 +++++++++++++++++++++++++++++++++++----------------
1 file changed, 59 insertions(+), 26 deletions(-)
diff --git a/CFiles/rand.c b/CFiles/rand.c
index 1c7306a..03d08db 100644
--- a/CFiles/rand.c
+++ b/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
+
+ //@ When used entirely with C.
+ getchar();
+ exit(1); // This exits the entire C program.
}
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
}
--
GitLab