学习在 Julia 中使用 Mamba 包

Learning to use the Mamba package in Julia

我正在尝试学习如何使用 Julia 中的 Mamba 包进行贝叶斯推理。虽然这个包很棒,但作为初学者,我发现文档中的信息有点匮乏。因此,我试图弄清楚如何实现一些非常简单的示例。

我试过的

我实现了一个对单变量正态分布的均值进行贝叶斯推理的示例。代码如下:

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)

这似乎工作得非常好(尽管可能有更好的编码方法?)。

我正在尝试做但不起作用。

在这个例子中,我试图对二元正态分布的均值进行贝叶斯推断。代码如下:

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)

您可能已经注意到,我认为必要的更改很少。但是,我在某处做错了,当我尝试 运行 时,我收到一个错误(类型之间的转换错误),这对我没有进一步的帮助。

感谢任何帮助。如果成功,我会考虑将这个简单示例贡献给其他新用户的 Mamba 文档。谢谢

附录:错误信息

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

正如我在您打开的 Mamba 问题上发布的那样:

这个问题是因为

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

是一个 2 x 30 维矩阵。你对 x 的随机节点进行编码的方式是

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

指定x是一个向量(1维的多维数组)。这就是 Stochastic 之后的 1 所表示的。它有助于用数学符号写出模型。因为 MvNormal 定义了向量而非矩阵的分布。也许您的模型类似于 X_1, ..., X_n iid MvNormal(mu, I) 在这种情况下您可以尝试类似

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)