In [1]:

    

using Mamba, Gadfly


    

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))


    

    




[Gadfly.Plot([Gadfly.Layer(nothing,Dict{Any,Any}(),Gadfly.StatisticElement[],Gadfly.Geom.LineGeometry(Gadfly.Stat.Identity(),false,2,symbol("")),nothing,0)],nothing,Data(
  x=[