Maximum likelihood optimization problem through Optim.jl

Specific Domains Optimization (Mathematical)
,
11

I think this is a numerical issue. If we consider c = exp(0.1) and d = exp(0.1) in your code (as given by your initial estimate), then your first call to to gamma_dens is (approximately) gamm_dens(0.0175, 1.1, 1.1, 185). The numbers that you get out of this are extremely small. The issue then comes from x^(a*t-1) - the exponent here is 203! So x^(a*t-1) = 0.0175^203 which may as well be zero (and Julia calls it zero, actually:

julia> x^(a*t-1)
0.0

So all your code is working, it’s just the numbers need some work.

You should use logpdf from Distributions so that you get stable results. Here’s a working version.

using SpecialFunctions
using Distributions, LinearAlgebra, Statistics
using Optim, NLSolversBase
det_R = [0.0175, 0.0055, 0.0059] # increments
det_w = [185, 163, 167] # corresponding time
function L(x, y)
    c = exp(x)
    d = exp(y)
    logpdf(Gamma(det_w[1] * c, d), det_R[1]) +
    logpdf(Gamma(det_w[2] * c, d), det_R[2]) +
    logpdf(Gamma(det_w[3] * c, d), det_R[3])
end
func = TwiceDifferentiable(x->-L(x[1], x[2]), [0.1, 0.1]; autodiff=:forward);
opt = optimize(func, [0.1, 0.1])
julia> opt.minimum
-12.159601199591297

julia> exp.(opt.minimizer)
2-element Vector{Float64}:
 0.02452153013876933
 0.0022884585316169555
2 Likes