#using Distributed
#@everywhere using MambaModels
using MambaModels
gr(size=(400,400))

# Data
line = Dict{Symbol, Any}()

howell1 = CSV.read(rel_path("..", "data", "Howell1.csv"), delim=';')
df = convert(DataFrame, howell1);

544 rows × 4 columns

Use only adults

df2 = filter(row -> row[:age] >= 18, df);
mean_weight = mean(df2[:weight])
df2[:weight_c] = convert(Vector{Float64}, df2[:weight]) .- mean_weight ;
line[:x] = convert(Array{Float64,1}, df2[:weight_c]);
line[:y] = convert(Array{Float64,1}, df2[:height]);
line[:xmat] = convert(Array{Float64,2}, [ones(length(line[:x])) line[:x]])

352×2 Array{Float64,2}:
 1.0    2.83512
 1.0   -8.50468
 1.0  -13.1256 
 1.0    8.05143
 1.0   -3.71361
 1.0   18.0021 
 1.0   -6.74701
 1.0   10.4895 
 1.0  -10.1206 
 1.0    9.49725
 ⋮             
 1.0    2.89182
 1.0   -5.58468
 1.0   -3.94041
 1.0   -4.16721
 1.0    2.04133
 1.0  -10.7443 
 1.0    7.17259
 1.0    9.07201
 1.0    7.54114

Model Specification

model = Model(
  y = Stochastic(1,
    (xmat, beta, s2) -> MvNormal(xmat * beta, sqrt(s2)),
    false
  ),
  beta = Stochastic(1, () -> MvNormal([178, 0], [sqrt(10000), sqrt(100)])),
  s2 = Stochastic(() -> Uniform(0, 50))
)

Object of type "Model"
-------------------------------------------------------------------------------
beta:
A monitored node of type "0-element ArrayStochastic{1}"
Float64[]
-------------------------------------------------------------------------------
y:
An unmonitored node of type "0-element ArrayStochastic{1}"
Float64[]
-------------------------------------------------------------------------------
s2:
A monitored node of type "ScalarStochastic"
NaN

Initial Values

inits = [
  Dict{Symbol, Any}(
    :y => line[:y],
    :beta => [rand(Normal(178, 100)), rand(Normal(0, 10))],
    :s2 => rand(Uniform(0, 50))
  )
  for i in 1:3
]

3-element Array{Dict{Symbol,Any},1}:
 Dict(:beta=>[200.815, 3.68352],:y=>[151.765, 139.7, 136.525, 156.845, 145.415, 163.83, 149.225, 168.91, 147.955, 165.1  …  156.21, 160.655, 146.05, 156.21, 152.4, 162.56, 142.875, 162.56, 156.21, 158.75],:s2=>30.1041)  
 Dict(:beta=>[203.827, -0.301006],:y=>[151.765, 139.7, 136.525, 156.845, 145.415, 163.83, 149.225, 168.91, 147.955, 165.1  …  156.21, 160.655, 146.05, 156.21, 152.4, 162.56, 142.875, 162.56, 156.21, 158.75],:s2=>49.33)  
 Dict(:beta=>[295.371, -0.447918],:y=>[151.765, 139.7, 136.525, 156.845, 145.415, 163.83, 149.225, 168.91, 147.955, 165.1  …  156.21, 160.655, 146.05, 156.21, 152.4, 162.56, 142.875, 162.56, 156.21, 158.75],:s2=>2.68406)

Tuning Parameters

scale1 = [0.5, 0.25]
summary1 = identity
eps1 = 0.5

scale2 = 0.5
summary2 = x -> [mean(x); sqrt(var(x))]
eps2 = 0.1

0.1

Define sampling scheme

scheme = [
  Mamba.NUTS([:beta]),
  Mamba.Slice([:s2], 10)
]

setsamplers!(model, scheme)

Object of type "Model"
-------------------------------------------------------------------------------
beta:
A monitored node of type "0-element ArrayStochastic{1}"
Float64[]
-------------------------------------------------------------------------------
y:
An unmonitored node of type "0-element ArrayStochastic{1}"
Float64[]
-------------------------------------------------------------------------------
s2:
A monitored node of type "ScalarStochastic"
NaN

MCMC Simulation

chn = mcmc(model, line, inits, 10000, burnin=1000, chains=3)

Object of type "ModelChains"

Iterations = 1001:10000
Thinning interval = 1
Chains = 1,2,3
Samples per chain = 9000

[25.4187 154.709 0.920336; 26.6117 154.77 0.98207; … ; 27.1471 154.283 0.867431; 25.4818 154.287 0.889286]

[25.374 154.785 0.893591; 25.7082 154.785 0.893591; … ; 26.8278 154.646 0.9003; 25.9587 154.446 0.900654]

[24.7338 154.574 0.900931; 25.2447 154.388 0.901082; … ; 25.9001 154.872 0.909671; 25.7867 154.872 0.909671]

Show draws summary

describe(chn)

Iterations = 1001:10000
Thinning interval = 1
Chains = 1,2,3
Samples per chain = 9000

Empirical Posterior Estimates:
            Mean          SD        Naive SE        MCSE        ESS
     s2  26.16130736 1.989228253 0.01210605763 0.01903204080 9000.0000
beta[1] 154.59957294 0.273740509 0.00166593169 0.00389110458 4949.1696
beta[2]   0.90468999 0.042338635 0.00025766473 0.00032517239 9000.0000

Quantiles:
            2.5%       25.0%        50.0%        75.0%        97.5%
     s2  22.5096692  24.7702763  26.07798872  27.44804314  30.29630489
beta[1] 154.0618462 154.4133301 154.60265035 154.78241808 155.13638960
beta[2]   0.8215179   0.8762166   0.90496442   0.93272244   0.98765192

Convert to MCMCChains.Chains object

chn2 = MCMCChains.Chains(chn.value, String.(chn.names))

Object of type Chains, with data of type 9000×3×3 Array{Float64,3}

Iterations        = 1:9000
Thinning interval = 1
Chains            = 1, 2, 3
Samples per chain = 9000
parameters        = beta[1], beta[2], s2

2-element Array{ChainDataFrame,1}

Summary Statistics
. Omitted printing of 1 columns
│ Row │ parameters │ mean    │ std       │ naive_se    │ mcse        │ ess     │
│     │ Symbol     │ Float64 │ Float64   │ Float64     │ Float64     │ Any     │
├─────┼────────────┼─────────┼───────────┼─────────────┼─────────────┼─────────┤
│ 1   │ beta[1]    │ 154.6   │ 0.273741  │ 0.00166593  │ 0.0038911   │ 5333.08 │
│ 2   │ beta[2]    │ 0.90469 │ 0.0423386 │ 0.000257665 │ 0.000325172 │ 17433.4 │
│ 3   │ s2         │ 26.1613 │ 1.98923   │ 0.0121061   │ 0.019032    │ 10857.3 │

Quantiles

│ Row │ parameters │ 2.5%     │ 25.0%    │ 50.0%    │ 75.0%    │ 97.5%    │
│     │ Symbol     │ Float64  │ Float64  │ Float64  │ Float64  │ Float64  │
├─────┼────────────┼──────────┼──────────┼──────────┼──────────┼──────────┤
│ 1   │ beta[1]    │ 154.062  │ 154.413  │ 154.603  │ 154.782  │ 155.136  │
│ 2   │ beta[2]    │ 0.821518 │ 0.876217 │ 0.904964 │ 0.932722 │ 0.987652 │
│ 3   │ s2         │ 22.5097  │ 24.7703  │ 26.078   │ 27.448   │ 30.2963  │

Describe the MCMCChains

MCMCChains.describe(chn2)

2-element Array{ChainDataFrame,1}

Summary Statistics
. Omitted printing of 1 columns
│ Row │ parameters │ mean    │ std       │ naive_se    │ mcse        │ ess     │
│     │ Symbol     │ Float64 │ Float64   │ Float64     │ Float64     │ Any     │
├─────┼────────────┼─────────┼───────────┼─────────────┼─────────────┼─────────┤
│ 1   │ beta[1]    │ 154.6   │ 0.273741  │ 0.00166593  │ 0.0038911   │ 5333.08 │
│ 2   │ beta[2]    │ 0.90469 │ 0.0423386 │ 0.000257665 │ 0.000325172 │ 17433.4 │
│ 3   │ s2         │ 26.1613 │ 1.98923   │ 0.0121061   │ 0.019032    │ 10857.3 │

Quantiles

│ Row │ parameters │ 2.5%     │ 25.0%    │ 50.0%    │ 75.0%    │ 97.5%    │
│     │ Symbol     │ Float64  │ Float64  │ Float64  │ Float64  │ Float64  │
├─────┼────────────┼──────────┼──────────┼──────────┼──────────┼──────────┤
│ 1   │ beta[1]    │ 154.062  │ 154.413  │ 154.603  │ 154.782  │ 155.136  │
│ 2   │ beta[2]    │ 0.821518 │ 0.876217 │ 0.904964 │ 0.932722 │ 0.987652 │
│ 3   │ s2         │ 22.5097  │ 24.7703  │ 26.078   │ 27.448   │ 30.2963  │

Plot chn2

MCMCChains.plot(chn2)

End of 04/m4.1m.jl

This page was generated using Literate.jl.