# I was greatly inspired by StatsWithJuliaBook
# https://statisticswithjulia.org/index.html
using Pkg
Pkg.activate("MAA204")
using StatsPlots, LaTeXStrings, Distributions, Random

  Activating environment at `C:\Users\david\Nextcloud\cours\MAA_204_Statistics\notebook_julia\MAA204\Project.toml`

Estimator

Estimator of $U\sim \mathcal{U}([0, a])$

a = 10 # true value (unkown)
N = 400
X = 10*rand(N) # Samples
# Various estimaros
a_1(X) = maximum(X) 
a_2(X) = 2*mean(X)
a_3(X) = 2*median(X)
a_4(X) = sqrt(12*var(X))

a_4 (generic function with 1 method)

plot([n -> a_1(X[1:n]), n-> a_2(X[1:n]), n-> a_3(X[1:n]), n-> a_4(X[1:n])],1:N, xlabel = "Number of sample n",legend=:bottom, label = ["Max(X)" "mean(X)" "median(X)" "sqrt(12 σ^2)"])

Law of Large Number

Example with Exponential distribution: the LLN applies $\mathbb{E}(|X|) = 1/\lambda$ $\Rightarrow$ the moving average converges

n = 10^6
λ = 10
data = rand(Exponential(1 / λ), n) # Julia convention is $\theta = 1/\lambda$
averages = accumulate(+, data) ./ collect(1:n)

plot(1:n, averages,
    legend = :none,
    xscale = :log10, xlims = (1, n), xlabel = "n", ylabel = "Running average")
hline!([1 / λ], c = :black, s = :dot)

Example with Cauchy distribution $f_X = \dfrac{1}{\pi}\dfrac{1}{1+x^2}$: the LLN does NOT apply because $\mathbb{E}(|X|)$ is infinite

$\Rightarrow$ the moving average does not converge

Random.seed!(808)
n = 10^6
U = rand(n)
data = tan.(U * pi .- pi / 2) # Generate Cauchy didtributed r.v. from uniform U with the inverse CDF.
averages = accumulate(+, data) ./ collect(1:n)

plot(1:n, averages,
    legend = :none,
    xscale = :log10, xlims = (1, n), xlabel = "n", ylabel = "Running average")
hline!([1 / λ], c = :black, s = :dot)

Central Limit theorem

n, N = 50, 10^5

dist1 = Uniform(1 - sqrt(3), 1 + sqrt(3))
dist2 = Exponential(1)
dist3 = Poisson(1)

data1 = [mean(rand(dist1, n)) for _ in 1:N]
data2 = [mean(rand(dist2, n)) for _ in 1:N]
data3 = [mean(rand(dist3, n)) for _ in 1:N]

plot(Normal(1, 1/sqrt(n)), label = "Normal(1,1/n)", lw = 2)
stephist!([data1 data2 data3], bins = 100,
    xlabel = L"x", ylabel = "p.d.f.", norm = :pdf,
    label = ["Average of Uniforms" "Average of Exponentials" "Average of Poisson"],
    normed = true, xlims = (0, 2), ylims = (0, 5))

Quantile

qq = 0.95
tau = quantile(Normal(), qq)
mu0, mu1, sd = 0, 20, 1
dist0, dist1 = Normal(mu0, sd), Normal(mu1, sd)
grid = -4:0.1:4
h0grid = tau:0.1:4
p1 = plot(grid, pdf.(dist0, grid), size = (800, 600),
    c = 1, label = L"f_X(x)", xlabel = L"x")
plot!(p1, h0grid, pdf.(dist0, h0grid),
    c = 1, fa = 0.2, fillrange = [0 1], label = "")
annotate!(p1, [(tau, -0.00, L"q_{\alpha}")])
p2 = plot(grid, cdf.(dist0, grid), size = (800, 600),
    c = 2, label = L"F_X(x)", xlabel = L"x")
vline!(p2, [tau], ls = :dash, label = :none)
hline!(p2, [qq], ls = :dash, c = 2, label = :none)
annotate!(p2, [(-4.5, qq - 0.01, L"\alpha")])
annotate!(p2, [(tau, -0.00, L"q_{\alpha}")])
plot(p1, p2, size = (950, 600))

Confidence interval

$t$-distribution example

Imagine the lifetime of some bubble light $X$ is given by a Normal distribution $X\sim \mathcal{N}(1950, 4200)$

With $n = 10$ observations with have the following $\textbf{exact}$ 95% confidence interval $$ I_{0.95} = \left[\bar{X}_{n} - t_{n-1,0.975}\dfrac{\hat{\sigma}_n}{\sqrt{n}}, \bar{X}_{n} + t_{n-1,0.975}\dfrac{\hat{\sigma}_n}{\sqrt{n}} \right] $$

n = 10 # Observation per experiment
N = 1000 # repetition of the experiment
true_mean = 1950
dist = Normal(true_mean, 65)

Normal{Float64}(μ=1950.0, σ=65.0)

X = rand(dist, n)
α = 0.05
t = quantile(TDist(n - 1), 1 - α / 2)
println("Empirical Mean = ", mean(X))
println("Empirical std = ", std(X))
println("Confidence interval  = [", mean(X) - t * std(X) / sqrt(n), " , ", mean(X) + t * std(X) / sqrt(n), "]")
plot(dist, label = "True distribution")
vline!([true_mean], s = :dot, c = :black, label = "True mean")
vline!([mean(X)], label = "Estimate mean")
vline!([mean(X) - t * std(X) / sqrt(n), mean(X) + t * std(X) / sqrt(n)], label = :none, c = 3, s = :dash)
xlabel!("Lifetime (hour)")

Empirical Mean = 1976.467080292447
Empirical std = 74.54915882840042
Confidence interval  = [1923.1378246902468 , 2029.7963358946472]

Bernoulli example

n = 1000 # Observation per experiment
N = 10000 # repetition of the experiment
true_p = 0.1
dist = Bernoulli(true_p)
X = rand(dist, N, n)
p̂ = vec(mean(X, dims = 2))
W = sqrt(n) * (p̂ .- true_p) ./ sqrt.(p̂ .* (1 .- p̂)) # asymptotically Normal(0,1)
plot(Normal(0, 1), label = "Normal(0,1)", lw = 2)
histogram!(W, label = "Distribution of p̂", norm = :pdf, alpha = 0.5)

Is it tight?

Always compare with Chebyshev which is valid for any r.v. (with second moment) and is not an approximation.

n = 10
println("q_{0.975} = ", quantile(Normal(), 0.975))
println("t_{n-1,0.975} = ", quantile(TDist(n-1), 0.975))
println("Chebyshev = ", 1 / (2 * sqrt(0.05)))

q_{0.975} = 1.9599639845400576
t_{n-1,0.975} = 2.262157162798205
Chebyshev = 2.23606797749979

α_range = 0.001:0.001:0.05
plot(α-> 1/(2*sqrt(α)), α_range, label="Chebyshev")
plot!(α-> quantile(Normal(),1-α/2), α_range, label = "Normal")
plot!(α-> quantile(TDist(n-1),1-α/2), α_range, label = "Student n=11")
xlabel!("α")