'Returns the model\'s Hessian matrix.
Calculates the Hessian matrix for the linear mixed effects
model with respect to the parameterization in which the
covariance matrix is represented directly (without square-root
transformation).
Parameters
params : MixedLMParams or array-like
The model parameters at which the Hessian is calculated.
If array-like, must contain the packed parameters in a
form that is compatible with this model instance.
Returns
hess : 2d ndarray
The Hessian matrix, evaluated at `params`.'

def hessian(self, params):

    if (type(params) is not MixedLMParams):
        params = MixedLMParams.from_packed(params, self.k_fe, self.k_re, use_sqrt=self.use_sqrt, has_fe=True)
    fe_params = params.fe_params
    vcomp = params.vcomp
    cov_re = params.cov_re
    if (self.k_re > 0):
        cov_re_inv = np.linalg.inv(cov_re)
    else:
        cov_re_inv = np.empty((0, 0))
    hess_fe = 0.0
    hess_re = np.zeros(((self.k_re2 + self.k_vc), (self.k_re2 + self.k_vc)))
    hess_fere = np.zeros(((self.k_re2 + self.k_vc), self.k_fe))
    fac = self.n_totobs
    if self.reml:
        fac -= self.exog.shape[1]
    rvir = 0.0
    xtvix = 0.0
    xtax = ([0.0] * (self.k_re2 + self.k_vc))
    m = (self.k_re2 + self.k_vc)
    B = np.zeros(m)
    D = np.zeros((m, m))
    F = [([0.0] * m) for k in range(m)]
    for (k, group) in enumerate(self.group_labels):
        vc_var = self._expand_vcomp(vcomp, group)
        exog = self.exog_li[k]
        (ex_r, ex2_r) = (self._aex_r[k], self._aex_r2[k])
        solver = _smw_solver(1.0, ex_r, ex2_r, cov_re_inv, (1 / vc_var))
        resid = self.endog_li[k]
        if (self.k_fe > 0):
            expval = np.dot(exog, fe_params)
            resid = (resid - expval)
        viexog = solver(exog)
        xtvix += np.dot(exog.T, viexog)
        vir = solver(resid)
        rvir += np.dot(resid, vir)
        for (jj1, matl1, matr1, vsl1, vsr1, sym1) in self._gen_dV_dPar(ex_r, solver, group):
            ul = _dot(viexog.T, matl1)
            ur = _dot(matr1.T, vir)
            hess_fere[jj1, :] += np.dot(ul, ur)
            if (not sym1):
                ul = _dot(viexog.T, matr1)
                ur = _dot(matl1.T, vir)
                hess_fere[jj1, :] += np.dot(ul, ur)
            if self.reml:
                ul = _dot(viexog.T, matl1)
                ur = (ul if sym1 else np.dot(viexog.T, matr1))
                ulr = _dot(ul, ur.T)
                xtax[jj1] += ulr
                if (not sym1):
                    xtax[jj1] += ulr.T
            ul = _dot(vir, matl1)
            ur = (ul if sym1 else _dot(vir, matr1))
            B[jj1] += (np.dot(ul, ur) * (1 if sym1 else 2))
            E = [(vsl1, matr1)]
            if (not sym1):
                E.append((vsr1, matl1))
            for (jj2, matl2, matr2, vsl2, vsr2, sym2) in self._gen_dV_dPar(ex_r, solver, group, jj1):
                re = sum([_multi_dot_three(matr2.T, x[0], x[1].T) for x in E])
                vt = (2 * _dot(_multi_dot_three(vir[None, :], matl2, re), vir[:, None]))
                if (not sym2):
                    le = sum([_multi_dot_three(matl2.T, x[0], x[1].T) for x in E])
                    vt += (2 * _dot(_multi_dot_three(vir[None, :], matr2, le), vir[:, None]))
                D[(jj1, jj2)] += vt
                if (jj1 != jj2):
                    D[(jj2, jj1)] += vt
                rt = (_dotsum(vsl2, re.T) / 2)
                if (not sym2):
                    rt += (_dotsum(vsr2, le.T) / 2)
                hess_re[(jj1, jj2)] += rt
                if (jj1 != jj2):
                    hess_re[(jj2, jj1)] += rt
                if self.reml:
                    ev = sum([_dot(x[0], _dot(x[1].T, viexog)) for x in E])
                    u1 = _dot(viexog.T, matl2)
                    u2 = _dot(matr2.T, ev)
                    um = np.dot(u1, u2)
                    F[jj1][jj2] += (um + um.T)
                    if (not sym2):
                        u1 = np.dot(viexog.T, matr2)
                        u2 = np.dot(matl2.T, ev)
                        um = np.dot(u1, u2)
                        F[jj1][jj2] += (um + um.T)
    hess_fe -= ((fac * xtvix) / rvir)
    hess_re = (hess_re - ((0.5 * fac) * ((D / rvir) - (np.outer(B, B) / (rvir ** 2)))))
    hess_fere = (((- fac) * hess_fere) / rvir)
    if self.reml:
        QL = [np.linalg.solve(xtvix, x) for x in xtax]
        for j1 in range((self.k_re2 + self.k_vc)):
            for j2 in range((j1 + 1)):
                a = _dotsum(QL[j1].T, QL[j2])
                a -= np.trace(np.linalg.solve(xtvix, F[j1][j2]))
                a *= 0.5
                hess_re[(j1, j2)] += a
                if (j1 > j2):
                    hess_re[(j2, j1)] += a
    m = ((self.k_fe + self.k_re2) + self.k_vc)
    hess = np.zeros((m, m))
    hess[0:self.k_fe, 0:self.k_fe] = hess_fe
    hess[0:self.k_fe, self.k_fe:] = hess_fere.T
    hess[self.k_fe:, 0:self.k_fe] = hess_fere
    hess[self.k_fe:, self.k_fe:] = hess_re
    return hess
