%% Simulation of continuous probability densities
%% 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)

%% 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


%% 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

%%
% 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)
%%
% which is approximately
double(b)

%%
% so the difference is approximately
double(b-a)
%%
% We could even use this to approximate $\pi$, since $b=2/\pi$.
piapprox = 2/a
%% 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



%%
% The actual area is $\pi/4$.
pi/4
%%
% The error is very small. This gives us an estimate for $\pi$.
piapprox = 4*a