0
$\begingroup$

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);

}
$\endgroup$
0
$\begingroup$

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:

Normal RV Histogram

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:

Dart Box

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.