using Mamba, Gadfly
Model and user-defined sampler specification
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
)
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)))))
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)
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[]
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"; }
line = Dict{Symbol, Any}(
:x => [1, 2, 3, 4, 5],
:y => [1, 3, 3, 3, 5]
)
line[:xmat] = [ones(5) line[:x]]
5x2 Array{Float64,2}: 1.0 1.0 1.0 2.0 1.0 3.0 1.0 4.0 1.0 5.0
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
]
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])
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]
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
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
show(plot(sim1))
[Gadfly.Plot([Gadfly.Layer(nothing,Dict{Any,Any}(),Gadfly.StatisticElement[],Gadfly.Geom.LineGeometry(Gadfly.Stat.Identity(),false,2,symbol("")),nothing,0)],nothing,Data( x=[