Hessian matrix for standard errors using Optim.jl example

Specific Domains Statistics
1

I’m looking at the maximum likelihood example on the Optim.jl page and trying it on a different likelihood function (truncated normal). I’m running into an issue where the covariance matrix returned using the Optim example method is not a valid covariance matrix. However, if I directly use the ForwardDiff package I get a valid covariance matrix, leaving me quite unsure what is going wrong with the below code. Any ideas as to what is going wrong?

The example:

using Optim, Distributions
import NLSolversBase, ForwardDiff

data = [2.59 0.93 0.70 1.83 1.46 2.04 1.31 1.55 1.19 0.16 2.09 1.32 0.18 0.12 2.17 1.52 0.84 1.85 2.29 0.94 0.01 0.89 2.76 3.15 0.52 1.91 2.11 0.81 2.46 1.30];
data = adjoint(data);

function nlog_lik(y,μ,σ2)
    n = length(y)
    llik = -n*(log(sqrt(2*π)) + log(sqrt(σ2)) + log(1 - cdf(Normal(),-μ/sqrt(σ2)))) - 1/(2*σ2)*sum((y .- μ).^2)
    nllik = -llik
    return nllik
end

nlog_lik_data = TwiceDifferentiable(vars -> nlog_lik(data, vars[1], vars[2]), [1.43,1.0]; autodiff=:forward);
opt = optimize(nlog_lik_data,[1.43,1.0],store_trace=true);
param = opt.minimizer;

# not a valid variance-covariance matrix
numerical_hessian = NLSolversBase.hessian!(nlog_lik_data,param);
@show var_cov_matrix = inv(numerical_hessian)  

# valid variance-covariance matrix
nlog_lik_data2 = (vars -> nlog_lik(data, vars[1], vars[2]))
numerical_hessian2 = ForwardDiff.hessian(nlog_lik_data2, param)
@show var_cov_matrix2 = inv(numerical_hessian2)

# should be the same?
var_cov_matrix = inv(numerical_hessian) = [-0.489987 1.71397; 0.360741 -0.882709] 
var_cov_matrix2 = inv(numerical_hessian2) = [0.0683089 -0.0587528; -0.0587528 0.143764]
2

Hello!

I wrote the original tutorial and I have amended it to the following code:

using Optim
using NLSolversBase
using ForwardDiff

n = 500                             # Number of observations
nvar = 2                            # Number of variables
β = ones(nvar) * 3.0                # True coefficients
x = [ones(n) randn(n, nvar - 1)]    # X matrix of explanatory variables plus constant
ε = randn(n) * 0.5                  # Error variance
y = x * β + ε                       # Generate Data

function Log_Likelihood(X, Y, β, log_σ)
    σ = exp(log_σ)
    llike = -n/2*log(2π) - n/2* log(σ^2) - (sum((Y - X * β).^2) / (2σ^2))
    llike = -llike
end

func = TwiceDifferentiable(vars -> Log_Likelihood(x, y, vars[1:nvar], vars[nvar + 1]), ones(nvar+1); autodiff=:forward)

opt = optimize(func, ones(nvar+1))

parameters = Optim.minimizer(opt)

parameters[nvar+1] = exp(parameters[nvar+1])

numerical_hessian = NLSolversBase.hessian!(func,parameters)

zz = vars -> Log_Likelihood(x, y, vars[1:nvar], vars[nvar + 1])

numerical_hessian2 = ForwardDiff.hessian(zz, parameters)

var_cov_matrix = inv(numerical_hessian)

var_cov_matrix2 = inv(numerical_hessian2)

β = parameters[1:nvar]

temp = diag(var_cov_matrix)
temp1 = temp[1:nvar]

t_stats = β./sqrt.(temp1)

When I run the above code, I get the same exact numerical Hessian using NLSolversBase.hessian! as I do with ForwardDiff.hessian so the issue may be with another part of your code (perhaps the likelihood is miscoded?).

You might also want to post your question in the Optimization forum as that is where I received help when I had an issue with optimization in general.

1 Like