'Fit a linear mixed model to the data.
Parameters
start_params: array-like or MixedLMParams
Starting values for the profile log-likelihood.  If not a
`MixedLMParams` instance, this should be an array
containing the packed parameters for the profile
log-likelihood, including the fixed effects
parameters.
reml : bool
If true, fit according to the REML likelihood, else
fit the standard likelihood using ML.
cov_pen : CovariancePenalty object
A penalty for the random effects covariance matrix
fe_pen : Penalty object
A penalty on the fixed effects
free : MixedLMParams object
If not `None`, this is a mask that allows parameters to be
held fixed at specified values.  A 1 indicates that the
correspondinig parameter is estimated, a 0 indicates that
it is fixed at its starting value.  Setting the `cov_re`
component to the identity matrix fits a model with
independent random effects.  Note that some optimization
methods do not respect this constraint (bfgs and lbfgs both
work).
full_output : bool
If true, attach iteration history to results
method : string
Optimization method.
Returns
A MixedLMResults instance.'

def fit(self, start_params=None, reml=True, niter_sa=0, do_cg=True, fe_pen=None, cov_pen=None, free=None, full_output=False, method='bfgs', **kwargs):

    _allowed_kwargs = ['gtol', 'maxiter']
    for x in kwargs.keys():
        if (x not in _allowed_kwargs):
            raise ValueError(('Argument    %s    not    allowed    for    MixedLM.fit' % x))
    if (method.lower() in ['newton', 'ncg']):
        raise ValueError(('method    %s    not    available    for    MixedLM' % method))
    self.reml = reml
    self.cov_pen = cov_pen
    self.fe_pen = fe_pen
    self._freepat = free
    if full_output:
        hist = []
    else:
        hist = None
    success = False
    if (start_params is None):
        params = MixedLMParams(self.k_fe, self.k_re, self.k_vc)
        params.fe_params = np.zeros(self.k_fe)
        params.cov_re = np.eye(self.k_re)
        params.vcomp = np.ones(self.k_vc)
    elif isinstance(start_params, MixedLMParams):
        params = start_params
    elif (len(start_params) == ((self.k_fe + self.k_re2) + self.k_vc)):
        params = MixedLMParams.from_packed(start_params, self.k_fe, self.k_re, self.use_sqrt, has_fe=True)
    elif (len(start_params) == (self.k_re2 + self.k_vc)):
        params = MixedLMParams.from_packed(start_params, self.k_fe, self.k_re, self.use_sqrt, has_fe=False)
    else:
        raise ValueError('invalid    start_params')
    if do_cg:
        kwargs['retall'] = (hist is not None)
        if ('disp' not in kwargs):
            kwargs['disp'] = False
        packed = params.get_packed(use_sqrt=self.use_sqrt, has_fe=False)
        for rep in range(5):
            rslt = super(MixedLM, self).fit(start_params=packed, skip_hessian=True, method=method, **kwargs)
            if rslt.mle_retvals['converged']:
                break
            packed = rslt.params
        params = np.atleast_1d(rslt.params)
        if (hist is not None):
            hist.append(rslt.mle_retvals)
    converged = rslt.mle_retvals['converged']
    if (not converged):
        msg = 'Gradient    optimization    failed.'
        warnings.warn(msg, ConvergenceWarning)
    params = MixedLMParams.from_packed(params, self.k_fe, self.k_re, use_sqrt=self.use_sqrt, has_fe=False)
    cov_re_unscaled = params.cov_re
    vcomp_unscaled = params.vcomp
    fe_params = self.get_fe_params(cov_re_unscaled, vcomp_unscaled)
    params.fe_params = fe_params
    scale = self.get_scale(fe_params, cov_re_unscaled, vcomp_unscaled)
    cov_re = (scale * cov_re_unscaled)
    vcomp = (scale * vcomp_unscaled)
    if (((self.k_re > 0) and (np.min(np.abs(np.diag(cov_re))) < 0.01)) or ((self.k_vc > 0) and (np.min(np.abs(vcomp)) < 0.01))):
        msg = 'The    MLE    may    be    on    the    boundary    of    the    parameter    space.'
        warnings.warn(msg, ConvergenceWarning)
    hess = self.hessian(params)
    hess_diag = np.diag(hess)
    if (free is not None):
        pcov = np.zeros_like(hess)
        pat = self._freepat.get_packed(use_sqrt=False, has_fe=True)
        ii = np.flatnonzero(pat)
        hess_diag = hess_diag[ii]
        if (len(ii) > 0):
            hess1 = hess[np.ix_(ii, ii)]
            pcov[np.ix_(ii, ii)] = np.linalg.inv((- hess1))
    else:
        pcov = np.linalg.inv((- hess))
    if np.any((hess_diag >= 0)):
        msg = 'The    Hessian    matrix    at    the    estimated    parameter    values    is    not    positive    definite.'
        warnings.warn(msg, ConvergenceWarning)
    params_packed = params.get_packed(use_sqrt=False, has_fe=True)
    results = MixedLMResults(self, params_packed, (pcov / scale))
    results.params_object = params
    results.fe_params = fe_params
    results.cov_re = cov_re
    results.vcomp = vcomp
    results.scale = scale
    results.cov_re_unscaled = cov_re_unscaled
    results.method = ('REML' if self.reml else 'ML')
    results.converged = converged
    results.hist = hist
    results.reml = self.reml
    results.cov_pen = self.cov_pen
    results.k_fe = self.k_fe
    results.k_re = self.k_re
    results.k_re2 = self.k_re2
    results.k_vc = self.k_vc
    results.use_sqrt = self.use_sqrt
    results.freepat = self._freepat
    return MixedLMResultsWrapper(results)
