"""
Limited dependent variable and qualitative variables.

Includes binary outcomes, count data, (ordered) ordinal data and limited
dependent variables.

General References
--------------------

A.C. Cameron and P.K. Trivedi.  `Regression Analysis of Count Data`.
    Cambridge, 1998

G.S. Madalla. `Limited-Dependent and Qualitative Variables in Econometrics`.
    Cambridge, 1983.

W. Greene. `Econometric Analysis`. Prentice Hall, 5th. edition. 2003.
"""

__all__ = [
    "CountModel",
    "GeneralizedPoisson",
    "Logit",
    "MNLogit",
    "NegativeBinomial",
    "NegativeBinomialP",
    "Poisson",
    "Probit",
]

from statsmodels.compat.pandas import Appender

import warnings

import numpy as np
from pandas import MultiIndex, get_dummies
from scipy import special, stats
from scipy.special import digamma, gammaln, loggamma, polygamma
from scipy.stats import nbinom

from statsmodels.base import _prediction_inference as pred
from statsmodels.base._constraints import fit_constrained_wrap
import statsmodels.base._parameter_inference as pinfer
from statsmodels.base.data import handle_data  # for mnlogit
from statsmodels.base.l1_slsqp import fit_l1_slsqp
import statsmodels.base.model as base
import statsmodels.base.wrapper as wrap
from statsmodels.distributions import genpoisson_p
import statsmodels.regression.linear_model as lm
from statsmodels.tools import data as data_tools, tools
from statsmodels.tools.decorators import cache_readonly
from statsmodels.tools.numdiff import approx_fprime_cs
from statsmodels.tools.sm_exceptions import (
    PerfectSeparationError,
    PerfectSeparationWarning,
    SpecificationWarning,
)

try:
    import cvxopt  # noqa:F401

    have_cvxopt = True
except ImportError:
    have_cvxopt = False


# TODO: When we eventually get user-settable precision, we need to change
#       this
FLOAT_EPS = np.finfo(float).eps

# Limit for exponentials to avoid overflow
EXP_UPPER_LIMIT = np.log(np.finfo(np.float64).max) - 1.0

# TODO: add options for the parameter covariance/variance
#       ie., OIM, EIM, and BHHH see Green 21.4

_discrete_models_docs = """
"""

_discrete_results_docs = """
    %(one_line_description)s

    Parameters
    ----------
    model : A DiscreteModel instance
    params : array_like
        The parameters of a fitted model.
    hessian : array_like
        The hessian of the fitted model.
    scale : float
        A scale parameter for the covariance matrix.

    Attributes
    ----------
    df_resid : float
        See model definition.
    df_model : float
        See model definition.
    llf : float
        Value of the loglikelihood
    %(extra_attr)s"""

_l1_results_attr = """    nnz_params : int
        The number of nonzero parameters in the model.  Train with
        trim_params == True or else numerical error will distort this.
    trimmed : bool array
        trimmed[i] == True if the ith parameter was trimmed from the model."""

_get_start_params_null_docs = """
Compute one-step moment estimator for null (constant-only) model

This is a preliminary estimator used as start_params.

Returns
-------
params : ndarray
    parameter estimate based one one-step moment matching

"""

_check_rank_doc = """
    check_rank : bool
        Check exog rank to determine model degrees of freedom. Default is
        True. Setting to False reduces model initialization time when
        exog.shape[1] is large.
    """


# helper for MNLogit (will be generally useful later)
def _numpy_to_dummies(endog):
    if endog.ndim == 2 and endog.dtype.kind not in ["S", "O"]:
        endog_dummies = endog
        ynames = range(endog.shape[1])
    else:
        dummies = get_dummies(endog, drop_first=False)
        ynames = {i: dummies.columns[i] for i in range(dummies.shape[1])}
        endog_dummies = np.asarray(dummies, dtype=float)

        return endog_dummies, ynames

    return endog_dummies, ynames


def _pandas_to_dummies(endog):
    if endog.ndim == 2:
        if endog.shape[1] == 1:
            yname = endog.columns[0]
            endog_dummies = get_dummies(endog.iloc[:, 0])
        else:  # assume already dummies
            yname = "y"
            endog_dummies = endog
    else:
        yname = endog.name
        if yname is None:
            yname = "y"
        endog_dummies = get_dummies(endog)
    ynames = endog_dummies.columns.tolist()

    return endog_dummies, ynames, yname


def _validate_l1_method(method):
    """
    As of 0.10.0, the supported values for `method` in `fit_regularized`
    are "l1" and "l1_cvxopt_cp".  If an invalid value is passed, raise
    with a helpful error message

    Parameters
    ----------
    method : str

    Raises
    ------
    ValueError
    """
    if method not in ["l1", "l1_cvxopt_cp"]:
        raise ValueError(
            "`method` = {method} is not supported, use either "
            '"l1" or "l1_cvxopt_cp"'.format(method=method)
        )


# Private Model Classes


class DiscreteModel(base.LikelihoodModel):
    """
    Abstract class for discrete choice models.

    This class does not do anything itself but lays out the methods and
    call signature expected of child classes in addition to those of
    statsmodels.model.LikelihoodModel.
    """

    def __init__(self, endog, exog, check_rank=True, **kwargs):
        self._check_rank = check_rank
        super().__init__(endog, exog, **kwargs)
        self.raise_on_perfect_prediction = False  # keep for backwards compat
        self.k_extra = 0

    def initialize(self):
        """
        Initialize is called by
        statsmodels.model.LikelihoodModel.__init__
        and should contain any preprocessing that needs to be done for a model.
        """
        if self._check_rank:
            # assumes constant
            rank = tools.matrix_rank(self.exog, method="qr")
        else:
            # If rank check is skipped, assume full
            rank = self.exog.shape[1]
        self.df_model = float(rank - 1)
        self.df_resid = float(self.exog.shape[0] - rank)

    def cdf(self, X):
        """
        The cumulative distribution function of the model.
        """
        raise NotImplementedError

    def pdf(self, X):
        """
        The probability density (mass) function of the model.
        """
        raise NotImplementedError

    def _check_perfect_pred(self, params, *args):
        endog = self.endog
        fittedvalues = self.predict(params)
        if np.allclose(fittedvalues - endog, 0):
            if self.raise_on_perfect_prediction:
                # backwards compatibility for attr raise_on_perfect_prediction
                msg = "Perfect separation detected, results not available"
                raise PerfectSeparationError(msg)
            else:
                msg = (
                    "Perfect separation or prediction detected, "
                    "parameter may not be identified"
                )
                warnings.warn(msg, PerfectSeparationWarning, stacklevel=2)

    @Appender(base.LikelihoodModel.fit.__doc__)
    def fit(
        self,
        start_params=None,
        method="newton",
        maxiter=35,
        full_output=1,
        disp=1,
        callback=None,
        **kwargs,
    ):
        """
        Fit the model using maximum likelihood.

        The rest of the docstring is from
        statsmodels.base.model.LikelihoodModel.fit
        """
        if callback is None:
            callback = self._check_perfect_pred
        else:
            pass  # TODO: make a function factory to have multiple call-backs

        mlefit = super().fit(
            start_params=start_params,
            method=method,
            maxiter=maxiter,
            full_output=full_output,
            disp=disp,
            callback=callback,
            **kwargs,
        )

        return mlefit  # It is up to subclasses to wrap results

    def fit_regularized(
        self,
        start_params=None,
        method="l1",
        maxiter="defined_by_method",
        full_output=1,
        disp=True,
        callback=None,
        alpha=0,
        trim_mode="auto",
        auto_trim_tol=0.01,
        size_trim_tol=1e-4,
        qc_tol=0.03,
        qc_verbose=False,
        **kwargs,
    ):
        """
        Fit the model using a regularized maximum likelihood.

        The regularization method AND the solver used is determined by the
        argument method.

        Parameters
        ----------
        start_params : array_like, optional
            Initial guess of the solution for the loglikelihood maximization.
            The default is an array of zeros.
        method : 'l1' or 'l1_cvxopt_cp'
            See notes for details.
        maxiter : {int, 'defined_by_method'}
            Maximum number of iterations to perform.
            If 'defined_by_method', then use method defaults (see notes).
        full_output : bool
            Set to True to have all available output in the Results object's
            mle_retvals attribute. The output is dependent on the solver.
            See LikelihoodModelResults notes section for more information.
        disp : bool
            Set to True to print convergence messages.
        fargs : tuple
            Extra arguments passed to the likelihood function, i.e.,
            loglike(x,*args).
        callback : callable callback(xk)
            Called after each iteration, as callback(xk), where xk is the
            current parameter vector.
        retall : bool
            Set to True to return list of solutions at each iteration.
            Available in Results object's mle_retvals attribute.
        alpha : non-negative scalar or numpy array (same size as parameters)
            The weight multiplying the l1 penalty term.
        trim_mode : 'auto, 'size', or 'off'
            If not 'off', trim (set to zero) parameters that would have been
            zero if the solver reached the theoretical minimum.
            If 'auto', trim params using the Theory above.
            If 'size', trim params if they have very small absolute value.
        size_trim_tol : float or 'auto' (default = 'auto')
            Tolerance used when trim_mode == 'size'.
        auto_trim_tol : float
            Tolerance used when trim_mode == 'auto'.
        qc_tol : float
            Print warning and do not allow auto trim when (ii) (above) is
            violated by this much.
        qc_verbose : bool
            If true, print out a full QC report upon failure.
        **kwargs
            Additional keyword arguments used when fitting the model.

        Returns
        -------
        Results
            A results instance.

        Notes
        -----
        Using 'l1_cvxopt_cp' requires the cvxopt module.

        Extra parameters are not penalized if alpha is given as a scalar.
        An example is the shape parameter in NegativeBinomial `nb1` and `nb2`.

        Optional arguments for the solvers (available in Results.mle_settings)::

            'l1'
                acc : float (default 1e-6)
                    Requested accuracy as used by slsqp
            'l1_cvxopt_cp'
                abstol : float
                    absolute accuracy (default: 1e-7).
                reltol : float
                    relative accuracy (default: 1e-6).
                feastol : float
                    tolerance for feasibility conditions (default: 1e-7).
                refinement : int
                    number of iterative refinement steps when solving KKT
                    equations (default: 1).

        Optimization methodology

        With :math:`L` the negative log likelihood, we solve the convex but
        non-smooth problem

        .. math:: \\min_\\beta L(\\beta) + \\sum_k\\alpha_k |\\beta_k|

        via the transformation to the smooth, convex, constrained problem
        in twice as many variables (adding the "added variables" :math:`u_k`)

        .. math:: \\min_{\\beta,u} L(\\beta) + \\sum_k\\alpha_k u_k,

        subject to

        .. math:: -u_k \\leq \\beta_k \\leq u_k.

        With :math:`\\partial_k L` the derivative of :math:`L` in the
        :math:`k^{th}` parameter direction, theory dictates that, at the
        minimum, exactly one of two conditions holds:

        (i) :math:`|\\partial_k L| = \\alpha_k`  and  :math:`\\beta_k \\neq 0`
        (ii) :math:`|\\partial_k L| \\leq \\alpha_k`  and  :math:`\\beta_k = 0`
        """
        _validate_l1_method(method)
        # Set attributes based on method
        cov_params_func = self.cov_params_func_l1

        # Bundle up extra kwargs for the dictionary kwargs.  These are
        # passed through super(...).fit() as kwargs and unpacked at
        # appropriate times
        alpha = np.array(alpha)
        assert alpha.min() >= 0
        try:
            kwargs["alpha"] = alpha
        except TypeError:
            kwargs = dict(alpha=alpha)
        kwargs["alpha_rescaled"] = kwargs["alpha"] / float(self.endog.shape[0])
        kwargs["trim_mode"] = trim_mode
        kwargs["size_trim_tol"] = size_trim_tol
        kwargs["auto_trim_tol"] = auto_trim_tol
        kwargs["qc_tol"] = qc_tol
        kwargs["qc_verbose"] = qc_verbose

        # Define default keyword arguments to be passed to super(...).fit()
        if maxiter == "defined_by_method":
            if method == "l1":
                maxiter = 1000
            elif method == "l1_cvxopt_cp":
                maxiter = 70

        # Parameters to pass to super(...).fit()
        # For the 'extra' parameters, pass all that are available,
        # even if we know (at this point) we will only use one.
        extra_fit_funcs = {"l1": fit_l1_slsqp}
        if have_cvxopt and method == "l1_cvxopt_cp":
            from statsmodels.base.l1_cvxopt import fit_l1_cvxopt_cp

            extra_fit_funcs["l1_cvxopt_cp"] = fit_l1_cvxopt_cp
        elif method.lower() == "l1_cvxopt_cp":
            raise ValueError(
                "Cannot use l1_cvxopt_cp as cvxopt "
                "was not found (install it, or use method='l1' instead)"
            )

        if callback is None:
            callback = self._check_perfect_pred
        else:
            pass  # make a function factory to have multiple call-backs

        mlefit = super().fit(
            start_params=start_params,
            method=method,
            maxiter=maxiter,
            full_output=full_output,
            disp=disp,
            callback=callback,
            extra_fit_funcs=extra_fit_funcs,
            cov_params_func=cov_params_func,
            **kwargs,
        )

        return mlefit  # up to subclasses to wrap results

    def cov_params_func_l1(self, likelihood_model, xopt, retvals):
        """
        Computes cov_params on a reduced parameter space
        corresponding to the nonzero parameters resulting from the
        l1 regularized fit.

        Returns a full cov_params matrix, with entries corresponding
        to zero'd values set to np.nan.
        """
        H = likelihood_model.hessian(xopt)
        trimmed = retvals["trimmed"]
        nz_idx = np.nonzero(~trimmed)[0]
        nnz_params = (~trimmed).sum()
        if nnz_params > 0:
            H_restricted = H[nz_idx[:, None], nz_idx]
            # Covariance estimate for the nonzero params
            H_restricted_inv = np.linalg.inv(-H_restricted)
        else:
            H_restricted_inv = np.zeros(0)

        cov_params = np.nan * np.ones(H.shape)
        cov_params[nz_idx[:, None], nz_idx] = H_restricted_inv

        return cov_params

    def predict(self, params, exog=None, which="mean", linear=None):
        """
        Predict response variable of a model given exogenous variables.
        """
        raise NotImplementedError

    def _derivative_exog(self, params, exog=None, dummy_idx=None, count_idx=None):
        """
        This should implement the derivative of the non-linear function
        """
        raise NotImplementedError

    def _derivative_exog_helper(
        self, margeff, params, exog, dummy_idx, count_idx, transform
    ):
        """
        Helper for _derivative_exog to wrap results appropriately
        """
        from .discrete_margins import _get_count_effects, _get_dummy_effects

        if count_idx is not None:
            margeff = _get_count_effects(
                margeff, exog, count_idx, transform, self, params
            )
        if dummy_idx is not None:
            margeff = _get_dummy_effects(
                margeff, exog, dummy_idx, transform, self, params
            )

        return margeff


class BinaryModel(DiscreteModel):
    _continuous_ok = False

    def __init__(self, endog, exog, offset=None, check_rank=True, **kwargs):
        # unconditional check, requires no extra kwargs added by subclasses
        self._check_kwargs(kwargs)
        super().__init__(endog, exog, offset=offset, check_rank=check_rank, **kwargs)
        if not issubclass(self.__class__, MultinomialModel):
            if not np.all((self.endog >= 0) & (self.endog <= 1)):
                raise ValueError("endog must be in the unit interval.")

        if offset is None:
            del self.offset

            if not self._continuous_ok and np.any(self.endog != np.round(self.endog)):
                raise ValueError("endog must be binary, either 0 or 1")

    def predict(self, params, exog=None, which="mean", linear=None, offset=None):
        """
        Predict response variable of a model given exogenous variables.

        Parameters
        ----------
        params : array_like
            Fitted parameters of the model.
        exog : array_like
            1d or 2d array of exogenous values.  If not supplied, the
            whole exog attribute of the model is used.
        which : {'mean', 'linear', 'var', 'prob'}, optional
            Statistic to predict. Default is 'mean'.

            - 'mean' returns the conditional expectation of endog E(y | x),
              i.e. exp of linear predictor.
            - 'linear' returns the linear predictor of the mean function.
            - 'var' returns the estimated variance of endog implied by the
              model.

            .. versionadded: 0.14

               ``which`` replaces and extends the deprecated ``linear``
               argument.

        linear : bool
            If True, returns the linear predicted values.  If False or None,
            then the statistic specified by ``which`` will be returned.

            .. deprecated: 0.14

               The ``linear` keyword is deprecated and will be removed,
               use ``which`` keyword instead.

        Returns
        -------
        array
            Fitted values at exog.
        """
        if linear is not None:
            msg = 'linear keyword is deprecated, use which="linear"'
            warnings.warn(msg, FutureWarning, stacklevel=2)
            if linear is True:
                which = "linear"

        # Use fit offset if appropriate
        if offset is None and exog is None and hasattr(self, "offset"):
            offset = self.offset
        elif offset is None:
            offset = 0.0

        if exog is None:
            exog = self.exog

        linpred = np.dot(exog, params) + offset

        if which == "mean":
            return self.cdf(linpred)
        elif which == "linear":
            return linpred
        if which == "var":
            mu = self.cdf(linpred)
            var_ = mu * (1 - mu)
            return var_
        else:
            raise ValueError(
                '`which` must be one of "mean", "linear", or "var"'
            )

    @Appender(DiscreteModel.fit_regularized.__doc__)
    def fit_regularized(
        self,
        start_params=None,
        method="l1",
        maxiter="defined_by_method",
        full_output=1,
        disp=1,
        callback=None,
        alpha=0,
        trim_mode="auto",
        auto_trim_tol=0.01,
        size_trim_tol=1e-4,
        qc_tol=0.03,
        **kwargs,
    ):

        _validate_l1_method(method)

        bnryfit = super().fit_regularized(
            start_params=start_params,
            method=method,
            maxiter=maxiter,
            full_output=full_output,
            disp=disp,
            callback=callback,
            alpha=alpha,
            trim_mode=trim_mode,
            auto_trim_tol=auto_trim_tol,
            size_trim_tol=size_trim_tol,
            qc_tol=qc_tol,
            **kwargs,
        )

        discretefit = L1BinaryResults(self, bnryfit)
        return L1BinaryResultsWrapper(discretefit)

    def fit_constrained(self, constraints, start_params=None, **fit_kwds):

        res = fit_constrained_wrap(self, constraints, start_params=None, **fit_kwds)
        return res

    fit_constrained.__doc__ = fit_constrained_wrap.__doc__

    def _derivative_predict(self, params, exog=None, transform="dydx", offset=None):
        """
        For computing marginal effects standard errors.

        This is used only in the case of discrete and count regressors to
        get the variance-covariance of the marginal effects. It returns
        [d F / d params] where F is the predict.

        Transform can be 'dydx' or 'eydx'. Checking is done in margeff
        computations for appropriate transform.
        """
        if exog is None:
            exog = self.exog
        linpred = self.predict(params, exog, offset=offset, which="linear")
        dF = self.pdf(linpred)[:, None] * exog
        if "ey" in transform:
            dF /= self.predict(params, exog, offset=offset)[:, None]
        return dF

    def _derivative_exog(
        self,
        params,
        exog=None,
        transform="dydx",
        dummy_idx=None,
        count_idx=None,
        offset=None,
    ):
        """
        For computing marginal effects returns dF(XB) / dX where F(.) is
        the predicted probabilities

        transform can be 'dydx', 'dyex', 'eydx', or 'eyex'.

        Not all of these make sense in the presence of discrete regressors,
        but checks are done in the results in get_margeff.
        """
        # Note: this form should be appropriate for
        #   group 1 probit, logit, logistic, cloglog, heckprob, xtprobit
        if exog is None:
            exog = self.exog

        linpred = self.predict(params, exog, offset=offset, which="linear")
        margeff = np.dot(self.pdf(linpred)[:, None], params[None, :])

        if "ex" in transform:
            margeff *= exog
        if "ey" in transform:
            margeff /= self.predict(params, exog)[:, None]

        return self._derivative_exog_helper(
            margeff, params, exog, dummy_idx, count_idx, transform
        )

    def _deriv_mean_dparams(self, params):
        """
        Derivative of the expected endog with respect to the parameters.

        Parameters
        ----------
        params : ndarray
            parameter at which score is evaluated

        Returns
        -------
        The value of the derivative of the expected endog with respect
        to the parameter vector.
        """
        link = self.link
        lin_pred = self.predict(params, which="linear")
        idl = link.inverse_deriv(lin_pred)
        dmat = self.exog * idl[:, None]
        return dmat

    def get_distribution(self, params, exog=None, offset=None):
        """Get frozen instance of distribution based on predicted parameters.

        Parameters
        ----------
        params : array_like
            The parameters of the model.
        exog : ndarray, optional
            Explanatory variables for the main count model.
            If ``exog`` is None, then the data from the model will be used.
        offset : ndarray, optional
            Offset is added to the linear predictor of the mean function with
            coefficient equal to 1.
            Default is zero if exog is not None, and the model offset if exog
            is None.
        exposure : ndarray, optional
            Log(exposure) is added to the linear predictor  of the mean
            function with coefficient equal to 1. If exposure is specified,
            then it will be logged by the method. The user does not need to
            log it first.
            Default is one if exog is is not None, and it is the model exposure
            if exog is None.

        Returns
        -------
        Instance of frozen scipy distribution.
        """
        mu = self.predict(params, exog=exog, offset=offset)
        # distr = stats.bernoulli(mu[:, None])
        distr = stats.bernoulli(mu)
        return distr


class MultinomialModel(BinaryModel):

    def _handle_data(self, endog, exog, missing, hasconst, **kwargs):
        if data_tools._is_using_ndarray_type(endog, None):
            endog_dummies, ynames = _numpy_to_dummies(endog)
            yname = "y"
        elif data_tools._is_using_pandas(endog, None):
            endog_dummies, ynames, yname = _pandas_to_dummies(endog)
        else:
            endog = np.asarray(endog)
            endog_dummies, ynames = _numpy_to_dummies(endog)
            yname = "y"

        if not isinstance(ynames, dict):
            ynames = dict(zip(range(endog_dummies.shape[1]), ynames))

        self._ynames_map = ynames
        data = handle_data(endog_dummies, exog, missing, hasconst, **kwargs)
        data.ynames = yname  # overwrite this to single endog name
        data.orig_endog = endog
        self.wendog = data.endog

        # repeating from upstream...
        for key in kwargs:
            if key in ["model_spec", "formula"]:  # leave attached to data
                continue
            try:
                setattr(self, key, data.__dict__.pop(key))
            except KeyError:
                pass
        return data

    def initialize(self):
        """
        Preprocesses the data for MNLogit.
        """
        super().initialize()
        # This is also a "whiten" method in other models (eg regression)
        self.endog = self.endog.argmax(1)  # turn it into an array of col idx
        self.J = self.wendog.shape[1]
        self.K = self.exog.shape[1]
        self.df_model *= self.J - 1  # for each J - 1 equation.
        self.df_resid = self.exog.shape[0] - self.df_model - (self.J - 1)

    def predict(self, params, exog=None, which="mean", linear=None):
        """
        Predict response variable of a model given exogenous variables.

        Parameters
        ----------
        params : array_like
            2d array of fitted parameters of the model. Should be in the
            order returned from the model.
        exog : array_like
            1d or 2d array of exogenous values.  If not supplied, the
            whole exog attribute of the model is used. If a 1d array is given
            it assumed to be 1 row of exogenous variables. If you only have
            one regressor and would like to do prediction, you must provide
            a 2d array with shape[1] == 1.
        which : {'mean', 'linear', 'var', 'prob'}, optional
            Statistic to predict. Default is 'mean'.

            - 'mean' returns the conditional expectation of endog E(y | x),
              i.e. exp of linear predictor.
            - 'linear' returns the linear predictor of the mean function.
            - 'var' returns the estimated variance of endog implied by the
              model.

            .. versionadded: 0.14

               ``which`` replaces and extends the deprecated ``linear``
               argument.

        linear : bool
            If True, returns the linear predicted values.  If False or None,
            then the statistic specified by ``which`` will be returned.

            .. deprecated: 0.14

               The ``linear` keyword is deprecated and will be removed,
               use ``which`` keyword instead.

        Notes
        -----
        Column 0 is the base case, the rest conform to the rows of params
        shifted up one for the base case.
        """
        if linear is not None:
            msg = 'linear keyword is deprecated, use which="linear"'
            warnings.warn(msg, FutureWarning, stacklevel=2)
            if linear is True:
                which = "linear"

        if exog is None:  # do here to accommodate user-given exog
            exog = self.exog
        if exog.ndim == 1:
            exog = exog[None]

        pred = super().predict(params, exog, which=which)
        if which == "linear":
            pred = np.column_stack((np.zeros(len(exog)), pred))
        return pred

    @Appender(DiscreteModel.fit.__doc__)
    def fit(
        self,
        start_params=None,
        method="newton",
        maxiter=35,
        full_output=1,
        disp=1,
        callback=None,
        **kwargs,
    ):
        if start_params is None:
            start_params = np.zeros(self.K * (self.J - 1))
        else:
            start_params = np.asarray(start_params)

        if callback is None:
            # placeholder until check_perfect_pred
            def callback(x, *args):
                return

        # skip calling super to handle results from LikelihoodModel
        mnfit = base.LikelihoodModel.fit(
            self,
            start_params=start_params,
            method=method,
            maxiter=maxiter,
            full_output=full_output,
            disp=disp,
            callback=callback,
            **kwargs,
        )
        mnfit.params = mnfit.params.reshape(self.K, -1, order="F")
        mnfit = MultinomialResults(self, mnfit)
        return MultinomialResultsWrapper(mnfit)

    @Appender(DiscreteModel.fit_regularized.__doc__)
    def fit_regularized(
        self,
        start_params=None,
        method="l1",
        maxiter="defined_by_method",
        full_output=1,
        disp=1,
        callback=None,
        alpha=0,
        trim_mode="auto",
        auto_trim_tol=0.01,
        size_trim_tol=1e-4,
        qc_tol=0.03,
        **kwargs,
    ):
        if start_params is None:
            start_params = np.zeros(self.K * (self.J - 1))
        else:
            start_params = np.asarray(start_params)
        mnfit = DiscreteModel.fit_regularized(
            self,
            start_params=start_params,
            method=method,
            maxiter=maxiter,
            full_output=full_output,
            disp=disp,
            callback=callback,
            alpha=alpha,
            trim_mode=trim_mode,
            auto_trim_tol=auto_trim_tol,
            size_trim_tol=size_trim_tol,
            qc_tol=qc_tol,
            **kwargs,
        )
        mnfit.params = mnfit.params.reshape(self.K, -1, order="F")
        mnfit = L1MultinomialResults(self, mnfit)
        return L1MultinomialResultsWrapper(mnfit)

    def _derivative_predict(self, params, exog=None, transform="dydx"):
        """
        For computing marginal effects standard errors.

        This is used only in the case of discrete and count regressors to
        get the variance-covariance of the marginal effects. It returns
        [d F / d params] where F is the predicted probabilities for each
        choice. dFdparams is of shape nobs x (J*K) x (J-1)*K.
        The zero derivatives for the base category are not included.

        Transform can be 'dydx' or 'eydx'. Checking is done in margeff
        computations for appropriate transform.
        """
        if exog is None:
            exog = self.exog
        if params.ndim == 1:  # will get flatted from approx_fprime
            params = params.reshape(self.K, self.J - 1, order="F")

        eXB = np.exp(np.dot(exog, params))
        sum_eXB = (1 + eXB.sum(1))[:, None]
        J = int(self.J)
        K = int(self.K)
        repeat_eXB = np.repeat(eXB, J, axis=1)
        X = np.tile(exog, J - 1)
        # this is the derivative wrt the base level
        F0 = -repeat_eXB * X / sum_eXB**2
        # this is the derivative wrt the other levels when
        # dF_j / dParams_j (ie., own equation)
        # NOTE: this computes too much, any easy way to cut down?
        F1 = eXB.T[:, :, None] * X * (sum_eXB - repeat_eXB) / (sum_eXB**2)
        F1 = F1.transpose((1, 0, 2))  # put the nobs index first

        # other equation index
        other_idx = ~np.kron(np.eye(J - 1), np.ones(K)).astype(bool)
        F1[:, other_idx] = (
            -eXB.T[:, :, None] * X * repeat_eXB / (sum_eXB**2)
        ).transpose((1, 0, 2))[:, other_idx]
        dFdX = np.concatenate((F0[:, None, :], F1), axis=1)

        if "ey" in transform:
            dFdX /= self.predict(params, exog)[:, :, None]
        return dFdX

    def _derivative_exog(
        self, params, exog=None, transform="dydx", dummy_idx=None, count_idx=None
    ):
        """
        For computing marginal effects returns dF(XB) / dX where F(.) is
