``````

In [1]:

``````

Model and user-defined sampler specification

``````

In [2]:

model = Model(

y = Stochastic(1,
(mu, s2) ->  MvNormal(mu, sqrt(s2)),
false
),

mu = Logical(1,
(xmat, beta) -> xmat * beta,
false
),

beta = Stochastic(1,
() -> MvNormal(2, sqrt(1000))
),

s2 = Stochastic(
() -> InverseGamma(0.001, 0.001)
)

)

Gibbs_beta = Sampler([:beta],
(beta, s2, xmat, y) ->
begin
beta_mean = mean(beta.distr)
beta_invcov = invcov(beta.distr)
Sigma = inv(xmat' * xmat / s2 + beta_invcov)
mu = Sigma * (xmat' * y / s2 + beta_invcov * beta_mean)
rand(MvNormal(mu, Sigma))
end
)

Gibbs_s2 = Sampler([:s2],
(mu, s2, y) ->
begin
a = length(y) / 2.0 + shape(s2.distr)
b = sumabs2(y - mu) / 2.0 + scale(s2.distr)
rand(InverseGamma(a, b))
end
)

``````
``````

Out[2]:

An object of type "Mamba.Sampler{Dict{Any,Any}}"
Sampling Block Nodes:
[:s2]

AST(:(\$(Expr(:lambda, Any[:(model::Model),:(block::Integer)], Any[Any[Any[:model,:Any,18],Any[:block,:Any,18],Any[:f,:Any,18]],Any[],0,Any[]], :(begin
model = (top(typeassert))(model,Mamba.Model)
block = (top(typeassert))(block,Mamba.Integer)
f = (anonymous function)
return f((Mamba.getindex)(model,:mu),(Mamba.getindex)(model,:s2),(Mamba.getindex)(model,:y))
end)))))

``````
``````

In [3]:

scheme1 = [NUTS(:beta),
Slice(:s2, 3.0)]

## No-U-Turn Sampling Scheme
scheme2 = [NUTS([:beta, :s2])]

## User-Defined Sampling Scheme
scheme3 = [Gibbs_beta, Gibbs_s2]

## Sampling Scheme Assignment
setsamplers!(model, scheme1)

``````
``````

Out[3]:

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

``````
``````

In [5]:

draw(model)

``````
``````

digraph MambaModel {
"y" [shape="ellipse", fillcolor="gray85", style="filled"];
"mu" [shape="diamond", fillcolor="gray85", style="filled"];
"mu" -> "y";
"s2" [shape="ellipse"];
"s2" -> "y";
"beta" [shape="ellipse"];
"beta" -> "mu";
"xmat" [shape="box", fillcolor="gray85", style="filled"];
"xmat" -> "mu";
}

``````
``````

In [7]:

line = Dict{Symbol, Any}(
:x => [1, 2, 3, 4, 5],
:y => [1, 3, 3, 3, 5]
)
line[:xmat] = [ones(5) line[:x]]

``````
``````

Out[7]:

5x2 Array{Float64,2}:
1.0  1.0
1.0  2.0
1.0  3.0
1.0  4.0
1.0  5.0

``````
``````

In [8]:

srand(123)

## Initial Values
inits = [
Dict{Symbol, Any}(
:y => line[:y],
:beta => rand(Normal(0, 1), 2),
:s2 => rand(Gamma(1, 1))
)
for i in 1:3
]

``````
``````

Out[8]:

3-element Array{Dict{Symbol,Any},1}:
Dict{Symbol,Any}(:y=>[1,3,3,3,5],:s2=>1.1830434911443408,:beta=>[1.1902678809862768,2.04817970778924])
Dict{Symbol,Any}(:y=>[1,3,3,3,5],:s2=>0.3057490088250573,:beta=>[1.142650902867199,0.45941562040708034])
Dict{Symbol,Any}(:y=>[1,3,3,3,5],:s2=>0.6108061237629374,:beta=>[-0.396679079295223,-0.6647125451916877])

``````
``````

In [9]:

sim1 = mcmc(model, line, inits, 10000, burnin=250, thin=2, chains=3);

``````
``````

MCMC Simulation of 10000 Iterations x 3 Chains...

Chain 1:   0% [0:31:50 of 0:31:52 remaining]
Chain 1:  10% [0:00:24 of 0:00:27 remaining]
Chain 1:  20% [0:00:13 of 0:00:16 remaining]
Chain 1:  30% [0:00:09 of 0:00:12 remaining]
Chain 1:  40% [0:00:06 of 0:00:10 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:08 remaining]
Chain 1:  80% [0:00:01 of 0:00:07 remaining]
Chain 1:  90% [0:00:01 of 0:00:07 remaining]
Chain 1: 100% [0:00:00 of 0:00:07 remaining]

Chain 2:   0% [0:00:05 of 0:00:05 remaining]
Chain 2:  10% [0:00:04 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:02 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:05 remaining]
Chain 2: 100% [0:00:00 of 0:00:05 remaining]

Chain 3:   0% [0:00:04 of 0:00:04 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:02 of 0:00:03 remaining]
Chain 3:  40% [0:00:02 of 0:00:03 remaining]
Chain 3:  50% [0:00:02 of 0:00:03 remaining]
Chain 3:  60% [0:00:01 of 0:00:03 remaining]
Chain 3:  70% [0:00:01 of 0:00:03 remaining]
Chain 3:  80% [0:00:01 of 0:00:03 remaining]
Chain 3:  90% [0:00:00 of 0:00:03 remaining]
Chain 3: 100% [0:00:00 of 0:00:03 remaining]

``````
``````

In [10]:

gewekediag(sim1) |> showall

``````
``````

Iterations = 252:10000
Thinning interval = 2
Chains = 1,2,3
Samples per chain = 4875

Geweke Diagnostic:
First Window Fraction = 0.1
Second Window Fraction = 0.5

Z-score p-value
beta[1]   1.237  0.2162
beta[2]  -1.568  0.1168
s2   1.710  0.0872

Z-score p-value
beta[1]  -1.457  0.1452
beta[2]   1.752  0.0797
s2  -1.428  0.1534

Z-score p-value
beta[1]   0.550  0.5824
beta[2]  -0.440  0.6597
s2   0.583  0.5596

``````
``````

In [11]:

describe(sim1)

``````
``````

Iterations = 252:10000
Thinning interval = 2
Chains = 1,2,3
Samples per chain = 4875

Empirical Posterior Estimates:
Mean       SD       Naive SE       MCSE       ESS
beta[1] 0.5971183 1.14894446 0.0095006014 0.016925598 4607.9743
beta[2] 0.8017036 0.34632566 0.0028637608 0.004793345 4875.0000
s2 1.2203777 2.00876760 0.0166104638 0.101798287  389.3843

Quantiles:
2.5%       25.0%       50.0%     75.0%     97.5%
beta[1] -1.74343373 0.026573102 0.59122696 1.1878720 2.8308472
beta[2]  0.12168742 0.628297573 0.80357822 0.9719441 1.5051573
s2  0.17091385 0.383671702 0.65371989 1.2206381 6.0313970

``````
``````

In [ ]:

show(plot(sim1))

``````
``````