Learning to use the Mamba package (also posted on stackoverflow)

Question

Learning to use the Mamba package (also posted on stackoverflow)

Opened this issue 8 years ago · 5 comments

Dear all,

I am sorry to be posting here, as far as I understand the forum here is intended for development issues mostly.

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 little difficult to follow. Hence I am trying to figure out how to 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.

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

Ps. The question has been previously posted on StackOverflow at: http://bit.ly/2d6CcY1

Answer 1 · 2016-09-27T14:43:07.000Z

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)

Answer 2 · 2016-09-28T09:40:26.000Z

Dear bdeonovic,

thanks very much for prompt reply.

Indeed I have problems understanding what the "1 right after Stochastic denotes." and I have to admit that I still do not understand it fully. What I currently understand is that, since x is 2x30 matrix what should come out of the Stochastic constructor for x is also a 2x30 object, is this right?

However: if this is the case, why doesn't the following code below work?
Notice that I have rearranged the data in x so that it now is Array{Array{Float64,1},1}.
I was hoping that the Stochastic constructor for x should be generating a result of the same type.
Evidently it is not as the code below does not work.
This probably means that I still fail to grasp something important.

Thanks again for all your time and patience.

using Mamba

## Model Specification

model = Model(

  x = Stochastic(1,
    (mu,N) ->
    MultivariateDistribution[
      begin
        MvNormal(mu, eye(2))
      end
      for i in 1:N
    ],
    false
  ),

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

)

## Data
N = 5
data = Dict{Symbol, Any}(
  :x => [randn(2)+13 for i=1:N],
  :N => N
)

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

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

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

Ps. the reason why I am insisting on re-writting things instead of using your working suggestion, is because as a next step I would like to do the same but instead of using a spherical multivariate gaussian like now, i.e. MvNormal(mu,eye(2)), I would like to use a gaussian with a full covariance matrix, i.e. MvNormal(mu,Sigma). This means that I will have to treat the two components of the mean mu jointly and not independently.

Answer 3 · 2016-09-28T17:00:49.000Z

You can get it to work with an arbitrary covariance matrix using the following code:

using Mamba

## Model Specification
model = Model(
  x = Stochastic(2,
    (mu, N, Sigma) ->
      MultivariateDistribution[
        begin
          MvNormal(mu, Sigma)
        end
        for i in 1:N
      ],
      false
  ),
  mu = Stochastic(1,
    () -> MvNormal(zeros(2), 1000.0),
    true
  )
)

## Data
data = Dict{Symbol, Any}(
  :x => randn(30,2)+13,  ##NOTE switch on dimensions!!!
  :P => 2,
  :N => 30,
  :Sigma => eye(2)  ##this can be anything, or you can put a prior on it to make it Stochastic (e.g. InverseWishart)
)
## 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)

Note I had to switch the dimensions on your x matrix because Mamba evaluates logpdf column wise on matrices.

A good example that this is similar to is birats (the beta parameter is similar to your x). The reason your code doesn't work is that the nodes have to be multidimensional arrays of Float64. In your most recent example you have Arrays of Arrays of Float64.

If you have any other questions please feel free to keep posting on the issues on the github. Brian and I have better visibility of posts on there (versus StackOverflow or the Julia Google group). But of course other users of Mamba will be able to respond to you on those other sites.

Answer 4 · 2016-09-29T09:35:14.000Z

Dear bdeonovic,

thanks very much for your reply indeed.

Despite your thorough answer, I am still unsure if I understand correctly how to
properly set the dimensions of the various quantities.

In an attempt to get a better understanding, I defined a function in the Julia command prompt:

f = (mu, N, Sigma) ->
      MultivariateDistribution[
        begin
          MvNormal(mu, Sigma)
        end
        for i in 1:N
      ]

This the function used in the Stochastic constructor for the data.

I then call this function with some arbitrary arguments and check the size of what is returned.
Unsuprisingly, when I execute:

size(f(zeros(2),30,eye(2))

I receive

(30,)

This coincides with the first dimension of data :x .
Conclusion 1: So this covers the first dimension of :x

I then also tried to define in the command prompt:

x = Stochastic(2,
    (mu, N, Sigma) ->
      MultivariateDistribution[
        begin
          MvNormal(mu, Sigma)
        end
        for i in 1:N
      ],
      false
  )

I then noticed that x is provided with a field called value whose type is Array{Float64,2}.
Conclusion 2: I guess that the 2 in Array{Float64,2} coincides with the second dimension of the 30x2 data x.

Are these two conclusions correct?

Thanks very much for your time!

ps. I had looked into the birats example before I posted here, but at the time it didn't help me. Now in light of our correspondence, I will have another look.

Answer 5 · 2016-10-02T02:22:18.000Z

If you look at your function f you will see that it basically creates an array of length N similar to how

Int64[5 for i in 1:10]

would create an array of length 10 that is filled with 5. This is why the dimension of the output of your function matches the first dimension of x.

Array{Float64,2} means this is a 2 dimensional array (i.e. matrix). A matrix with 4 rows and 5 columns would also be a Array{Float64, 2}. We need to specify x as Stochastic(2,...) because your data xis stored as a matrix.