m4.4m · MambaModels.jl

#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

	height	weight	age	male
	Float64⍰	Float64⍰	Float64⍰	Int64⍰
1	151.765	47.8256	63.0	1
2	139.7	36.4858	63.0	0
3	136.525	31.8648	65.0	0
4	156.845	53.0419	41.0	1
5	145.415	41.2769	51.0	0
6	163.83	62.9926	35.0	1
7	149.225	38.2435	32.0	0
8	168.91	55.48	27.0	1
9	147.955	34.8699	19.0	0
10	165.1	54.4877	54.0	1
11	154.305	49.8951	47.0	0
12	151.13	41.2202	66.0	1
13	144.78	36.0322	73.0	0
14	149.9	47.7	20.0	0
15	150.495	33.8493	65.3	0
16	163.195	48.5627	36.0	1
17	157.48	42.3258	44.0	1
18	143.942	38.3569	31.0	0
19	121.92	19.6179	12.0	1
20	105.41	13.948	8.0	0
21	86.36	10.4893	6.5	0
22	161.29	48.9879	39.0	1
23	156.21	42.7227	29.0	0
24	129.54	23.5868	13.0	1
⋮	⋮	⋮	⋮	⋮

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=>[234.716, 1.74853],: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=>23.6892)
 Dict(:beta=>[274.993, -10.9435],: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=>35.776)
 Dict(:beta=>[195.061, 7.02901],: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=>38.2237)

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)

MCMC Simulation of 10000 Iterations x 3 Chains...

Chain 1:   0% [0:45:36 of 0:45:39 remaining]
Chain 1:  10% [0:00:29 of 0:00:32 remaining]
Chain 1:  20% [0:00:16 of 0:00:20 remaining]
Chain 1:  30% [0:00:10 of 0:00:15 remaining]
Chain 1:  40% [0:00:07 of 0:00:12 remaining]
Chain 1:  50% [0:00:05 of 0:00:11 remaining]
Chain 1:  60% [0:00:04 of 0:00:10 remaining]
Chain 1:  70% [0:00:03 of 0:00:09 remaining]
Chain 1:  80% [0:00:02 of 0:00:08 remaining]
Chain 1:  90% [0:00:01 of 0:00:08 remaining]
Chain 1: 100% [0:00:00 of 0:00:08 remaining]

Chain 2:   0% [0:00:03 of 0:00:03 remaining]
Chain 2:  10% [0:00:03 of 0:00:04 remaining]
Chain 2:  20% [0:00:03 of 0:00:04 remaining]
Chain 2:  30% [0:00:03 of 0:00:04 remaining]
Chain 2:  40% [0:00:03 of 0:00:04 remaining]
Chain 2:  50% [0:00:02 of 0:00:04 remaining]
Chain 2:  60% [0:00:02 of 0:00:04 remaining]
Chain 2:  70% [0:00:01 of 0:00:04 remaining]
Chain 2:  80% [0:00:01 of 0:00:04 remaining]
Chain 2:  90% [0:00:00 of 0:00:04 remaining]
Chain 2: 100% [0:00:00 of 0:00:04 remaining]

Chain 3:   0% [0:00:02 of 0:00:02 remaining]
Chain 3:  10% [0:00:03 of 0:00:04 remaining]
Chain 3:  20% [0:00:03 of 0:00:04 remaining]
Chain 3:  30% [0:00:03 of 0:00:04 remaining]
Chain 3:  40% [0:00:02 of 0:00:04 remaining]
Chain 3:  50% [0:00:02 of 0:00:04 remaining]
Chain 3:  60% [0:00:02 of 0:00:04 remaining]
Chain 3:  70% [0:00:01 of 0:00:04 remaining]
Chain 3:  80% [0:00:01 of 0:00:04 remaining]
Chain 3:  90% [0:00:00 of 0:00:04 remaining]
Chain 3: 100% [0:00:00 of 0:00:04 remaining]

Object of type "ModelChains"

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

[24.5338 154.186 0.896712; 29.2975 154.186 0.896712; … ; 28.2548 154.422 0.858305; 28.0598 154.954 0.964935]

[25.9371 154.413 0.875213; 26.4889 154.546 0.89657; … ; 24.0353 154.949 0.889058; 24.2107 154.489 0.898575]

[27.9343 154.776 0.831786; 28.6161 154.335 0.960829; … ; 23.2665 154.636 0.884034; 24.3545 154.406 0.914905]

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.15372585 1.992994994 0.01212898128 0.02070715765 9000.000
beta[1] 154.59176159 0.273999633 0.00166750866 0.00370996861 5454.557
beta[2]   0.90505901 0.042259614 0.00025718382 0.00030447388 9000.000

Quantiles:
            2.5%        25.0%        50.0%        75.0%       97.5%
     s2  22.5331015  24.76312186  26.07240120  27.43844975  30.3452136
beta[1] 154.0449654 154.40813312 154.59342969 154.77637257 155.1237737
beta[2]   0.8216441   0.87691367   0.90513733   0.93351717   0.9878548

Convert to MCMCChains.Chains object

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

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

Log evidence      = 0.0
Iterations        = 1:9000
Thinning interval = 1
Chains            = Chain1, Chain2, Chain3
Samples per chain = 9000
parameters        = s2, beta[1], beta[2]

parameters
          Mean     SD   Naive SE  MCSE     ESS
beta[1] 154.5918 0.2740   0.0017 0.0037 5454.557
beta[2]   0.9051 0.0423   0.0003 0.0003 9000.000
     s2  26.1537 1.9930   0.0121 0.0207 9000.000

Describe the MCMCChains

MCMCChains.describe(chn2)

Log evidence      = 0.0
Iterations        = 1:9000
Thinning interval = 1
Chains            = Chain1, Chain2, Chain3
Samples per chain = 9000
parameters        = s2, beta[1], beta[2]

Empirical Posterior Estimates:
====================================================
parameters
          Mean     SD   Naive SE  MCSE     ESS
beta[1] 154.5918 0.2740   0.0017 0.0037 5454.557
beta[2]   0.9051 0.0423   0.0003 0.0003 9000.000
     s2  26.1537 1.9930   0.0121 0.0207 9000.000

Quantiles:
====================================================
parameters
          2.5%     25.0%    50.0%    75.0%    97.5%
beta[1] 153.3728 154.4081 154.5934 154.7764 155.5619
beta[2]   0.7348   0.8769   0.9051   0.9335   1.0686
     s2  19.5225  24.7631  26.0724  27.4384  36.2415

Plot chn2

MCMCChains.plot(chn2)

End of 04/m4.1m.jl

This page was generated using Literate.jl.