I’m trying to optimize this log-likelihood function:
-log(σ^2) / 2 - ϵ^2 / 2 + log( 2Φ(a * ϵ) )
Here ϵ depends on data and some more parameters.
I’m trying to optimize this in two ways:
- Put the sum of these likelihoods (over the data) right into
@NLobjective - Put the likelihood into a function, register it and put the sum over invocations of this function into
@NLobjective.
I expected the results (parameter estimates) to be the same, but they’re different, even though the likelihood is exactly the same in both cases.
The code looks like this:
begin
import Pkg
Pkg.activate(temp=true)
Pkg.add([
"CSV", "DataFrames", "Ipopt", "JuMP", "SpecialFunctions"
], io=devnull)
Pkg.status()
end
using CSV, DataFrames
using JuMP
import Ipopt
using SpecialFunctions: erf
const AV = AbstractVector{T} where T;
"Normal CDF"
Φ(x::Real) = (1 + erf(x / sqrt(2))) / 2
# 1. Load data
df = CSV.read("data_sample.csv", DataFrame)
# 2. Define log-likelihood for one observation
logLik_SF(ϵ::Real, a::Real, σ::Real) = -log(σ^2) / 2 - ϵ^2 / 2 + log( 2Φ(a * ϵ) )
# 3. Define function for estimation
function estimate_SF(Y::AV{<:Real}, L::AV{<:Real}, K::AV{<:Real}; optimizer=Ipopt.Optimizer, ll_sum::Bool=true)
N = length(Y)
y = log.(Y)
l = log.(L)
k = log.(K)
model = Model(optimizer)
set_optimizer_attribute(model, MOI.Silent(), true)
@variable(model, α, start=0.5) # elasticity
@variable(model, a, start=0) # skewed normal shape
@variable(model, σ ≥ 0) # skewed normal scale
ϵ = [
@NLexpression(
model,
(y[i] - (α * l[i] + (1 - α) * k[i])) / σ
)
for i ∈ eachindex(y)
]
if ll_sum
register(model, :Φ, 1, Φ, autodiff=true)
@NLobjective(
model, Max,
# Pasted the content of `logLik_SF`
# and replaced `ϵ` with `ϵ[i]`
sum(
-log(σ^2) / 2 - ϵ[i]^2 / 2 + log( 2Φ(a * ϵ[i]) )
for i ∈ eachindex(ϵ)
)
)
else
register(model, :logLik_SF, 3, logLik_SF, autodiff=true)
@NLobjective(
model, Max,
# Just call `logLik_SF`
sum(
logLik_SF(ϵ[i], a, σ)
for i ∈ eachindex(ϵ)
)
)
end
optimize!(model)
model
end
println("Estimate using hard-coded expression:")
model = estimate_SF(
df[!, :Y], df[!, :L], df[!, :C],
ll_sum=true
)
a, α, σ = value(model[:a]), value(model[:α]), value(model[:σ])
@show a α σ
println("\nEstimate using function call:")
model = estimate_SF(
df[!, :Y], df[!, :L], df[!, :C],
ll_sum=false
)
a, α, σ = value(model[:a]), value(model[:α]), value(model[:σ])
@show a α σ
The data come from data_sample.csv: Some data - Pastebin.com.
I’m running the code like this: julia-1.7 report.jl
Output:
Status `/private/var/folders/ys/3h0gnqns4b98zb66_vl_m35m0000gn/T/jl_vWtB7k/Project.toml`
[336ed68f] CSV v0.9.11
[a93c6f00] DataFrames v1.3.0
[b6b21f68] Ipopt v0.7.0
[4076af6c] JuMP v0.21.5
[276daf66] SpecialFunctions v2.0.0
Estimate using hard-coded expression:
a = -0.5226232660699477
α = 0.1438052406305529
σ = 0.2727175471168534
Estimate using function call:
a = 1.5968811212735676
α = 0.18443033301359052
σ = 0.3456179651592732
As you can see, I’m getting different results depending on whether I use -log(σ^2) / 2 - ϵ[i]^2 / 2 + log( 2Φ(a * ϵ[i]) ) or logLik_SF(ϵ[i], a, σ) in the sum in the objective function. …but they’re literally the same exact function!
What could be the problem here?