# Maximum likelihood optimization problem through Optim.jl

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, x, x, det_w)) +
log(gamma_dens(det_R, x, x, det_w)) +
log(gamma_dens(det_R, x, x, det_w))
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!

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) and exp(x) 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

Thank you very much for the reply.

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, c, d, det_w)) +
log(gamma_dens(det_R, c, d, det_w)) +
log(gamma_dens(det_R, c, d, det_w))
end
func = TwiceDifferentiable((x, y)->L(x, y), [0.1, 0.1]; autodiff=:forward);
opt = optimize(func, [0.1, 0.1])


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, x), [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

Or just x -> L(x,x) ?

1 Like

Oops, missed that - you’re right.

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.

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

func = TwiceDifferentiable(x->-L(x, x), [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.)

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

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 * c, d), det_R) +
logpdf(Gamma(det_w * c, d), det_R) +
logpdf(Gamma(det_w * c, d), det_R)
end
func = TwiceDifferentiable(x->-L(x, x), [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

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, x, x, det_w)) +
log(gamma_dens(det_R, x, x, det_w)) +
log(gamma_dens(det_R, x, x, det_w))
end

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

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

lower_bounds = [0.0, 0.0]; upper_bounds = [Inf, Inf];
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

Thank you very much for providing this solution.

1 Like

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.

I checked the solution offered by @legola18 , 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.