pseudorandom.cpp

随机数类,产生随机数

 */
double PseudoRandom::fdist(const double nu1, const double nu2)
{

	double Y1 =  gamma (nu1 / 2, 2.0);
	double Y2 =  gamma (nu2 / 2, 2.0);

	double f = (Y1 * nu2) / (Y2 * nu1);

	return f;
}

/**** The T Distribution ****/
/* The t-distribution arises in statistics. 
 * If Y_1 has a normal distribution and Y_2 has a chi-squared distribution with \nu degrees of freedom then the ratio,
 * X = { Y_1 \over \sqrt{Y_2 / \nu} }
 * has a t-distribution t(x;\nu) with \nu degrees of freedom. 
 * The t-distribution has the form
 * p(x) dx = (Gamma((nu + 1)/2)/(sqrt(pi nu) Gamma(nu/2)) * (1 + (x^2)/nu)^-((nu + 1)/2) dx
 * The method used here is the one described in Knuth 
 */
double PseudoRandom::tdist (const double nu)
{
	if (nu <= 2)
    {
		double Y1 = gaussian (0, 1.0);
		double Y2 = chisq (nu);

		double t = Y1 / sqrt (Y2 / nu);

		return t;
    }
	else
    {
		double Y1, Y2, Z, t;
		do
        {
			Y1 = gaussian (0, 1.0);
			Y2 = exponential (1 / (nu/2 - 1));

			Z = Y1 * Y1 / (nu - 2);
        }
		while (1 - Z < 0 || exp (-Y2 - Z) > (1 - Z));

      /* Note that there is a typo in Knuth's formula, the line below
       * is taken from the original paper of Marsaglia, Mathematics of
       * Computation, 34 (1980), p 234-256 
       */

		t = Y1 / sqrt ((1 - 2 / nu) * (1 - Z));
		return t;
    }
}

/**** The Beta Distribution ****/
/* The beta distribution has the form
 * p(x) dx = (Gamma(a + b)/(Gamma(a) Gamma(b))) x^(a-1) (1-x)^(b-1) dx
 * The method used here is the one described in Knuth 
 */
double PseudoRandom::beta(const double a, const double b)
{
	double x1 = gamma (a, 1.0);
	double x2 = gamma (b, 1.0);
	
	return x1 / (x1 + x2);
}

/**** The Logistic Distribution ****/
/* The logistic distribution has the form,
 * p(x) dx = (1/a) exp(-x/a) / (1 + exp(-x/a))^2 dx
 * for -infty < x < infty 
 */
double PseudoRandom::logistic(const double a)
{
	double x, z;

	do
    {
		x = randDblExc();
    }
	while (x == 1);

	z = log (x / (1 - x));

	return a * z;
}

/**** The Pareto Distribution ****/
/* The Pareto distribution has the form,
 * p(x) dx = (a/b) / (x/b)^(a+1) dx     for x >= b
 */
double PseudoRandom::pareto(double a, const double b)
{
	double x = randDblExc();

	double z = pow (x, -1 / a);

	return b * z;
}

/**** The Weibull Distribution ****/
/* The Weibull distribution has the form,
 * p(x) dx = (b/a) (x/a)^(b-1) exp(-(x/a)^b) dx
 */
double PseudoRandom::weibull(const double a, const double b)
{
	double x = randDblExc();

	double z = pow (-log (x), 1 / b);

	return a * z;
}

/**** The Type I Gumbel Distribution ****/
/* The Type I Gumbel distribution has the form,
 * p(x) dx = a b exp(-(b exp(-ax) + ax)) dx
 * for -\infty < x < \infty.
 */
double PseudoRandom::gumbel1(const double a, const double b)
{
	double x = randDblExc();

	double z = (log(b) - log(-log(x))) / a;

	return z;
}

/**** The Type II Gumbel Distribution ****/
/* The Type II Gumbel distribution has the form,
 * p(x) dx = b a x^-(a+1) exp(-b x^-a)) dx
 * for 0 < x < \infty.
 */
double PseudoRandom::gumbel2(const double a, const double b)
{
	double x = randDblExc();

	double z = pow(-b / log(x), 1/a);

	return z;
}

/**** The Dirichlet Probability Distribution ****/
/* The Dirichlet probability distribution of order K-1 is 
 * p(\theta_1,...,\theta_K) d\theta_1 ... d\theta_K = (1/Z) \prod_i=1,K \theta_i^{alpha_i - 1} \delta(1 -\sum_i=1,K \theta_i)
 * The normalization factor Z can be expressed in terms of gamma functions:
 * Z = {\prod_i=1,K \Gamma(\alpha_i)} / {\Gamma( \sum_i=1,K \alpha_i)}  
 * The K constants, \alpha_1,...,\alpha_K, must be positive. The K parameters, 
 * \theta_1,...,\theta_K are nonnegative and sum to 1.
 * The random variates are generated by sampling K values from gamma
 * distributions with parameters a=\alpha_i, b=1, and renormalizing. 
 * See A.M. Law, W.D. Kelton, Simulation Modeling and Analysis (1991).
 * Gavin E. Crooks <[e-mail]> (2002)
 */
void PseudoRandom::dirichlet(const unsigned K, const double alpha[], double theta[])
{
	unsigned int i;
	double norm = 0.0;

	for (i = 0; i < K; i++)
    {
		theta[i] = gamma (alpha[i], 1.0);
    }

	for (i = 0; i < K; i++)
    {
		norm += theta[i];
    }

	for (i = 0; i < K; i++)
    {
		theta[i] /= norm;
    }
}

/**** The Poisson Distribution ****/
/* The poisson distribution has the form
 * p(n) = (mu^n / n!) exp(-mu) 
 * for n = 0, 1, 2, ... . 
 * The method used here is the one from Knuth.
 */
unsigned int PseudoRandom::poisson(double mu)
{
	double emu;
  	double prod = 1.0;
  	unsigned int k = 0;
  	
  	while (mu > 10)
    {
    	unsigned int m = mu * (7.0 / 8.0);

      	double X = gamma_int (m);
      	
      	if (X >= mu)
        {
        	return k + binomial (mu / X, m - 1);
        }
        else
        {
        	k += m;
          	mu -= X; 
        }
    }

  // This following method works well when mu is smal

  	emu = exp (-mu);

  	do
  	{
  		prod *= randExc();
      	k++;
    }
    while (prod > emu);
    
    return k - 1;
}

/**** The Bernoulli Distribution ****/
/* The bernoulli distribution has the form,
 * prob(0) = 1-p, prob(1) = p
 */
unsigned int PseudoRandom::bernoulli(double p)
{
	double u = randExc();
	
	if (u < p)
	{
		return 1 ;
	}
	else
	{
		return 0 ;
	}
}

/**** The Binomial Distribution ****/
/* The binomial distribution has the form,
 * prob(k) =  n!/(k!(n-k)!) *  p^k (1-p)^(n-k) for k = 0, 1, ..., n
 * This is the algorithm from Knuth 
 */
unsigned int PseudoRandom::binomial(double p, unsigned int n)
{
	unsigned int i, a, b, k = 0;
	
	while (n > 10)        // This parameter is tunable 
	{
		double X;
		a = 1 + (n / 2);
		b = 1 + n - a;
		
		X = beta ((double) a, (double) b);
		
		if (X >= p)
		{
			n = a - 1;
			p /= X;
		}
		else
		{
			k += a;
			n = b - 1;
			p = (p - X) / (1 - X);
		}
	}
	
	for (i = 0; i < n; i++)
	{
		double u = randExc();
		if (u < p)
			{
				k++;
			}
	}
	
	return k;
}

/**** The Multinomial Distribution ****/
/* The multinomial distribution has the form
 *
 *                                    N!           n_1  n_2      n_K
 * prob(n_1, n_2, ... n_K) = -------------------- p_1  p_2  ... p_K
 *                           (n_1! n_2! ... n_K!) 
 *
 * where n_1, n_2, ... n_K are nonnegative integers, sum_{k=1,K} n_k = N,
 * and p = (p_1, p_2, ..., p_K) is a probability distribution. 
 * Random variates are generated using the conditional binomial method.
 * This scales well with N and does not require a setup step.
 * Ref: 
 * C.S. David, The computer generation of multinomial random variates,
 * Comp. Stat. Data Anal. 16 (1993) 205-217
 */
void PseudoRandom::multinomial(const unsigned int K, const unsigned int n, const double p[], unsigned int d[])
{
	unsigned int i;
	double norm = 0.0;
	double sum_p = 0.0;

	unsigned int sum_n = 0;

  /* p[k] may contain non-negative weights that do not sum to 1.0.
   * Even a probability distribution will not exactly sum to 1.0
   * due to rounding errors. 
   */

	for (i = 0; i < K; i++)
    {
		norm += p[i];
    }

	for (i = 0; i < K; i++)
    {
		if (p[i] > 0.0)
        {
			d[i] = binomial (p[i] / (norm - sum_p), n - sum_n);
        }
		else
        {
			d[i] = 0;
        }
		sum_p += p[i];
		sum_n += d[i];
    }

}
/**** The Negative Binomial Distribution ****/
/* The negative binomial distribution has the form,
 * prob(k) =  Gamma(n + k)/(Gamma(n) Gamma(k + 1))  p^n (1-p)^k 
 * for k = 0, 1, ... . Note that n does not have to be an integer.
 * This is the Leger's algorithm (given in the answers in Knuth) 
 */
unsigned int PseudoRandom::negative_binomial(double p, double n)
{
	double X = gamma (n, 1.0) ;
	unsigned int k = poisson (X*(1-p)/p) ;
	return k ;
}

/**** The Pascal Distribution ****/
/* The Pascal distribution is a negative binomial with valued integer n
 * prob(k) =  (n - 1 + k)!/(n!(k - 1)!) *  p^n (1-p)^k for k = 0, 1, ..., n
 */
unsigned int PseudoRandom::pascal_(double p, unsigned int n)
{
	/* This is a separate interface for the pascal distribution so that
	 * it can be optimized differently from the negative binomial in future.
	 * e.g. if n < 10 it might be faster to generate the Pascal
	 * distributions as the sum of geometric variates directly.  
	 */
     
	unsigned int k = negative_binomial (p, (double) n);
  	return k;
}

/**** The Geometric Distribution ****/
/* Geometric distribution (bernoulli trial with probability p) 
 * prob(k) =  p (1 - p)^(k-1) for n = 1, 2, 3, ...
 * It gives the distribution of "waiting times" for an event that occurs with probability p. 
 */
unsigned int PseudoRandom::geometric(const double p)
{
	double u = randDblExc();
	
	unsigned int k;
	
	if (p == 1)
	{
		k = 1;
	}
	else
	{
		k = log (u) / log (1 - p) + 1;
	}
	
	return k;
}

/**** The Hypergeometric Distribution ****/
/* The hypergeometric distribution has the form,
 * prob(k) =  choose(n1,t) choose(n2, t-k) / choose(n1+n2,t)
 * where choose(a,b) = a!/(b!(a-b)!) 
 * n1 + n2 is the total population (tagged plus untagged)
 * n1 is the tagged population
 * t is the number of samples taken (without replacement)
 * k is the number of tagged samples found
 */

unsigned int PseudoRandom::hypergeometric(unsigned int n1, unsigned int n2, unsigned int t)
{
	const unsigned int n = n1 + n2;

	unsigned int i = 0;
	unsigned int a = n1;
	unsigned int b = n1 + n2;
	unsigned int k = 0;

	if (t > n)
    {
		t = n ;
    }

	if (t < n / 2) 
    {
		for (i = 0 ; i < t ; i++)
        {
			double u = randExc();
          
			if (b * u < a)
            {
				k++ ;
				if (k == n1)
					return k ;
              a-- ;
            }
			b-- ;
        }
		return k;
    }
	else
    {
		for (i = 0 ; i < n - t ; i++)
        {
			double u = randExc();
          
			if (b * u < a)
            {
				k++ ;
				if (k == n1)
					return n1 - k ;
				a-- ;
            }
			b-- ;
        }
		return n1 - k;
    }
}

/**** The Logarithmic Distribution ****/
/* Logarithmic distribution 
 * prob(n) =   p^n / (n log(1/(1-p)) for n = 1, 2, 3, ...
 * We use Kemp's second accelerated generator, from Luc Devroye's book
 * on "Non-Uniform Random Variate Generation", Springer 
 */
unsigned int PseudoRandom::logarithmic(const double p)
{
	double c = log (1-p) ;

	double v = randDblExc();
  
	if (v >= p)
    {
		return 1 ;
    }
	else
    {
		double u = randDblExc();      
		double q = 1 - exp (c * u);

		if (v <= q*q)
        {
			double x = 1 + log(v)/log(q) ;
			return x ;
        }
		else if (v <= q)
        {
			return 2;
        }