Calling mcmc() from within a function
arzwa opened this issue · 6 comments
Hi,
I'm very interested in using Mamba for a current project. I have enjoyed the tutorial and some examples so far.
However I have noticed some rather unexpected errors I believe when calling mcmc() from within another function (also briefly mentioned in this closed issue #105 (comment)). For example, calling the blr() function below gives the error ERROR: MethodError: no method matching (::getfield(Main, Symbol("##127#128")))(::Model) The applicable method may be too new: running in world age 25145, while current world is 25151.
Example function:
function blr()
# Model specification
model = Model(
y = Stochastic(1, (μ, σ²) -> MvNormal(μ, √σ²), false),
μ = Logical(1, (X, β) -> X * β, false),
β = Stochastic(1, () -> MvNormal(2, √1000)),
σ² = Stochastic(() -> InverseGamma(0.001, 0.001))
)
# Sampling scheme
scheme1 = [NUTS(:β), Slice(:σ², 3.0)]
setsamplers!(model, scheme1)
# Data
line = Dict{Symbol, Any}(:x => [1, 2, 3, 4, 5], :y => [1, 3, 3, 3, 5])
line[:X] = [ones(5) line[:x]]
# Initial values
inits = [Dict{Symbol, Any}(:y => line[:y], :β => rand(Normal(0, 1), 2),
:σ² => rand(Gamma(1, 1))) for i in 1:3]
# MCMC
sim1 = mcmc(model, line, inits, 10000, burnin=250, thin=2, chains=3)
return sim1
end The following code runs however without any issues:
function blr()
# Model specification
model = Model(
y = Stochastic(1, (μ, σ²) -> MvNormal(μ, √σ²), false),
μ = Logical(1, (X, β) -> X * β, false),
β = Stochastic(1, () -> MvNormal(2, √1000)),
σ² = Stochastic(() -> InverseGamma(0.001, 0.001))
)
# Sampling scheme
scheme1 = [NUTS(:β), Slice(:σ², 3.0)]
setsamplers!(model, scheme1)
# Data
line = Dict{Symbol, Any}(:x => [1, 2, 3, 4, 5], :y => [1, 3, 3, 3, 5])
line[:X] = [ones(5) line[:x]]
# Initial values
inits = [Dict{Symbol, Any}(:y => line[:y], :β => rand(Normal(0, 1), 2),
:σ² => rand(Gamma(1, 1))) for i in 1:3]
# MCMC
return model, line, inits
end
model, line, inits = blr()
mcmc(model, line, inits, 10000, burnin=250, thin=2, chains=3)So am I missing some details here or is this a bug?
Thanks in advance,
Arthur
I think this seems to be a bug. I have recreated with the following:
julia> versioninfo()
Julia Version 1.0.1
Commit 0d713926f8 (2018-09-29 19:05 UTC)
Platform Info:
OS: Windows (x86_64-w64-mingw32)
CPU: Intel(R) Core(TM) i7-8550U CPU @ 1.80GHz
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-6.0.0 (ORCJIT, skylake)
Environment:
JULIA_EDITOR = "C:\Users\mark\AppData\Local\atom\app-1.32.2\atom.exe" -a
JULIA_NUM_THREADS = 4
installed package status:
(v1.0) pkg> status
Status `C:\Users\mark\.julia\environments\v1.0\Project.toml`
[324073ee] Astro v0.0.0 #master (https://github.com/cormullion/Astro.jl)
[c52e3926] Atom v0.7.10
[336ed68f] CSV v0.4.2
[5d742f6a] CSVFiles v0.10.0
[159f3aea] Cairo v0.5.6
[a93c6f00] DataFrames v0.14.1
[31c24e10] Distributions v0.16.4
[38e38edf] GLM v1.0.1
[c91e804a] Gadfly v1.0.0
[09f84164] HypothesisTests v0.8.0
[7073ff75] IJulia v1.14.1
[e5e0dc1b] Juno v0.5.3
[b964fa9f] LaTeXStrings v1.0.3
[5424a776] Mamba v0.12.0
[86f7a689] NamedArrays v0.9.1
[47be7bcc] ORCA v0.2.0
[f0f68f2c] PlotlyJS v0.12.0
[91a5bcdd] Plots v0.21.0
[438e738f] PyCall v1.18.5
[d330b81b] PyPlot v2.6.3
[612083be] Queryverse v0.1.0
[6f49c342] RCall v0.13.0
[ce6b1742] RDatasets v0.6.1
[79098fc4] Rmath v0.5.0
[60ddc479] StatPlots v0.8.2
[fce5fe82] Turing v0.5.1
[c17dfb99] WinRPM v0.4.2
[8ba89e20] Distributed
[37e2e46d] LinearAlgebra
Example obtained from docs (https://mambajl.readthedocs.io/en/latest/examples/seeds.html), works when run outside a function with a minor modification to fully qualify the use of I (identity matrix) i.e. changed I to LinearAlgebra.I:
using Distributed
@everywhere using Mamba
import LinearAlgebra
function testymcmc()
## Data
seeds = Dict{Symbol, Any}(
:r => [10, 23, 23, 26, 17, 5, 53, 55, 32, 46, 10, 8, 10, 8, 23, 0, 3, 22, 15,
32, 3],
:n => [39, 62, 81, 51, 39, 6, 74, 72, 51, 79, 13, 16, 30, 28, 45, 4, 12, 41,
30, 51, 7],
:x1 => [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
:x2 => [0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1]
)
seeds[:N] = length(seeds[:r])
## Model Specification
model = Model(
r = Stochastic(1,
(alpha0, alpha1, x1, alpha2, x2, alpha12, b, n, N) ->
UnivariateDistribution[(
p = invlogit(alpha0 + alpha1 * x1[i] + alpha2 * x2[i] +
alpha12 * x1[i] * x2[i] + b[i]);
Binomial(n[i], p)) for i in 1:N
],
false
),
b = Stochastic(1,
s2 -> Normal(0, sqrt(s2)),
false
),
alpha0 = Stochastic(
() -> Normal(0, 1000)
),
alpha1 = Stochastic(
() -> Normal(0, 1000)
),
alpha2 = Stochastic(
() -> Normal(0, 1000)
),
alpha12 = Stochastic(
() -> Normal(0, 1000)
),
s2 = Stochastic(
() -> InverseGamma(0.001, 0.001)
)
)
## Initial Values
inits = [
Dict(:r => seeds[:r], :alpha0 => 0, :alpha1 => 0, :alpha2 => 0,
:alpha12 => 0, :s2 => 0.01, :b => zeros(seeds[:N])),
Dict(:r => seeds[:r], :alpha0 => 0, :alpha1 => 0, :alpha2 => 0,
:alpha12 => 0, :s2 => 1, :b => zeros(seeds[:N]))
]
## Sampling Scheme
scheme = [AMM([:alpha0, :alpha1, :alpha2, :alpha12], 0.01 * Matrix{Float64}(LinearAlgebra.I, 4, 4)),
AMWG(:b, 0.01),
AMWG(:s2, 0.1)]
setsamplers!(model, scheme)
## MCMC Simulations
sim = mcmc(model, seeds, inits, 12500, burnin=2500, thin=2, chains=2)
end
out = testymcmc()
The error and call stack:
MethodError: no method matching (::getfield(Main, Symbol("##123#124")))(::Model)
The applicable method may be too new: running in world age 25222, while current world is 25232.
Closest candidates are:
#123(::Model) at none:0 (method too new to be called from this world context.)
setinits!(::ScalarStochastic, ::Model, ::Int64) at dependent.jl:163
setinits!(::Model, ::Dict{Symbol,Any}) at initialization.jl:11
setinits!(::Model, ::Array{Dict{Symbol,Any},1}) at initialization.jl:24
#mcmc#23(::Int64, ::Int64, ::Int64, ::Bool, ::Function, ::Model, ::Dict{Symbol,Any}, ::Array{Dict{Symbol,Any},1}, ::Int64) at mcmc.jl:30
(::getfield(Mamba, Symbol("#kw##mcmc")))(::NamedTuple{(:burnin, :thin, :chains),Tuple{Int64,Int64,Int64}}, ::typeof(mcmc), ::Model, ::Dict{Symbol,Any}, ::Array{Dict{Symbol,Any},1}, ::Int64) at none:0
testymcmc() at julia_tmp4.jl:76
top-level scope at none:0
include_string(::Module, ::String, ::String) at loading.jl:1005
include_string(::Module, ::String, ::String, ::Int64) at eval.jl:30
(::getfield(Atom, Symbol("##114#119")){String,Int64,String})() at eval.jl:94
withpath(::getfield(Atom, Symbol("##114#119")){String,Int64,String}, ::String) at utils.jl:30
withpath at eval.jl:46 [inlined]
#113 at eval.jl:93 [inlined]
with_logstate(::getfield(Atom, Symbol("##113#118")){String,Int64,String}, ::Base.CoreLogging.LogState) at logging.jl:397
with_logger at logging.jl:493 [inlined]
#112 at eval.jl:92 [inlined]
hideprompt(::getfield(Atom, Symbol("##112#117")){String,Int64,String}) at repl.jl:85
macro expansion at eval.jl:91 [inlined]
macro expansion at dynamic.jl:24 [inlined]
(::getfield(Atom, Symbol("##111#116")))(::Dict{String,Any}) at eval.jl:86
handlemsg(::Dict{String,Any}, ::Dict{String,Any}) at comm.jl:164
(::getfield(Atom, Symbol("##19#21")){Array{Any,1}})() at task.jl:259
Please can someone take a look?
Also see: #105
Further to above the problems seems to relate solely to calling the Mamba.mcmc function. I can encapsulate everything else within a function then return the elements in a Dict that can be used to call mcmc outside the function.
It is a known problem with no current solution. Mamba needs to restructure how it handles and passes around anonymous functions.
Thanks Benjamin. I will take a look to see what is involved.
Probably a big task. Other people in other packages have had similar issues: JuliaLang/julia#21356
Just a tip for people who have this problem and need a quick fix:
As a way to get around the problem, while the bug is not fixed, put the model definition outside the function and reference it. Also in the model definition make sure that all variables with input data are named, so that you can change the input values inside the function.
For me this is a much more flexible way to still be able to structure my program, as the model definition itself doesn't change, just the input data.