Random numbers with arbitrary probability distribution using rand() only

5
votes

I am looking for a simple algorithm how to generate pseudo-random floating point numbers using only ANSI rand() function but with arbitrary probability distribution. For a simple uniform distribution I use following code:

x = (float)rand() / (float)RAND_MAX;

Of course it is not very accurate, but enough for my needs. I need also other distributions like logistic and gaussian. Ideally I have to define an arbitrary pdf using a simple vector of finite length, e.g. for logistic pdf this vector may look like: 

logistic_pdf = {0., 0.26894, 0.33924, 0.41742, 0.5, 0.58257, 0.66075, 1.};

and for uniform (using same dimensionality 8):

uniform_pdf = {0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125};

This is only an idea. But I am not sure how to implement it efficiently using rand()->{0...RAND_MAX} only.

5
calgorithmmathrandomprobability
I think there may be more to this than you might have guessed. I recommend you consult Simulation by Sheldon Ross as it has nice coverage on how you generate random numbers for different distributions in both the discrete and continuous cases. - frankc

5 Answers

3
votes

there is no simple algorithm to do arbitrary complex things. you have to find the inverse probability integral transform for each of your 'arbirary' distributions.

2
votes

There is a linear time initialization constant time sampling algorithm for arbitrary discrete distributions:

http://web.eecs.utk.edu/~vose/Publications/random.pdf

Check it out! It is quite clever and not hard to implement at all.

2
votes

Your question doesn't have a simple answer. Luc Devroye took over 800 pages to address it very comprehensively in his text "Non-Uniform Random Variate Generation".

1
votes

Take a look at this link. Here an example for Poisson distribution:

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

int Poisson ( double ev ) {
      int         n = 0;      // counter of iterations 
      double      em;         // e^(-ev), where v is the expected value
      double      x;          // pseudorandom number

      em = exp (-ev);
      x = rand() / (double) RAND_MAX;     // check your C compiler docs
                                          // for the correct constant name
      while (x > em) { 
            n++;
            x *=  rand() / (double) RAND_MAX;
            }
      return n;
      } 

main () {
  int i;
  for (i = 0; i< 1000; i++) {
    printf("new Poisson value: %d\n", Poisson(.133333) );
    }
  }
0
votes

You're going to have to do some research because the method for drawing psuedo-random values from each of those distributions will differ. A reasonable place to start is Wikipedia. They have methods for generating values from a normal / Gaussian distribution and from a logistic distribution. Others which may be of interest are the exponential distribution, the beta distribution, and the gamma distribution.

Alternatively, if you have some source data you are interested in replicating, you can create a histogram of that data you have and generate a CDF from that data. Then simply generate X~U(0,1), and determine which bin in your histogram this corresponds to, scaling linearly between the upper and lower bounds of your bin. This is the essence of the inverse probability integral transform method stafan mentions.