In [98]:
using Distributions
using Gadfly
In [99]:
T = 1000; # Number of darts launched (simulated)
R = 1; # Circle radius
theta = linspace(0,pi,100);
x = R*cos(theta);
y1 = R*sin(theta);
y2 = -y1;
# Launch darts
points = [rand(Uniform(-R,R), T) rand(Uniform(-R,R), T)];
plot(
layer(x=x, y=y1, Geom.line),
layer(x=x, y=y2, Geom.line),
layer(x=points[:,1], y=points[:,2], Geom.point)
)
Out[99]:
In [100]:
# Count the number of darts inside the circle to estimate the area
insideCircle(x, y) = ((x^2 + y^2) < R^2) ? 1 : 0
N = reduce(+, mapslices(p -> insideCircle(p[1], p[2]), points, 2))
# Estimate Pi
pi_hat = 4*(N/T)
Out[100]:
We would like to approximate the integral $I$ defined by:
$$ I=\int_a^b \frac{1}{\sqrt{2\pi}} \exp^{-x^2/2} dx. $$With $a = -1.96$ and $b = 1.96$ (95% confidence interval).
We rewrite $$ I = \int_a^b g(x)\,f(x)\,dx \\ h(x) = \frac{1}{\sqrt{2\pi}} \exp^{-x^2/2} \\ g(x) = h(x)\,(b-a) \\ f(x) = \frac{\mathbb{1}_{[a,b]}}{b-a} $$
(Théorème de transfert)
$$ I = \mathbb{E}_f[g(X)] = \int_a^b g(x)\,f(x)\,dx = \int_a^b h(x)\,(b-a)\,f(x)\,dx \\ = \frac{(b-a)}{(b-a)} \int_a^b h(x)\,dx = \int_a^b h(x)\,dx $$$\{x_i\}_{1\lt i\lt N} \sim \mathcal{U}_{[a,b]}$, for large N:
$$ \hat{I} = \frac{1}{N}\sum_{i=1}^{N}g(x_i) $$
In [101]:
T = 100;
a = -1.96;
b = 1.96;
x1 = linspace(a, b, T);
x2 = linspace(-5, 5, T);
h(x) = pdf(Normal(0, 1), x);
g(x) = h(x)*(b-a);
# Plot the integral to compute (note: only between a and b)
plot(x=x2, y=h(x2), Geom.line)
Out[101]:
In [102]:
X = rand(Uniform(a, b), T);
G = g(X);
# Estimate I by E_f[g(X)]
I_hat = (1/T) * sum(G)
Out[102]:
In [103]:
for T = [100,1000,10000,100000]
X = rand(Uniform(a, b), T);
G = g(X);
I_hat = (1/T) * sum(G)
println("Estimate of I for T = ", T, ": ", I_hat)
end
In [ ]: