Miximum Likelihood - using Optim package

Question

Dear users of the language julia. I have a problem when using the optimize function of the Optim package. What is the error of the code below?

using Optim
using Distributions

rng = MersenneTwister(1234);

d = Weibull(1,1)
x = rand(d,1000)


function pdf_weibull(x, lambda, k)  
    k/lambda * (x/lambda).^(k-1) * exp((-x/lambda)^k)

end


function obj(x::Vector, lambda, k)
    soma = 0
    for i in x
        soma = soma + log(pdf_weibull(i,lambda,k))
    end
    -soma
end

obj(x, pars) = obj(x, pars...)

optimize(vars -> obj(x, vars...), [1.0,1.0])

Output

julia> optimize(vars -> obj(x, vars...), [1.0,1.0])
ERROR: DomainError:
Exponentiation yielding a complex result requires a complex argument.
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.
Stacktrace:
 [1] nan_dom_err at ./math.jl:300 [inlined]
 [2] ^ at ./math.jl:699 [inlined]
 [3] (::##2#4)(::Float64, ::Float64, ::Float64) at ./<missing>:0
 [4] pdf_weibull(::Float64, ::Float64, ::Float64) at ./REPL[6]:2
 [5] obj(::Array{Float64,1}, ::Float64, ::Float64) at ./REPL[7]:4
 [6] (::##5#6)(::Array{Float64,1}) at ./REPL[11]:1
 [7] value(::NLSolversBase.NonDifferentiable{Float64,Array{Float64,1},Val{false}}, ::Array{Float64,1}) at /home/pedro/.julia/v0.6/NLSolversBase/src/interface.jl:19
 [8] initial_state(::Optim.NelderMead{Optim.AffineSimplexer,Optim.AdaptiveParameters}, ::Optim.Options{Float64,Void}, ::NLSolversBase.NonDifferentiable{Float64,Array{Float64,1},Val{false}}, ::Array{Float64,1}) at /home/pedro/.julia/v0.6/Optim/src/multivariate/solvers/zeroth_order/nelder_mead.jl:139
 [9] optimize(::NLSolversBase.NonDifferentiable{Float64,Array{Float64,1},Val{false}}, ::Array{Float64,1}, ::Optim.NelderMead{Optim.AffineSimplexer,Optim.AdaptiveParameters}, ::Optim.Options{Float64,Void}) at /home/pedro/.julia/v0.6/Optim/src/multivariate/optimize/optimize.jl:25
 [10] #optimize#151(::Array{Any,1}, ::Function, ::Tuple{##5#6}, ::Array{Float64,1}) at /home/pedro/.julia/v0.6/Optim/src/multivariate/optimize/interface.jl:62
 [11] #optimize#148(::Array{Any,1}, ::Function, ::Function, ::Array{Float64,1}) at /home/pedro/.julia/v0.6/Optim/src/multivariate/optimize/interface.jl:52
 [12] optimize(::Function, ::Array{Float64,1}) at /home/pedro/.julia/v0.6/Optim/src/multivariate/optimize/interface.jl:52
 [13] macro expansion at ./REPL.jl:97 [inlined]
 [14] (::Base.REPL.##1#2{Base.REPL.REPLBackend})() at ./event.jl:73

It is a simple problem to obtain the maximum likelihood estimates of the parameters that index the weibull distribution.

Best regards.

What happens when you run this? Do you get an error message? If so, please tell us what it is. Do you get unexpected results? If so, what are they and what do you expect to get? — divibisan
I expect the maximum likelihood estimates which in this case will estimate lambda and k by approximately 1. — Pedro Rafael
Have you tried making the changes that the error message suggests? If so, show us what you've tried — divibisan
I have not tried because I just do not understand why the first error message refers to the complex numbers. Could you explain? — Pedro Rafael

Bogumił Kamiński Bogumił Kamiński · Accepted Answer · 2018-06-01T06:40:42

The reason of your problem is that your definition of pdf_weibull is incorrect. Here is a corrected definition:

function pdf_weibull(x, lambda, k)
    k/lambda * (x/lambda)^(k-1) * exp(-(x/lambda)^k)
end

Note that I have moved - sign in exp part of the expression. If you change this all will work as expected.

Now - why does Julia complain with DomainError. The reason is that because of the error in your code you try to calculate the value of something like (-1.0)^0.5. In Julia ^ is implemented in type stable way. This means, in particular, that when it is passed Float64 as both arguments it guarantees to return Float64 or throw an error. Clearly (-1.0)^0.5 cannot be computed in real domain - that is why an error is thrown. If you passed (-1+0im)^0.5 then there would be no error as we are passing a complex number to ^ so the result can also be complex number in a type stable way.

Miximum Likelihood - using Optim package

1 Answers