randompackage.c

lavarnd是一个随机数发生器,要求 2.4.9-31 kernel以上

/* Random Number Generator for non-uniform distributions */
/* Authors : Weili Chen, Zixuan Ma                       */
/* Anyone can use it for any purposes                    */

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>

/* FUNCTION DECLARATIONS */

/* Discrete random variable generators */
int bernoulli(double p);
int binomial(double p, int n);
int negativeBinomial(double p, int r); 
int poisson(int lambda);
int geometric(double p);

/* Continuous random variable generators */
double uniform();
double exponential(double lambda);
double weibull(double k, double lambda);
double normal(double mu, double sigma);
double lognormal(double mu, double sigma);
double chisquare(int dof);
double t(int dof);
double F(int m, int n);
double erlang(int k, double rate);

/* Helper functions for generators */
int P(int m, int n);
int C(int m, int n);
int factor(int n);
double normalhelper();


/* 
 * Generates a Bernoulli random variable 
 * p is the probability of sucess
 * This generator uses the Accept/Reject method
 */
int bernoulli(double p) {
  double u = uniform();

  if(u < p)
    return 1;
  else 
    return 0;
}

/* 
 * Generates a Binomial random variable 
 * p is the probability of sucess
 * n is the number of trials
 * This generator uses the Accept/Reject method
 */
int binomial(double p, int n)
{
  if( n < 0 ) {
    return -1;
  }
  
  double u = uniform();
  double sum = 0;
  int x = 0;
  double temp = C(0, n) * pow(p, 0) * pow(1 - p, n);
  for(; sum < u; x++) {
    sum += temp;
    temp = temp * p / (1 - p) * (n - x) / (x + 1);
  }
  return x - 1;
}


/* 
 * Generates a Negative-Binomial random variable 
 * p is the probability of sucess
 * r is the (number of successes + 1)
 * This function assumes r is an integer
 * This generator uses the Accept/Reject method
 */
int negativeBinomial(double p, int r) 
{
  if( r < 1 ) {
    return -1;
  }
  
  double u = uniform();
  double sum = 0;
  int x = 0;
  for(; sum < u; x ++) {
    sum += C(r - 1, x) * pow(p, r) * pow(1 - p, x);;
  }
  return x - 1;
}


/* 
 * Generates a Poisson random variable 
 * lambda is the mean of the distibution
 * This generator uses the Accept/Reject method
 */
int poisson(int lambda)
{
  if( lambda < 1 ) {
    return -1;
  }
  int x = 0;
  double u = uniform();
  double sum = 0;
  double temp = exp(-lambda) * pow(lambda, 0) / factor(0);
  for(; sum < u; x++) {
    sum += temp;
    temp = temp * lambda / (x + 1);
  }
  return x - 1;
}


/* 
 * Generates a Geometric random variable 
 * p is the probability of success
 * This generator uses the Accept/Reject method
 */
int geometric(double p)
{
  double u = uniform();
  double sum = 0;
  int x = 0;
  for(; sum < u; x ++) {
    sum += p * pow(1 - p, x);
  }
  return x - 1;
}


/* 
 * Generates a Uniform random variable between 0 and 1
 * This generator is part of the C library
 */
double uniform()
{
  int r = rand();
  return (double)r / RAND_MAX;
}


/* 
 * Generates an Exponential random variable 
 * lambda is the rate parameter
 * This generator uses the Inverse Transform  method
 */
double exponential(double lambda)
{
  return -log(uniform()) / lambda;
}


/* 
 * Generates a two-parameter Weibull random variable 
 * k is the shape parameter
 * lambda is the scale parameter
 * This generator uses the Inverse Transform  method
 */
double weibull(double k, double lambda)
{
  return pow(-log(uniform()), 1/k) * lambda;
}


/* 
 * Generates a Normal random variable 
 * mu is the mean of the distribution
 * sigma is the stdev of the distribution
 * This generator uses the Covolution method on N(0,1)
 */
double normal(double mu, double sigma)
{
  return sigma * normalhelper() + mu;
}


/* 
 * Generates a Log-normal random variable 
 * mu is the mean of the log of the variable
 * sigma is the stdev of the log of the variable
 * This generator uses the Convolution method
 */
double lognormal(double mu, double sigma) {
  return exp(sigma * normalhelper() + mu);
}


/* 
 * Generates a Chi-squared random variable 
 * dof is the degree of freedom
 * This generator uses the Convolution method
 */
double chisquare(int dof)
{
  if( dof < 1 ) {
    return 0;
  }
  
  int i = 0;
  double chi = 0;
  for(i = 0; i < dof; i ++) {
    chi += pow(normalhelper(), 2);
  }
  return chi;
}


/* 
 * Generates a t-distribution random variable 
 * dof is the degree of freedom
 * This generator uses the Convolution method
 */
double t(int dof) {
  return normalhelper() / sqrt(chisquare(dof) / dof);
}


/* 
 * Generates a F random variable 
 * m is the degree of freedom of first chi-square distribution
 * n is the degree of freedom of second chi-square distribution
 * This generator uses the Convolution method
 */
double F(int m, int n){
  return (chisquare(m) / m) /(chisquare(n) / n);
}


/* 
 * Generates an Erlang random variable 
 * k is the shape parameter
 * rate is the rate parameter
 * This generator uses the Composition method
 */
double erlang(int k, double rate) {
  
  if (k < 0) return 0;

  int i = 0;
  double erl = 0;
 
  for (i = 0; i < k; i++) {
    erl += -log(uniform());
  }
  
  return erl / rate;
}


/*
 * MAIN FUNCTION
 * Main function only for testing only.
 * It generates a large number of variables for statisical anaylsis.
 * It can be commented out if not needed.
 */
int main(int argc, char **argv)
{
  srand(time(0));
  
  /* Simulate 100000 variables. */
  int numIterations = 100000;

  /* INDICATE IF IT IS A DISCRETE DISTRIBUTION */
  int discrete = 1;
 
  double a[numIterations];
  double sum = 0;
  double min = 0;
  double max = 0;
  int i = 0;
  int j = 0;
  
  /* Simulating N(0,1) variables */
  for(i = 0; i < numIterations; i ++ ) {
    
    /* TRY DIFFERENT DISTRIBUTIONS HERE */
    a[i] = poisson(20);
    sum += a[i];
    
    if (a[i] < min) min = a[i]; 
    if (a[i] > max) max = a[i]; 
  }

  /* Simple statistics printout */
  double mean = sum / numIterations, var = 0;
  
  for(i = 0; i < numIterations; i ++ ) {
    var += (a[i] - mean) * (a[i] - mean);
  }
  
  printf("Mean  = % .3f\n", mean);
  printf("Stdev = % .3f\n", sqrt(var / numIterations));  
  printf("Max   = % .3f\n", max);
  printf("Min   = % .3f\n", min);
  printf("\n");
  
  
  /* FREQUENCY DISTRIBUTION PRINTOUT */

  /* Distributed in 20 buckets for printout. */
  int numBuckets = 40;
  if (discrete) numBuckets = max - min;

  /* Height of distribution for printout. */
  int maxStars = 30;

  int maxBucket = 0;
  double varPerStar = 0;

  printf("Frequency Distribution\n");
  printf("----------------------\n");
  
  int buckets[numBuckets];
  double bucketsize = (max - min) / numBuckets;
  
  for(i = 0; i < numBuckets; i ++ ) {
    buckets[i] = 0;  
  }
  
  for (i = 0; i < numIterations; i++){
    buckets[(int) floor((a[i] - min) / bucketsize)]++;
  }
  
  for(i = 0; i < numBuckets; i ++ ) {
    if (buckets[i] > maxBucket) maxBucket = buckets[i];  
  }
  
  varPerStar = maxBucket / maxStars;
  
  for(i = 0; i < numBuckets; i ++ ) {
    printf("% 10.2f to % 10.2f : ", i * bucketsize + min, (i + 1)* bucketsize + min); 
    
    int bucketStars = ((int) floor(buckets[i] / varPerStar));

    for (j = 0; j < bucketStars; j++){
      printf("*");
    }
    
    for (j = 0; j < maxStars - bucketStars; j++){
      printf(" ");
    }
    
    printf("  %d \n", buckets[i]);   
  }
}



/* 
 * HELPER FUNCTIONS
 * The functions below are helper functions for the random variable generators.
 */


/* Generates a N(0, 1) variable
 * This generator uses the Accept/Reject method
 */
double normalhelper()
{
  double c = sqrt( 2 * M_E / M_PI );
  double t = exponential( 1 );
  
  while( t > (sqrt(2 / M_PI) * exp(t - t * t / 2) / c)) {
    t = exponential( 1 );
  }
  
  if(rand() % 2 == 0 ) t = -1 * t;
  
  return 2 * t;
}


/* Calculates m permutate n */
int P(int m, int n)
{
  int ret = 1, i = 0;
  if( m > n || m < 0 || n < 0 ) {
    return 0;
  }
  
  for(i = n; i > n - m; i --) {
    ret *= i;
  }
  return ret;
}


/* Calculates n factorial */
int factor(int n)
{
  int i = 1;
  int prod = 1;
  
  if( n < 0 ) {
    return 0;
  }
  else if( n == 0 ) {
    return 1;
  }
  else {
    for (i = 1; i <= n; i++)
      prod = prod * i;

    return prod;
  }
}

/* Calculates m choose n */
int C(int m, int n)
{
  int ret = 1;
  if( m > n ) {
    return 0;
  }
  
  return P(m, n) / factor(m);
}