4

I am trying to learn how to use the Mamba package in Julia for doing Bayesian inference. Though the package is great, as a beginner I find the documentation a bit scarce in information. Hence I am trying to figure out how implement some very simple examples.

What I have tried

I implemented an example for doing Bayesian inference for the mean of a univariate normal distribution. The code follows:

using Mamba

## Model Specification

model = Model(

  x = Stochastic(1,
    mu -> Normal(mu, 2.0),
    false
  ),

  mu = Stochastic(
    () -> Normal(0.0, 1000.0),
    true
  )

)

## Data
data = Dict{Symbol, Any}(
  :x => randn(30)*2+13
)

## Initial Values
inits = [
  Dict{Symbol, Any}(
    :x => data[:x],
    :mu => randn()*1
  )
]

## Sampling Scheme Assignment
scheme1 = NUTS([:mu])
setsamplers!(model, [scheme1])

sim1 = mcmc(model, data, inits, 10000, burnin=250, thin=2, chains=1);
describe(sim1)

This seems to work absolutely fine (though there may be better ways perhaps of coding this?).

What I am trying to do and doesn't work.

In this example, I am trying to do bayesian inference for the mean of a bivariate normal distribution. The code follows:

using Mamba

## Model Specification

model = Model(

  x = Stochastic(1,
    mu -> MvNormal(mu, eye(2)),
    false
  ),

  mu = Stochastic(1,
    () -> MvNormal(zeros(2), 1000.0),
    true
  )

)

## Data
data = Dict{Symbol, Any}(
  :x => randn(2,30)+13
)

## Initial Values
inits = [
  Dict{Symbol, Any}(
    :x => data[:x],
    :mu => randn(2)*1
  )
]

## Sampling Scheme Assignment
scheme1 = NUTS([:mu])
setsamplers!(model, [scheme1])

sim1 = mcmc(model, data, inits, 10000, burnin=250, thin=2, chains=1);
describe(sim1)

As you may notice, the changes that I suppose are necessary are minimal. However, I am doing somewhere something wrong and when I attempt to run this I get an error (a conversion between types error) which does not help me further.

Any help appreciated. If this works out, I will consider contribute this simple example to the Mamba documentation for other new users. Thanks.

Addendum: the error message

ERROR: MethodError: Cannot `convert` an object of type Array{Float64,2} to an object of type Array{Float64,1}
This may have arisen from a call to the constructor Array{Float64,1}(...),
since type constructors fall back to convert methods.
 in setinits!(::Mamba.ArrayStochastic{1}, ::Mamba.Model, ::Array{Float64,2}) at /lhome/lgiannins/.julia/v0.5/Mamba/src/model/dependent.jl:164
 in setinits!(::Mamba.Model, ::Dict{Symbol,Any}) at /lhome/lgiannins/.julia/v0.5/Mamba/src/model/initialization.jl:11
 in setinits!(::Mamba.Model, ::Array{Dict{Symbol,Any},1}) at /lhome/lgiannins/.julia/v0.5/Mamba/src/model/initialization.jl:24
 in #mcmc#29(::Int64, ::Int64, ::Int64, ::Bool, ::Function, ::Mamba.Model, ::Dict{Symbol,Any}, ::Array{Dict{Symbol,Any},1}, ::Int64) at /lhome/lgiannins/.julia/v0.5/Mamba/src/model/mcmc.jl:29
 in (::Mamba.#kw##mcmc)(::Array{Any,1}, ::Mamba.#mcmc, ::Mamba.Model, ::Dict{Symbol,Any}, ::Array{Dict{Symbol,Any},1}, ::Int64) at ./<missing>:0
user1438310
  • 787
  • 4
  • 14
  • Please add the traceback (or at least the error message) so that someone can help you find what is going on. – Alexander Morley Sep 24 '16 at 11:03
  • good point. Updated post above. – user1438310 Sep 26 '16 at 19:23
  • given the traceback I would guess that one of your entries has the wrong dimensions? – Alexander Morley Sep 26 '16 at 21:57
  • I tried of course changing the declaration "Stochastic(1,..." to "Stochastic(2,..." but to no avail. – user1438310 Sep 27 '16 at 12:32
  • Ok, now I am thinking that the problem must at how I specify the data x. Unfortunately, in the Mamba documentation I could't find an example where multivariate data are specified as the outcome of a distribution. – user1438310 Sep 27 '16 at 13:01
  • Have you looked at this? https://mambajl.readthedocs.io/en/latest/examples/birats.html – jfish003 Sep 27 '16 at 14:03
  • Yes, thanks, I have looked at this, but I don't know how to adapt it to the my simple example above. The thing is I don't quite understand the meaning of the single integer value for the Stochastic constructors entirely.... – user1438310 Sep 27 '16 at 14:13

1 Answers1

2

As I posted on the Mamba issue you openned:

The issue is because

data[:x]
2x30 Array{Float64,2}:

is a matrix of dimension 2 x 30. The way you coded up the stochastic node for x is

 x = Stochastic(1,
    mu -> MvNormal(mu, eye(2)),
    false
  ),

which specifies that x is a vector (multidimensional array of dimension 1). That's what the 1 right after Stochastic denotes. It helps to write out the model in math notation. Because the MvNormal defines a distribution on a vector, not a matrix. Perhaps your model is something like X_1, ..., X_n iid MvNormal(mu, I) in which case you can try something like

using Mamba

## Model Specification

model = Model(
  x = Stochastic(2,
    (mu, N, P) ->
      UnivariateDistribution[
      begin
        Normal(mu[i], 1)
      end
      for i in 1:P, j in 1:N
    ],
    false
  ),
  mu = Stochastic(1,
    () -> MvNormal(zeros(2), 1000.0),
    true
  )
)

## Data
data = Dict{Symbol, Any}(
:x => randn(2,30)+13,
:P => 2,
:N => 30
)
## Initial Values
inits = [
  Dict{Symbol, Any}(
    :x => data[:x],
    :mu => randn(2)*1
  )
]

## Sampling Scheme Assignment
scheme1 = NUTS([:mu])
setsamplers!(model, [scheme1])

sim1 = mcmc(model, data, inits, 10000, burnin=250, thin=2, chains=1);
describe(sim1)
bdeonovic
  • 4,130
  • 7
  • 40
  • 70