pseudorandom.cpp

随机数类,产生随机数

		else
        {
			return 1 ;
        }
    }
}

/**** The Multivariate Gaussian Distribution ****/
/*Random vectors from a multivariate Normal distribution
 *CALL:  r = mnormrnd(n,mu,S)
 * r = matrix of random numbers from the multivariate normal
 *     distribution with mean mu and covariance matrix S.
 * n = number of sample vectors
 * method = 'svd'  Singular value decomp.  (stable, quite fast)
 * S must be a symmetric, semi-positive definite matrix with size equal
 * to the length of mu.
 * Example
 * mu = [0 5]; S = [1 0.45; 0.45 0.25];
 * r = mnormrnd(100,mu,S)
 */
Matrix<double> PseudoRandom::mnormrnd(int cases, Matrix<double> & mu, Matrix<double> & sigma)
{
	int m;
	m = sigma.RowNum();
	
    Matrix<double> w(m,1),v(m,m),U(m,m),T(m,m),r(m,cases),randn(m,cases),mu_matrix(m,cases);
	// Cholesky decomposition
	U = sigma.Chold();
    T = U.Tran();
	
	for (int i=1;i<=m;i++){
		for (int j=1;j<=cases;j++){
			randn.SetValue(i,j,gaussian(0.0,1.0));
			mu_matrix.SetValue(i,j,mu(i,1));
		}
	}
	//cout <<mu_matrix<<endl;
	//cout <<prod(T,randn)<<endl;
	r = T*randn + mu_matrix;
	//cout<<r<<endl;
	return r;
}
/**** The Truncated Multivariate Gaussian Distribution ****/
/*Procedure for drawing from truncated multivariate normal based on
Geweke's code. i.e draws from
     xdraw is N(amu,sigma) subject to a < d*x < b
     where N(.,) in the n-variate Normal, a and b are nx1
Note that d is restricted to being a nonsingular nxn matrix
la and lb are nx1 vectors set to one if no upper/lower bounds
kstep=order of Gibbs within the constraint rows
*/
Matrix<double> PseudoRandom::tnorm_rnd(int n, Matrix<double> & amu, Matrix<double> & sigma, Matrix<double> & a, 
	                         Matrix<double> & b,Matrix<double> & la, Matrix<double> & lb, Matrix<double> & d,Matrix<int> & kstep)
{	
	int niter=10;
	double t1,aa;
	int m = amu.RowNum();
	//transform to work in terms of z=d*x
	Matrix<double> z(m,1), anu(m,1), a1(m,1), b1(m,1), h(m,1), xdraw(m,1);
	Matrix<int> indx(n,1);
	Matrix<double> d_copy(m,m), dinv(m,m), tau(m,m), tau_(m,m), tau_copy(m,m), tauinv(m,m), c(m,m), xdraw_matrix(m,n);
	for (int col=1;col<=n;col++){
		
		d_copy = d;
		dinv = d_copy.Inv();
		anu = d * amu;
		
		tau_ = d*sigma;
		tau = tau_*(d.Tran());
	
		tau_copy = tau;
		tauinv = tau_copy.Inv();
	
		a1 = a - anu;
		b1 = b - anu;
		for (int i=1; i<=m; i++)
		{
			aa = tauinv(i,i);
			h.SetValue(i,1,1/sqrt(aa));
			for (int j=1; j<=m; j++)
			{
				c.SetValue(i,j,-tauinv(i,j)/aa);
			}
		}
		
		for (int initer=1; initer<=niter; initer++)
		{
			for (int k=1; k<=m; k++)
			{
				int i = kstep(k,1);
				aa = 0;
				for (int j=1; j<=m; j++)
				{
					if (i != j)
					{
						aa = aa + c(i,j) * z(j,1);
					}
				}
				
				if (la(i,1)==1)
				{
					t1=normrt_rnd(0,1,(b1(i,1)-aa)/h(i,1));
				}
				else 
				if (lb(i,1)==1)
				{
					t1=normlt_rnd(0,1,(a1(i,1)-aa)/h(i,1));
				}
				else
				{
					t1=normt_rnd(0,1,(a1(i,1)-aa)/h(i,1),(b1(i,1)-aa)/h(i,1));
				}
				z.SetValue(i,1,aa + h(i,1)*t1);
				}
		}
		
		//Transform back to x
		xdraw = dinv*z;
		for (int row=1; row<=m; row++)
		{
			xdraw_matrix.SetValue(row,col,xdraw(row,1) + amu(row,1));
		}
	}
		return xdraw_matrix;
}


/**** The Gaussian Distribution PDF****/
double PseudoRandom::gaussian_pdf(const double x, const double mu, const double sigma)
{
	double u;
	double p;
	
	u = (x - mu) / fabs(sigma);
	p = (1 / (sqrt (2 * M_PI) * fabs (sigma))) * exp (-u * u / 2);
	
	return p;
}

/**** The Bivariate Gaussian Distribution PDF****/
double PseudoRandom::bivariate_gaussian_pdf(const double x, const double y, const double sigma_x, const double sigma_y, const double rho)
{
  double u = x / sigma_x ;
  double v = y / sigma_y ;
  double c = 1 - rho*rho ;
  double p = (1 / (2 * M_PI * sigma_x * sigma_y * sqrt(c))) 
    * exp (-(u * u - 2 * rho * u * v + v * v) / (2 * c));
  return p;
}

/**** The Exponential Distribution PDF****/
double PseudoRandom::exponential_pdf(const double x, const double mu)
{
  if (x < 0)
    {
      return 0 ;
    }
  else
    {
      double p = exp (-x/mu)/mu;
      
      return p;
    }
}

/**** The Uniform Distribution PDF****/
double PseudoRandom::flat_pdf(double x, const double a, const double b)
{
  if (x < b && x >= a)
    {
      return 1 / (b - a);
    }
  else
    {
      return 0;
    }
}

//export template <class T>
//T PseudoRandom::mean(const T data[], const unsigned int stride, const unsigned int size)
//{
//	/* Compute the arithmetic mean of a dataset using the recurrence relation 
//     mean_(n) = mean(n-1) + (data[n] - mean(n-1))/(n+1)   */
//
//	long double mean = 0;
//	unsigned int i;
//
//	for (i = 0; i < size; i++)
//    {
//		mean += (data[i * stride] - mean) / (i + 1);
//    }
//
//  return mean;
//}
	

/*** Others  ***/

double PseudoRandom::gamma_int(const unsigned int a)
{
	if (a < 12)
	{
		unsigned int i;
      	double prod = 1;
      	
      	for (i = 0; i < a; i++)
        {
        	prod *= randDblExc();
        }
        
        /* Note: for 12 iterations we are safe against underflow, since
         * the smallest positive random number is O(2^-32). This means
         * the smallest possible product is 2^(-12*32) = 10^-116 which
         * is within the range of double precision. 
         */
         
         return -log (prod);
	}
	else
	{
		return gamma_large((double) a);
	}
}

double PseudoRandom::gamma_large(const double a)
{
	/* Works only if a > 1, and is most efficient if a is large
	 * This algorithm, reported in Knuth, is attributed to Ahrens.  A
	 * faster one, we are told, can be found in: J. H. Ahrens and
	 * U. Dieter, Computing 12 (1974) 223-246.  
	 */
     
	double sqa, x, y, v;
	sqa = sqrt (2 * a - 1);
  	do
  	{
  		do
  		{
			y = tan (M_PI * randExc());
			x = sqa * y + a - 1;
  		}
  		while (x <= 0);
  		v = randExc();
  	}
  	while (v > (1 + y * y) * exp ((a - 1) * log (x / (a - 1)) - sqa * y));
  		
  	return x;
}

double PseudoRandom::gamma_frac(const double a)
{
	// This is exercise 16 from Knuth; see page 135, and the solution is on page 551.
     
	double p, q, x, u, v;
	p = M_E / (a + M_E);
	do
	{
		u = randExc();
		v = randDblExc();
		
		if (u < p)
		{
			x = exp ((1 / a) * log (v));
			q = exp (-x);
		}
		else
		{
			 x = 1 - log (v);
			 q = exp ((a - 1) * log (x));
		}
	}
	while (randExc() >= q);
	
	return x;
}

double PseudoRandom::normal_01_cdf ( double x )
{
  double a1 = 0.398942280444E+00;
  double a2 = 0.399903438504E+00;
  double a3 = 5.75885480458E+00;
  double a4 = 29.8213557808E+00;
  double a5 = 2.62433121679E+00;
  double a6 = 48.6959930692E+00;
  double a7 = 5.92885724438E+00;
  double b0 = 0.398942280385E+00;
  double b1 = 3.8052E-08;
  double b2 = 1.00000615302E+00;
  double b3 = 3.98064794E-04;
  double b4 = 1.98615381364E+00;
  double b5 = 0.151679116635E+00;
  double b6 = 5.29330324926E+00;
  double b7 = 4.8385912808E+00;
  double b8 = 15.1508972451E+00;
  double b9 = 0.742380924027E+00;
  double b10 = 30.789933034E+00;
  double b11 = 3.99019417011E+00;
  double cdf;
  double q;
  double y;
//
//  |X| <= 1.28.
//
  if ( fabs ( x ) <= 1.28 )
  {
    y = 0.5 * x * x;

    q = 0.5 - fabs ( x ) * ( a1 - a2 * y / ( y + a3 - a4 / ( y + a5 
      + a6 / ( y + a7 ) ) ) );
//
//  1.28 < |X| <= 12.7
//
  }
  else if ( fabs ( x ) <= 12.7 )
  {
    y = 0.5 * x * x;

    q = exp ( - y ) * b0 / ( fabs ( x ) - b1 
      + b2 / ( fabs ( x ) + b3 
      + b4 / ( fabs ( x ) - b5 
      + b6 / ( fabs ( x ) + b7 
      - b8 / ( fabs ( x ) + b9 
      + b10 / ( fabs ( x ) + b11 ) ) ) ) ) );
//
//  12.7 < |X|
//
  }
  else
  {
    q = 0.0;
  }
//
//  Take account of negative X.
//
  if ( x < 0.0 )
  {
    cdf = q;
  }
  else
  {
    cdf = 1.0 - q;
  }

  return cdf;
}

double PseudoRandom::normal_01_cdf_inv ( double p )
{
  double a[8] = { 
    3.3871328727963666080,     1.3314166789178437745e+2,
    1.9715909503065514427e+3,  1.3731693765509461125e+4,
    4.5921953931549871457e+4,  6.7265770927008700853e+4,
    3.3430575583588128105e+4,  2.5090809287301226727e+3 };
  double b[8] = {
    1.0,                       4.2313330701600911252e+1,
    6.8718700749205790830e+2,  5.3941960214247511077e+3,
    2.1213794301586595867e+4,  3.9307895800092710610e+4,
    2.8729085735721942674e+4,  5.2264952788528545610e+3 };
  double c[8] = {
    1.42343711074968357734,     4.63033784615654529590,
    5.76949722146069140550,     3.64784832476320460504,
    1.27045825245236838258,     2.41780725177450611770e-1,
    2.27238449892691845833e-2,  7.74545014278341407640e-4 };
  double const1 = 0.180625;
  double const2 = 1.6;
  double d[8] = {
    1.0,                        2.05319162663775882187,
    1.67638483018380384940,     6.89767334985100004550e-1,
    1.48103976427480074590e-1,  1.51986665636164571966e-2,
    5.47593808499534494600e-4,  1.05075007164441684324e-9 };
  double e[8] = {
    6.65790464350110377720,     5.46378491116411436990,
    1.78482653991729133580,     2.96560571828504891230e-1,
    2.65321895265761230930e-2,  1.24266094738807843860e-3,
    2.71155556874348757815e-5,  2.01033439929228813265e-7 };
  double f[8] = {
    1.0,                        5.99832206555887937690e-1,
    1.36929880922735805310e-1,  1.48753612908506148525e-2,
    7.86869131145613259100e-4,  1.84631831751005468180e-5,
    1.42151175831644588870e-7,  2.04426310338993978564e-15 };
  double q;
  double r;
  double split1 = 0.425;
  double split2 = 5.0;
  double value;

  if ( p <= 0.0 )
  {
    value = -d_huge ( );
    return value;
  }

  if ( 1.0 <= p )
  {
    value = d_huge ( );
    return value;
  }

  q = p - 0.5;

  if ( fabs ( q ) <= split1 )
  {
    r = const1 - q * q;
    value = q * dpoly_value ( 8, a, r ) / dpoly_value ( 8, b, r );
  }
  else
  {
    if ( q < 0.0 )
    {
      r = p;
    }
    else
    {
      r = 1.0 - p;
    }

    if ( r <= 0.0 )
    {
      value = -1.0;
      exit ( 1 );
    }

    r = sqrt ( -log ( r ) );

    if ( r <= split2 )
    {
      r = r - const2;
      value = dpoly_value ( 8, c, r ) / dpoly_value ( 8, d, r ); 
     }
     else
     {
       r = r - split2;
       value = dpoly_value ( 8, e, r ) / dpoly_value ( 8, f, r );
    }

    if ( q < 0.0 )
    {
      value = -value;
    }

  }

  return value;
}

double PseudoRandom::d_huge ( void )
{
  return HUGE_VAL;
}

double PseudoRandom::dpoly_value ( int n, double a[], double x )
{
	{
  int i;
  double value;

  value = 0.0;

  for ( i = n-1; 0 <= i; i-- )
  {
    value = value * x + a[i];
  }

  return value;
}
}