Thread: How to get a random number that follows a particular distribution?

OK, a PDF (probability density function) is a function, such that the closed integral from x=a to x=b is equal to the probability that a random variable which follows this distribution will be between a and b.

One example is the standard normal distribution (gaussian with mean = 0, std deviation =1): x will be between -1.96 and +1.96 95% of the time, because the integral from -1.96 to 1.96 of the PDF is 0.95.

And I do see what you're saying about i(F(x)) -- F(x) is the cumulative density function (CDF), and can be derived from the PDF by simple integration. I hadn't thought about the inverse of this, but I believe you are correct, if you distribute a value between 0 and 1 uniformly, this variable should have the same properties as the original distribution.

Thanks! Now, as I think I can integrate and invert the functions by hand, i should be able to use this.

Well, after looking at it, I *can't* integrate by hand, as there IS no explicit expression for the CDF of a gaussian distribution.

However, being as the inverse-CDF of a gaussian distro is highly desirable, there's a very nifty algorithm which I found on the web which approximates the inverse CDF with incredibly low error (something like +- 10^-18 or somesuch).

With that, I can at least generate gaussian variables easily now, something I couldn't do this morning, and need to do before next week friday.

If (X,Y) is a random vector uniformly distributed within the unit circle , ie. it has a density f such that
f(x,y)=1 if x*x+y*y<1
=0 otherwise
and S=X*X+Y*Y then
r1=X*sqrt(-2*log(S)/S)
r2=Y*sqrt(-2*log(S)/S)
then r1 and r2 are independent and normally distributed with mean 0 and variance 1.
so a possible function could be

#include <math.h>
#include <time.h>
#include <stdlib.h>
double gennor(void)
{
static int i=0;
double x,y,s;
static double r1,r2;
if(i==1)
{
i=0;
return r2;
}
do
{
x = 2.*((double) rand()/(RAND_MAX+1.0))-1. ;
y = 2.*((double) rand()/(RAND_MAX+1.0))-1. ;
s=x*x+y*y;
} while(s>1.);
r1=x*sqrt(-2.*log(s)/s);
r2=y*sqrt(-2.*log(s)/s);
i=0;
return r1;
}

int main(void)
{
double d;
int i;
srand( time(0) );
for(i=0;i<1000;i++)
printf("%lf\n",gennor());
return 0;
}

there is a slight mistake above
the density above is actually
f(x,y) = 1/pi if x*x+y*y<1
= 0 ow
the program is ok