[playground] Julia Turing.jl : Bayesian Cognitive Modeling - Comparing Gaussian means

May 30, 2023 ·playground·julia·turing·bayesian

bayesian cognitive modeling chapter 8 with julia programming language

Github: https://github.com/quangtiencs/bayesian-cognitive-modeling-with-turing.jl

Bayesian Cognitive Modeling is one of the classical books for Bayesian Inference. The old version used WinBUGS/JAG software as the main implementation. You can find other implementations, such as Stan and PyMC, in the below link. I reimplemented these source codes with Julia Programming Language & Turing library in this tutorial.

WinBUGS/JAGS (official) & Stan: https://bayesmodels.com/
PyMC: https://github.com/pymc-devs/pymc-resources/tree/main/BCM

using Pkg
using Logging
using DynamicPPL, Turing
using Zygote, ReverseDiff
using StatsPlots, Random
using LaTeXStrings
using CSV
using DataFrames
using SpecialFunctions
using LinearAlgebra
using FillArrays
using CSV, DataFrames
using LogExpFunctions
using KernelDensity

format=:png

:png

Random.seed!(6)

TaskLocalRNG()

8.1 One-sample comparison

\delta \sim \text{Cauchy} (0, 1)

\sigma \sim \text{Cauchy} (0, 1)_{\mathcal I(0,∞)}

\mu = \delta\sigma

x_{i} \sim \text{Gaussian}(\mu,1/\sigma^2)

winter = [
    -0.05, 0.41, 0.17, -0.13, 0.00, -0.05, 0.00, 0.17, 0.29, 0.04, 0.21, 0.08, 0.37,
    0.17, 0.08, -0.04, -0.04, 0.04, -0.13, -0.12, 0.04, 0.21, 0.17, 0.17, 0.17,
    0.33, 0.04, 0.04, 0.04, 0.00, 0.21, 0.13, 0.25, -0.05, 0.29, 0.42, -0.05, 0.12,
    0.04, 0.25, 0.12
]

summer = [
    0.00, 0.38, -0.12, 0.12, 0.25, 0.12, 0.13, 0.37, 0.00, 0.50, 0.00, 0.00, -0.13,
    -0.37, -0.25, -0.12, 0.50, 0.25, 0.13, 0.25, 0.25, 0.38, 0.25, 0.12, 0.00, 0.00,
    0.00, 0.00, 0.25, 0.13, -0.25, -0.38, -0.13, -0.25, 0.00, 0.00, -0.12, 0.25,
    0.00, 0.50, 0.00
]
x = winter - summer  # allowed because it is a within-subjects design
x = x / std(x);

@model function OneSampleComparison(x)
    delta ~ Cauchy(0,1)
    sigma ~ truncated(Cauchy(0,1), lower=0)

    mu = delta * sigma
    
    for i in eachindex(x)
        x[i] ~ Normal(mu, sigma)
    end
end

iterations=5_000
burnin=2500

model_one_sample_comparison = OneSampleComparison(x)
chain = sample(model_one_sample_comparison, NUTS(), iterations, burnin=burnin)

┌ Info: Found initial step size
└   ϵ = 0.05
Sampling: 100%|█████████████████████████████████████████| Time: 0:00:00





Chains MCMC chain (5000×14×1 Array{Float64, 3}):

Iterations        = 1001:1:6000
Number of chains  = 1
Samples per chain = 5000
Wall duration     = 8.43 seconds
Compute duration  = 8.43 seconds
parameters        = delta, sigma
internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size

Summary Statistics
  parameters      mean       std      mcse    ess_bulk    ess_tail      rhat   ⋯
      Symbol   Float64   Float64   Float64     Float64     Float64   Float64   ⋯

       delta    0.1199    0.1508    0.0021   5048.5898   3389.8513    1.0000   ⋯
       sigma    1.0069    0.1123    0.0016   5440.4428   3822.7463    1.0000   ⋯
                                                                1 column omitted

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5% 
      Symbol   Float64   Float64   Float64   Float64   Float64 

       delta   -0.1707    0.0183    0.1148    0.2244    0.4147
       sigma    0.8106    0.9269    0.9968    1.0770    1.2567

plot(chain, size=(800, 500),
    left_margin=10Plots.mm, bottom_margin=10Plots.mm, fmt=format)

plot(-3:0.1:3, [pdf(Cauchy(0, 1), e) for e in -3:0.1:3], 
    size=(500,300), fmt=format, linewidth=2, label="Prior")
k = kde(chain[:delta] |> vec)
plot!(-3:0.05:3, pdf(k, -3:0.05:3), linewidth=2, ls=:dash, label="Posterior")

plot!([0, 0], [pdf(Cauchy(0, 1), 0), pdf(k, 0)], linewidth=2, label=false)
scatter!([0, 0], [pdf(Cauchy(0, 1), 0), pdf(k, 0)], label=false)

xlabel!("Delta")
ylabel!("Density")
xlims!(-3, 3)

bayes_factor_1 = pdf(k, 0) / pdf(Cauchy(0, 1), 0)

6.028442134709729

8.2 Order-restricted one-sample comparison

\delta \sim \text{Cauchy} (0, 1)_{\mathcal I(-∞,0)}

\sigma \sim \text{Cauchy} (0, 1)_{\mathcal I(0,∞)}

\mu = \delta\sigma

x_{i} \sim \text{Gaussian}(\mu,1/\sigma^2)

winter = [
    -0.05, 0.41, 0.17, -0.13, 0.00, -0.05, 0.00, 0.17, 0.29, 0.04, 0.21, 0.08, 0.37,
    0.17, 0.08, -0.04, -0.04, 0.04, -0.13, -0.12, 0.04, 0.21, 0.17, 0.17, 0.17,
    0.33, 0.04, 0.04, 0.04, 0.00, 0.21, 0.13, 0.25, -0.05, 0.29, 0.42, -0.05, 0.12,
    0.04, 0.25, 0.12
]

summer = [
    0.00, 0.38, -0.12, 0.12, 0.25, 0.12, 0.13, 0.37, 0.00, 0.50, 0.00, 0.00, -0.13,
    -0.37, -0.25, -0.12, 0.50, 0.25, 0.13, 0.25, 0.25, 0.38, 0.25, 0.12, 0.00, 0.00,
    0.00, 0.00, 0.25, 0.13, -0.25, -0.38, -0.13, -0.25, 0.00, 0.00, -0.12, 0.25,
    0.00, 0.50, 0.00
]
x = winter - summer
x = x / std(x);

@model function OrderRestrictedOneSampleComparison(x)
    delta ~ truncated(Cauchy(0, 1), upper=0)
    sigma ~ truncated(Cauchy(0, 1), lower=0)

    mu = delta * sigma

    for i in eachindex(x)
        x[i] ~ Normal(mu, sigma)
    end
end

iterations = 5_000
burnin = 2500

model_order_restricted_one_sample_comparison = OrderRestrictedOneSampleComparison(x)
chain = sample(model_order_restricted_one_sample_comparison, NUTS(), iterations, burnin=burnin)

┌ Info: Found initial step size
└   ϵ = 0.05
Sampling: 100%|█████████████████████████████████████████| Time: 0:00:00





Chains MCMC chain (5000×14×1 Array{Float64, 3}):

Iterations        = 1001:1:6000
Number of chains  = 1
Samples per chain = 5000
Wall duration     = 2.7 seconds
Compute duration  = 2.7 seconds
parameters        = delta, sigma
internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size

Summary Statistics
  parameters      mean       std      mcse    ess_bulk    ess_tail      rhat   ⋯
      Symbol   Float64   Float64   Float64     Float64     Float64   Float64   ⋯

       delta   -0.0882    0.0724    0.0016   1910.6339   1633.7101    1.0010   ⋯
       sigma    1.0177    0.1144    0.0022   2863.8006   2769.0360    0.9998   ⋯
                                                                1 column omitted

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5% 
      Symbol   Float64   Float64   Float64   Float64   Float64 

       delta   -0.2688   -0.1261   -0.0705   -0.0319   -0.0030
       sigma    0.8225    0.9365    1.0082    1.0885    1.2658

plot(chain, size=(800, 500),
    left_margin=10Plots.mm, bottom_margin=10Plots.mm, fmt=format)

plot(-3:0.1:0, [pdf(Cauchy(0, 1), e) for e in -3:0.1:0], 
    size=(500,300), fmt=format, linewidth=2, label="Prior")

k = kde(chain[:delta] |> vec)
plot!(-3:0.01:0, pdf(k, -3:0.01:0), linewidth=2, ls=:dash, label="Posterior")

plot!([0, 0], [pdf(Cauchy(0, 1), 0), pdf(k, 0)], linewidth=2, label=false)
scatter!([0, 0], [pdf(Cauchy(0, 1), 0), pdf(k, 0)], label=false)

xlabel!("Delta")
ylabel!("Density")
xlims!(-2, 0.2)

bayes_factor_2 = pdf(k, 0) / pdf(Cauchy(0, 1), 0)

12.430334525019

8.3 Two-sample comparison

\delta \sim \text{Cauchy} (0, 1)

\mu \sim \text{Cauchy} (0, 1)

\sigma \sim \text{Cauchy} (0, 1)_{\mathcal I(0,∞)}

\alpha = \delta\sigma

x_{i} \sim \text{Gaussian}(\mu+\frac{\alpha}{2},1/\sigma^2)

y_{i} \sim \text{Gaussian}(\mu-\frac{\alpha}{2},1/\sigma^2)

x = [70, 80, 79, 83, 77, 75, 84, 78, 75, 75, 78, 82, 74, 81, 72, 70, 75, 72, 76, 77]
y = [56, 80, 63, 62, 67, 71, 68, 76, 79, 67, 76, 74, 67, 70, 62, 65, 72, 72, 69, 71]

n1 = length(x)
n2 = length(y)

y = y .- mean(x)
y = y ./ std(x)
x = (x .- mean(x)) ./ std(x);

@model function TwoSampleComparison(x, y)
    delta ~ Cauchy(0, 1)
    mu ~ Cauchy(0, 1)
    sigma ~ truncated(Cauchy(0, 1), lower=0)
    alpha = delta * sigma
    
    for i in eachindex(x)
        x[i] ~ Normal(mu + alpha / 2., sigma)
        y[i] ~ Normal(mu - alpha / 2., sigma)
    end
end

iterations = 5_000
burnin = 2500

model_two_sample_comparison = TwoSampleComparison(x, y)
chain = sample(model_two_sample_comparison, NUTS(), iterations, burnin=burnin)

┌ Info: Found initial step size
└   ϵ = 0.2
Sampling: 100%|█████████████████████████████████████████| Time: 0:00:00





Chains MCMC chain (5000×15×1 Array{Float64, 3}):

Iterations        = 1001:1:6000
Number of chains  = 1
Samples per chain = 5000
Wall duration     = 3.55 seconds
Compute duration  = 3.55 seconds
parameters        = delta, mu, sigma
internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size

Summary Statistics
  parameters      mean       std      mcse    ess_bulk    ess_tail      rhat   ⋯
      Symbol   Float64   Float64   Float64     Float64     Float64   Float64   ⋯

       delta    1.3101    0.3555    0.0050   4974.9642   2996.4651    1.0002   ⋯
          mu   -0.8608    0.2047    0.0026   6454.6084   3856.2363    1.0000   ⋯
       sigma    1.3021    0.1505    0.0022   4500.3307   3901.1558    0.9999   ⋯
                                                                1 column omitted

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5% 
      Symbol   Float64   Float64   Float64   Float64   Float64 

       delta    0.6234    1.0708    1.3141    1.5506    2.0071
          mu   -1.2664   -0.9922   -0.8625   -0.7328   -0.4527
       sigma    1.0437    1.2006    1.2892    1.3922    1.6319

plot(chain, size=(800,800), left_margin=10Plots.mm, bottom_margin=10Plots.mm, fmt=format)

plot(-3:0.1:3, [pdf(Cauchy(0, 1), e) for e in -3:0.1:3], 
    size=(500,300), fmt=format, linewidth=2, label="Prior")

k = kde(chain[:delta] |> vec)
plot!(-3:0.01:3, pdf(k, -3:0.01:3), linewidth=2, ls=:dash, label="Posterior")

plot!([0, 0], [pdf(Cauchy(0, 1), 0), pdf(k, 0)], linewidth=2, label=false)
scatter!([0, 0], [pdf(Cauchy(0, 1), 0), pdf(k, 0)], label=false)

xlabel!("Delta")
ylabel!("Density")
xlims!(-3, 3)

bayes_factor_3 = pdf(Cauchy(0, 1), 0) / pdf(k, 0)

503.9036133433892

Pkg.status()

Status `~/quangtiencs_projects/bayesian-cognitive-modeling-with-turing.jl/Project.toml`
  [336ed68f] CSV v0.10.11
  [a93c6f00] DataFrames v1.6.1
⌃ [366bfd00] DynamicPPL v0.23.0
  [7073ff75] IJulia v1.24.2
  [5ab0869b] KernelDensity v0.6.7
  [b964fa9f] LaTeXStrings v1.3.0
  [2ab3a3ac] LogExpFunctions v0.3.24
  [91a5bcdd] Plots v1.38.17
  [37e2e3b7] ReverseDiff v1.15.0
  [276daf66] SpecialFunctions v2.3.0
  [f3b207a7] StatsPlots v0.15.6
  [fce5fe82] Turing v0.28.1
  [e88e6eb3] Zygote v0.6.62
Info Packages marked with ⌃ have new versions available and may be upgradable.