Non-uniform distributed random numbers

Creating non-uniform distrubed random numbers is not straightforward. Usally programming languages like C provide functions like rand() which returns an integer random number between 0 and RAND_MAX or drand48() with returns a double random number between 0 and 1. The numbers returned by these functions are uniform distributed meaning that if you run them very often you get each number by the same amount.
If you want random numbers in other ranges, these can be rescaled easily. For example, if you need integer random numbers between 1 and 10 you simply could get a random number x from rand() and apply the function $$f(x) = 1 + 9 \cdot x$$ to it.

But what if you want to distribute them in a non-uniform way? Assume we want to distribute them in a power-law $$x^n$$ between $$x_\mathrm{min}$$ and $$x_\mathrm{max}$$. The probability distribution function P(x) would be given as
$P(x) = \begin{cases}\frac{1}{N} x^n & x_\mathrm{min} \leq x \leq x_\mathrm{max} \\ 0 & \mathrm{otherwise}\end{cases}$
where N is a constant to assure that the integral over all propabilities is 1. If $$n \ne -1$$, N can be calculated as
$N =\int_{x_\mathrm{min}}^{x_\mathrm{max}} P(x) \, dx = \frac{x_\mathrm{max}^{n+1}-x_\mathrm{min}^{n+1}}{n+1}$

To get random numbers which follow this probability distribution function we can use the inverse transformation method. Therefore we need the inverse of the cumulative distribution function. The cumulative distribution function F(x) is given as
$F(x) = \int_{-\infty}^x P(x')\,dx' = \frac{1}{N} \frac{x^{n+1}-x_\mathrm{min}^{n+1}}{n+1} = \frac{x^{n+1}-x_\mathrm{min}^{n+1}}{x_\mathrm{max}^{n+1}-x_\mathrm{min}^{n+1}}$
Now we can solve this for $$F^{-1}(x)$$ and get
$F^{-1}(x) = \left(x \left( x_\mathrm{max}^{n+1}-x_\mathrm{min}^{n+1} \right) + x_\mathrm{min}^{n+1} \right)^{\frac{1}{n+1}}$

The same calculus can be done for $$n=1$$ and we receive
$F^{-1}(x) = x_\mathrm{min} \cdot \mathrm{e}^{x \ln \frac{x_\mathrm{max}}{x_\mathrm{min}}}$

We can apply these functions to our uniform distributed random numbers between 0 and 1 and receive random numbers in the desired range. The following figures display the distribution we get for $$n=1$$ with $$x_\mathrm{min}=1$$ and $$x_\mathrm{max}=10$$ for different sample sizes (N=100,10000,1000000). We binned the random numbers into bins of 0.1 width and rescaled them by the total number of random numbers generated. For comparison, we plotted the probability distribution function on top of it.

The sample data was created using the following example code:

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

const unsigned int N = 100;
const double MIN = 1;
const double MAX = 10;
const double EXPONENT = 2;

double f(double x) {
if (EXPONENT == -1) {
return exp(log(MAX-MIN)*x)*MIN;
} else {
return exp(log(x*(-pow(MIN,EXPONENT+1)+pow(MAX,EXPONENT+1))+pow(MIN,EXPONENT+1))/(EXPONENT+1));
}
}

int main(int argc, char** argv) {
// provide initial random seed
srand(time(NULL));

printf("# i\tuniform\tnon-uniform\n");
for (unsigned int i = 0; i < N; ++i) {
double x = drand48();

printf("%i\t%.18lg\t%.18lg\n", i, x, f(x));
}

exit(EXIT_FAILURE);
}

This entry was posted in Software and tagged by twam. Bookmark the permalink.