Monte Carlo
Lecture 11
March 3, 2024
MC Example: Finding \(\pi\)
Finding \(\pi\) by sampling random values from the unit square and computing the fraction in the unit circle. This is an example of Monte Carlo integration.
\[\frac{\text{Area of Circle}}{\text{Area of Square}} = \frac{\pi}{4}\]
Code
Logging.disable_logging(Logging.Info)
function circleShape(r)
θ = LinRange(0, 2 * π, 500)
r * sin.(θ), r * cos.(θ)
end
nsamp = 3000
unif = Uniform(-1, 1)
x = rand(unif, (nsamp, 2))
l = mapslices(v -> sum(v.^2), x, dims=2)
in_circ = l .< 1
pi_est = [4 * mean(in_circ[1:i]) for i in 1:nsamp]
plt1 = plot(
1,
xlim = (-1, 1),
ylim = (-1, 1),
legend = false,
markersize = 4,
framestyle = :origin,
tickfontsize=16,
grid=:false
)
plt2 = plot(
1,
xlim = (1, nsamp),
ylim = (3, 3.5),
legend = :false,
linewidth=3,
color=:black,
tickfontsize=16,
guidefontsize=16,
xlabel="Iteration",
ylabel="Estimate",
right_margin=5mm
)
hline!(plt2, [π], color=:red, linestyle=:dash)
plt = plot(plt1, plt2, layout=Plots.grid(2, 1, heights=[2/3, 1/3]), size=(600, 500))
plot!(plt, circleShape(1), linecolor=:blue, lw=1, aspectratio=1, subplot=1)
mc_anim = @animate for i = 1:nsamp
if l[i] < 1
scatter!(plt[1], Tuple(x[i, :]), color=:blue, markershape=:x, subplot=1)
else
scatter!(plt[1], Tuple(x[i, :]), color=:red, markershape=:x, subplot=1)
end
push!(plt, 2, i, pi_est[i])
end every 100
gif(mc_anim, "figures/mc_pi.gif", fps=3)
MC Example: Dice
What is the probability of rolling 4 dice for a total of 19?
Can simulate dice rolls and find the frequency of 19s among the samples.
Code
function dice_roll_repeated(n_trials, n_dice)
dice_dist = DiscreteUniform(1, 6)
roll_results = zeros(n_trials)
for i=1:n_trials
roll_results[i] = sum(rand(dice_dist, n_dice))
end
return roll_results
end
nsamp = 10000
# roll four dice 10000 times
rolls = dice_roll_repeated(nsamp, 4)
# calculate probability of 19
sum(rolls .== 19) / length(rolls)
# initialize storage for frequencies by sample length
avg_freq = zeros(length(rolls))
std_freq = zeros(length(rolls))
# compute average frequencies of 19
avg_freq[1] = (rolls[1] == 19)
count = 1
for i=2:length(rolls)
avg_freq[i] = (avg_freq[i-1] * (i-1) + (rolls[i] == 19)) / i
std_freq[i] = 1/sqrt(i-1) * std(rolls[1:i] .== 19)
end
plt = plot(
1,
xlim = (1, nsamp),
ylim = (0, 0.1),
legend = :false,
tickfontsize=16,
guidefontsize=16,
xlabel="Iteration",
ylabel="Estimate",
right_margin=8mm,
color=:black,
linewidth=3,
size=(600, 400)
)
hline!(plt, [0.0432], color="red",
linestyle=:dash)
mc_anim = @animate for i = 1:nsamp
push!(plt, 1, i, avg_freq[i])
end every 100
gif(mc_anim, "figures/mc_dice.gif", fps=10)
Airshed Model
Goal: Find the probability of exceeding the 1-hour SO2 average exposure concentration standard, which is 0.14 ppm.
\[\mathbb{P}[\text{SO}_2(\theta) > 0.14] = \int \mathbb{I}(\text{SO}_2(\theta) > 0.14) p(\theta) d\theta\]
Airshed Model
Figure 4: Illustration of the airshed, including notation.
\[\frac{dC}{dt} = \frac{u}{L} C_\text{in} + \frac{S-D}{WHL} - \left(\frac{u}{L} + k\right)C\]