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))

        # Blocks for the fixed and random effects parameters.
        hess_fe = 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.
        xtvix = 0.
        xtax = [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.] * 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., ex_r, ex2_r, cov_re_inv, 1 / vc_var)

            # The residuals
            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)

                # V^{-1} * dV/d_theta
                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

        # Put the blocks together to get the Hessian.
        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
Returns the model's Hessian matrix.
dataset/ETHPy150Open statsmodels/statsmodels/statsmodels/regression/mixed_linear_model.py/MixedLM.hessian