Simulation of continuous probability densities

Contents

Plot of random points in the unit square

I'll plot 1000 random points in the unit using the command scatter to plot the points. If I give the command scatter(x,y) the points will be fairly large. So I'll give a third argument, designed to make them large enough to see but not too large.

n=10^3
x=rand(1,n);
y=rand(1,n);
scatter(x,y,3)

n =

        1000

Plot of random points and a curve

I'll plot again, adding the curve $y=\sin(x)$ to the plot. I'll use the plot command to plot it, after creating a vector X of points at which to plot it.

n=10^3
x=rand(1,n);
y=rand(1,n);
scatter(x,y,3)
hold on
X = 0:0.001:1;
plot(X,sin(pi*X),'r','LineWidth',2)
hold off

n =

        1000

Monte Carlo simulation

This is used to estimate an integral. Here I use it to estimate the integral of sin(pi*x), 0<x<1.

n=10^6
x=rand(1,n);
y=rand(1,n);
a = sum(y<sin(pi*x))/n

n =

     1000000


a =

    0.6370

I'll compare this with the actual integral. I want to calculate that symbolically, so I first introduce a symbolic variable, t.

syms t
b = int(sin(sym('pi')*t),t,0,1)

 
b =
 
                                      2
                                     ----
                                      pi

which is approximately

double(b)

ans =

    0.6366

so the difference is approximately

double(b-a)

ans =

  -3.4723e-04

We could even use this to approximate $\pi$, since $b=2/\pi$.

piapprox = 2/a

piapprox =

    3.1399

Approximating $\mathbf \pi$ using Monte Carlo simulation

Let E be the circle of radius $1/2$, center $(1/2,1/2)$, $E=\{(x,y):(x-1/2)^2+(y-1/2)^2<1/4$. I use Monte Carlo simulation to estimate the area.

Important note

Because x and y are n component vectors, I have to use .^2 instead of ^2 to square each component. Similarly, to multiply or divide vectors component by component, you have to use .* or ./

n=10^6
x=rand(1,n);
y=rand(1,n);
a = sum((x-1/2).^2+(y-1/2).^2<1/4)/n

n =

     1000000


a =

    0.7851

The actual area is $\pi/4$.

pi/4

ans =

    0.7854

The error is very small. This gives us an estimate for $\pi$.

piapprox = 4*a

piapprox =

    3.1403


Published with MATLAB® R2015a