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!

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.

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.

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.

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

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.

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.