# Maximum likelihood optimization problem through Optim.jl
**URL:** https://discourse.julialang.org/t/maximum-likelihood-optimization-problem-through-optim-jl/84096
**Category:** Optimization (Mathematical)
**Tags:** optim, mle
**Created:** [July 12, 2022, 9:38am UTC](https://discourse.julialang.org/t/maximum-likelihood-optimization-problem-through-optim-jl/84096 "2022-07-12T09:38:55Z")
**Posts on this page:** 15
**Page:** 1
### Author:  [@Jian\_ZUO](https://discourse.julialang.org/u/Jian_ZUO)
#### Post date: [July 12, 2022, 9:38am UTC](https://discourse.julialang.org/t/maximum-likelihood-optimization-problem-through-optim-jl/84096/1 "2022-07-12T09:38:55Z")
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:
```julia
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!
---
### Author:  [@DanielVandH](https://discourse.julialang.org/u/DanielVandH)
#### Post date: [July 12, 2022, 10:15am UTC](https://discourse.julialang.org/t/maximum-likelihood-optimization-problem-through-optim-jl/84096/2 "2022-07-12T10:15:25Z")
> `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](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](https://julianlsolvers.github.io/Optim.jl/stable/#examples/generated/maxlikenlm/) for an actual Optim example for finding MLEs too that might be of interest to you.
---
### Author:  [@Jian\_ZUO](https://discourse.julialang.org/u/Jian_ZUO)
#### Post date: [July 12, 2022, 10:41am UTC](https://discourse.julialang.org/t/maximum-likelihood-optimization-problem-through-optim-jl/84096/3 "2022-07-12T10:41:18Z")
Thank you very much for the reply.
I will try your suggestions.
---
### Author:  [@Jian\_ZUO](https://discourse.julialang.org/u/Jian_ZUO)
#### Post date: [July 12, 2022, 11:39am UTC](https://discourse.julialang.org/t/maximum-likelihood-optimization-problem-through-optim-jl/84096/4 "2022-07-12T11:39:42Z")
I have changed my solution based on the example in Optim, still not working:
```julia
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!
---
### Author:  [@DanielVandH](https://discourse.julialang.org/u/DanielVandH)
#### Post date: [July 12, 2022, 11:50am UTC](https://discourse.julialang.org/t/maximum-likelihood-optimization-problem-through-optim-jl/84096/5 "2022-07-12T11:50:40Z")
See the line
```julia
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
```julia
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
```julia
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.
---
### Author:  [@jd-foster](https://discourse.julialang.org/u/jd-foster)
#### Post date: [July 12, 2022, 11:56am UTC](https://discourse.julialang.org/t/maximum-likelihood-optimization-problem-through-optim-jl/84096/6 "2022-07-12T11:56:01Z")
> [@DanielVandH](#):
>
> `(x, y)->L(x[1], x[2])`
Or just `x -> L(x[1],x[2])` ?
---
### Author:  [@DanielVandH](https://discourse.julialang.org/u/DanielVandH)
#### Post date: [July 12, 2022, 11:57am UTC](https://discourse.julialang.org/t/maximum-likelihood-optimization-problem-through-optim-jl/84096/7 "2022-07-12T11:57:53Z")
Oops, missed that - you’re right.
---
### Author:  [@Jian\_ZUO](https://discourse.julialang.org/u/Jian_ZUO)
#### Post date: [July 12, 2022, 12:02pm UTC](https://discourse.julialang.org/t/maximum-likelihood-optimization-problem-through-optim-jl/84096/8 "2022-07-12T12:02:58Z")
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.
---
### Author:  [@DanielVandH](https://discourse.julialang.org/u/DanielVandH)
#### Post date: [July 12, 2022, 12:03pm UTC](https://discourse.julialang.org/t/maximum-likelihood-optimization-problem-through-optim-jl/84096/9 "2022-07-12T12:03:44Z")
Does it work if you multiply by `-1` first? Try it with
```julia
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`.)
---
### Author:  [@Jian\_ZUO](https://discourse.julialang.org/u/Jian_ZUO)
#### Post date: [July 12, 2022, 12:08pm UTC](https://discourse.julialang.org/t/maximum-likelihood-optimization-problem-through-optim-jl/84096/10 "2022-07-12T12:08:58Z")
No. Same as before.
I mean add -, the objective value after optimization is Inf.
The optimized results are still the same as initial values.
---
### Author:  [@DanielVandH](https://discourse.julialang.org/u/DanielVandH)
#### Post date: [July 12, 2022, 12:17pm UTC](https://discourse.julialang.org/t/maximum-likelihood-optimization-problem-through-optim-jl/84096/11 "2022-07-12T12:17:21Z")
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
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.
```julia
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
```
---
### Author:  [@jd-foster](https://discourse.julialang.org/u/jd-foster)
#### Post date: [July 12, 2022, 12:28pm UTC](https://discourse.julialang.org/t/maximum-likelihood-optimization-problem-through-optim-jl/84096/12 "2022-07-12T12:28:10Z")
I don’t know the theory behind your problem, but some adjustments to your original formulation gives a result:
```julia
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.
---
### Author:  [@Jian\_ZUO](https://discourse.julialang.org/u/Jian_ZUO)
#### Post date: [July 12, 2022, 12:33pm UTC](https://discourse.julialang.org/t/maximum-likelihood-optimization-problem-through-optim-jl/84096/13 "2022-07-12T12:33:59Z")
Thank you very much for providing this solution.
---
### Author:  [@Jian\_ZUO](https://discourse.julialang.org/u/Jian_ZUO)
#### Post date: [July 12, 2022, 12:36pm UTC](https://discourse.julialang.org/t/maximum-likelihood-optimization-problem-through-optim-jl/84096/14 "2022-07-12T12:36:49Z")
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.
---
### Author:  [@Jian\_ZUO](https://discourse.julialang.org/u/Jian_ZUO)
#### Post date: [July 12, 2022, 1:11pm UTC](https://discourse.julialang.org/t/maximum-likelihood-optimization-problem-through-optim-jl/84096/15 "2022-07-12T13:11:50Z")
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.