Максимальное правдоподобие - использование пакета Optim - PullRequest
Новая
       9

Максимальное правдоподобие - использование пакета Optim

0 голосов
/ 01 июня 2018

Уважаемые пользователи языка Юлия .У меня проблема при использовании функции optimize пакета Optim .В чем заключается ошибка приведенного ниже кода? 

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])

Выход 

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

Это простая задача - получить оценки максимального правдоподобия параметров, которые индексируютраспределение Вейбулла.

С уважением.

1 Ответ

0 голосов
/ 01 июня 2018

Причина вашей проблемы в том, что ваше определение pdf_weibull неверно.Вот исправленное определение: 

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

Обратите внимание, что я переместил - знак exp части выражения.Если вы измените это, все будет работать как положено.

Теперь - почему Джулия жалуется на DomainError.Причина в том, что из-за ошибки в вашем коде вы пытаетесь вычислить значение что-то вроде (-1.0)^0.5.У Юлии ^ реализован типоустойчивый способ.Это, в частности, означает, что при передаче Float64 в качестве обоих аргументов он гарантирует возврат Float64 или выдачу ошибки.Ясно, что (-1.0)^0.5 не может быть вычислено в реальном домене, поэтому выдается ошибка.Если вы передадите (-1+0im)^0.5, то ошибки не будет, поскольку мы передаем комплексное число в ^, поэтому результатом также может быть комплексное число в стабильном виде.

...