Hi all. I’m working on an EM algorithm which will necessarily call many optimization solves (using Optim.jl). That means that I’ll need to optimize various functions while holding other values constant. The obvious way to do this, to me, would be to make a closure. You can see in the function laplace_urchin I do this to optimize nlyfbs_urchin with respect to b but holding all else constant. However, when running laplace_urchin I get an UndefVarError for the θ′ variable, which is an argument to laplace_urchin. I guess this is because at the time the closure is created (when laplace_urchin is defined), it looks in the global scope to find θ′ and then can’t find it when it’s called even though it is a local argument to the function scope. Is my diagnosis of the problem right?
I’m looking for help to both understand better how the closure mechanism causes this problem, and also recommendations to fix my code. I attached the MWE below as well as the stacktrace. Thanks!
using DataFrames
using Distributions
using LinearAlgebra
using Optim, LineSearches
using DifferentiationInterface
import ForwardDiff
const ad_sys = AutoForwardDiff()
urchin = DataFrame(
age = [1.0,6.0,27.0],
vol = [0.1,1.5,35.0]
)
# index into a single vector for b (random effects)
const log_g_ix = 1:nrow(urchin)
const log_p_ix = range(start=nrow(urchin)+1, length=nrow(urchin))
# index into a single vector for θ (fixed effects)
const log_ω_ix = 1
const μ_g_ix = 2
const log_σ_g_ix = 3
const μ_p_ix = 4
const log_σ_p_ix = 5
const log_σ_ix = 6
# initial values for fixed/random effects
th = [
-3.0,
-0.3,
-1.5,
0.15,
-1.5,
-1.37
]
n = nrow(urchin)
b = [fill(th[2],n); fill(th[4],n)]
function model_urchin_vol(ω, g, p, a)
aₘ = log(p / (g*ω))/g
return ifelse(a < aₘ, ω*exp(g*a), p/g + p*(a-aₘ))
end
function lfyb_urchin(b, θ, urchin)
# extract fixed effects
σ = exp(θ[log_σ_ix])
σ_g = exp(θ[log_σ_g_ix])
σ_p = exp(θ[log_σ_p_ix])
ω = exp(θ[log_ω_ix])
# extract random effects
log_g = b[log_g_ix]
log_p = b[log_p_ix]
# calculate the loglikelihood of y and b
logll = 0.0
for i in axes(urchin,1)
v = model_urchin_vol(ω, exp(log_g[i]), exp(log_p[i]), urchin[i, :age])
logll += logpdf(Normal(sqrt(v), σ), sqrt(urchin[i, :vol]))
logll += logpdf(Normal(θ[μ_g_ix], σ_g), log_g[i])
logll += logpdf(Normal(θ[μ_p_ix], σ_p), log_p[i])
end
return logll
end
function lyfbs_urchin(b, s, θ, θ′, urchin)
return s * lfyb_urchin(b, θ, urchin) + lfyb_urchin(b, θ′, urchin)
end
function nlyfbs_urchin(b, s, θ, θ′, urchin)
return -lyfbs_urchin(b, s, θ, θ′, urchin)
end
prep_g_nlfybs = prepare_gradient(nlyfbs_urchin, ad_sys, zero(b), Constant(0.0), Constant(th), Constant(th), Constant(urchin))
g_nlfybs!(G, b) = gradient!(nlyfbs_urchin, G, prep_g_nlfybs, ad_sys, b, Constant(0.1), Constant(th), Constant(th), Constant(urchin))
function laplace_urchin(b, s, θ, θ′, urchin)
nlfybs_b_hat = optimize(
b -> nlyfbs_urchin(b, s, θ, Θ′, urchin),
g_nlfybs!,
b,
LBFGS(;alphaguess=InitialStatic(scaled=true), linesearch=BackTracking())
)
b_hat = Optim.minimizer(nlfybs_b_hat)
H = hessian(lyfbs_urchin, prep_sp_h_lfybs, sp_ad_sys, b_hat, Constant(s), Constant(θ), Constant(θ′), Constant(urchin))
return logdet(-H)
end
laplace_urchin(b, 0.5, th, th+rand(length(th)), urchin)
Stacktrace
ERROR: UndefVarError: `Θ′` not defined in `Main`
Suggestion: check for spelling errors or missing imports.
Stacktrace:
[1] (::var"#11#12"{Float64, Vector{Float64}, DataFrame})(b::Vector{Float64})
@ Main ~/Desktop/git/MessyNapkin/Statistics/Urchin/bug.jl:79
[2] (::NLSolversBase.var"#fg!#8"{var"#11#12"{…}, typeof(g_nlfybs!)})(gx::Vector{Float64}, x::Vector{Float64})
@ NLSolversBase ~/.julia/packages/NLSolversBase/totsP/src/objective_types/abstract.jl:14
[3] value_gradient!!(obj::OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, x::Vector{Float64})
@ NLSolversBase ~/.julia/packages/NLSolversBase/totsP/src/interface.jl:82
[4] initial_state(method::LBFGS{…}, options::Optim.Options{…}, d::OnceDifferentiable{…}, initial_x::Vector{…})
@ Optim ~/.julia/packages/Optim/8dE7C/src/multivariate/solvers/first_order/l_bfgs.jl:168
[5] optimize
@ ~/.julia/packages/Optim/8dE7C/src/multivariate/optimize/optimize.jl:43 [inlined]
[6] #optimize#79
@ ~/.julia/packages/Optim/8dE7C/src/multivariate/optimize/interface.jl:301 [inlined]
[7] optimize (repeats 2 times)
@ ~/.julia/packages/Optim/8dE7C/src/multivariate/optimize/interface.jl:289 [inlined]
[8] laplace_urchin(b::Vector{Float64}, s::Float64, θ::Vector{Float64}, θ′::Vector{Float64}, urchin::DataFrame)
@ Main ~/Desktop/git/MessyNapkin/Statistics/Urchin/bug.jl:78
[9] top-level scope
@ ~/Desktop/git/MessyNapkin/Statistics/Urchin/bug.jl:89
Some type information was truncated. Use `show(err)` to see complete types.
. Thanks. EM algorithm is: