Hessians to compute variances¶

We have our point estimates, but now we want variances of the estimates
We can approximate them (under some asymptotic conditions) using the Hessian at the (negative) MLE minimum
Think of the the Hessian as a matrix representing the curvature at the minimum point. If we move slightly away from that minimum, how much change do we see in the likelihood?
- high curvature (think of a tight peak) => low variance of estimate.
- low curvature (think of a broad hill) => high variance of estimate.

To extract the parameter estimates, we invert the Hessian at the MLE, $\hat{\theta}$, and take the diagonal:

$$\text{Var}(\hat{\theta}_i) = \text{diag}(H(\hat{\theta})^{-1})_i $$

The Hessian is the matrix of all second derivatives. autograd has this built in:

In [1]:

from autograd import numpy as np
from autograd import elementwise_grad, value_and_grad, hessian
from scipy.optimize import minimize

In [2]:

T = (np.random.exponential(size=1000)/1.5) ** 2.3
E = np.random.binomial(1, 0.95, size=1000)

In [3]:

# seen all this...
def cumulative_hazard(params, t):
    lambda_, rho_ = params
    return (t / lambda_) ** rho_

hazard = elementwise_grad(cumulative_hazard, argnum=1)

def log_hazard(params, t):
    return np.log(hazard(params, t))

def log_likelihood(params, t, e):
    return np.sum(e * log_hazard(params, t)) - np.sum(cumulative_hazard(params, t))

def negative_log_likelihood(params, t, e):
    return -log_likelihood(params, t, e)

from autograd import value_and_grad

results = minimize(
        value_and_grad(negative_log_likelihood), 
        x0 = np.array([1.0, 1.0]),
        method=None, 
        args=(T, E),
        jac=True,
        bounds=((0.00001, None), (0.00001, None)))

print(results)

      fun: 249.99413038558657
 hess_inv: <2x2 LbfgsInvHessProduct with dtype=float64>
      jac: array([-0.00046244,  0.00022406])
  message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
     nfev: 11
      nit: 7
   status: 0
  success: True
        x: array([0.46459167, 0.44953271])

In [4]:

estimates_ = results.x

# Note: this will produce a matrix
hessian(negative_log_likelihood)(estimates_, T, E)

Out[4]:

array([[ 887.54177234, -758.67183105],
       [-758.67183105, 8295.5360413 ]])

In [5]:

H = hessian(negative_log_likelihood)(estimates_, T, E)
variance_ = np.diag(np.linalg.inv(H))
print(variance_)

[0.00122226 0.00013077]

In [6]:

std_error = np.sqrt(variance_)
print(std_error)

[0.03496082 0.01143546]

In [7]:

print("lambda_  %.3f (±%.2f)" % (estimates_[0], std_error[0]))
print("rho_     %.3f (±%.2f)" % (estimates_[1], std_error[1]))

lambda_  0.465 (±0.03)
rho_     0.450 (±0.01)

From here you can construct confidence intervals (CIs), and such. Let's move to Part 6, which is where you connect these abstract parameters to business logic.