In [1]:
using Mamba
In [2]:
# Data
schools_dat = Dict{Symbol, Any}(
    :J => 8,
    :y => [28,  8, -3,  7, -1,  1, 18, 12],
    :sigma => [15, 10, 16, 11,  9, 11, 10, 18])
Out[2]:
Dict{Symbol,Any} with 3 entries:
  :J     => 8
  :y     => [28, 8, -3, 7, -1, 1, 18, 12]
  :sigma => [15, 10, 16, 11, 9, 11, 10, 18]
In [3]:
# Model
model = Model(
    y = Stochastic(1,
        (theta, sigma)->MvNormal(theta, sigma),
        false
    ),
    theta = Logical(1,(mu, tau, eta)->mu + tau*eta,false),
    eta = Stochastic(1, ()-> MvNormal(8, 1.0)),
    mu  = Stochastic(()-> Normal(0, 100)),
    tau = Stochastic(() -> Rayleigh(100))
)
Out[3]:
Object of type "Mamba.Model"
-------------------------------------------------------------------------------
y:
An unmonitored node of type "0-element Mamba.ArrayStochastic{1}"
Float64[]
-------------------------------------------------------------------------------
eta:
A monitored node of type "0-element Mamba.ArrayStochastic{1}"
Float64[]
-------------------------------------------------------------------------------
mu:
A monitored node of type "Mamba.ScalarStochastic"
NaN
-------------------------------------------------------------------------------
theta:
An unmonitored node of type "0-element Mamba.ArrayLogical{1}"
Float64[]
-------------------------------------------------------------------------------
tau:
A monitored node of type "Mamba.ScalarStochastic"
NaN
In [4]:
# Initial Values
inits = [
    Dict(:y=>schools_dat[:y],
         :eta=>rand(Normal(0,1),8),
         :mu =>rand(Normal(0,1)),
         :tau => rand(Rayleigh(1))
        )for i in 1:4
]
Out[4]:
4-element Array{Dict{Symbol,Any},1}:
 Dict{Symbol,Any}(Pair{Symbol,Any}(:y, [28, 8, -3, 7, -1, 1, 18, 12]),Pair{Symbol,Any}(:eta, [0.468247, 0.466959, 1.7241, 1.63449, 0.466135, 0.575605, -0.100683, -0.160592]),Pair{Symbol,Any}(:mu, 0.867688),Pair{Symbol,Any}(:tau, 2.25625))      
 Dict{Symbol,Any}(Pair{Symbol,Any}(:y, [28, 8, -3, 7, -1, 1, 18, 12]),Pair{Symbol,Any}(:eta, [-1.00713, 1.24793, 0.262724, -0.446687, 1.1828, -0.961785, 0.148227, -0.995843]),Pair{Symbol,Any}(:mu, 0.389412),Pair{Symbol,Any}(:tau, 0.365725))    
 Dict{Symbol,Any}(Pair{Symbol,Any}(:y, [28, 8, -3, 7, -1, 1, 18, 12]),Pair{Symbol,Any}(:eta, [-0.272355, -1.39426, 1.51095, -0.352887, -0.0123941, 1.31394, -0.583011, 0.545557]),Pair{Symbol,Any}(:mu, -0.0108867),Pair{Symbol,Any}(:tau, 1.55082))
 Dict{Symbol,Any}(Pair{Symbol,Any}(:y, [28, 8, -3, 7, -1, 1, 18, 12]),Pair{Symbol,Any}(:eta, [-1.25073, -0.709787, 0.397141, 0.848035, 0.874183, 1.48974, 2.59394, -0.640735]),Pair{Symbol,Any}(:mu, -1.67913),Pair{Symbol,Any}(:tau, 0.561097))    
In [5]:
## Sampling Scheme
scheme = [NUTS([:mu,:eta]),Slice(:tau,3)]
setsamplers!(model, scheme)
Out[5]:
Object of type "Mamba.Model"
-------------------------------------------------------------------------------
y:
An unmonitored node of type "0-element Mamba.ArrayStochastic{1}"
Float64[]
-------------------------------------------------------------------------------
eta:
A monitored node of type "0-element Mamba.ArrayStochastic{1}"
Float64[]
-------------------------------------------------------------------------------
mu:
A monitored node of type "Mamba.ScalarStochastic"
NaN
-------------------------------------------------------------------------------
theta:
An unmonitored node of type "0-element Mamba.ArrayLogical{1}"
Float64[]
-------------------------------------------------------------------------------
tau:
A monitored node of type "Mamba.ScalarStochastic"
NaN
In [6]:
## MCMC Simulations
sim = mcmc(model, schools_dat, inits, 1000,burnin=500,thin=1, chains=4)
describe(sim)
MCMC Simulation of 1000 Iterations x 4 Chains...

Chain 1:   1% [0:05:59 of 0:06:03 remaining]
Chain 1:  10% [0:00:34 of 0:00:38 remaining]
Chain 1:  20% [0:00:16 of 0:00:20 remaining]
Chain 1:  30% [0:00:10 of 0:00:14 remaining]
Chain 1:  40% [0:00:07 of 0:00:11 remaining]
Chain 1:  50% [0:00:05 of 0:00:09 remaining]
Chain 1:  60% [0:00:03 of 0:00:08 remaining]
Chain 1:  70% [0:00:02 of 0:00:07 remaining]
Chain 1:  80% [0:00:01 of 0:00:07 remaining]
Chain 1:  90% [0:00:01 of 0:00:06 remaining]
Chain 1: 100% [0:00:00 of 0:00:06 remaining]

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

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

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

Iterations = 501:1000
Thinning interval = 1
Chains = 1,2,3,4
Samples per chain = 500

Empirical Posterior Estimates:
           Mean         SD       Naive SE      MCSE        ESS   
   tau 11.973184934 7.11380920 0.159069609 1.355458557  27.544313
    mu  5.383857841 8.33140697 0.186295923 1.518361677  30.108300
eta[1]  0.689808179 0.86688345 0.019384103 0.082986222 109.121277
eta[2]  0.053703914 0.72766689 0.016271126 0.052481258 192.245755
eta[3] -0.292673712 0.83395925 0.018647896 0.048122119 300.331012
eta[4]  0.067606021 0.78853122 0.017632094 0.070188465 126.213641
eta[5] -0.222602894 0.89875205 0.020096707 0.125052049  51.653311
eta[6] -0.189599057 0.73919283 0.016528854 0.055794386 175.523185
eta[7]  0.558052649 0.72465191 0.016203709 0.060600547 142.990042
eta[8]  0.109216120 0.82952584 0.018548762 0.043986389 355.650334

Quantiles:
           2.5%        25.0%        50.0%       75.0%       97.5%  
   tau   1.5725862  6.351361252 10.413302141 16.97363892 27.2000793
    mu -10.4120857 -0.067713166  6.905032150 10.84400497 18.8816636
eta[1]  -1.1728667  0.127502592  0.761832699  1.45184594  2.2259377
eta[2]  -1.5722271 -0.330428430  0.165527304  0.38331406  1.5034164
eta[3]  -1.9701499 -0.838346301 -0.226005597  0.22082322  1.4088399
eta[4]  -1.5492714 -0.408090584  0.188031876  0.42874867  1.6458610
eta[5]  -1.9629868 -0.956440458 -0.235636020  0.61545328  1.1176237
eta[6]  -1.7834951 -0.641866111 -0.042730689  0.17197250  1.1906110
eta[7]  -1.0552292  0.125584701  0.651236802  1.03628935  1.8405283
eta[8]  -1.7019302 -0.265362431  0.008813590  0.61733141  1.7965113

In [7]:
p=plot(sim)
draw(p)
Value -2 -1 0 1 2 3 4 0 1 2 3 4 Density eta[1] Iteration 500 600 700 800 900 1000 -4 -2 0 2 4 Value eta[1] Value -20 -10 0 10 20 30 0.00 0.05 0.10 0.15 0.20 Density mu Iteration 500 600 700 800 900 1000 -20 -10 0 10 20 30 Value mu Value -10 0 10 20 30 40 0.00 0.05 0.10 0.15 Density tau Iteration 500 600 700 800 900 1000 0 10 20 30 40 Value tau
Press ENTER to draw next plot
STDIN> 
Value -3 -2 -1 0 1 2 3 0 1 2 3 4 Density eta[4] Iteration 500 600 700 800 900 1000 -3 -2 -1 0 1 2 3 Value eta[4] Value -3 -2 -1 0 1 2 3 0 1 2 3 Density eta[3] Iteration 500 600 700 800 900 1000 -4 -2 0 2 4 Value eta[3] Value -3 -2 -1 0 1 2 3 0 1 2 3 Density eta[2] Iteration 500 600 700 800 900 1000 -3 -2 -1 0 1 2 3 Value eta[2]
Press ENTER to draw next plot
STDIN> 
Value -2 -1 0 1 2 3 0 5 10 15 20 25 Density eta[7] Iteration 500 600 700 800 900 1000 -3 -2 -1 0 1 2 3 Value eta[7] Value -3 -2 -1 0 1 2 3 0 1 2 3 4 Density eta[6] Iteration 500 600 700 800 900 1000 -3 -2 -1 0 1 2 3 Value eta[6] Value -3 -2 -1 0 1 2 0.0 0.5 1.0 1.5 2.0 Density eta[5] Iteration 500 600 700 800 900 1000 -3 -2 -1 0 1 2 3 Value eta[5]
Press ENTER to draw next plot
STDIN> 
Value -3 -2 -1 0 1 2 3 0 1 2 3 4 Density eta[8] Iteration 500 600 700 800 900 1000 -3 -2 -1 0 1 2 3 Value eta[8]