'Returns the estimated error variance based on given estimates
of the slopes and random effects covariance matrix.
Parameters
fe_params : array-like
The regression slope estimates
cov_re : 2d array-like
Estimate of the random effects covariance matrix
vcomp : array-like
Estimate of the variance components
Returns
scale : float
The estimated error variance.'

def get_scale(self, fe_params, cov_re, vcomp):

    try:
        cov_re_inv = np.linalg.inv(cov_re)
    except np.linalg.LinAlgError:
        cov_re_inv = None
    qf = 0.0
    for (group_ix, group) in enumerate(self.group_labels):
        vc_var = self._expand_vcomp(vcomp, group)
        exog = self.exog_li[group_ix]
        (ex_r, ex2_r) = (self._aex_r[group_ix], self._aex_r2[group_ix])
        solver = _smw_solver(1.0, ex_r, ex2_r, cov_re_inv, (1 / vc_var))
        resid = self.endog_li[group_ix]
        if (self.k_fe > 0):
            expval = np.dot(exog, fe_params)
            resid = (resid - expval)
        mat = solver(resid)
        qf += np.dot(resid, mat)
    if self.reml:
        qf /= (self.n_totobs - self.k_fe)
    else:
        qf /= self.n_totobs
    return qf
