I like to show you how the restriction method can be implemented in EMT. Again, it would be faster to do this in C, but the EMT method is fast enough if it is done in the right way.

The restriction method uses the generator for one distribution to produce samples for another distribution. It can be applied, if we have a generator for a distribution with density f and want a generator for the density g, and if there is a constant c>0 such that

\(g(x) \le c f(x)\)

for all x. Our example is the Cauchy distribution with density

\(f(x) = \dfrac{1}{\pi(1+x^2)}\)

which has an easy generator. Written as an EMT function, we can generate a matrix of nxm random values with

function prand(n,m) := -1/tan(pi*random(n,m))

The reason for this is that the cumulative distribution is simply

\(F(x) = \dfrac{1}{\pi} \arctan(x) + \dfrac{1}{2}\)

as you can check by taking the derivative and the value at minus infinity. With some trigonometric identities you get the inverse in the form

\(x = – \dfrac{1}{\tan(\pi x)}\)

Now we want a generator for

\(g(x) = \dfrac{c}{1+|x|^3}\)

The constant c can be computed with

>&1/(2*integrate(1/(1+x^3),x,0,inf))

                                   3/2
                                  3
                                  ----
                                  4 pi

But it is the beauty of the rejection method that c does not matter. In fact the cumulative distribution can also be computed exactly. But for the inverse, we would need numerical methods. So we completely avoid this road.

We implement the rejection method with the following code.

>function p1randr(n,m) ...
$k=n*m;
$v=[];
$repeat
$  x=prand(1,2*k);
$  a=random(1,2*k);
$  i=nonzeros(2*a/(1+x^2)<1/(1+abs(x)^3)); 
$  v=v|x[i]; 
$  until length(v)>k;
$end;
$return redim(v[1:k],n,m);
$endfunction

This generates twice as much values as needed and accepts the ones with

\(\dfrac{2a}{1+x^2} < \dfrac{1}{1+|x|^3}.\)

where the value a is random in [0,1].

Why does this work? To see this, we look at the plots of g and cf. In our case this looks like the following plot.

with the curves

\(\dfrac{1}{1+|x|^3}, \dfrac{2}{1+x^2}.\)

Now we can check that (x,acf(x)) is in fact a random point under the red curve if x is distributed with respect to the distribution belonging to the red function. Since we accept only the points under the black curve, we get random points under the black curve. Thus the accepted values of x are distributed with the desired distribution.

Here is a simulation, which shows the rejected points in red.

>plot2d("2/(1+x^2)",-5,5,color=red); plot2d("1/(1+abs(x)^3)",>add):
>function prand(n,m) := -1/tan(pi*random(n,m))
>x=prand(1,10000); y=random(size(x))*2/(1+x^2);
>c=ones(size(x))*black;
>c[nonzeros(y>1/(1+abs(x)^3))]=red;
>plot2d(x,y,>points,color=c,style=".",>add):