julia - Julia 中的最大似然法

Question

我正在尝试使用最大似然估计 Julia 中的正态线性模型。根据 Optim 文档中关于不更改的值，我使用以下代码通过拦截和匿名函数来模拟该过程:

using Optim

nobs = 500
nvar = 1
β = ones(nvar)*3.0
x = [ones(nobs) randn(nobs,nvar-1)]
ε = randn(nobs)*0.5
y = x*β + ε

function LL_anon(X, Y, β, σ)
  -(-length(X)*log(2π)/2 - length(X)*log(σ) - (sum((Y - X*β).^2) / (2σ^2)))
end
LL_anon(X,Y, pars) = LL_anon(X,Y, pars...)

res2 = optimize(vars -> LL_anon(x,y, vars...), [1.0,1.0]) # Start values: β=1.0, σ=1.0

这实际上恢复了参数，我收到了以下输出:

 * Algorithm: Nelder-Mead
 * Starting Point: [1.0,1.0]
 * Minimizer: [2.980587812647935,0.5108406803949835]
 * Minimum: 3.736217e+02
 * Iterations: 47
 * Convergence: true
   *  √(Σ(yᵢ-ȳ)²)/n < 1.0e-08: true
   * Reached Maximum Number of Iterations: false
 * Objective Calls: 92

但是，当我尝试设置 nvar = 2，即一个截距加上一个额外的协变量时，我收到以下错误消息:

MethodError: no method matching LL_anon(::Array{Float64,2}, ::Array{Float64,1}, ::Float64, ::Float64, ::Float64)
Closest candidates are:
  LL_anon(::Any, ::Any, ::Any, ::Any) at In[297]:2
  LL_anon(::Array{Float64,1}, ::Array{Float64,1}, ::Any, ::Any) at In[113]:2
  LL_anon(::Any, ::Any, ::Any) at In[297]:4
  ...

Stacktrace:
 [1] (::##245#246)(::Array{Float64,1}) at .\In[299]:1
 [2] value!!(::NLSolversBase.NonDifferentiable{Float64,Array{Float64,1},Val{false}}, ::Array{Float64,1}) at C:\Users\dolacomb\.julia\v0.6\NLSolversBase\src\interface.jl:9
 [3] initial_state(::Optim.NelderMead{Optim.AffineSimplexer,Optim.AdaptiveParameters}, ::Optim.Options{Float64,Void}, ::NLSolversBase.NonDifferentiable{Float64,Array{Float64,1},Val{false}}, ::Array{Float64,1}) at C:\Users\dolacomb\.julia\v0.6\Optim\src\multivariate/solvers/zeroth_order\nelder_mead.jl:136
 [4] optimize(::NLSolversBase.NonDifferentiable{Float64,Array{Float64,1},Val{false}}, ::Array{Float64,1}, ::Optim.NelderMead{Optim.AffineSimplexer,Optim.AdaptiveParameters}, ::Optim.Options{Float64,Void}) at C:\Users\dolacomb\.julia\v0.6\Optim\src\multivariate/optimize\optimize.jl:25
 [5] #optimize#151(::Array{Any,1}, ::Function, ::Tuple{##245#246}, ::Array{Float64,1}) at C:\Users\dolacomb\.julia\v0.6\Optim\src\multivariate/optimize\interface.jl:62
 [6] #optimize#148(::Array{Any,1}, ::Function, ::Function, ::Array{Float64,1}) at C:\Users\dolacomb\.julia\v0.6\Optim\src\multivariate/optimize\interface.jl:52
 [7] optimize(::Function, ::Array{Float64,1}) at C:\Users\dolacomb\.julia\v0.6\Optim\src\multivariate/optimize\interface.jl:52

我不确定为什么添加一个额外的变量会导致这个问题，但它似乎是一个类型不稳定的问题。

第二个问题是，当我使用我原来的工作示例并将起始值设置为 [2.0,2.0] 时，我收到以下错误消息:

log will only return a complex result if called with a complex argument. Try log(complex(x)).

Stacktrace:
 [1] nan_dom_err at .\math.jl:300 [inlined]
 [2] log at .\math.jl:419 [inlined]
 [3] LL_anon(::Array{Float64,2}, ::Array{Float64,1}, ::Float64, ::Float64) at .\In[302]:2
 [4] (::##251#252)(::Array{Float64,1}) at .\In[304]:1
 [5] value(::NLSolversBase.NonDifferentiable{Float64,Array{Float64,1},Val{false}}, ::Array{Float64,1}) at C:\Users\dolacomb\.julia\v0.6\NLSolversBase\src\interface.jl:19
 [6] update_state!(::NLSolversBase.NonDifferentiable{Float64,Array{Float64,1},Val{false}}, ::Optim.NelderMeadState{Array{Float64,1},Float64,Array{Float64,1}}, ::Optim.NelderMead{Optim.AffineSimplexer,Optim.AdaptiveParameters}) at C:\Users\dolacomb\.julia\v0.6\Optim\src\multivariate/solvers/zeroth_order\nelder_mead.jl:193
 [7] optimize(::NLSolversBase.NonDifferentiable{Float64,Array{Float64,1},Val{false}}, ::Array{Float64,1}, ::Optim.NelderMead{Optim.AffineSimplexer,Optim.AdaptiveParameters}, ::Optim.Options{Float64,Void}, ::Optim.NelderMeadState{Array{Float64,1},Float64,Array{Float64,1}}) at C:\Users\dolacomb\.julia\v0.6\Optim\src\multivariate/optimize\optimize.jl:51
 [8] optimize(::NLSolversBase.NonDifferentiable{Float64,Array{Float64,1},Val{false}}, ::Array{Float64,1}, ::Optim.NelderMead{Optim.AffineSimplexer,Optim.AdaptiveParameters}, ::Optim.Options{Float64,Void}) at C:\Users\dolacomb\.julia\v0.6\Optim\src\multivariate/optimize\optimize.jl:25
 [9] #optimize#151(::Array{Any,1}, ::Function, ::Tuple{##251#252}, ::Array{Float64,1}) at C:\Users\dolacomb\.julia\v0.6\Optim\src\multivariate/optimize\interface.jl:62

同样，我不确定为什么会发生这种情况，并且由于起始值非常重要，我想知道如何解决这个问题并且它们与真实值相差不远。

如有任何帮助，我们将不胜感激!

Answer 1

飞溅是导致问题的原因。例如。它将 [1, 2, 3] 转换为三个参数，而您的函数只接受两个参数。

使用以下调用:

res2 = optimize(vars -> LL_anon(x,y, vars[1:end-1], vars[end]), [1.0,1.0,1.0])

您可以从代码中删除以下行

LL_anon(X,Y, pars) = LL_anon(X,Y, pars...)

或替换为:

LL_anon(X,Y, pars) = LL_anon(X,Y, pars[1:end-1], pars[end])

但优化例程不会调用它，除非您将调用更改为:

res2 = optimize(vars -> LL_anon(x,y, vars), [1.0,1.0,1.0])

最后 - 为了获得这段代码的良好性能，我建议将它全部包装在一个函数中。

编辑:现在我看到第二个问题。原因是 σ 在优化过程中可能变为负值，然后 log(σ) 失败。在这种情况下，最简单的做法是像这样获取 log(abs(σ))):

function LL_anon(X, Y, β, σ)
    -(-length(X)*log(2π)/2 - length(X)*log(abs(σ)) - (sum((Y - X*β).^2) / (2σ^2)))
end

当然，您必须将 σ 的绝对值作为解决方案，因为您可能会从优化程序中获得负值。

一种更简洁的方法是优化例如log(σ) 不是 σ 这样的:

function LL_anon(X, Y, β, logσ)
    -(-length(X)*log(2π)/2 - length(X)*logσ - (sum((Y - X*β).^2) / (2(exp(logσ))^2)))
end

但是你必须在优化完成后使用 exp(logσ) 来恢复 σ。