# Different results when registering the log-likelihood function vs pasting its contents

**URL:** <https://discourse.julialang.org/t/different-results-when-registering-the-log-likelihood-function-vs-pasting-its-contents/73214>\
**Category:** Optimization (Mathematical)\
**Tags:** question, jump, ipopt\
**Created:** [December 16, 2021, 4:28pm UTC](https://discourse.julialang.org/t/different-results-when-registering-the-log-likelihood-function-vs-pasting-its-contents/73214 "2021-12-16T16:28:20Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![ForceBru](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/forcebru/32/21389_2.png) [@ForceBru](https://discourse.julialang.org/u/ForceBru)\
**Post date:** [December 16, 2021, 4:28pm UTC](https://discourse.julialang.org/t/different-results-when-registering-the-log-likelihood-function-vs-pasting-its-contents/73214/1 "2021-12-16T16:28:20Z")

</div>

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:

1. Put the sum of these likelihoods (over the data) right into `@NLobjective`
2. 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:

```julia
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](https://pastebin.com/6Nety2ZN).

I’m running the code like this: `julia-1.7 report.jl`

Output:

```julia
      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?

---

<div class="post-metadata">

**Author:** ![ForceBru](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/forcebru/32/21389_2.png) [@ForceBru](https://discourse.julialang.org/u/ForceBru)\
**Post date:** [December 16, 2021, 4:37pm UTC](https://discourse.julialang.org/t/different-results-when-registering-the-log-likelihood-function-vs-pasting-its-contents/73214/2 "2021-12-16T16:37:38Z")

</div>

Random observation (that _doesn’t_ solve the issue):

- installing packages like in the post installs:
  - “old” Ipopt v0.7.0
  - “old” JuMP v0.21.5
  - new SpecialFunctions v2.0.0

- installing Ipopt and JuMP first, and then SecialFunctions installs:
  - _much_ newer Ipopt v0.9.1
  - new JuMP v0.22.1
  - old SpecialFunctions v1.8.1

Multiple equilibria in package management strike again! The problem happens on all of these versions, though…

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [December 16, 2021, 7:17pm UTC](https://discourse.julialang.org/t/different-results-when-registering-the-log-likelihood-function-vs-pasting-its-contents/73214/3 "2021-12-16T19:17:51Z")

</div>

If you look at the Ipopt log of the two problems, they are not equivalent. In the `register` case, we disable second-order derivatives: [Feature Request: Second Derivatives for User Defined Functions · Issue #1198 · jump-dev/JuMP.jl · GitHub](https://github.com/jump-dev/JuMP.jl/issues/1198)

> Multiple equilibria in package management strike again!

Again, this is expected. The package manager does’t know whether you want the latest JuMP or the latest SpecialFunctions. Use `] add JuMP@0.22` to force it to choose JuMP.

---

<div class="post-metadata">

**Author:** ![ForceBru](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/forcebru/32/21389_2.png) [@ForceBru](https://discourse.julialang.org/u/ForceBru)\
**Post date:** [December 16, 2021, 7:27pm UTC](https://discourse.julialang.org/t/different-results-when-registering-the-log-likelihood-function-vs-pasting-its-contents/73214/4 "2021-12-16T19:27:37Z")

</div>

> [@odow](#):
>
> this is expected

Yeah, I get this now - I’m just still kinda frustrated with this, but that’s beside the point.

> [@odow](#):
>
> In the `register` case, we disable second-order derivatives

Okay, got it, thanks!