# Identifying and understanding algorithm of random number generation

Studying different ways to generate random numbers according to a distribution and the below algorithm describes the "box method". A search on Google led to the Box-Mueller method. Are they related? Also, what would be a simple implementation of this algorithm for $$f(x)=\sin{(x)}$$ on $$[0,\pi]$$? Is it correct that $$y$$ should be scaled to $$[0,1]$$?

We generate two random numbers $$x$$ and $$y$$. We scale $$x$$ so that it gives a random point in the restricted range we want to generate random numbers in. Now we scale $$y$$ so it matches the range from 0 to the maximum value of the function $$f(x)$$ we want to generate. Now we accept $$x$$ if $$y < f(x)$$ and reject $$x$$ otherwise.

Implementation example:

#include <math.h>

int main(){

double x = (double)rand();
double y = (double)rand();

double xs = fmod(x,pi);
double ys = fmod(y,1);

if(ys < sin(xs)) return sin(xs);

}


First of all $$f(x)$$ needs to be a density, so $$\sin(x)$$ doesn't do the job. Let's sample from standard normal distribution, i.e. $$f(x)=\frac{1}{\sqrt{2\pi}}e^{-x^2/2}$$ in the range $$[-3,3]$$, where most of the samples reside.

A note for your C code: You shouldn't return the number from main function. Either write to a file or print out for plotting (or directly plot if you have suitable libraries). Most random samplers give you uniform RVs in $$[0,1]$$, so let's assume we have such $$X,Y$$. You can achieve that in C by dividing the generated number with INT_MAX macro.

First we scale $$X$$ to our new range, i.e. $$X\leftarrow6X-3$$. Now, our $$X$$ is uniform in the range $$[-3,3]$$. Also, $$\max f(x)=1/\sqrt{2\pi}$$, and therefore let $$Y\leftarrow \sqrt{2\pi}Y$$. Then, we compare $$Y$$ and $$f(X)$$ and accept the scaled $$X$$ if $$Y$$ is smaller. Here is an example C code:

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

int main()
{
srand(time(NULL));

int N = 100000;
for(int i = 0; i < N; i++){
double x = ((double) rand()) / INT_MAX;
double y = ((double) rand()) / INT_MAX;

x = 6 * x - 3;
y = y / sqrt(2*M_PI);

if( y < 1/sqrt(2*M_PI) * exp(-x*x/2) )
printf("%f\n", x);
}

return 0;
}


If you redirect the standard output to a file, via > opearator while running, you'll see a histogram as follows:

Now, why does it work?

By sampling $$X$$ within $$[a,b]$$ (i.e. the interested range) and $$Y$$ within $$[0,f_{max}]$$, you create a box in the x-y plane, and sample random points on it as follows, just like you're throwing darts to a rectangular window:

Then, you accept the points which are under the curve, i.e. $$y, which actually makes up the histogram for you, i.e. the higher $$f(x)$$ in the region, the higher the probability of acceptance.