'Fit a model in which the fixed effects parameters are
penalized.  The dependence parameters are held fixed at their
estimated values in the unpenalized model.
Parameters
method : string of Penalty object
Method for regularization.  If a string, must be \'l1\'.
alpha : array-like
Scalar or vector of penalty weights.  If a scalar, the
same weight is applied to all coefficients; if a vector,
it contains a weight for each coefficient.  If method is a
Penalty object, the weights are scaled by alpha.  For L1
regularization, the weights are used directly.
ceps : positive real scalar
Fixed effects parameters smaller than this value
in magnitude are treaded as being zero.
ptol : positive real scalar
Convergence occurs when the sup norm difference
between successive values of `fe_params` is less than
`ptol`.
maxit : integer
The maximum number of iterations.
fit_kwargs : keywords
Additional keyword arguments passed to fit.
Returns
A MixedLMResults instance containing the results.
Notes
The covariance structure is not updated as the fixed effects
parameters are varied.
The algorithm used here for L1 regularization is a"shooting"
or cyclic coordinate descent algorithm.
If method is \'l1\', then `fe_pen` and `cov_pen` are used to
obtain the covariance structure, but are ignored during the
L1-penalized fitting.
References
Friedman, J. H., Hastie, T. and Tibshirani, R. Regularized
Paths for Generalized Linear Models via Coordinate
Descent. Journal of Statistical Software, 33(1) (2008)
http://www.jstatsoft.org/v33/i01/paper
http://statweb.stanford.edu/~tibs/stat315a/Supplements/fuse.pdf'

def fit_regularized(self, start_params=None, method='l1', alpha=0, ceps=0.0001, ptol=1e-06, maxit=200, **fit_kwargs):

    if (isinstance(method, string_types) and (method.lower() != 'l1')):
        raise ValueError('Invalid    regularization    method')
    if isinstance(method, Penalty):
        method.alpha = alpha
        fit_kwargs.update({'fe_pen': method})
        return self.fit(**fit_kwargs)
    if np.isscalar(alpha):
        alpha = (alpha * np.ones(self.k_fe, dtype=np.float64))
    mdf = self.fit(**fit_kwargs)
    fe_params = mdf.fe_params
    cov_re = mdf.cov_re
    vcomp = mdf.vcomp
    scale = mdf.scale
    try:
        cov_re_inv = np.linalg.inv(cov_re)
    except np.linalg.LinAlgError:
        cov_re_inv = None
    for itr in range(maxit):
        fe_params_s = fe_params.copy()
        for j in range(self.k_fe):
            if (abs(fe_params[j]) < ceps):
                continue
            fe_params[j] = 0.0
            expval = np.dot(self.exog, fe_params)
            resid_all = (self.endog - expval)
            (a, b) = (0.0, 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])
                resid = resid_all[self.row_indices[group]]
                solver = _smw_solver(scale, ex_r, ex2_r, cov_re_inv, (1 / vc_var))
                x = exog[:, j]
                u = solver(x)
                a += np.dot(u, x)
                b -= (2 * np.dot(u, resid))
            pwt1 = alpha[j]
            if (b > pwt1):
                fe_params[j] = ((- (b - pwt1)) / (2 * a))
            elif (b < (- pwt1)):
                fe_params[j] = ((- (b + pwt1)) / (2 * a))
        if (np.abs((fe_params_s - fe_params)).max() < ptol):
            break
    params_prof = mdf.params.copy()
    params_prof[0:self.k_fe] = fe_params
    scale = self.get_scale(fe_params, mdf.cov_re_unscaled, mdf.vcomp)
    hess = self.hessian(params_prof)
    pcov = (np.nan * np.ones_like(hess))
    ii = (np.abs(params_prof) > ceps)
    ii[self.k_fe:] = True
    ii = np.flatnonzero(ii)
    hess1 = hess[ii, :][:, ii]
    pcov[np.ix_(ii, ii)] = np.linalg.inv((- hess1))
    params_object = MixedLMParams.from_components(fe_params, cov_re=cov_re)
    results = MixedLMResults(self, params_prof, (pcov / scale))
    results.params_object = params_object
    results.fe_params = fe_params
    results.cov_re = cov_re
    results.scale = scale
    results.cov_re_unscaled = mdf.cov_re_unscaled
    results.method = mdf.method
    results.converged = True
    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
    return MixedLMResultsWrapper(results)
