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

Median

function my_median(x::Vector)
    n = length(x)
    if isodd(n)
        return sort(x)[(n+1)÷2]
    else
        #! This is a convention that can be changed
        # takes the smallest of the two middle points -> this definition works with qualitative variables
        return sort(x)[n÷2] 
    end
end

my_median (generic function with 1 method)

my_median(["A","F","B","C","D","E"])

"C"

my_median([1, 5, 6, 8, 0.2, π, -1.2,1])

1.0

Quantile

function my_quantile(x::Vector, α::Real)
    n = length(x)
    if !isinteger(n*α) # If it is NOT an integer (! denote the negation here)
        return sort(x)[ceil(Int, n*α)]
    else 
        #! Here it works only for reals
        return 1/2*(sort(x)[Int(n*α)] + sort(x)[Int(n*α) + 1] ) # some convention
    end
end

my_quantile (generic function with 1 method)

my_quantile(["A","F","B","C","D","E","F","D","A"], 0.99)

"F"

my_quantile([1, 5, 6, 8, 0.2, π, -1.2,1], 0.5)

2.0707963267948966

N = 100
distribution = Exponential()
X = rand(distribution, N)
α_range = 1e-3:1e-3:1-1e-3
plot(α_range, quantile(distribution, α_range), label = "Theoritical quantile")
plot!(α_range, quantile(X, α_range), label = "Empirical quantile")
xlabel!("α", legend = :topleft)

Summary of one data set

histogram(X, label="Some data")

summarystats(X)

Summary Stats:
Length:         100
Missing Count:  0
Mean:           0.940679
Minimum:        0.004529
1st Quartile:   0.289284
Median:         0.696518
3rd Quartile:   1.301311
Maximum:        3.828397

Comparing two data sets

Laplace distribution $ f_{\mathrm{Laplace}}(x) = \dfrac{\lambda}{2}e^{-\lambda |x|} $

μ = 20
d₁ = Normal(μ, μ)
d₂ = Laplace(μ, μ/1.5)
 
n = 1000
X = rand(d₁,n)
Y = rand(d₂,n);
println("E(X) = ", mean(X), ", σ₁ = ", std(X))
println("E(Y) = ", mean(Y), ", σ₂ = ", std(Y))

E(X) = 21.03795030904557, σ₁ = 20.607354180304306
E(Y) = 20.939345572216812, σ₂ = 18.944859691781925

histogram([X,Y], alpha = 0.5, label=["Normal distribution Data" "Laplace distribution Data"], yscale =:identity)

qqnorm(X, ms=3, msw=0, label="Normal distribution Data")
qqnorm!(Y, ms=3, msw=0, label="Laplace distribution Data",
        xlabel="Normal Theoretical Quantiles",
        ylabel="Data Empirical Quantiles", legend=:topleft)

Another example

μ = 20
d₁ = Normal(-2,sqrt(5/3))
d₂ = TDist(5) # Student distribution
 
n = 10_000
X = rand(d₁,n)
Y = rand(d₂,n);
println("E(X) = ", mean(X), ", σ₁ = ", std(X))
println("E(Y) = ", mean(Y), ", σ₂ = ", std(Y))

E(X) = -2.0011132997932584, σ₁ = 1.2829128841947814
E(Y) = -0.010827771686638463, σ₂ = 1.2553641434714244

histogram([X,Y], alpha = 0.5, label=["Normal distribution Data" "Student's t-distribution Data"])

qqnorm(X, ms=3, msw=0, label="Normal Data")
qqnorm!(Y,  ms=3, msw=0, label="Student's t-distribution Data",
        xlabel="Normal Theoretical Quantiles",
        ylabel="Data Empirical Quantiles", legend=:topleft)

qqplot(X, Y, ms=3, msw=0, c = 3, label=:none)
xlabel!("Normal Theoretical Quantiles")
ylabel!("Student's t-distribution Data")