In [1]:

using Plots

n, L = 3, 10^6
X = rand(n, L)
S = vec(sum(X, dims=1))
plot(title="X_i ∼ Uniform(0, 1),  n = $n")
histogram!(S; norm=true, alpha=0.3, bin=201, label="X_1+…+X_n")

Out[1]:

In [2]:

using Plots

n, L = 5, 10^6
X = rand(n, L)
S = vec(sum(X, dims=1))
plot(title="X_i ∼ Uniform(0, 1),  n = $n")
histogram!(S; norm=true, alpha=0.3, bin=201, label="X_1+…+X_n")

Out[2]:

In [3]:

using Plots
f(x,μ,σ²) = 1/√(2π*σ²)*exp(-(x-μ)^2/(2σ²))

n, L = 10, 10^6
X = rand(n, L)
S = vec(sum(X, dims=1))
plot(title="X_i ∼ Uniform(0, 1),  n = $n")
histogram!(S; norm=true, alpha=0.3, bin=201, label="X_1+…+X_n")
plot!(x -> f(x, n/2, n/12), extrema(S)...; label="normal distribution", lw=2, ls=:dash)

Out[3]:

In [4]:

using Distributions, StatsPlots
default(titlefontsize=10)

function plotCLT(; diststr="Gamma(0.1, 10)", n=10, L=10^5, bin=101)
    dist = eval(Meta.parse(diststr))
    X = rand(dist, n, L)
    S = vec(sum(X; dims=1))
    sk = round(skewness(dist), digits=3)
    plot(title="n = $n,  X_i ∼ $diststr,  skewness = $sk")
    histogram!(S; norm=true, alpha=0.3, bin, label="X_1+…+X_n")
    plot!(Normal(n*mean(dist), √(n*var(dist))); label="normal distribution", lw=2, ls=:dash)
end

plotCLT()

Out[4]:

In [5]:

plotCLT(n=100)

Out[5]:

In [6]:

plotCLT(n=1000)

Out[6]:

In [7]:

histogram(rand(Gamma(0.1, 10), 1000); norm=true, alpha=0.3, label="sample of Gamma(0.1, 10)", xlim=(-0.4, 8))
plot!(Gamma(0.1, 10), 0.01, 8; label="Gamma(0.1, 10)", lw=2)

Out[7]:

In [8]:

plotCLT(diststr="Bernoulli(0.3)", n=100, bin=-0.5:60.5, L=10^6)

Out[8]:

In [9]:

plotCLT(diststr="Bernoulli(0.03)", n=100, bin=-0.5:10.5, L=10^6)

Out[9]:

In [10]:

plotCLT(diststr="Bernoulli(0.03)", n=100, bin=-0.5:10.5, L=10^6)
plot!(Poisson(3); label="Poisson distribution", color=:blue)

Out[10]:

In [11]:

using QuadGK
myskewness(d) = quadgk(x -> (x - mean(d))^3*pdf(d, x), extrema(d)...)[1]/std(d)^3
@show skewness(Gamma(0.1, 10))
@show myskewness(Gamma(0.1, 10));

skewness(Gamma(0.1, 10)) = 6.324555320336758
myskewness(Gamma(0.1, 10)) = 6.324555547644784

In [12]:

using Distributions, StatsPlots

MixNormal(a, b) = MixtureModel([Normal(), Normal(b, 1)], [1-a, a])

function Distributions.skewness(d::MixtureModel)
    quadgk(x -> (x - mean(d))^3*pdf(d, x), extrema(d)...)[1]/std(d)^3
end

function Distributions.kurtosis(d::MixtureModel)
    quadgk(x -> (x - mean(d))^4*pdf(d, x), extrema(d)...)[1]/std(d)^4
end

@show std(MixNormal(0.05, 50))
@show skewness(MixNormal(0.05, 50))
@show kurtosis(MixNormal(0.05, 50))

histogram(rand(MixNormal(0.05, 50), 10^4); norm=true, alpha=0.3, 
    label="sample of MixNormal(0.05, 50)", bin=200)
plot!(x -> pdf(MixNormal(0.05, 50), x), -4, 54; 
    label="MixNormal(0.05, 50)", lw=2)

std(MixNormal(0.05, 50)) = 10.943034314119645
skewness(MixNormal(0.05, 50)) = -0.0167645562973185
kurtosis(MixNormal(0.05, 50)) = 0.005270853944982805

Out[12]:

In [13]:

plotCLT(diststr="MixNormal(0.05, 50)", n=30)

Out[13]:

In [14]:

plotCLT(diststr="MixNormal(0.05, 50)", n=100)

Out[14]:

In [15]:

plotCLT(diststr="MixNormal(0.05, 50)", n=300)

Out[15]:

In [ ]: