# I was greatly inspired by StatsWithJuliaBook
# https://statisticswithjulia.org/index.html
using StatsPlots, Distributions, Random, LaTeXStrings

The Monty Hall problem

Random.seed!(1);

# Play the Monty Hall game with
# false: I keep the same door
# true: I change door
function montyHall(switchPolicy)
    nDoor = 3
    prize, choice = rand(1:nDoor), rand(1:nDoor)
    if prize == choice
        revealed = rand(setdiff(1:nDoor,choice))
    else
        revealed = rand(setdiff(1:nDoor,[prize,choice]))
    end

    if switchPolicy
        choice = setdiff(1:nDoor,[revealed,choice])[1]
    end
    return choice == prize
end

function montyHall(switchPolicy, nDoor)
    # n Doors!
    # We can change after n-2 doors have been opened
    doors = collect(1:nDoor)
    prize, choice = rand(doors), rand(doors)
    
    # Not the smartest implementation but it is clear
    while length(doors)>2
        if prize == choice
            revealed = rand(setdiff(doors, choice))
        else
            revealed = rand(setdiff(doors,[prize,choice]))
        end
        setdiff!(doors, revealed)
    end

    if switchPolicy
        choice = setdiff(doors,choice)[1]
    end
    return choice == prize
end

montyHall (generic function with 2 methods)

N = 10^3
println("Success probability with policy I (stay): ", 
	sum([montyHall(false) for _ in 1:N])/N)
println("Success probability with policy II (switch): ", 
	sum([montyHall(true) for _ in 1:N])/N)

Success probability with policy I (stay): 0.315
Success probability with policy II (switch): 0.652

# 100 doors, open 98, do Zyou still want to keep your original choice?
N = 10^3
println("Success probability with policy I (stay): ", 
	sum([montyHall(false,100) for _ in 1:N])/N)
println("Success probability with policy II (switch): ", 
	sum([montyHall(true,100) for _ in 1:N])/N)

Success probability with policy I (stay): 0.011
Success probability with policy II (switch): 0.985

Plot distribution (discrete and continuous)

pDiscrete = [0.25, 0.25, 0.5]
xGridD = 0:2
pContinuous(x) = 3/4*(1 - x^2)
xGridC = -1:0.01:1
pContinuous2(x) = x < 0 ? x+1 : 1-x

pContinuous2 (generic function with 1 method)

p1 = plot(xGridD, line=:stem, pDiscrete, marker=:circle, c=:blue, ms=6, msw=0, ylabel = "PMF")
p2 = plot(xGridC, pContinuous.(xGridC), c=:blue, ylabel = "PDF (Density)")
p3 = plot(xGridC, pContinuous2.(xGridC), c=:blue, ylabel = "PDF (Density)");
plot(p1, p2, p3, layout=(1,3), legend=false, ylims=(0,1.1), xlabel=L"x", size=(800, 600))

Transformation

N = 10000
X = rand(N) # Uniform on [0,1]
Y = X.^2 # Takes the square of every X[i]
histogram([X,Y], alpha = 0.5, label = ["X" "Y"], norm = :pdf, xlabel = L"x")
plot!(x-> 1, 0:0.01:1, c = 1, lw = 3, label = L"f_X(x)")
plot!(x-> 1/(2*sqrt(x)), 0:0.01:1, c = 2, lw = 3, label = L"f_Y(x)")

Chebychev (Markov) bound

BienayméChebychev(d,a) = var(d)/a^2

BienayméChebychev (generic function with 1 method)

dc = Normal(0,0.25)
xr = 0:0.01:5
plot(x-> 2(1-cdf(dc,x)),xr, label=L"P(|Z-0|\geq a)", c=3)
plot!(x-> BienayméChebychev(dc,x), xr, label="Bienaymé Chebychev")
title!("$dc")
ylims!(0,1)
xlabel!(L"a")

Are theses bounds sharp? i.e. is the bound reached by some distributions?

Exemple: $\mathbb{P}(X=-1) = \mathbb{P}(X=-1) = 1/(2k^2)$ and $ \mathbb{P}(X=0) = 1-1/k^2$ for $k\geq 1$

k = 2 # k≥1
xs = [-1.,0,1] # Values of X
ps = [1/(2*k^2),1-1/k^2, 1/(2*k^2)] # Associated probability
d = DiscreteNonParametric(xs, ps)

println("E(X) = ", mean(d)) # 0
println("Var(X) = ", var(d)) # 1/k^2

xr = 0:0.01:5
plot(x-> 2(1-cdf(d,x)),xr, label=L"P(|X-0|\geq a)")
plot!(x-> BienayméChebychev(d,x), xr, label="Bienaymé Chebychev")
plot!(x-> 2(1-cdf(dc,x)),xr, label=L"P(|Z-0|\geq a)", c=3)
ylims!(0,1)
xlabel!(L"a")

E(X) = 0.0
Var(X) = 0.25