学习在Julia中使用Mamba包

2019年7月20日 215次阅读

我正在尝试学习如何使用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的多维数组).这就是随机指数之后的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)
点赞