using CairoMakie
using Distributions
using Random
Random.seed!(123);

Let’s start with a toy problem of a multivariate normal distribution of $X$ and $Y$, where $$\begin{gathered} \left[\begin{array}{r}X\\ Y\end{array}\right]\sim N(\left[\begin{array}{r}{\mu }_X\\ {\mu }_Y\end{array}\right],\mathrm{\Sigma }) \\ \mathrm{\Sigma }\sim \left(\begin{array}{ll}{\sigma }_X^2& {\sigma }_X{\sigma }_Y\rho \\ {\sigma }_X{\sigma }_Y\rho & {\sigma }_Y^2\end{array}\right) \end{gathered}$$

const N = 100_000
const μ = [0, 0]
const Σ = [1 0.8; 0.8 1]

const mvnormal = MvNormal(μ, Σ)

FullNormal(
dim: 2
μ: [0.0, 0.0]
Σ: [1.0 0.8; 0.8 1.0]
)

FullNormal(
dim: 2
μ: [0.0, 0.0]
Σ: [1.0 0.8; 0.8 1.0])

ErrorException: syntax: missing comma or ) in argument list

x = -3:0.01:3
y = -3:0.01:3
dens_mvnormal = [pdf(mvnormal, [i, j]) for i in x, j in y]
f, ax, c = contourf(x, y, dens_mvnormal; axis=(; xlabel=L"X", ylabel=L"Y"))
Colorbar(f[1, 2], c)

Colorbar()

Metropolis and Metropolis-Hastings

The Metropolis algorithm widely uses a proposal distribution $J_t({\theta }^{\ast })$ and t to define the next value of the distribution $P^{\ast }({\theta }^{\ast }|data)$ After that the Wilfred Keith Hastings proposed a new algorithm called Metropolis-Hastings algorithm.

Reference

• https://storopoli.github.io/Bayesian-Julia/pages/05_MCMC/

using LinearAlgebra: eigvals, eigvecs
#source: https://discourse.julialang.org/t/plot-ellipse-in-makie/82814/4
function getellipsepoints(cx, cy, rx, ry, θ)
    t = range(0, 2 * pi; length=100)
    ellipse_x_r = @. rx * cos(t)
    ellipse_y_r = @. ry * sin(t)
    R = [cos(θ) sin(θ); -sin(θ) cos(θ)]
    r_ellipse = [ellipse_x_r ellipse_y_r] * R
    x = @. cx + r_ellipse[:, 1]
    y = @. cy + r_ellipse[:, 2]
    return (x, y)
end
function getellipsepoints(μ, Σ; confidence=0.95)
    quant = sqrt(quantile(Chisq(2), confidence))
    cx = μ[1]
    cy = μ[2]

    egvs = eigvals(Σ)
    if egvs[1] > egvs[2]
        idxmax = 1
        largestegv = egvs[1]
        smallesttegv = egvs[2]
    else
        idxmax = 2
        largestegv = egvs[2]
        smallesttegv = egvs[1]
    end

    rx = quant * sqrt(largestegv)
    ry = quant * sqrt(smallesttegv)

    eigvecmax = eigvecs(Σ)[:, idxmax]
    θ = atan(eigvecmax[2] / eigvecmax[1])
    if θ < 0
        θ += 2 * π
    end

    return getellipsepoints(cx, cy, rx, ry, θ)
end

f, ax, l = lines(
    getellipsepoints(μ, Σ; confidence=0.9)...;
    label="90% HPD",
    linewidth=2,
    axis=(; limits=(-3, 3, -3, 3), xlabel=L"\theta_1", ylabel=L"\theta_2")
)
axislegend(ax)
record(f, joinpath(@OUTPUT, "met_anim.gif"); framerate=5) do frame
    for i in 1:100
        scatter!(ax, (X_met[i, 1], X_met[i, 2]); color=(:red, 0.5))
        linesegments!(X_met[i:(i+1), 1], X_met[i:(i+1), 2]; color=(:green, 0.5))
        recordframe!(frame)
    end
end;

LoadError: LoadError: UndefVarError: @OUTPUT not defined
in expression starting at /Users/a182501/quarto/program-chunk/posts/mcmc-julia.ipynb:48

onst warmup = 1_000

f, ax, l = lines(
    getellipsepoints(μ, Σ; confidence=0.9)...;
    label="90% HPD",
    linewidth=2,
    axis=(; limits=(-3, 3, -3, 3), xlabel=L"\theta_1", ylabel=L"\theta_2"),
)
axislegend(ax)
scatter!(
    ax,
    X_met[warmup:(warmup + 1_000), 1],
    X_met[warmup:(warmup + 1_000), 2];
    color=(:red, 0.3),
)

f, ax, l = lines(
    getellipsepoints(μ, Σ; confidence=0.9)...;
    label="90% HPD",
    linewidth=2,
    axis=(; limits=(-3, 3, -3, 3), xlabel=L"\theta_1", ylabel=L"\theta_2"),
)
axislegend(ax)
scatter!(ax, X_met[warmup:end, 1], X_met[warmup:end, 2]; color=(:red, 0.3))