'Evaluate the (profile) log-likelihood of the linear mixed
effects model.
Parameters
params : MixedLMParams, or array-like.
The parameter value.  If array-like, must be a packed
parameter vector containing only the covariance
parameters.
profile_fe : boolean
If True, replace the provided value of `fe_params` with
the GLS estimates.
Returns
The log-likelihood value at `params`.
Notes
The scale parameter `scale` is always profiled out of the
log-likelihood.  In addition, if `profile_fe` is true the
fixed effects parameters are also profiled out.'

def loglike(self, params, profile_fe=True):

    if (type(params) is not MixedLMParams):
        params = MixedLMParams.from_packed(params, self.k_fe, self.k_re, self.use_sqrt, has_fe=False)
    cov_re = params.cov_re
    vcomp = params.vcomp
    if profile_fe:
        fe_params = self.get_fe_params(cov_re, vcomp)
    else:
        fe_params = params.fe_params
    if (self.k_re > 0):
        try:
            cov_re_inv = np.linalg.inv(cov_re)
        except np.linalg.LinAlgError:
            cov_re_inv = None
        (_, cov_re_logdet) = np.linalg.slogdet(cov_re)
    else:
        cov_re_inv = np.zeros((0, 0))
        cov_re_logdet = 0
    expval = np.dot(self.exog, fe_params)
    resid_all = (self.endog - expval)
    likeval = 0.0
    if ((self.cov_pen is not None) and (self.k_re > 0)):
        likeval -= self.cov_pen.func(cov_re, cov_re_inv)
    if (self.fe_pen is not None):
        likeval -= self.fe_pen.func(fe_params)
    (xvx, qf) = (0.0, 0.0)
    for (k, group) in enumerate(self.group_labels):
        vc_var = self._expand_vcomp(vcomp, group)
        cov_aug_logdet = (cov_re_logdet + np.sum(np.log(vc_var)))
        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 = resid_all[self.row_indices[group]]
        ld = _smw_logdet(1.0, ex_r, ex2_r, cov_re_inv, (1 / vc_var), cov_aug_logdet)
        likeval -= (ld / 2.0)
        u = solver(resid)
        qf += np.dot(resid, u)
        if self.reml:
            mat = solver(exog)
            xvx += np.dot(exog.T, mat)
    if self.reml:
        likeval -= (((self.n_totobs - self.k_fe) * np.log(qf)) / 2.0)
        (_, ld) = np.linalg.slogdet(xvx)
        likeval -= (ld / 2.0)
        likeval -= (((self.n_totobs - self.k_fe) * np.log((2 * np.pi))) / 2.0)
        likeval += (((self.n_totobs - self.k_fe) * np.log((self.n_totobs - self.k_fe))) / 2.0)
        likeval -= ((self.n_totobs - self.k_fe) / 2.0)
    else:
        likeval -= ((self.n_totobs * np.log(qf)) / 2.0)
        likeval -= ((self.n_totobs * np.log((2 * np.pi))) / 2.0)
        likeval += ((self.n_totobs * np.log(self.n_totobs)) / 2.0)
        likeval -= (self.n_totobs / 2.0)
    return likeval
