Mamba.jl icon indicating copy to clipboard operation
Mamba.jl copied to clipboard

Incompatible Initial Value?

Open jordiabante opened this issue 6 years ago • 2 comments

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!

jordiabante avatar Jul 06 '18 13:07 jordiabante

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)

brian-j-smith avatar Jul 06 '18 15:07 brian-j-smith

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

jordiabante avatar Jul 06 '18 23:07 jordiabante