Юлия - Оптим - Градиент в наблюдении - PullRequest
Новая
       9

Юлия - Оптим - Градиент в наблюдении

0 голосов
/ 29 августа 2018

Я занимаюсь разработкой специальной многочленной логистической модели с использованием Юлии. Он отлично работает (хотя я уверен, что это можно улучшить!)

Я написал функцию правдоподобия и использую Optim для оценки параметров и стандартных ошибок.

Я хотел бы сейчас разработать некоторые надежные оценки. Исходя из R, я бы использовал пакет сэндвич. Я не нашел никакого эквивалента в Юлии. Так что я мог бы разработать что-то конкретное, я думаю.

Для этого мне понадобится значение градиента для каждого наблюдения (строки). Я не нахожу способ сделать это с помощью Optim. (Когда я использую gradient!(func, x), я получаю сумму градиентов по строкам, а это не то, что я ищу)

Есть ли способ сделать это, используя OnceDifferentiable или TwiceDifferentiable?

В качестве альтернативы, есть ли пакет, эквивалентный R Sandwich, который бы ускользнул от моих исследований в Google?

Код, который я разработал до сих пор: 

LLIK_MNL = function(init::Array{Float64})
    b = init
    u1 = X1*b
    u2 = X2*b
    u3 = X3*b
    umax = max.(u1, u2, u3)
    umin = min.(u1, u2, u3)
    ucensor = (umax + umin)/2

    exu1 = exp.(u1 - ucensor)
    exu2 = exp.(u2 - ucensor)
    exu3 = exp.(u3 - ucensor)
    sexu = exu1 .+ exu2 .+ exu3 

    Pr=(choice1 .* exu1 + choice2 .* exu2 + choice3 .* exu3) ./ sexu
    LL = sum(log.(Pr))
    return -LL
end

func = TwiceDifferentiable(var -> LLIK_MNL(var), beta_ini)
opt = optimize(func, beta_ini)

est_MNL = Optim.minimizer(opt)
numerical_hessian = hessian!(func, est_MNL)
var_cov_matrix = inv(numerical_hessian)
temp = diag(var_cov_matrix)
t_stats = est_MNL ./ sqrt.(temp)
pp = 2 * cdf.(Normal(), -abs.(t_stats))
hcat(est_MNL, sqrt.(temp), t_stats, round.(pp, 4))
...