Maximum likelihood optimization problem through Optim.jl

Specific Domains Optimization (Mathematical)
,
1

Hi,
I am trying to solve a likelihood function in Optim as follows:
I have some increments which are gamma-distributed (Ga(a*t, β)):
det_x = [0.0175, 0.0055, 0.0059] # increments
det_t = [185, 163, 167] # corresponding time
I want to estimate parameters a, and b from the above data.
I used the following program:

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
## density function for Gamma(a*t, β)
gamma_dens(x, a, β, t) = (x^(a*t-1)*exp(-x/β))/(β^(a*t)*gamma(a*t))
function L(x)
    log(gamma_dens(det_R[1], x[1], x[2], det_w[1])) + 
    log(gamma_dens(det_R[2], x[1], x[2], det_w[2])) +
    log(gamma_dens(det_R[3], x[1], x[2], det_w[3])) 
end
x0 = [1e-3, 1e-3]
opt = optimize(L, x0)

But I got error message saying:
Output exceeds the [size limit](command:workbench.action.openSettings?[). Open the full output data [in a text editor](command:workbench.action.openLargeOutput?81dfbb88-756a-4313-986f-379f044b325b)
DomainError with -0.0245: Exponentiation yielding a complex result requires a complex argument. Replace x^y with (x+0im)^y, Complex(x)^y, or similar.
How can I solve this optimization problem?
Thank you very much!

2

DomainError with -0.0245: Exponentiation yielding a complex result requires a complex argument. Replace x^y

Your x values are going negative, giving you complex values. You can define some bounds on your parameter values (remembering that \alpha, \beta > 0): https://julianlsolvers.github.io/Optim.jl/stable/#user/minimization/#box-constrained-optimization.

Another nice way would be, in your L(x) function, use exp(x[1]) and exp(x[2]) so that you don’t have to use any bounds, taking \alpha = \exp(x_1) and \beta = \exp(x_2) at the end.

See also Optim.jl for an actual Optim example for finding MLEs too that might be of interest to you.

2 Likes
3

Thank you very much for the reply.
I will try your suggestions.

4

I have changed my solution based on the example in Optim, still not working:

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
## density function for Gamma(a*t, β)
gamma_dens(x, a, β, t) = (x^(a*t-1)*exp(-x/β))/(β^(a*t)*gamma(a*t))
function L(x, y)
    c = exp(x)
    d = exp(y)
    log(gamma_dens(det_R[1], c, d, det_w[1])) + 
    log(gamma_dens(det_R[2], c, d, det_w[2])) +
    log(gamma_dens(det_R[3], c, d, det_w[3])) 
end
func = TwiceDifferentiable((x, y)->L(x, y), [0.1, 0.1]; autodiff=:forward);
opt = optimize(func, [0.1, 0.1])

Thank you in advance!

5

See the line

func = TwiceDifferentiable((x, y)->L(x, y), [0.1, 0.1]; autodiff=:forward)

Optim requires that your functions have a single input, so you should write this as

func = TwiceDifferentiable(x->L(x[1], x[2]), [0.1, 0.1]; autodiff=:forward)

I believe it should work then. I believe you got this from the line

func = TwiceDifferentiable(vars -> Log_Likelihood(x, y, vars[1:nvar], vars[nvar + 1]),
                           ones(nvar+1); autodiff=:forward);

in the example I linked - notice that x and y are just data in their example, but the actual parameters being optimised are those in vars - the function just takes in a single input.

1 Like
6

Or just x -> L(x[1],x[2]) ?

1 Like
7

Oops, missed that - you’re right.

8

Thank for the reply and explanation.
It works now. But I got an output of objective as -Inf and Optim.minimize(opt) are the same as
the initial values I have set.
I need to check the objective function further.

9

Does it work if you multiply by -1 first? Try it with

func = TwiceDifferentiable(x->-L(x[1], x[2]), [0.1, 0.1]; autodiff=:forward)

(Optim minimises by default, so this will make the maximum of L be given by the minimum of -L.)

10

No. Same as before.
I mean add -, the objective value after optimization is Inf.
The optimized results are still the same as initial values.

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
12

I don’t know the theory behind your problem, but some adjustments to your original formulation gives a result:

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
## density function for Gamma(a*t, β)
gamma_dens(x, a, β, t) = (x^(a*t-1)*exp(-x/β))/(β^(a*t)*gamma(a*t))
function L(x)
    log(gamma_dens(det_R[1], x[1], x[2], det_w[1])) + 
    log(gamma_dens(det_R[2], x[1], x[2], det_w[2])) +
    log(gamma_dens(det_R[3], x[1], x[2], det_w[3])) 
end

x0 = [1e-3, 1e-3]

func = TwiceDifferentiable(L, x0; autodiff=:forward);

lower_bounds = [0.0, 0.0]; upper_bounds = [Inf, Inf];
inner_optimizer = GradientDescent()
opt = optimize(func.f, func.df, lower_bounds, upper_bounds, x0, Fminbox(inner_optimizer))

x_sol = opt.minimizer

But it’s hard to see the result as meaningful since it seems to highly depend on your initial point x0, as is characteristic of (nonconvex) nonlinear optimisation problems.

2 Likes
13

Thank you very much for providing this solution.

1 Like
14

Thank you very much for checking this. Your remarks are really important for my problem.
For my problem, I just try to fit some observations (data) into Gamma distribution to estimate
the shape and scale parameter.
I will carefully check the solution.

15

I checked the solution offered by @DanielVandH , it gives a stable
solution around (0.039, 0.0014). I tested different initial values.
I also run your solution, it gives different values for different x_0 as
you said. It is weird. I see in your solution you used a Boxed contained
solver.