# I was greatly inspired by StatsWithJuliaBook
# https://statisticswithjulia.org/index.html
using Plots, LaTeXStrings, LinearAlgebra, Distributions; pyplot()
# Plots.PyPlot.rc("text", usetex = "true")
using StatsBase
using Measures

Plotting multidimensional densities

Consider the r.v. $(X,Y)$ with density $f_{X,Y}(x,y) = 9/8(4x+y)\sqrt{(1-x)(1-y)}$

delta = 0.01
grid = 0:delta:1
f(x,y) = 9/8*(4x+y)*sqrt((1-x)*(1-y))
z = [f(x,y) for y in grid, x in grid];

# Need to work on the 3D plot display
p1 = surface(grid, grid, z, 
	c=cgrad([:blue, :red]), la=1, camera=(60,50),
	ylabel="y", zlabel=L"f(x,y)")
p2 = contourf(grid, grid, z, 
	c=cgrad([:blue, :red]))
p2 = contour!(grid, grid, z, xlims=(0,1), ylims=(0,1), ylabel="y", ratio=:equal)

plot(p1, p2, size=(800, 400), xlabel="x", margin=10mm)

Normalization

densityIntegral = sum(z)*delta^2
println("2-dimensional Riemann sum over density: ", densityIntegral)

2-dimensional Riemann sum over density: 1.0063787264382458

$ \mathbb{P}(X<Y)$ = $\int_{x<y} f(x,y) dx dy$

probB = sum([sum([f(x,y)*delta for y in x:delta:1])*delta for x in grid])
println("2-dimensional Riemann sum to evaluate probability: ", probB)

2-dimensional Riemann sum to evaluate probability: 0.3932640388868346

$\operatorname{Cor}(X,Y) = 0$ does not imply Independence.

N = 100
X = rand(Uniform(-1,1), N)
f(x) = 1/2*(3*x^2 - 1)
Y = f.(X)
println("Cor(X,Y) = ", cor(X,Y))
scatter(X,Y,xlabel = "X", ylabel= "Y", label = :none)

Cor(X,Y) = -0.0679111373269616

Same mean, same variance, same covariance BUT different distribution

N = 10^5

SigY = [ 6 4 ;
         4 9]
muY = [15 ; 
       20]
A = cholesky(SigY).L

rngGens = [()->rand(Normal()), 
           ()->rand(Uniform(-sqrt(3),sqrt(3))),
           ()->rand(Exponential())-1]

rv(rg) = A*[rg(),rg()] + muY
    
data = [[rv(r) for _ in 1:N] for r in rngGens]

stats(data) = begin
    data1, data2 = first.(data),last.(data)
    println(round(mean(data1),digits=2), "\t",round(mean(data2),digits=2),"\t",
            round(var(data1),digits=2), "\t", round(var(data2),digits=2), "\t",
            round(cov(data1,data2),digits=2))
end

println("Mean1\tMean2\tVar1\tVar2\tCov")
for d in data
    stats(d)
end

scatter(first.(data[1]), last.(data[1]),  ms=1, msw=0, label="Normal", alpha = 0.75)
scatter!(first.(data[2]), last.(data[2]),  ms=1, msw=0, label="Uniform", alpha = 0.75)
scatter!(first.(data[3]),last.(data[3]), ms=1,msw=0,label="Exponential",
	xlims=(0,40), ylims=(0,40), legend=:bottomright, ratio=:equal,
    xlabel=L"X_1", ylabel=L"X_2", alpha = 0.5)

Mean1	Mean2	Var1	Var2	Cov
15.01	20.0	6.0	9.01	4.01
15.0	19.99	6.0	9.01	4.01
14.99	20.0	5.93	8.93	3.95