search

brian-j-smith / Mamba.jl

Markov chain Monte Carlo (MCMC) for Bayesian analysis in julia
Other
253 stars 52 forks source link

Incompatible Initial Value? #141

Open jordiabante opened 6 years ago

jordiabante commented 6 years ago

Following an example on the tutorial page, I came up with the code below. 

## Dependencies
using Mamba

## Data
obs = Dict{Symbol, Any}(
    :m1 => 10,
    :m2 => 7,
    :m3 => 3,
    :m4 => 2,
    :m5 => 1
)

## Model
model = Model(
    # alpha's are the same in all cases (known to be small)
    alpha = Stochastic(
        () -> Beta(1,10)
    ),
    # beta's are specific to each point (no previous knowledge)
    beta1 = Stochastic(
        () -> Beta(1,1)
    ),
    beta2 = Stochastic(
        () -> Beta(1,1)
    ),
    beta3 = Stochastic(
        () -> Beta(1,1)
    ),
    beta4 = Stochastic(
        () -> Beta(1,1)
    ),
    # Observations are binomial samples from corresponding hidden N
    m1 = Stochastic(1,
        (n1,alpha) -> Binomial(n1,alpha),
        false
    ),
    m2 = Stochastic(1,
        (n2,alpha) -> Binomial(n2,alpha),
        false
    ),
    m3 = Stochastic(1,
        (n3,alpha) -> Binomial(n3,alpha),
        false
    ),
    m4 = Stochastic(1,
        (n4,alpha) -> Binomial(n4,alpha),
        false
    ),
    m5 = Stochastic(1,
        (n5,alpha) -> Binomial(n5,alpha),
        false
    ),
    # Hidden N's are binomial samples from previous N
    n1 = Stochastic(
        () -> Poisson(100)
    ),
    n2 = Stochastic(1,
        (n1,beta1) -> Binomial(n1,beta1)
    ),
    n3 = Stochastic(1,
        (n2,beta2) -> Binomial(n2,beta2)
    ),
    n4 = Stochastic(1,
        (n3,beta3) -> Binomial(n3,beta3)
    ),
    n5 = Stochastic(1,
        (n4,beta4) -> Binomial(n4,beta4)
    )
)

## Initial Values
inits = [
  Dict(:m1 => obs[:m1], :m2 => obs[:m2], :m3 => obs[:m3], :m4 => obs[:m4],
       :m5 => obs[:m5], :n1 => 100, :n2 => 70, :n3 => 50, :n4 => 40,
       :n5 => 30, :alpha => 0.1, :beta1 => 0.5, :beta2 => 0.5, :beta3 => 0.5,
       :beta4 => 0.5)
]

## Sampling Scheme
scheme = [NUTS([:n1,:n2,:n3,:n4,:n5,:alpha,:beta1,:beta2,:beta3,:beta4])]
setsamplers!(model, scheme)

## MCMC Simulations
sim = mcmc(model, obs, inits, 1000, burnin=250, thin=2, chains=1)

and I get the following error:

ERROR: ArgumentError: incompatible initial value for node : n2
Stacktrace:
 [1] setinits!(::Mamba.ArrayStochastic{1}, ::Mamba.Model, ::Int64) at /Applications/JuliaPro-0.6.3.1.app/Contents/Resources/pkgs-0.6.3.1/v0.6/Mamba/src/model/dependent.jl:177
 [2] setinits!(::Mamba.Model, ::Dict{Symbol,Any}) at /Applications/JuliaPro-0.6.3.1.app/Contents/Resources/pkgs-0.6.3.1/v0.6/Mamba/src/model/initialization.jl:11
 [3] setinits!(::Mamba.Model, ::Array{Dict{Symbol,Any},1}) at /Applications/JuliaPro-0.6.3.1.app/Contents/Resources/pkgs-0.6.3.1/v0.6/Mamba/src/model/initialization.jl:24
 [4] #mcmc#29(::Int64, ::Int64, ::Int64, ::Bool, ::Function, ::Mamba.Model, ::Dict{Symbol,Any}, ::Array{Dict{Symbol,Any},1}, ::Int64) at /Applications/JuliaPro-0.6.3.1.app/Contents/Resources/pkgs-0.6.3
.1/v0.6/Mamba/src/model/mcmc.jl:29
 [5] (::Mamba.#kw##mcmc)(::Array{Any,1}, ::Mamba.#mcmc, ::Mamba.Model, ::Dict{Symbol,Any}, ::Array{Dict{Symbol,Any},1}, ::Int64) at ./<missing>:0
 [6] macro expansion at /Applications/JuliaPro-0.6.3.1.app/Contents/Resources/pkgs-0.6.3.1/v0.6/Atom/src/repl.jl:118 [inlined]
 [7] anonymous at ./<missing>:?

I tried tracing it back, but since I am not that familiar yet with the package, I didn't find what's wrong with it.

Thanks!

brian-j-smith commented 6 years ago

A couple tips. First, model nodes that are scalar, which all seem to be in this example, should be specified with calls to Stochastic or Logical without a dimension as the first argument. The alpha, betas, and n1 are specified correctly in this regard; n2-n5 and m1-m5 are not. The dimension argument sets the dimension for nodes that are arrays of model parameters. Second, the NUTS sampler is appropriate for sampling continuous nodes (alpha and betas) but not for discrete ones (n1-n5). The DGS sampler can be used for discrete nodes with finite support (n2-n5) but not for a discrete node with infinite support (n1 ~ Poisson(100)). For the latter, a continuous distribution could be used as an approximation as shown in the modified code below.

Hope that helps.

## Dependencies
using Mamba

## Data
obs = Dict{Symbol, Any}(
    :m1 => 10,
    :m2 => 7,
    :m3 => 3,
    :m4 => 2,
    :m5 => 1
)

## Model
model = Model(
    # alpha's are the same in all cases (known to be small)
    alpha = Stochastic(
        () -> Beta(1,10)
    ),
    # beta's are specific to each point (no previous knowledge)
    beta1 = Stochastic(
        () -> Beta(1,1)
    ),
    beta2 = Stochastic(
        () -> Beta(1,1)
    ),
    beta3 = Stochastic(
        () -> Beta(1,1)
    ),
    beta4 = Stochastic(
        () -> Beta(1,1)
    ),
    # Observations are binomial samples from corresponding hidden N
    m1 = Stochastic(
        (n1,alpha) -> Binomial(n1,alpha),
        false
    ),
    m2 = Stochastic(
        (n2,alpha) -> Binomial(n2,alpha),
        false
    ),
    m3 = Stochastic(
        (n3,alpha) -> Binomial(n3,alpha),
        false
    ),
    m4 = Stochastic(
        (n4,alpha) -> Binomial(n4,alpha),
        false
    ),
    m5 = Stochastic(
        (n5,alpha) -> Binomial(n5,alpha),
        false
    ),
    # Hidden N's are binomial samples from previous N
    n1_real = Stochastic(
        () -> Gamma(100, 1)
    ),
    n1 = Logical(
        n1_real -> round(Int, n1_real)
    ),
    n2 = Stochastic(
        (n1,beta1) -> Binomial(n1,beta1)
    ),
    n3 = Stochastic(
        (n2,beta2) -> Binomial(n2,beta2)
    ),
    n4 = Stochastic(
        (n3,beta3) -> Binomial(n3,beta3)
    ),
    n5 = Stochastic(
        (n4,beta4) -> Binomial(n4,beta4)
    )
)

## Initial Values
inits = [
  Dict(:m1 => obs[:m1], :m2 => obs[:m2], :m3 => obs[:m3], :m4 => obs[:m4],
       :m5 => obs[:m5], :n1_real => 100, :n2 => 70, :n3 => 50, :n4 => 40,
       :n5 => 30, :alpha => 0.1, :beta1 => 0.5, :beta2 => 0.5, :beta3 => 0.5,
       :beta4 => 0.5)
]

## Sampling Scheme
scheme = [NUTS(:n1_real), DGS([:n2, :n3, :n4, :n5]),
          NUTS([:alpha, :beta1, :beta2, :beta3, :beta4])]
setsamplers!(model, scheme)

## MCMC Simulations
sim = mcmc(model, obs, inits, 1000, burnin=250, thin=2, chains=1)
jordiabante commented 6 years ago

Thank you so much for your feedback, this was very helpful.