"""
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__ = ["Poisson", "Logit", "Probit", "MNLogit", "NegativeBinomial",
           "GeneralizedPoisson", "NegativeBinomialP", "CountModel"]

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.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.base._constraints import fit_constrained_wrap
import statsmodels.base._parameter_inference as pinfer
from statsmodels.base import _prediction_inference as pred
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, category=PerfectSeparationWarning)

    @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:
            delattr(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)
            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.

        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('Only `which` is "mean", "linear" or "var" are'
                             ' available.')

    @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 ['design_info', '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)
            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
            callback = lambda x, *args : None
        # 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
        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.

        For Multinomial models the marginal effects are

        P[j] * (params[j] - sum_k P[k]*params[k])

        It is returned unshaped, so that each row contains each of the J
        equations. This makes it easier to take derivatives of this for
        standard errors. If you want average marginal effects you can do
        margeff.reshape(nobs, K, J, order='F).mean(0) and the marginal effects
        for choice J are in column J
        """
        J = int(self.J)  # number of alternative choices
        K = int(self.K)  # number of variables
        # Note: this form should be appropriate for
        #   group 1 probit, logit, logistic, cloglog, heckprob, xtprobit
        if exog is None:
            exog = self.exog
        if params.ndim == 1:  # will get flatted from approx_fprime
            params = params.reshape(K, J-1, order='F')

        zeroparams = np.c_[np.zeros(K), params]  # add base in

        cdf = self.cdf(np.dot(exog, params))

        # TODO: meaningful interpretation for `iterm`?
        iterm = np.array([cdf[:, [i]] * zeroparams[:, i]
                          for i in range(int(J))]).sum(0)

        margeff = np.array([cdf[:, [j]] * (zeroparams[:, j] - iterm)
                            for j in range(J)])

        # swap the axes to make sure margeff are in order nobs, K, J
        margeff = np.transpose(margeff, (1, 2, 0))

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

        margeff = self._derivative_exog_helper(margeff, params, exog,
                                               dummy_idx, count_idx, transform)
        return margeff.reshape(len(exog), -1, order='F')

    def get_distribution(self, params, exog=None, offset=None):
        """get frozen instance of distribution
        """
        raise NotImplementedError


class CountModel(DiscreteModel):
    def __init__(self, endog, exog, offset=None, exposure=None, missing='none',
                 check_rank=True, **kwargs):
        self._check_kwargs(kwargs)
        super().__init__(endog, exog, check_rank, missing=missing,
                         offset=offset, exposure=exposure, **kwargs)
        if exposure is not None:
            self.exposure = np.asarray(self.exposure)
            self.exposure = np.log(self.exposure)
        if offset is not None:
            self.offset = np.asarray(self.offset)
        self._check_inputs(self.offset, self.exposure, self.endog)
        if offset is None:
            delattr(self, 'offset')
        if exposure is None:
            delattr(self, 'exposure')

        # promote dtype to float64 if needed
        dt = np.promote_types(self.endog.dtype, np.float64)
        self.endog = np.asarray(self.endog, dt)
        dt = np.promote_types(self.exog.dtype, np.float64)
        self.exog = np.asarray(self.exog, dt)


    def _check_inputs(self, offset, exposure, endog):
        if offset is not None and offset.shape[0] != endog.shape[0]:
            raise ValueError("offset is not the same length as endog")

        if exposure is not None and exposure.shape[0] != endog.shape[0]:
            raise ValueError("exposure is not the same length as endog")

    def _get_init_kwds(self):
        # this is a temporary fixup because exposure has been transformed
        # see #1609
        kwds = super()._get_init_kwds()
        if 'exposure' in kwds and kwds['exposure'] is not None:
            kwds['exposure'] = np.exp(kwds['exposure'])
        return kwds

    def _get_predict_arrays(self, exog=None, offset=None, exposure=None):

        # convert extras if not None
        if exposure is not None:
            exposure = np.log(np.asarray(exposure))
        if offset is not None:
            offset = np.asarray(offset)

        # get defaults
        if exog is None:
            # prediction is in-sample
            exog = self.exog
            if exposure is None:
                exposure = getattr(self, 'exposure', 0)
            if offset is None:
                offset = getattr(self, 'offset', 0)
        else:
            # user specified
            exog = np.asarray(exog)
            if exposure is None:
                exposure = 0
            if offset is None:
                offset = 0

        return exog, offset, exposure

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

        Parameters
        ----------
        params : array_like
            Model parameters
        exog : array_like, optional
            Design / exogenous data. Is exog is None, model exog is used.
        exposure : array_like, optional
            Log(exposure) is added to the linear prediction with
            coefficient equal to 1. If exposure is not provided and exog
            is None, uses the model's exposure if present.  If not, uses
            0 as the default value.
        offset : array_like, optional
            Offset is added to the linear prediction with coefficient
            equal to 1. If offset is not provided and exog
            is None, uses the model's offset if present.  If not, uses
            0 as the default value.
        which : 'mean', 'linear', 'var', 'prob' (optional)
            Statitistic 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' variance of endog implied by the likelihood model
            - 'prob' predicted probabilities for counts.

            .. 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
        -----
        If exposure is specified, then it will be logged by the method.
        The user does not need to log it first.
        """
        if linear is not None:
            msg = 'linear keyword is deprecated, use which="linear"'
            warnings.warn(msg, FutureWarning)
            if linear is True:
                which = "linear"

        # the following is copied from GLM predict (without family/link check)
        # 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.

        # Use fit exposure if appropriate
        if exposure is None and exog is None and hasattr(self, 'exposure'):
            # Already logged
            exposure = self.exposure
        elif exposure is None:
            exposure = 0.
        else:
            exposure = np.log(exposure)

        if exog is None:
            exog = self.exog

        fitted = np.dot(exog, params[:exog.shape[1]])
        linpred = fitted + exposure + offset
        if which == "mean":
            return np.exp(linpred)
        elif which.startswith("lin"):
            return linpred
        else:
            raise ValueError('keyword which has to be "mean" and "linear"')

    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 predict.

        Transform can be 'dydx' or 'eydx'. Checking is done in margeff
        computations for appropriate transform.
        """
        if exog is None:
            exog = self.exog
        #NOTE: this handles offset and exposure
        dF = self.predict(params, exog)[:,None] * exog
        if 'ey' in transform:
            dF /= self.predict(params, exog)[:,None]
        return dF

    def _derivative_exog(self, params, exog=None, transform="dydx",
                         dummy_idx=None, count_idx=None):
        """
        For computing marginal effects. These are the marginal effects
        d F(XB) / dX
        For the Poisson model F(XB) is the predicted counts rather than
        the 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.
        """
        # group 3 poisson, nbreg, zip, zinb
        if exog is None:
            exog = self.exog
        k_extra = getattr(self, 'k_extra', 0)
        params_exog = params if k_extra == 0 else params[:-k_extra]
        margeff = self.predict(params, exog)[:,None] * params_exog[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.
        """
        from statsmodels.genmod.families import links
        link = links.Log()
        lin_pred = self.predict(params, which="linear")
        idl = link.inverse_deriv(lin_pred)
        dmat = self.exog * idl[:, None]
        if self.k_extra > 0:
            dmat_extra = np.zeros((dmat.shape[0], self.k_extra))
            dmat = np.column_stack((dmat, dmat_extra))
        return dmat


    @Appender(DiscreteModel.fit.__doc__)
    def fit(self, start_params=None, method='newton', maxiter=35,
            full_output=1, disp=1, callback=None, **kwargs):
        cntfit = super().fit(start_params=start_params,
                             method=method,
                             maxiter=maxiter,
                             full_output=full_output,
                             disp=disp,
                             callback=callback,
                             **kwargs)
        discretefit = CountResults(self, cntfit)
        return CountResultsWrapper(discretefit)

    @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)

        cntfit = 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 = L1CountResults(self, cntfit)
        return L1CountResultsWrapper(discretefit)


# Public Model Classes


class Poisson(CountModel):
    __doc__ = """
    Poisson Model

    %(params)s
    %(extra_params)s

    Attributes
    ----------
    endog : ndarray
        A reference to the endogenous response variable
    exog : ndarray
        A reference to the exogenous design.
    """ % {'params': base._model_params_doc,
           'extra_params':
           """offset : array_like
        Offset is added to the linear prediction with coefficient equal to 1.
    exposure : array_like
        Log(exposure) is added to the linear prediction with coefficient
        equal to 1.
        """ + base._missing_param_doc + _check_rank_doc}

    @cache_readonly
    def family(self):
        from statsmodels.genmod import families
        return families.Poisson()

    def cdf(self, X):
        """
        Poisson model cumulative distribution function

        Parameters
        ----------
        X : array_like
            `X` is the linear predictor of the model.  See notes.

        Returns
        -------
        The value of the Poisson CDF at each point.

        Notes
        -----
        The CDF is defined as

        .. math:: \\exp\\left(-\\lambda\\right)\\sum_{i=0}^{y}\\frac{\\lambda^{i}}{i!}

        where :math:`\\lambda` assumes the loglinear model. I.e.,

        .. math:: \\ln\\lambda_{i}=X\\beta

        The parameter `X` is :math:`X\\beta` in the above formula.
        """
        y = self.endog
        return stats.poisson.cdf(y, np.exp(X))

    def pdf(self, X):
        """
        Poisson model probability mass function

        Parameters
        ----------
        X : array_like
            `X` is the linear predictor of the model.  See notes.

        Returns
        -------
        pdf : ndarray
            The value of the Poisson probability mass function, PMF, for each
            point of X.

        Notes
        -----
        The PMF is defined as

        .. math:: \\frac{e^{-\\lambda_{i}}\\lambda_{i}^{y_{i}}}{y_{i}!}

        where :math:`\\lambda` assumes the loglinear model. I.e.,

        .. math:: \\ln\\lambda_{i}=x_{i}\\beta

        The parameter `X` is :math:`x_{i}\\beta` in the above formula.
        """
        y = self.endog
        return np.exp(stats.poisson.logpmf(y, np.exp(X)))

    def loglike(self, params):
        """
        Loglikelihood of Poisson model

        Parameters
        ----------
        params : array_like
            The parameters of the model.

        Returns
        -------
        loglike : float
            The log-likelihood function of the model evaluated at `params`.
            See notes.

        Notes
        -----
        .. math:: \\ln L=\\sum_{i=1}^{n}\\left[-\\lambda_{i}+y_{i}x_{i}^{\\prime}\\beta-\\ln y_{i}!\\right]
        """
        offset = getattr(self, "offset", 0)
        exposure = getattr(self, "exposure", 0)
        XB = np.dot(self.exog, params) + offset + exposure
        endog = self.endog
        return np.sum(
            -np.exp(np.clip(XB, None, EXP_UPPER_LIMIT))
            + endog * XB
            - gammaln(endog + 1)
        )

    def loglikeobs(self, params):
        """
        Loglikelihood for observations of Poisson model

        Parameters
        ----------
        params : array_like
            The parameters of the model.

        Returns
        -------
        loglike : array_like
            The log likelihood for each observation of the model evaluated
            at `params`. See Notes

        Notes
        -----
        .. math:: \\ln L_{i}=\\left[-\\lambda_{i}+y_{i}x_{i}^{\\prime}\\beta-\\ln y_{i}!\\right]

        for observations :math:`i=1,...,n`
        """
        offset = getattr(self, "offset", 0)
        exposure = getattr(self, "exposure", 0)
        XB = np.dot(self.exog, params) + offset + exposure
        endog = self.endog
        #np.sum(stats.poisson.logpmf(endog, np.exp(XB)))
        return -np.exp(XB) +  endog*XB - gammaln(endog+1)

    @Appender(_get_start_params_null_docs)
    def _get_start_params_null(self):
        offset = getattr(self, "offset", 0)
        exposure = getattr(self, "exposure", 0)
        const = (self.endog / np.exp(offset + exposure)).mean()
        params = [np.log(const)]
        return params

    @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 and self.data.const_idx is not None:
            # k_params or k_exog not available?
            start_params = 0.001 * np.ones(self.exog.shape[1])
            start_params[self.data.const_idx] = self._get_start_params_null()[0]

        kwds = {}
        if kwargs.get('cov_type') is not None:
            kwds['cov_type'] = kwargs.get('cov_type')
            kwds['cov_kwds'] = kwargs.get('cov_kwds', {})

        cntfit = super(CountModel, self).fit(start_params=start_params,
                                             method=method,
                                             maxiter=maxiter,
                                             full_output=full_output,
                                             disp=disp,
                                             callback=callback,
                                             **kwargs)

        discretefit = PoissonResults(self, cntfit, **kwds)
        return PoissonResultsWrapper(discretefit)

    @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)

        cntfit = super(CountModel, self).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 = L1PoissonResults(self, cntfit)
        return L1PoissonResultsWrapper(discretefit)

    def fit_constrained(self, constraints, start_params=None, **fit_kwds):
        """fit the model subject to linear equality constraints

        The constraints are of the form   `R params = q`
        where R is the constraint_matrix and q is the vector of
        constraint_values.

        The estimation creates a new model with transformed design matrix,
        exog, and converts the results back to the original parameterization.

        Parameters
        ----------
        constraints : formula expression or tuple
            If it is a tuple, then the constraint needs to be given by two
            arrays (constraint_matrix, constraint_value), i.e. (R, q).
            Otherwise, the constraints can be given as strings or list of
            strings.
            see t_test for details
        start_params : None or array_like
            starting values for the optimization. `start_params` needs to be
            given in the original parameter space and are internally
            transformed.
        **fit_kwds : keyword arguments
            fit_kwds are used in the optimization of the transformed model.

        Returns
        -------
        results : Results instance
        """

        #constraints = (R, q)
        # TODO: temporary trailing underscore to not overwrite the monkey
        #       patched version
        # TODO: decide whether to move the imports
        from patsy import DesignInfo
        from statsmodels.base._constraints import (fit_constrained,
                                                   LinearConstraints)

        # same pattern as in base.LikelihoodModel.t_test
        lc = DesignInfo(self.exog_names).linear_constraint(constraints)
        R, q = lc.coefs, lc.constants

        # TODO: add start_params option, need access to tranformation
        #       fit_constrained needs to do the transformation
        params, cov, res_constr = fit_constrained(self, R, q,
                                                  start_params=start_params,
                                                  fit_kwds=fit_kwds)
        #create dummy results Instance, TODO: wire up properly
        res = self.fit(maxiter=0, method='nm', disp=0,
                       warn_convergence=False) # we get a wrapper back
        res.mle_retvals['fcall'] = res_constr.mle_retvals.get('fcall', np.nan)
        res.mle_retvals['iterations'] = res_constr.mle_retvals.get(
                                                        'iterations', np.nan)
        res.mle_retvals['converged'] = res_constr.mle_retvals['converged']
        res._results.params = params
        res._results.cov_params_default = cov
        cov_type = fit_kwds.get('cov_type', 'nonrobust')
        if cov_type != 'nonrobust':
            res._results.normalized_cov_params = cov # assume scale=1
        else:
            res._results.normalized_cov_params = None
        k_constr = len(q)
        res._results.df_resid += k_constr
        res._results.df_model -= k_constr
        res._results.constraints = LinearConstraints.from_patsy(lc)
        res._results.k_constr = k_constr
        res._results.results_constrained = res_constr
        return res

    def score(self, params):
        """
        Poisson model score (gradient) vector of the log-likelihood

        Parameters
        ----------
        params : array_like
            The parameters of the model

        Returns
        -------
        score : ndarray, 1-D
            The score vector of the model, i.e. the first derivative of the
            loglikelihood function, evaluated at `params`

        Notes
        -----
        .. math:: \\frac{\\partial\\ln L}{\\partial\\beta}=\\sum_{i=1}^{n}\\left(y_{i}-\\lambda_{i}\\right)x_{i}

        where the loglinear model is assumed

        .. math:: \\ln\\lambda_{i}=x_{i}\\beta
        """
        offset = getattr(self, "offset", 0)
        exposure = getattr(self, "exposure", 0)
        X = self.exog
        L = np.exp(np.dot(X,params) + offset + exposure)
        return np.dot(self.endog - L, X)

    def score_obs(self, params):
        """
        Poisson model Jacobian of the log-likelihood for each observation

        Parameters
        ----------
        params : array_like
            The parameters of the model

        Returns
        -------
        score : array_like
            The score vector (nobs, k_vars) of the model evaluated at `params`

        Notes
        -----
        .. math:: \\frac{\\partial\\ln L_{i}}{\\partial\\beta}=\\left(y_{i}-\\lambda_{i}\\right)x_{i}

        for observations :math:`i=1,...,n`

        where the loglinear model is assumed

        .. math:: \\ln\\lambda_{i}=x_{i}\\beta
        """
        offset = getattr(self, "offset", 0)
        exposure = getattr(self, "exposure", 0)
        X = self.exog
        L = np.exp(np.dot(X,params) + offset + exposure)
        return (self.endog - L)[:,None] * X

    def score_factor(self, params):
        """
        Poisson model score_factor for each observation

        Parameters
        ----------
        params : array_like
            The parameters of the model

        Returns
        -------
        score : array_like
            The score factor (nobs, ) of the model evaluated at `params`

        Notes
        -----
        .. math:: \\frac{\\partial\\ln L_{i}}{\\partial\\beta}=\\left(y_{i}-\\lambda_{i}\\right)

        for observations :math:`i=1,...,n`

        where the loglinear model is assumed

        .. math:: \\ln\\lambda_{i}=x_{i}\\beta
        """
        offset = getattr(self, "offset", 0)
        exposure = getattr(self, "exposure", 0)
        X = self.exog
        L = np.exp(np.dot(X,params) + offset + exposure)
        return (self.endog - L)

    def hessian(self, params):
        """
        Poisson model Hessian matrix of the loglikelihood

        Parameters
        ----------
        params : array_like
            The parameters of the model

        Returns
        -------
        hess : ndarray, (k_vars, k_vars)
            The Hessian, second derivative of loglikelihood function,
            evaluated at `params`

        Notes
        -----
        .. math:: \\frac{\\partial^{2}\\ln L}{\\partial\\beta\\partial\\beta^{\\prime}}=-\\sum_{i=1}^{n}\\lambda_{i}x_{i}x_{i}^{\\prime}

        where the loglinear model is assumed

        .. math:: \\ln\\lambda_{i}=x_{i}\\beta
        """
        offset = getattr(self, "offset", 0)
        exposure = getattr(self, "exposure", 0)
        X = self.exog
        L = np.exp(np.dot(X,params) + exposure + offset)
        return -np.dot(L*X.T, X)

    def hessian_factor(self, params):
        """
        Poisson model Hessian factor

        Parameters
        ----------
        params : array_like
            The parameters of the model

        Returns
        -------
        hess : ndarray, (nobs,)
            The Hessian factor, second derivative of loglikelihood function
            with respect to the linear predictor evaluated at `params`

        Notes
        -----
        .. math:: \\frac{\\partial^{2}\\ln L}{\\partial\\beta\\partial\\beta^{\\prime}}=-\\sum_{i=1}^{n}\\lambda_{i}

        where the loglinear model is assumed

        .. math:: \\ln\\lambda_{i}=x_{i}\\beta
        """
        offset = getattr(self, "offset", 0)
        exposure = getattr(self, "exposure", 0)
        X = self.exog
        L = np.exp(np.dot(X,params) + exposure + offset)
        return -L

    def _deriv_score_obs_dendog(self, params, scale=None):
        """derivative of score_obs w.r.t. endog

        Parameters
        ----------
        params : ndarray
            parameter at which score is evaluated
        scale : None or float
            If scale is None, then the default scale will be calculated.
            Default scale is defined by `self.scaletype` and set in fit.
            If scale is not None, then it is used as a fixed scale.

        Returns
        -------
        derivative : ndarray_2d
            The derivative of the score_obs with respect to endog. This
            can is given by `score_factor0[:, None] * exog` where
            `score_factor0` is the score_factor without the residual.
        """
        return self.exog

    def predict(self, params, exog=None, exposure=None, offset=None,
                which='mean', linear=None, y_values=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, optional
            1d or 2d array of exogenous values.  If not supplied, then the
            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.
        offset : array_like, optional
            Offset is added to the linear predictor with coefficient equal
            to 1.
            Default is zero if exog is not None, and the model offset if exog
            is None.
        exposure : array_like, optional
            Log(exposure) is added to the linear prediction with coefficient
            equal to 1.
            Default is one if exog is is not None, and is the model exposure
            if exog is None.
        which : 'mean', 'linear', 'var', 'prob' (optional)
            Statitistic 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.
            - 'prob' return probabilities for counts from 0 to max(endog) or
              for y_values if those are provided.

            .. versionadded: 0.14

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

        linear : bool
            The ``linear` keyword is deprecated and will be removed,
            use ``which`` keyword instead.
            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.

        y_values : array_like
            Values of the random variable endog at which pmf is evaluated.
            Only used if ``which="prob"``
        """
        # Note docstring is reused by other count models

        if linear is not None:
            msg = 'linear keyword is deprecated, use which="linear"'
            warnings.warn(msg, FutureWarning)
            if linear is True:
                which = "linear"

        if which.startswith("lin"):
            which = "linear"
        if which in ["mean", "linear"]:
            return super().predict(params, exog=exog, exposure=exposure,
                                   offset=offset,
                                   which=which, linear=linear)
        # TODO: add full set of which
        elif which == "var":
            mu = self.predict(params, exog=exog,
                              exposure=exposure, offset=offset,
                              )
            return mu
        elif which == "prob":
            if y_values is not None:
                y_values = np.atleast_2d(y_values)
            else:
                y_values = np.atleast_2d(
                    np.arange(0, np.max(self.endog) + 1))
            mu = self.predict(params, exog=exog,
                              exposure=exposure, offset=offset,
                              )[:, None]
            # uses broadcasting
            return stats.poisson._pmf(y_values, mu)
        else:
            raise ValueError('Value of the `which` option is not recognized')

    def _prob_nonzero(self, mu, params=None):
        """Probability that count is not zero

        internal use in Censored model, will be refactored or removed
        """
        prob_nz = - np.expm1(-mu)
        return prob_nz

    def _var(self, mu, params=None):
        """variance implied by the distribution

        internal use, will be refactored or removed
        """
        return mu

    def get_distribution(self, params, exog=None, exposure=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 subclass.
        """
        mu = self.predict(params, exog=exog, exposure=exposure, offset=offset)
        distr = stats.poisson(mu)
        return distr


class GeneralizedPoisson(CountModel):
    __doc__ = """
    Generalized Poisson Model

    %(params)s
    %(extra_params)s

    Attributes
    ----------
    endog : ndarray
        A reference to the endogenous response variable
    exog : ndarray
        A reference to the exogenous design.
    """ % {'params': base._model_params_doc,
           'extra_params':
               """
    p : scalar
        P denotes parameterizations for GP regression. p=1 for GP-1 and
        p=2 for GP-2. Default is p=1.
    offset : array_like
        Offset is added to the linear prediction with coefficient equal to 1.
    exposure : array_like
        Log(exposure) is added to the linear prediction with coefficient
        equal to 1.""" + base._missing_param_doc + _check_rank_doc}

    def __init__(self, endog, exog, p=1, offset=None,
                 exposure=None, missing='none', check_rank=True, **kwargs):
        super().__init__(endog,
                         exog,
                         offset=offset,
                         exposure=exposure,
                         missing=missing,
                         check_rank=check_rank,
                         **kwargs)
        self.parameterization = p - 1
        self.exog_names.append('alpha')
        self.k_extra = 1
        self._transparams = False

    def _get_init_kwds(self):
        kwds = super()._get_init_kwds()
        kwds['p'] = self.parameterization + 1
        return kwds

    def _get_exogs(self):
        return (self.exog, None)

    def loglike(self, params):
        """
        Loglikelihood of Generalized Poisson model

        Parameters
        ----------
        params : array_like
            The parameters of the model.

        Returns
        -------
        loglike : float
            The log-likelihood function of the model evaluated at `params`.
            See notes.

        Notes
        -----
        .. math:: \\ln L=\\sum_{i=1}^{n}\\left[\\mu_{i}+(y_{i}-1)*ln(\\mu_{i}+
            \\alpha*\\mu_{i}^{p-1}*y_{i})-y_{i}*ln(1+\\alpha*\\mu_{i}^{p-1})-
            ln(y_{i}!)-\\frac{\\mu_{i}+\\alpha*\\mu_{i}^{p-1}*y_{i}}{1+\\alpha*
            \\mu_{i}^{p-1}}\\right]
        """
        return np.sum(self.loglikeobs(params))

    def loglikeobs(self, params):
        """
        Loglikelihood for observations of Generalized Poisson model

        Parameters
        ----------
        params : array_like
            The parameters of the model.

        Returns
        -------
        loglike : ndarray
            The log likelihood for each observation of the model evaluated
            at `params`. See Notes

        Notes
        -----
        .. math:: \\ln L=\\sum_{i=1}^{n}\\left[\\mu_{i}+(y_{i}-1)*ln(\\mu_{i}+
            \\alpha*\\mu_{i}^{p-1}*y_{i})-y_{i}*ln(1+\\alpha*\\mu_{i}^{p-1})-
            ln(y_{i}!)-\\frac{\\mu_{i}+\\alpha*\\mu_{i}^{p-1}*y_{i}}{1+\\alpha*
            \\mu_{i}^{p-1}}\\right]

        for observations :math:`i=1,...,n`
        """
        if self._transparams:
            alpha = np.exp(params[-1])
        else:
            alpha = params[-1]
        params = params[:-1]
        p = self.parameterization
        endog = self.endog
        mu = self.predict(params)
        mu_p = np.power(mu, p)
        a1 = 1 + alpha * mu_p
        a2 = mu + (a1 - 1) * endog
        a1 = np.maximum(1e-20, a1)
        a2 = np.maximum(1e-20, a2)
        return (np.log(mu) + (endog - 1) * np.log(a2) - endog *
                np.log(a1) - gammaln(endog + 1) - a2 / a1)

    @Appender(_get_start_params_null_docs)
    def _get_start_params_null(self):
        offset = getattr(self, "offset", 0)
        exposure = getattr(self, "exposure", 0)

        const = (self.endog / np.exp(offset + exposure)).mean()
        params = [np.log(const)]
        mu = const * np.exp(offset + exposure)
        resid = self.endog - mu
        a = self._estimate_dispersion(mu, resid, df_resid=resid.shape[0] - 1)
        params.append(a)

        return np.array(params)

    def _estimate_dispersion(self, mu, resid, df_resid=None):
        q = self.parameterization
        if df_resid is None:
            df_resid = resid.shape[0]
        a = ((np.abs(resid) / np.sqrt(mu) - 1) * mu**(-q)).sum() / df_resid
        return a


    @Appender(
        """
        use_transparams : bool
            This parameter enable internal transformation to impose
            non-negativity. True to enable. Default is False.
            use_transparams=True imposes the no underdispersion (alpha > 0)
            constraint. In case use_transparams=True and method="newton" or
            "ncg" transformation is ignored.
        """)
    @Appender(DiscreteModel.fit.__doc__)
    def fit(self, start_params=None, method='bfgs', maxiter=35,
            full_output=1, disp=1, callback=None, use_transparams=False,
            cov_type='nonrobust', cov_kwds=None, use_t=None, optim_kwds_prelim=None,
            **kwargs):
        if use_transparams and method not in ['newton', 'ncg']:
            self._transparams = True
        else:
            if use_transparams:
                warnings.warn('Parameter "use_transparams" is ignored',
                              RuntimeWarning)
            self._transparams = False

        if start_params is None:
            offset = getattr(self, "offset", 0) + getattr(self, "exposure", 0)
            if np.size(offset) == 1 and offset == 0:
                offset = None
            kwds_prelim = {'disp': 0, 'skip_hessian': True,
                           'warn_convergence': False}
            if optim_kwds_prelim is not None:
                kwds_prelim.update(optim_kwds_prelim)
            mod_poi = Poisson(self.endog, self.exog, offset=offset)
            with warnings.catch_warnings():
                warnings.simplefilter("always")
                res_poi = mod_poi.fit(**kwds_prelim)
            start_params = res_poi.params
            a = self._estimate_dispersion(res_poi.predict(), res_poi.resid,
                                          df_resid=res_poi.df_resid)
            start_params = np.append(start_params, max(-0.1, a))

        if callback is None:
            # work around perfect separation callback #3895
            callback = lambda *x: x

        mlefit = super().fit(start_params=start_params,
                             maxiter=maxiter,
                             method=method,
                             disp=disp,
                             full_output=full_output,
                             callback=callback,
                             **kwargs)
        if optim_kwds_prelim is not None:
            mlefit.mle_settings["optim_kwds_prelim"] = optim_kwds_prelim
        if use_transparams and method not in ["newton", "ncg"]:
            self._transparams = False
            mlefit._results.params[-1] = np.exp(mlefit._results.params[-1])

        gpfit = GeneralizedPoissonResults(self, mlefit._results)
        result = GeneralizedPoissonResultsWrapper(gpfit)

        if cov_kwds is None:
            cov_kwds = {}

        result._get_robustcov_results(cov_type=cov_type,
                                      use_self=True, use_t=use_t, **cov_kwds)
        return result

    @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)

        if np.size(alpha) == 1 and alpha != 0:
            k_params = self.exog.shape[1] + self.k_extra
            alpha = alpha * np.ones(k_params)
            alpha[-1] = 0

        alpha_p = alpha[:-1] if (self.k_extra and np.size(alpha) > 1) else alpha
        self._transparams = False
        if start_params is None:
            offset = getattr(self, "offset", 0) + getattr(self, "exposure", 0)
            if np.size(offset) == 1 and offset == 0:
                offset = None
            mod_poi = Poisson(self.endog, self.exog, offset=offset)
            with warnings.catch_warnings():
                warnings.simplefilter("always")
                start_params = mod_poi.fit_regularized(
                    start_params=start_params, method=method, maxiter=maxiter,
                    full_output=full_output, disp=0, callback=callback,
                    alpha=alpha_p, trim_mode=trim_mode,
                    auto_trim_tol=auto_trim_tol, size_trim_tol=size_trim_tol,
                    qc_tol=qc_tol, **kwargs).params
            start_params = np.append(start_params, 0.1)

        cntfit = super(CountModel, self).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 = L1GeneralizedPoissonResults(self, cntfit)
        return L1GeneralizedPoissonResultsWrapper(discretefit)

    def score_obs(self, params):
        if self._transparams:
            alpha = np.exp(params[-1])
        else:
            alpha = params[-1]

        params = params[:-1]
        p = self.parameterization
        exog = self.exog
        y = self.endog[:,None]
        mu = self.predict(params)[:,None]
        mu_p = np.power(mu, p)
        a1 = 1 + alpha * mu_p
        a2 = mu + alpha * mu_p * y
        a3 = alpha * p * mu ** (p - 1)
        a4 = a3 * y
        dmudb = mu * exog

        dalpha = (mu_p * (y * ((y - 1) / a2 - 2 / a1) + a2 / a1**2))
        dparams = dmudb * (-a4 / a1 +
                           a3 * a2 / (a1 ** 2) +
                           (1 + a4) * ((y - 1) / a2 - 1 / a1) +
                           1 / mu)

        return np.concatenate((dparams, np.atleast_2d(dalpha)),
                              axis=1)

    def score(self, params):
        score = np.sum(self.score_obs(params), axis=0)
        if self._transparams:
            score[-1] == score[-1] ** 2
            return score
        else:
            return score

    def score_factor(self, params, endog=None):
        params = np.asarray(params)
        if self._transparams:
            alpha = np.exp(params[-1])
        else:
            alpha = params[-1]

        params = params[:-1]
        p = self.parameterization
        y = self.endog if endog is None else endog

        mu = self.predict(params)
        mu_p = np.power(mu, p)
        a1 = 1 + alpha * mu_p
        a2 = mu + alpha * mu_p * y
        a3 = alpha * p * mu ** (p - 1)
        a4 = a3 * y
        dmudb = mu

        dalpha = (mu_p * (y * ((y - 1) / a2 - 2 / a1) + a2 / a1**2))
        dparams = dmudb * (-a4 / a1 +
                           a3 * a2 / (a1 ** 2) +
                           (1 + a4) * ((y - 1) / a2 - 1 / a1) +
                           1 / mu)

        return dparams, dalpha

    def _score_p(self, params):
        """
        Generalized Poisson model derivative of the log-likelihood by p-parameter

        Parameters
        ----------
        params : array_like
            The parameters of the model

        Returns
        -------
        dldp : float
            dldp is first derivative of the loglikelihood function,
        evaluated at `p-parameter`.
        """
        if self._transparams:
            alpha = np.exp(params[-1])
        else:
            alpha = params[-1]
        params = params[:-1]
        p = self.parameterization
        y = self.endog[:,None]
        mu = self.predict(params)[:,None]
        mu_p = np.power(mu, p)
        a1 = 1 + alpha * mu_p
        a2 = mu + alpha * mu_p * y

        dp = np.sum(np.log(mu) * ((a2 - mu) * ((y - 1) / a2 - 2 / a1) +
                                   (a1 - 1) * a2 / a1 ** 2))
        return dp

    def hessian(self, params):
        """
        Generalized Poisson model Hessian matrix of the loglikelihood

        Parameters
        ----------
        params : array_like
            The parameters of the model

        Returns
        -------
        hess : ndarray, (k_vars, k_vars)
            The Hessian, second derivative of loglikelihood function,
            evaluated at `params`
        """
        if self._transparams:
            alpha = np.exp(params[-1])
        else:
            alpha = params[-1]

        params = params[:-1]
        p = self.parameterization
        exog = self.exog
        y = self.endog[:,None]
        mu = self.predict(params)[:,None]
        mu_p = np.power(mu, p)
        a1 = 1 + alpha * mu_p
        a2 = mu + alpha * mu_p * y
        a3 = alpha * p * mu ** (p - 1)
        a4 = a3 * y
        a5 = p * mu ** (p - 1)
        dmudb = mu * exog

        # for dl/dparams dparams
        dim = exog.shape[1]
        hess_arr = np.empty((dim+1,dim+1))

        for i in range(dim):
            for j in range(i + 1):
                hess_val = np.sum(mu * exog[:,i,None] * exog[:,j,None] *
                    (mu * (a3 * a4 / a1**2 -
                           2 * a3**2 * a2 / a1**3 +
                           2 * a3 * (a4 + 1) / a1**2 -
                           a4 * p / (mu * a1) +
                           a3 * p * a2 / (mu * a1**2) +
                           (y - 1) * a4 * (p - 1) / (a2 * mu) -
                           (y - 1) * (1 + a4)**2 / a2**2 -
                           a4 * (p - 1) / (a1 * mu)) +
                     ((y - 1) * (1 + a4) / a2 -
                      (1 + a4) / a1)), axis=0)
                hess_arr[i, j] = np.squeeze(hess_val)
        tri_idx = np.triu_indices(dim, k=1)
        hess_arr[tri_idx] = hess_arr.T[tri_idx]

        # for dl/dparams dalpha
        dldpda = np.sum((2 * a4 * mu_p / a1**2 -
                         2 * a3 * mu_p * a2 / a1**3 -
                         mu_p * y * (y - 1) * (1 + a4) / a2**2 +
                         mu_p * (1 + a4) / a1**2 +
                         a5 * y * (y - 1) / a2 -
                         2 * a5 * y / a1 +
                         a5 * a2 / a1**2) * dmudb,
                        axis=0)

        hess_arr[-1,:-1] = dldpda
        hess_arr[:-1,-1] = dldpda

        # for dl/dalpha dalpha
        dldada = mu_p**2 * (3 * y / a1**2 -
                            (y / a2)**2. * (y - 1) -
                            2 * a2 / a1**3)

        hess_arr[-1,-1] = dldada.sum()

        return hess_arr

    def hessian_factor(self, params):
        """
        Generalized Poisson model Hessian matrix of the loglikelihood

        Parameters
        ----------
        params : array-like
            The parameters of the model

        Returns
        -------
        hess : ndarray, (nobs, 3)
            The Hessian factor, second derivative of loglikelihood function
            with respect to linear predictor and dispersion parameter
            evaluated at `params`
            The first column contains the second derivative w.r.t. linpred,
            the second column contains the cross derivative, and the
            third column contains the second derivative w.r.t. the dispersion
            parameter.

        """
        params = np.asarray(params)
        if self._transparams:
            alpha = np.exp(params[-1])
        else:
            alpha = params[-1]

        params = params[:-1]
        p = self.parameterization
        y = self.endog
        mu = self.predict(params)
        mu_p = np.power(mu, p)
        a1 = 1 + alpha * mu_p
        a2 = mu + alpha * mu_p * y
        a3 = alpha * p * mu ** (p - 1)
        a4 = a3 * y
        a5 = p * mu ** (p - 1)
        dmudb = mu

        dbb = mu * (
             mu * (a3 * a4 / a1**2 -
                   2 * a3**2 * a2 / a1**3 +
                   2 * a3 * (a4 + 1) / a1**2 -
                   a4 * p / (mu * a1) +
                   a3 * p * a2 / (mu * a1**2) +
                   a4 / (mu * a1) -
                   a3 * a2 / (mu * a1**2) +
                   (y - 1) * a4 * (p - 1) / (a2 * mu) -
                   (y - 1) * (1 + a4)**2 / a2**2 -
                   a4 * (p - 1) / (a1 * mu) -
                   1 / mu**2) +
             (-a4 / a1 +
              a3 * a2 / a1**2 +
              (y - 1) * (1 + a4) / a2 -
              (1 + a4) / a1 +
              1 / mu))

        # for dl/dlinpred dalpha
        dba = ((2 * a4 * mu_p / a1**2 -
                         2 * a3 * mu_p * a2 / a1**3 -
                         mu_p * y * (y - 1) * (1 + a4) / a2**2 +
                         mu_p * (1 + a4) / a1**2 +
                         a5 * y * (y - 1) / a2 -
                         2 * a5 * y / a1 +
                         a5 * a2 / a1**2) * dmudb)

        # for dl/dalpha dalpha
        daa = mu_p**2 * (3 * y / a1**2 -
                            (y / a2)**2. * (y - 1) -
                            2 * a2 / a1**3)

        return dbb, dba, daa

    @Appender(Poisson.predict.__doc__)
    def predict(self, params, exog=None, exposure=None, offset=None,
                which='mean', y_values=None):

        if exog is None:
            exog = self.exog

        if exposure is None:
            exposure = getattr(self, 'exposure', 0)
        elif exposure != 0:
            exposure = np.log(exposure)

        if offset is None:
            offset = getattr(self, 'offset', 0)

        fitted = np.dot(exog, params[:exog.shape[1]])
        linpred = fitted + exposure + offset

        if which == 'mean':
            return np.exp(linpred)
        elif which == 'linear':
            return linpred
        elif which == 'var':
            mean = np.exp(linpred)
            alpha = params[-1]
            pm1 = self.parameterization  # `p - 1` in GPP
            var_ = mean * (1 + alpha * mean**pm1)**2
            return var_
        elif which == 'prob':
            if y_values is None:
                y_values = np.atleast_2d(np.arange(0, np.max(self.endog)+1))
            mu = self.predict(params, exog=exog, exposure=exposure,
                              offset=offset)[:, None]
            return genpoisson_p.pmf(y_values, mu, params[-1],
                                    self.parameterization + 1)
        else:
            raise ValueError('keyword \'which\' not recognized')

    def _deriv_score_obs_dendog(self, params):
        """derivative of score_obs w.r.t. endog

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

        Returns
        -------
        derivative : ndarray_2d
            The derivative of the score_obs with respect to endog.
        """
        # code duplication with NegativeBinomialP
        from statsmodels.tools.numdiff import _approx_fprime_cs_scalar

        def f(y):
            if y.ndim == 2 and y.shape[1] == 1:
                y = y[:, 0]
            sf = self.score_factor(params, endog=y)
            return np.column_stack(sf)

        dsf = _approx_fprime_cs_scalar(self.endog[:, None], f)
        # deriv is 2d vector
        d1 = dsf[:, :1] * self.exog
        d2 = dsf[:, 1:2]

        return np.column_stack((d1, d2))

    def _var(self, mu, params=None):
        """variance implied by the distribution

        internal use, will be refactored or removed
        """
        alpha = params[-1]
        pm1 = self.parameterization  # `p-1` in GPP
        var_ = mu * (1 + alpha * mu**pm1)**2
        return var_

    def _prob_nonzero(self, mu, params):
        """Probability that count is not zero

        internal use in Censored model, will be refactored or removed
        """
        alpha = params[-1]
        pm1 = self.parameterization  # p-1 in GPP
        prob_zero = np.exp(- mu / (1 + alpha * mu**pm1))
        prob_nz = 1 - prob_zero
        return prob_nz

    @Appender(Poisson.get_distribution.__doc__)
    def get_distribution(self, params, exog=None, exposure=None, offset=None):
        """get frozen instance of distribution
        """
        mu = self.predict(params, exog=exog, exposure=exposure, offset=offset)
        p = self.parameterization + 1
        # distr = genpoisson_p(mu[:, None], params[-1], p)
        distr = genpoisson_p(mu, params[-1], p)
        return distr


class Logit(BinaryModel):
    __doc__ = """
    Logit Model

    {params}
    offset : array_like
        Offset is added to the linear prediction with coefficient equal to 1.
    {extra_params}

    Attributes
    ----------
    endog : ndarray
        A reference to the endogenous response variable
    exog : ndarray
        A reference to the exogenous design.
    """.format(params=base._model_params_doc,
           extra_params=base._missing_param_doc + _check_rank_doc)

    _continuous_ok = True

    @cache_readonly
    def link(self):
        from statsmodels.genmod.families import links
        link = links.Logit()
        return link

    def cdf(self, X):
        """
        The logistic cumulative distribution function

        Parameters
        ----------
        X : array_like
            `X` is the linear predictor of the logit model.  See notes.

        Returns
        -------
        1/(1 + exp(-X))

        Notes
        -----
        In the logit model,

        .. math:: \\Lambda\\left(x^{\\prime}\\beta\\right)=
                  \\text{Prob}\\left(Y=1|x\\right)=
                  \\frac{e^{x^{\\prime}\\beta}}{1+e^{x^{\\prime}\\beta}}
        """
        X = np.asarray(X)
        return 1/(1+np.exp(-X))

    def pdf(self, X):
        """
        The logistic probability density function

        Parameters
        ----------
        X : array_like
            `X` is the linear predictor of the logit model.  See notes.

        Returns
        -------
        pdf : ndarray
            The value of the Logit probability mass function, PMF, for each
            point of X. ``np.exp(-x)/(1+np.exp(-X))**2``

        Notes
        -----
        In the logit model,

        .. math:: \\lambda\\left(x^{\\prime}\\beta\\right)=\\frac{e^{-x^{\\prime}\\beta}}{\\left(1+e^{-x^{\\prime}\\beta}\\right)^{2}}
        """
        X = np.asarray(X)
        return np.exp(-X)/(1+np.exp(-X))**2

    @cache_readonly
    def family(self):
        from statsmodels.genmod import families
        return families.Binomial()

    def loglike(self, params):
        """
        Log-likelihood of logit model.

        Parameters
        ----------
        params : array_like
            The parameters of the logit model.

        Returns
        -------
        loglike : float
            The log-likelihood function of the model evaluated at `params`.
            See notes.

        Notes
        -----
        .. math::

           \\ln L=\\sum_{i}\\ln\\Lambda
           \\left(q_{i}x_{i}^{\\prime}\\beta\\right)

        Where :math:`q=2y-1`. This simplification comes from the fact that the
        logistic distribution is symmetric.
        """
        q = 2*self.endog - 1
        linpred = self.predict(params, which="linear")
        return np.sum(np.log(self.cdf(q * linpred)))

    def loglikeobs(self, params):
        """
        Log-likelihood of logit model for each observation.

        Parameters
        ----------
        params : array_like
            The parameters of the logit model.

        Returns
        -------
        loglike : ndarray
            The log likelihood for each observation of the model evaluated
            at `params`. See Notes

        Notes
        -----
        .. math::

           \\ln L=\\sum_{i}\\ln\\Lambda
           \\left(q_{i}x_{i}^{\\prime}\\beta\\right)

        for observations :math:`i=1,...,n`

        where :math:`q=2y-1`. This simplification comes from the fact that the
        logistic distribution is symmetric.
        """
        q = 2*self.endog - 1
        linpred = self.predict(params, which="linear")
        return np.log(self.cdf(q * linpred))

    def score(self, params):
        """
        Logit model score (gradient) vector of the log-likelihood

        Parameters
        ----------
        params : array_like
            The parameters of the model

        Returns
        -------
        score : ndarray, 1-D
            The score vector of the model, i.e. the first derivative of the
            loglikelihood function, evaluated at `params`

        Notes
        -----
        .. math:: \\frac{\\partial\\ln L}{\\partial\\beta}=\\sum_{i=1}^{n}\\left(y_{i}-\\Lambda_{i}\\right)x_{i}
        """

        y = self.endog
        X = self.exog
        fitted = self.predict(params)
        return np.dot(y - fitted, X)

    def score_obs(self, params):
        """
        Logit model Jacobian of the log-likelihood for each observation

        Parameters
        ----------
        params : array_like
            The parameters of the model

        Returns
        -------
        jac : array_like
            The derivative of the loglikelihood for each observation evaluated
            at `params`.

        Notes
        -----
        .. math:: \\frac{\\partial\\ln L_{i}}{\\partial\\beta}=\\left(y_{i}-\\Lambda_{i}\\right)x_{i}

        for observations :math:`i=1,...,n`
        """

        y = self.endog
        X = self.exog
        fitted = self.predict(params)
        return (y - fitted)[:,None] * X

    def score_factor(self, params):
        """
        Logit model derivative of the log-likelihood with respect to linpred.

        Parameters
        ----------
        params : array_like
            The parameters of the model

        Returns
        -------
        score_factor : array_like
            The derivative of the loglikelihood for each observation evaluated
            at `params`.

        Notes
        -----
        .. math:: \\frac{\\partial\\ln L_{i}}{\\partial\\beta}=\\left(y_{i}-\\lambda_{i}\\right)

        for observations :math:`i=1,...,n`

        where the loglinear model is assumed

        .. math:: \\ln\\lambda_{i}=x_{i}\\beta
        """
        y = self.endog
        fitted = self.predict(params)
        return (y - fitted)

    def hessian(self, params):
        """
        Logit model Hessian matrix of the log-likelihood

        Parameters
        ----------
        params : array_like
            The parameters of the model

        Returns
        -------
        hess : ndarray, (k_vars, k_vars)
            The Hessian, second derivative of loglikelihood function,
            evaluated at `params`

        Notes
        -----
        .. math:: \\frac{\\partial^{2}\\ln L}{\\partial\\beta\\partial\\beta^{\\prime}}=-\\sum_{i}\\Lambda_{i}\\left(1-\\Lambda_{i}\\right)x_{i}x_{i}^{\\prime}
        """
        X = self.exog
        L = self.predict(params)
        return -np.dot(L*(1-L)*X.T,X)

    def hessian_factor(self, params):
        """
        Logit model Hessian factor

        Parameters
        ----------
        params : array_like
            The parameters of the model

        Returns
        -------
        hess : ndarray, (nobs,)
            The Hessian factor, second derivative of loglikelihood function
            with respect to the linear predictor evaluated at `params`
        """
        L = self.predict(params)
        return -L * (1 - L)

    @Appender(DiscreteModel.fit.__doc__)
    def fit(self, start_params=None, method='newton', maxiter=35,
            full_output=1, disp=1, callback=None, **kwargs):
        bnryfit = super().fit(start_params=start_params,
                              method=method,
                              maxiter=maxiter,
                              full_output=full_output,
                              disp=disp,
                              callback=callback,
                              **kwargs)

        discretefit = LogitResults(self, bnryfit)
        return BinaryResultsWrapper(discretefit)

    def _deriv_score_obs_dendog(self, params):
        """derivative of score_obs w.r.t. endog

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

        Returns
        -------
        derivative : ndarray_2d
            The derivative of the score_obs with respect to endog. This
            can is given by `score_factor0[:, None] * exog` where
            `score_factor0` is the score_factor without the residual.
        """
        return self.exog


class Probit(BinaryModel):
    __doc__ = """
    Probit Model

    {params}
    offset : array_like
        Offset is added to the linear prediction with coefficient equal to 1.
    {extra_params}

    Attributes
    ----------
    endog : ndarray
        A reference to the endogenous response variable
    exog : ndarray
        A reference to the exogenous design.
    """.format(params=base._model_params_doc,
           extra_params=base._missing_param_doc + _check_rank_doc)

    @cache_readonly
    def link(self):
        from statsmodels.genmod.families import links
        link = links.Probit()
        return link

    def cdf(self, X):
        """
        Probit (Normal) cumulative distribution function

        Parameters
        ----------
        X : array_like
            The linear predictor of the model (XB).

        Returns
        -------
        cdf : ndarray
            The cdf evaluated at `X`.

        Notes
        -----
        This function is just an alias for scipy.stats.norm.cdf
        """
        return stats.norm._cdf(X)

    def pdf(self, X):
        """
        Probit (Normal) probability density function

        Parameters
        ----------
        X : array_like
            The linear predictor of the model (XB).

        Returns
        -------
        pdf : ndarray
            The value of the normal density function for each point of X.

        Notes
        -----
        This function is just an alias for scipy.stats.norm.pdf
        """
        X = np.asarray(X)
        return stats.norm._pdf(X)


    def loglike(self, params):
        """
        Log-likelihood of probit model (i.e., the normal distribution).

        Parameters
        ----------
        params : array_like
            The parameters of the model.

        Returns
        -------
        loglike : float
            The log-likelihood function of the model evaluated at `params`.
            See notes.

        Notes
        -----
        .. math:: \\ln L=\\sum_{i}\\ln\\Phi\\left(q_{i}x_{i}^{\\prime}\\beta\\right)

        Where :math:`q=2y-1`. This simplification comes from the fact that the
        normal distribution is symmetric.
        """

        q = 2*self.endog - 1
        linpred = self.predict(params, which="linear")
        return np.sum(np.log(np.clip(self.cdf(q * linpred), FLOAT_EPS, 1)))

    def loglikeobs(self, params):
        """
        Log-likelihood of probit model for each observation

        Parameters
        ----------
        params : array_like
            The parameters of the model.

        Returns
        -------
        loglike : array_like
            The log likelihood for each observation of the model evaluated
            at `params`. See Notes

        Notes
        -----
        .. math:: \\ln L_{i}=\\ln\\Phi\\left(q_{i}x_{i}^{\\prime}\\beta\\right)

        for observations :math:`i=1,...,n`

        where :math:`q=2y-1`. This simplification comes from the fact that the
        normal distribution is symmetric.
        """

        q = 2*self.endog - 1
        linpred = self.predict(params, which="linear")
        return np.log(np.clip(self.cdf(q*linpred), FLOAT_EPS, 1))


    def score(self, params):
        """
        Probit model score (gradient) vector

        Parameters
        ----------
        params : array_like
            The parameters of the model

        Returns
        -------
        score : ndarray, 1-D
            The score vector of the model, i.e. the first derivative of the
            loglikelihood function, evaluated at `params`

        Notes
        -----
        .. math:: \\frac{\\partial\\ln L}{\\partial\\beta}=\\sum_{i=1}^{n}\\left[\\frac{q_{i}\\phi\\left(q_{i}x_{i}^{\\prime}\\beta\\right)}{\\Phi\\left(q_{i}x_{i}^{\\prime}\\beta\\right)}\\right]x_{i}

        Where :math:`q=2y-1`. This simplification comes from the fact that the
        normal distribution is symmetric.
        """
        y = self.endog
        X = self.exog
        XB = self.predict(params, which="linear")
        q = 2*y - 1
        # clip to get rid of invalid divide complaint
        L = q*self.pdf(q*XB)/np.clip(self.cdf(q*XB), FLOAT_EPS, 1 - FLOAT_EPS)
        return np.dot(L,X)

    def score_obs(self, params):
        """
        Probit model Jacobian for each observation

        Parameters
        ----------
        params : array_like
            The parameters of the model

        Returns
        -------
        jac : array_like
            The derivative of the loglikelihood for each observation evaluated
            at `params`.

        Notes
        -----
        .. math:: \\frac{\\partial\\ln L_{i}}{\\partial\\beta}=\\left[\\frac{q_{i}\\phi\\left(q_{i}x_{i}^{\\prime}\\beta\\right)}{\\Phi\\left(q_{i}x_{i}^{\\prime}\\beta\\right)}\\right]x_{i}

        for observations :math:`i=1,...,n`

        Where :math:`q=2y-1`. This simplification comes from the fact that the
        normal distribution is symmetric.
        """
        y = self.endog
        X = self.exog
        XB = self.predict(params, which="linear")
        q = 2*y - 1
        # clip to get rid of invalid divide complaint
        L = q*self.pdf(q*XB)/np.clip(self.cdf(q*XB), FLOAT_EPS, 1 - FLOAT_EPS)
        return L[:,None] * X

    def score_factor(self, params):
        """
        Probit model Jacobian for each observation

        Parameters
        ----------
        params : array-like
            The parameters of the model

        Returns
        -------
        score_factor : array_like (nobs,)
            The derivative of the loglikelihood function for each observation
            with respect to linear predictor evaluated at `params`

        Notes
        -----
        .. math:: \\frac{\\partial\\ln L_{i}}{\\partial\\beta}=\\left[\\frac{q_{i}\\phi\\left(q_{i}x_{i}^{\\prime}\\beta\\right)}{\\Phi\\left(q_{i}x_{i}^{\\prime}\\beta\\right)}\\right]x_{i}

        for observations :math:`i=1,...,n`

        Where :math:`q=2y-1`. This simplification comes from the fact that the
        normal distribution is symmetric.
        """
        y = self.endog
        XB = self.predict(params, which="linear")
        q = 2*y - 1
        # clip to get rid of invalid divide complaint
        L = q*self.pdf(q*XB)/np.clip(self.cdf(q*XB), FLOAT_EPS, 1 - FLOAT_EPS)
        return L


    def hessian(self, params):
        """
        Probit model Hessian matrix of the log-likelihood

        Parameters
        ----------
        params : array_like
            The parameters of the model

        Returns
        -------
        hess : ndarray, (k_vars, k_vars)
            The Hessian, second derivative of loglikelihood function,
            evaluated at `params`

        Notes
        -----
        .. math:: \\frac{\\partial^{2}\\ln L}{\\partial\\beta\\partial\\beta^{\\prime}}=-\\lambda_{i}\\left(\\lambda_{i}+x_{i}^{\\prime}\\beta\\right)x_{i}x_{i}^{\\prime}

        where

        .. math:: \\lambda_{i}=\\frac{q_{i}\\phi\\left(q_{i}x_{i}^{\\prime}\\beta\\right)}{\\Phi\\left(q_{i}x_{i}^{\\prime}\\beta\\right)}

        and :math:`q=2y-1`
        """
        X = self.exog
        XB = self.predict(params, which="linear")
        q = 2*self.endog - 1
        L = q*self.pdf(q*XB)/self.cdf(q*XB)
        return np.dot(-L*(L+XB)*X.T,X)

    def hessian_factor(self, params):
        """
        Probit model Hessian factor of the log-likelihood

        Parameters
        ----------
        params : array-like
            The parameters of the model

        Returns
        -------
        hess : ndarray, (nobs,)
            The Hessian factor, second derivative of loglikelihood function
            with respect to linear predictor evaluated at `params`

        Notes
        -----
        .. math:: \\frac{\\partial^{2}\\ln L}{\\partial\\beta\\partial\\beta^{\\prime}}=-\\lambda_{i}\\left(\\lambda_{i}+x_{i}^{\\prime}\\beta\\right)x_{i}x_{i}^{\\prime}

        where

        .. math:: \\lambda_{i}=\\frac{q_{i}\\phi\\left(q_{i}x_{i}^{\\prime}\\beta\\right)}{\\Phi\\left(q_{i}x_{i}^{\\prime}\\beta\\right)}

        and :math:`q=2y-1`
        """
        XB = self.predict(params, which="linear")
        q = 2 * self.endog - 1
        L = q * self.pdf(q * XB) / self.cdf(q * XB)
        return -L * (L + XB)

    @Appender(DiscreteModel.fit.__doc__)
    def fit(self, start_params=None, method='newton', maxiter=35,
            full_output=1, disp=1, callback=None, **kwargs):
        bnryfit = super().fit(start_params=start_params,
                              method=method,
                              maxiter=maxiter,
                              full_output=full_output,
                              disp=disp,
                              callback=callback,
                              **kwargs)
        discretefit = ProbitResults(self, bnryfit)
        return BinaryResultsWrapper(discretefit)

    def _deriv_score_obs_dendog(self, params):
        """derivative of score_obs w.r.t. endog

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

        Returns
        -------
        derivative : ndarray_2d
            The derivative of the score_obs with respect to endog. This
            can is given by `score_factor0[:, None] * exog` where
            `score_factor0` is the score_factor without the residual.
        """

        linpred = self.predict(params, which="linear")

        pdf_ = self.pdf(linpred)
        # clip to get rid of invalid divide complaint
        cdf_ = np.clip(self.cdf(linpred), FLOAT_EPS, 1 - FLOAT_EPS)
        deriv = pdf_ / cdf_ / (1 - cdf_)  # deriv factor
        return deriv[:, None] * self.exog


class MNLogit(MultinomialModel):
    __doc__ = """
    Multinomial Logit Model

    Parameters
    ----------
    endog : array_like
        `endog` is an 1-d vector of the endogenous response.  `endog` can
        contain strings, ints, or floats or may be a pandas Categorical Series.
        Note that if it contains strings, every distinct string will be a
        category.  No stripping of whitespace is done.
    exog : array_like
        A nobs x k array where `nobs` is the number of observations and `k`
        is the number of regressors. An intercept is not included by default
        and should be added by the user. See `statsmodels.tools.add_constant`.
    {extra_params}

    Attributes
    ----------
    endog : ndarray
        A reference to the endogenous response variable
    exog : ndarray
        A reference to the exogenous design.
    J : float
        The number of choices for the endogenous variable. Note that this
        is zero-indexed.
    K : float
        The actual number of parameters for the exogenous design.  Includes
        the constant if the design has one.
    names : dict
        A dictionary mapping the column number in `wendog` to the variables
        in `endog`.
    wendog : ndarray
        An n x j array where j is the number of unique categories in `endog`.
        Each column of j is a dummy variable indicating the category of
        each observation. See `names` for a dictionary mapping each column to
        its category.

    Notes
    -----
    See developer notes for further information on `MNLogit` internals.
    """.format(extra_params=base._missing_param_doc + _check_rank_doc)

    def __init__(self, endog, exog, check_rank=True, **kwargs):
        super().__init__(endog, exog, check_rank=check_rank, **kwargs)

        # Override cov_names since multivariate model
        yname = self.endog_names
        ynames = self._ynames_map
        ynames = MultinomialResults._maybe_convert_ynames_int(ynames)
        # use range below to ensure sortedness
        ynames = [ynames[key] for key in range(int(self.J))]
        idx = MultiIndex.from_product((ynames[1:], self.data.xnames),
                                      names=(yname, None))
        self.data.cov_names = idx

    def pdf(self, eXB):
        """
        NotImplemented
        """
        raise NotImplementedError

    def cdf(self, X):
        """
        Multinomial logit cumulative distribution function.

        Parameters
        ----------
        X : ndarray
            The linear predictor of the model XB.

        Returns
        -------
        cdf : ndarray
            The cdf evaluated at `X`.

        Notes
        -----
        In the multinomial logit model.
        .. math:: \\frac{\\exp\\left(\\beta_{j}^{\\prime}x_{i}\\right)}{\\sum_{k=0}^{J}\\exp\\left(\\beta_{k}^{\\prime}x_{i}\\right)}
        """
        eXB = np.column_stack((np.ones(len(X)), np.exp(X)))
        return eXB/eXB.sum(1)[:,None]

    def loglike(self, params):
        """
        Log-likelihood of the multinomial logit model.

        Parameters
        ----------
        params : array_like
            The parameters of the multinomial logit model.

        Returns
        -------
        loglike : float
            The log-likelihood function of the model evaluated at `params`.
            See notes.

        Notes
        -----
        .. math::

           \\ln L=\\sum_{i=1}^{n}\\sum_{j=0}^{J}d_{ij}\\ln
           \\left(\\frac{\\exp\\left(\\beta_{j}^{\\prime}x_{i}\\right)}
           {\\sum_{k=0}^{J}
           \\exp\\left(\\beta_{k}^{\\prime}x_{i}\\right)}\\right)

        where :math:`d_{ij}=1` if individual `i` chose alternative `j` and 0
        if not.
        """
        params = params.reshape(self.K, -1, order='F')
        d = self.wendog
        logprob = np.log(self.cdf(np.dot(self.exog,params)))
        return np.sum(d * logprob)

    def loglikeobs(self, params):
        """
        Log-likelihood of the multinomial logit model for each observation.

        Parameters
        ----------
        params : array_like
            The parameters of the multinomial logit model.

        Returns
        -------
        loglike : array_like
            The log likelihood for each observation of the model evaluated
            at `params`. See Notes

        Notes
        -----
        .. math::

           \\ln L_{i}=\\sum_{j=0}^{J}d_{ij}\\ln
           \\left(\\frac{\\exp\\left(\\beta_{j}^{\\prime}x_{i}\\right)}
           {\\sum_{k=0}^{J}
           \\exp\\left(\\beta_{k}^{\\prime}x_{i}\\right)}\\right)

        for observations :math:`i=1,...,n`

        where :math:`d_{ij}=1` if individual `i` chose alternative `j` and 0
        if not.
        """
        params = params.reshape(self.K, -1, order='F')
        d = self.wendog
        logprob = np.log(self.cdf(np.dot(self.exog,params)))
        return d * logprob

    def score(self, params):
        """
        Score matrix for multinomial logit model log-likelihood

        Parameters
        ----------
        params : ndarray
            The parameters of the multinomial logit model.

        Returns
        -------
        score : ndarray, (K * (J-1),)
            The 2-d score vector, i.e. the first derivative of the
            loglikelihood function, of the multinomial logit model evaluated at
            `params`.

        Notes
        -----
        .. math:: \\frac{\\partial\\ln L}{\\partial\\beta_{j}}=\\sum_{i}\\left(d_{ij}-\\frac{\\exp\\left(\\beta_{j}^{\\prime}x_{i}\\right)}{\\sum_{k=0}^{J}\\exp\\left(\\beta_{k}^{\\prime}x_{i}\\right)}\\right)x_{i}

        for :math:`j=1,...,J`

        In the multinomial model the score matrix is K x J-1 but is returned
        as a flattened array to work with the solvers.
        """
        params = params.reshape(self.K, -1, order='F')
        firstterm = self.wendog[:,1:] - self.cdf(np.dot(self.exog,
                                                  params))[:,1:]
        #NOTE: might need to switch terms if params is reshaped
        return np.dot(firstterm.T, self.exog).flatten()

    def loglike_and_score(self, params):
        """
        Returns log likelihood and score, efficiently reusing calculations.

        Note that both of these returned quantities will need to be negated
        before being minimized by the maximum likelihood fitting machinery.
        """
        params = params.reshape(self.K, -1, order='F')
        cdf_dot_exog_params = self.cdf(np.dot(self.exog, params))
        loglike_value = np.sum(self.wendog * np.log(cdf_dot_exog_params))
        firstterm = self.wendog[:, 1:] - cdf_dot_exog_params[:, 1:]
        score_array = np.dot(firstterm.T, self.exog).flatten()
        return loglike_value, score_array

    def score_obs(self, params):
        """
        Jacobian matrix for multinomial logit model log-likelihood

        Parameters
        ----------
        params : ndarray
            The parameters of the multinomial logit model.

        Returns
        -------
        jac : array_like
            The derivative of the loglikelihood for each observation evaluated
            at `params` .

        Notes
        -----
        .. math:: \\frac{\\partial\\ln L_{i}}{\\partial\\beta_{j}}=\\left(d_{ij}-\\frac{\\exp\\left(\\beta_{j}^{\\prime}x_{i}\\right)}{\\sum_{k=0}^{J}\\exp\\left(\\beta_{k}^{\\prime}x_{i}\\right)}\\right)x_{i}

        for :math:`j=1,...,J`, for observations :math:`i=1,...,n`

        In the multinomial model the score vector is K x (J-1) but is returned
        as a flattened array. The Jacobian has the observations in rows and
        the flattened array of derivatives in columns.
        """
        params = params.reshape(self.K, -1, order='F')
        firstterm = self.wendog[:,1:] - self.cdf(np.dot(self.exog,
                                                  params))[:,1:]
        #NOTE: might need to switch terms if params is reshaped
        return (firstterm[:,:,None] * self.exog[:,None,:]).reshape(self.exog.shape[0], -1)

    def hessian(self, params):
        """
        Multinomial logit Hessian matrix of the log-likelihood

        Parameters
        ----------
        params : array_like
            The parameters of the model

        Returns
        -------
        hess : ndarray, (J*K, J*K)
            The Hessian, second derivative of loglikelihood function with
            respect to the flattened parameters, evaluated at `params`

        Notes
        -----
        .. math:: \\frac{\\partial^{2}\\ln L}{\\partial\\beta_{j}\\partial\\beta_{l}}=-\\sum_{i=1}^{n}\\frac{\\exp\\left(\\beta_{j}^{\\prime}x_{i}\\right)}{\\sum_{k=0}^{J}\\exp\\left(\\beta_{k}^{\\prime}x_{i}\\right)}\\left[\\boldsymbol{1}\\left(j=l\\right)-\\frac{\\exp\\left(\\beta_{l}^{\\prime}x_{i}\\right)}{\\sum_{k=0}^{J}\\exp\\left(\\beta_{k}^{\\prime}x_{i}\\right)}\\right]x_{i}x_{l}^{\\prime}

        where
        :math:`\\boldsymbol{1}\\left(j=l\\right)` equals 1 if `j` = `l` and 0
        otherwise.

        The actual Hessian matrix has J**2 * K x K elements. Our Hessian
        is reshaped to be square (J*K, J*K) so that the solvers can use it.

        This implementation does not take advantage of the symmetry of
        the Hessian and could probably be refactored for speed.
        """
        params = params.reshape(self.K, -1, order='F')
        X = self.exog
        pr = self.cdf(np.dot(X,params))
        partials = []
        J = self.J
        K = self.K
        for i in range(J-1):
            for j in range(J-1): # this loop assumes we drop the first col.
                if i == j:
                    partials.append(\
                        -np.dot(((pr[:,i+1]*(1-pr[:,j+1]))[:,None]*X).T,X))
                else:
                    partials.append(-np.dot(((pr[:,i+1]*-pr[:,j+1])[:,None]*X).T,X))
        H = np.array(partials)
        # the developer's notes on multinomial should clear this math up
        H = np.transpose(H.reshape(J-1, J-1, K, K), (0, 2, 1, 3)).reshape((J-1)*K, (J-1)*K)
        return H


#TODO: Weibull can replaced by a survival analsysis function
# like stat's streg (The cox model as well)
#class Weibull(DiscreteModel):
#    """
#    Binary choice Weibull model
#
#    Notes
#    ------
#    This is unfinished and untested.
#    """
##TODO: add analytic hessian for Weibull
#    def initialize(self):
#        pass
#
#    def cdf(self, X):
#        """
#        Gumbell (Log Weibull) cumulative distribution function
#        """
##        return np.exp(-np.exp(-X))
#        return stats.gumbel_r.cdf(X)
#        # these two are equivalent.
#        # Greene table and discussion is incorrect.
#
#    def pdf(self, X):
#        """
#        Gumbell (LogWeibull) probability distribution function
#        """
#        return stats.gumbel_r.pdf(X)
#
#    def loglike(self, params):
#        """
#        Loglikelihood of Weibull distribution
#        """
#        X = self.exog
#        cdf = self.cdf(np.dot(X,params))
#        y = self.endog
#        return np.sum(y*np.log(cdf) + (1-y)*np.log(1-cdf))
#
#    def score(self, params):
#        y = self.endog
#        X = self.exog
#        F = self.cdf(np.dot(X,params))
#        f = self.pdf(np.dot(X,params))
#        term = (y*f/F + (1 - y)*-f/(1-F))
#        return np.dot(term,X)
#
#    def hessian(self, params):
#        hess = nd.Jacobian(self.score)
#        return hess(params)
#
#    def fit(self, start_params=None, method='newton', maxiter=35, tol=1e-08):
## The example had problems with all zero start values, Hessian = 0
#        if start_params is None:
#            start_params = OLS(self.endog, self.exog).fit().params
#        mlefit = super(Weibull, self).fit(start_params=start_params,
#                method=method, maxiter=maxiter, tol=tol)
#        return mlefit
#


class NegativeBinomial(CountModel):
    __doc__ = """
    Negative Binomial Model

    %(params)s
    %(extra_params)s

    Attributes
    ----------
    endog : ndarray
        A reference to the endogenous response variable
    exog : ndarray
        A reference to the exogenous design.

    References
    ----------
    Greene, W. 2008. "Functional forms for the negative binomial model
        for count data". Economics Letters. Volume 99, Number 3, pp.585-590.
    Hilbe, J.M. 2011. "Negative binomial regression". Cambridge University
        Press.
    """ % {'params': base._model_params_doc,
           'extra_params':
           """loglike_method : str
        Log-likelihood type. 'nb2','nb1', or 'geometric'.
        Fitted value :math:`\\mu`
        Heterogeneity parameter :math:`\\alpha`

        - nb2: Variance equal to :math:`\\mu + \\alpha\\mu^2` (most common)
        - nb1: Variance equal to :math:`\\mu + \\alpha\\mu`
        - geometric: Variance equal to :math:`\\mu + \\mu^2`
    offset : array_like
        Offset is added to the linear prediction with coefficient equal to 1.
    exposure : array_like
        Log(exposure) is added to the linear prediction with coefficient
        equal to 1.
    """ + base._missing_param_doc + _check_rank_doc}

    def __init__(self, endog, exog, loglike_method='nb2', offset=None,
                 exposure=None, missing='none', check_rank=True, **kwargs):
        super().__init__(endog,
                         exog,
                         offset=offset,
                         exposure=exposure,
                         missing=missing,
                         check_rank=check_rank,
                         **kwargs)
        self.loglike_method = loglike_method
        self._initialize()
        if loglike_method in ['nb2', 'nb1']:
            self.exog_names.append('alpha')
            self.k_extra = 1
        else:
            self.k_extra = 0
        # store keys for extras if we need to recreate model instance
        # we need to append keys that do not go to super
        self._init_keys.append('loglike_method')

    def _initialize(self):
        if self.loglike_method == 'nb2':
            self.hessian = self._hessian_nb2
            self.score = self._score_nbin
            self.loglikeobs = self._ll_nb2
            self._transparams = True  # transform lnalpha -> alpha in fit
        elif self.loglike_method == 'nb1':
            self.hessian = self._hessian_nb1
            self.score = self._score_nb1
            self.loglikeobs = self._ll_nb1
            self._transparams = True  # transform lnalpha -> alpha in fit
        elif self.loglike_method == 'geometric':
            self.hessian = self._hessian_geom
            self.score = self._score_geom
            self.loglikeobs = self._ll_geometric
        else:
            raise ValueError('Likelihood type must "nb1", "nb2" '
                             'or "geometric"')

    # Workaround to pickle instance methods
    def __getstate__(self):
        odict = self.__dict__.copy()  # copy the dict since we change it
        del odict['hessian']
        del odict['score']
        del odict['loglikeobs']
        return odict

    def __setstate__(self, indict):
        self.__dict__.update(indict)
        self._initialize()

    def _ll_nbin(self, params, alpha, Q=0):
        if np.any(np.iscomplex(params)) or np.iscomplex(alpha):
            gamma_ln = loggamma
        else:
            gamma_ln = gammaln
        endog = self.endog
        mu = self.predict(params)
        size = 1/alpha * mu**Q
        prob = size/(size+mu)
        coeff = (gamma_ln(size+endog) - gamma_ln(endog+1) -
                 gamma_ln(size))
        llf = coeff + size*np.log(prob) + endog*np.log(1-prob)
        return llf

    def _ll_nb2(self, params):
        if self._transparams:  # got lnalpha during fit
            alpha = np.exp(params[-1])
        else:
            alpha = params[-1]
        return self._ll_nbin(params[:-1], alpha, Q=0)

    def _ll_nb1(self, params):
        if self._transparams:  # got lnalpha during fit
            alpha = np.exp(params[-1])
        else:
            alpha = params[-1]
        return self._ll_nbin(params[:-1], alpha, Q=1)

    def _ll_geometric(self, params):
        # we give alpha of 1 because it's actually log(alpha) where alpha=0
        return self._ll_nbin(params, 1, 0)

    def loglike(self, params):
        r"""
        Loglikelihood for negative binomial model

        Parameters
        ----------
        params : array_like
            The parameters of the model. If `loglike_method` is nb1 or
            nb2, then the ancillary parameter is expected to be the
            last element.

        Returns
        -------
        llf : float
            The loglikelihood value at `params`

        Notes
        -----
        Following notation in Greene (2008), with negative binomial
        heterogeneity parameter :math:`\alpha`:

        .. math::

           \lambda_i &= exp(X\beta) \\
           \theta &= 1 / \alpha \\
           g_i &= \theta \lambda_i^Q \\
           w_i &= g_i/(g_i + \lambda_i) \\
           r_i &= \theta / (\theta+\lambda_i) \\
           ln \mathcal{L}_i &= ln \Gamma(y_i+g_i) - ln \Gamma(1+y_i) + g_iln (r_i) + y_i ln(1-r_i)

        where :math`Q=0` for NB2 and geometric and :math:`Q=1` for NB1.
        For the geometric, :math:`\alpha=0` as well.
        """
        llf = np.sum(self.loglikeobs(params))
        return llf

    def _score_geom(self, params):
        exog = self.exog
        y = self.endog[:, None]
        mu = self.predict(params)[:, None]
        dparams = exog * (y-mu)/(mu+1)
        return dparams.sum(0)

    def _score_nbin(self, params, Q=0):
        """
        Score vector for NB2 model
        """
        if self._transparams: # lnalpha came in during fit
            alpha = np.exp(params[-1])
        else:
            alpha = params[-1]
        params = params[:-1]
        exog = self.exog
        y = self.endog[:,None]
        mu = self.predict(params)[:,None]
        a1 = 1/alpha * mu**Q
        prob = a1 / (a1 + mu)  # a1 aka "size" in _ll_nbin
        if Q == 1:  # nb1
            # Q == 1 --> a1 = mu / alpha --> prob = 1 / (alpha + 1)
            dgpart = digamma(y + a1) - digamma(a1)
            dparams = exog * a1 * (np.log(prob) +
                       dgpart)
            dalpha = ((alpha * (y - mu * np.log(prob) -
                              mu*(dgpart + 1)) -
                       mu * (np.log(prob) +
                           dgpart))/
                       (alpha**2*(alpha + 1))).sum()

        elif Q == 0:  # nb2
            dgpart = digamma(y + a1) - digamma(a1)
            dparams = exog*a1 * (y-mu)/(mu+a1)
            da1 = -alpha**-2
            dalpha = (dgpart + np.log(a1)
                        - np.log(a1+mu) - (y-mu)/(a1+mu)).sum() * da1

        #multiply above by constant outside sum to reduce rounding error
        if self._transparams:
            return np.r_[dparams.sum(0), dalpha*alpha]
        else:
            return np.r_[dparams.sum(0), dalpha]

    def _score_nb1(self, params):
        return self._score_nbin(params, Q=1)

    def _hessian_geom(self, params):
        exog = self.exog
        y = self.endog[:,None]
        mu = self.predict(params)[:,None]

        # for dl/dparams dparams
        dim = exog.shape[1]
        hess_arr = np.empty((dim, dim))
        const_arr = mu*(1+y)/(mu+1)**2
        for i in range(dim):
            for j in range(dim):
                if j > i:
                    continue
                hess_arr[i,j] = np.squeeze(
                    np.sum(-exog[:,i,None] * exog[:,j,None] * const_arr,
                           axis=0
                           )
                )
        tri_idx = np.triu_indices(dim, k=1)
        hess_arr[tri_idx] = hess_arr.T[tri_idx]
        return hess_arr


    def _hessian_nb1(self, params):
        """
        Hessian of NB1 model.
        """
        if self._transparams: # lnalpha came in during fit
            alpha = np.exp(params[-1])
        else:
            alpha = params[-1]

        params = params[:-1]
        exog = self.exog
        y = self.endog[:,None]
        mu = self.predict(params)[:,None]

        a1 = mu/alpha
        dgpart = digamma(y + a1) - digamma(a1)
        prob = 1 / (1 + alpha)  # equiv: a1 / (a1 + mu)

        # for dl/dparams dparams
        dim = exog.shape[1]
        hess_arr = np.empty((dim+1,dim+1))
        #const_arr = a1*mu*(a1+y)/(mu+a1)**2
        # not all of dparams
        dparams = exog / alpha * (np.log(prob) +
                                  dgpart)

        dmudb = exog*mu
        xmu_alpha = exog * a1
        trigamma = (special.polygamma(1, a1 + y) -
                    special.polygamma(1, a1))
        for i in range(dim):
            for j in range(dim):
                if j > i:
                    continue
                hess_arr[i,j] = np.squeeze(
                    np.sum(
                        dparams[:,i,None] * dmudb[:,j,None] +
                        xmu_alpha[:,i,None] * xmu_alpha[:,j,None] * trigamma,
                        axis=0
                    )
                )
        tri_idx = np.triu_indices(dim, k=1)
        hess_arr[tri_idx] = hess_arr.T[tri_idx]

        # for dl/dparams dalpha
        # da1 = -alpha**-2
        dldpda = np.sum(-a1 * dparams + exog * a1 *
                        (-trigamma*mu/alpha**2 - prob), axis=0)

        hess_arr[-1,:-1] = dldpda
        hess_arr[:-1,-1] = dldpda

        log_alpha = np.log(prob)
        alpha3 = alpha**3
        alpha2 = alpha**2
        mu2 = mu**2
        dada = ((alpha3*mu*(2*log_alpha + 2*dgpart + 3) -
                 2*alpha3*y +
                 4*alpha2*mu*(log_alpha + dgpart) +
                 alpha2 * (2*mu - y) +
                 2*alpha*mu2*trigamma + mu2 * trigamma + alpha2 * mu2 * trigamma +
                 2*alpha*mu*(log_alpha + dgpart)
                 )/(alpha**4*(alpha2 + 2*alpha + 1)))
        hess_arr[-1,-1] = dada.sum()

        return hess_arr

    def _hessian_nb2(self, params):
        """
        Hessian of NB2 model.
        """
        if self._transparams: # lnalpha came in during fit
            alpha = np.exp(params[-1])
        else:
            alpha = params[-1]
        a1 = 1/alpha
        params = params[:-1]

        exog = self.exog
        y = self.endog[:,None]
        mu = self.predict(params)[:,None]
        prob = a1 / (a1 + mu)
        dgpart = digamma(a1 + y) - digamma(a1)

        # for dl/dparams dparams
        dim = exog.shape[1]
        hess_arr = np.empty((dim+1,dim+1))
        const_arr = a1*mu*(a1+y)/(mu+a1)**2
        for i in range(dim):
            for j in range(dim):
                if j > i:
                    continue
                hess_arr[i,j] = np.sum(-exog[:,i,None] * exog[:,j,None] *
                                       const_arr, axis=0).squeeze()
        tri_idx = np.triu_indices(dim, k=1)
        hess_arr[tri_idx] = hess_arr.T[tri_idx]

        # for dl/dparams dalpha
        da1 = -alpha**-2
        dldpda = -np.sum(mu*exog*(y-mu)*a1**2/(mu+a1)**2 , axis=0)
        hess_arr[-1,:-1] = dldpda
        hess_arr[:-1,-1] = dldpda

        # for dl/dalpha dalpha
        #NOTE: polygamma(1,x) is the trigamma function
        da2 = 2*alpha**-3
        dalpha = da1 * (dgpart +
                    np.log(prob) - (y - mu)/(a1+mu))
        dada = (da2 * dalpha/da1 + da1**2 * (special.polygamma(1, a1+y) -
                    special.polygamma(1, a1) + 1/a1 - 1/(a1 + mu) +
                    (y - mu)/(mu + a1)**2)).sum()
        hess_arr[-1,-1] = dada

        return hess_arr

    #TODO: replace this with analytic where is it used?
    def score_obs(self, params):
        sc = approx_fprime_cs(params, self.loglikeobs)
        return sc

    @Appender(Poisson.predict.__doc__)
    def predict(self, params, exog=None, exposure=None, offset=None,
                which='mean', linear=None, y_values=None):

        if linear is not None:
            msg = 'linear keyword is deprecated, use which="linear"'
            warnings.warn(msg, FutureWarning)
            if linear is True:
                which = "linear"

        # avoid duplicate computation for get-distribution
        if which == "prob":
            distr = self.get_distribution(
                params,
                exog=exog,
                exposure=exposure,
                offset=offset
                )
            if y_values is None:
                y_values = np.arange(0, np.max(self.endog) + 1)
            else:
                y_values = np.asarray(y_values)

            assert y_values.ndim == 1
            y_values = y_values[..., None]
            return distr.pmf(y_values).T

        exog, offset, exposure = self._get_predict_arrays(
            exog=exog,
            offset=offset,
            exposure=exposure
            )

        fitted = np.dot(exog, params[:exog.shape[1]])
        linpred = fitted + exposure + offset
        if which == "mean":
            return np.exp(linpred)
        elif which.startswith("lin"):
            return linpred
        elif which == "var":
            mu = np.exp(linpred)
            if self.loglike_method == 'geometric':
                var_ = mu * (1 + mu)
            else:
                if self.loglike_method == 'nb2':
                    p = 2
                elif self.loglike_method == 'nb1':
                    p = 1
                alpha = params[-1]
                var_ = mu * (1 + alpha * mu**(p - 1))
            return var_
        else:
            raise ValueError('keyword which has to be "mean" and "linear"')

    @Appender(_get_start_params_null_docs)
    def _get_start_params_null(self):
        offset = getattr(self, "offset", 0)
        exposure = getattr(self, "exposure", 0)
        const = (self.endog / np.exp(offset + exposure)).mean()
        params = [np.log(const)]
        mu = const * np.exp(offset + exposure)
        resid = self.endog - mu
        a = self._estimate_dispersion(mu, resid, df_resid=resid.shape[0] - 1)
        params.append(a)
        return np.array(params)

    def _estimate_dispersion(self, mu, resid, df_resid=None):
        if df_resid is None:
            df_resid = resid.shape[0]
        if self.loglike_method == 'nb2':
            #params.append(np.linalg.pinv(mu[:,None]).dot(resid**2 / mu - 1))
            a = ((resid**2 / mu - 1) / mu).sum() / df_resid
        else: #self.loglike_method == 'nb1':
            a = (resid**2 / mu - 1).sum() / df_resid
        return a

    def fit(self, start_params=None, method='bfgs', maxiter=35,
            full_output=1, disp=1, callback=None,
            cov_type='nonrobust', cov_kwds=None, use_t=None,
            optim_kwds_prelim=None, **kwargs):

        # Note: do not let super handle robust covariance because it has
        # transformed params
        self._transparams = False # always define attribute
        if self.loglike_method.startswith('nb') and method not in ['newton',
                                                                   'ncg']:
            self._transparams = True # in case same Model instance is refit
        elif self.loglike_method.startswith('nb'): # method is newton/ncg
            self._transparams = False # because we need to step in alpha space

        if start_params is None:
            # Use poisson fit as first guess.
            #TODO, Warning: this assumes exposure is logged
            offset = getattr(self, "offset", 0) + getattr(self, "exposure", 0)
            if np.size(offset) == 1 and offset == 0:
                offset = None
            kwds_prelim = {'disp': 0, 'skip_hessian': True, 'warn_convergence': False}
            if optim_kwds_prelim is not None:
                kwds_prelim.update(optim_kwds_prelim)
            mod_poi = Poisson(self.endog, self.exog, offset=offset)
            with warnings.catch_warnings():
                warnings.simplefilter("always")
                res_poi = mod_poi.fit(**kwds_prelim)
            start_params = res_poi.params
            if self.loglike_method.startswith('nb'):
                a = self._estimate_dispersion(res_poi.predict(), res_poi.resid,
                                              df_resid=res_poi.df_resid)
                start_params = np.append(start_params, max(0.05, a))
        else:
            if self._transparams is True:
                # transform user provided start_params dispersion, see #3918
                start_params = np.array(start_params, copy=True)
                start_params[-1] = np.log(start_params[-1])

        if callback is None:
            # work around perfect separation callback #3895
            callback = lambda *x: x

        mlefit = super().fit(start_params=start_params,
                             maxiter=maxiter, method=method, disp=disp,
                             full_output=full_output, callback=callback,
                             **kwargs)
        if optim_kwds_prelim is not None:
            mlefit.mle_settings["optim_kwds_prelim"] = optim_kwds_prelim
        # TODO: Fix NBin _check_perfect_pred
        if self.loglike_method.startswith('nb'):
            # mlefit is a wrapped counts results
            self._transparams = False # do not need to transform anymore now
            # change from lnalpha to alpha
            if method not in ["newton", "ncg"]:
                mlefit._results.params[-1] = np.exp(mlefit._results.params[-1])

            nbinfit = NegativeBinomialResults(self, mlefit._results)
            result = NegativeBinomialResultsWrapper(nbinfit)
        else:
            result = mlefit

        if cov_kwds is None:
            cov_kwds = {}  #TODO: make this unnecessary ?
        result._get_robustcov_results(cov_type=cov_type, use_self=True, use_t=use_t, **cov_kwds)
        return result


    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)

        if self.loglike_method.startswith('nb') and (np.size(alpha) == 1 and
                                                     alpha != 0):
            # do not penalize alpha if alpha is scalar
            k_params = self.exog.shape[1] + self.k_extra
            alpha = alpha * np.ones(k_params)
            alpha[-1] = 0

        # alpha for regularized poisson to get starting values
        alpha_p = alpha[:-1] if (self.k_extra and np.size(alpha) > 1) else alpha

        self._transparams = False
        if start_params is None:
            # Use poisson fit as first guess.
            #TODO, Warning: this assumes exposure is logged
            offset = getattr(self, "offset", 0) + getattr(self, "exposure", 0)
            if np.size(offset) == 1 and offset == 0:
                offset = None
            mod_poi = Poisson(self.endog, self.exog, offset=offset)
            with warnings.catch_warnings():
                warnings.simplefilter("always")
                start_params = mod_poi.fit_regularized(
                    start_params=start_params, method=method, maxiter=maxiter,
                    full_output=full_output, disp=0, callback=callback,
                    alpha=alpha_p, trim_mode=trim_mode,
                    auto_trim_tol=auto_trim_tol, size_trim_tol=size_trim_tol,
                    qc_tol=qc_tol, **kwargs).params
            if self.loglike_method.startswith('nb'):
                start_params = np.append(start_params, 0.1)

        cntfit = super(CountModel, self).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 = L1NegativeBinomialResults(self, cntfit)
        return L1NegativeBinomialResultsWrapper(discretefit)

    @Appender(Poisson.get_distribution.__doc__)
    def get_distribution(self, params, exog=None, exposure=None, offset=None):
        """get frozen instance of distribution
        """
        mu = self.predict(params, exog=exog, exposure=exposure, offset=offset)
        if self.loglike_method == 'geometric':
            # distr = stats.geom(1 / (1 + mu[:, None]), loc=-1)
            distr = stats.geom(1 / (1 + mu), loc=-1)
        else:
            if self.loglike_method == 'nb2':
                p = 2
            elif self.loglike_method == 'nb1':
                p = 1

            alpha = params[-1]
            q = 2 - p
            size = 1. / alpha * mu**q
            prob = size / (size + mu)
            # distr = nbinom(size[:, None], prob[:, None])
            distr = nbinom(size, prob)

        return distr


class NegativeBinomialP(CountModel):
    __doc__ = """
    Generalized Negative Binomial (NB-P) Model

    %(params)s
    %(extra_params)s

    Attributes
    ----------
    endog : ndarray
        A reference to the endogenous response variable
    exog : ndarray
        A reference to the exogenous design.
    p : scalar
        P denotes parameterizations for NB-P regression. p=1 for NB-1 and
        p=2 for NB-2. Default is p=1.
    """ % {'params': base._model_params_doc,
           'extra_params':
               """p : scalar
        P denotes parameterizations for NB regression. p=1 for NB-1 and
        p=2 for NB-2. Default is p=2.
    offset : array_like
        Offset is added to the linear prediction with coefficient equal to 1.
    exposure : array_like
        Log(exposure) is added to the linear prediction with coefficient
        equal to 1.
        """ + base._missing_param_doc + _check_rank_doc}

    def __init__(self, endog, exog, p=2, offset=None,
                 exposure=None, missing='none', check_rank=True,
                 **kwargs):
        super().__init__(endog,
                         exog,
                         offset=offset,
                         exposure=exposure,
                         missing=missing,
                         check_rank=check_rank,
                         **kwargs)
        self.parameterization = p
        self.exog_names.append('alpha')
        self.k_extra = 1
        self._transparams = False

    def _get_init_kwds(self):
        kwds = super()._get_init_kwds()
        kwds['p'] = self.parameterization
        return kwds

    def _get_exogs(self):
        return (self.exog, None)

    def loglike(self, params):
        """
        Loglikelihood of Generalized Negative Binomial (NB-P) model

        Parameters
        ----------
        params : array_like
            The parameters of the model.

        Returns
        -------
        loglike : float
            The log-likelihood function of the model evaluated at `params`.
            See notes.
        """
        return np.sum(self.loglikeobs(params))

    def loglikeobs(self, params):
        """
        Loglikelihood for observations of Generalized Negative Binomial (NB-P) model

        Parameters
        ----------
        params : array_like
            The parameters of the model.

        Returns
        -------
        loglike : ndarray
            The log likelihood for each observation of the model evaluated
            at `params`. See Notes
        """
        if self._transparams:
            alpha = np.exp(params[-1])
        else:
            alpha = params[-1]

        params = params[:-1]
        p = self.parameterization
        y = self.endog

        mu = self.predict(params)
        mu_p = mu**(2 - p)
        a1 = mu_p / alpha
        a2 = mu + a1

        llf = (gammaln(y + a1) - gammaln(y + 1) - gammaln(a1) +
               a1 * np.log(a1) + y * np.log(mu) -
               (y + a1) * np.log(a2))

        return llf

    def score_obs(self, params):
        """
        Generalized Negative Binomial (NB-P) model score (gradient) vector of the log-likelihood for each observations.

        Parameters
        ----------
        params : array_like
            The parameters of the model

        Returns
        -------
        score : ndarray, 1-D
            The score vector of the model, i.e. the first derivative of the
            loglikelihood function, evaluated at `params`
        """
        if self._transparams:
            alpha = np.exp(params[-1])
        else:
            alpha = params[-1]

        params = params[:-1]
        p = 2 - self.parameterization
        y = self.endog

        mu = self.predict(params)
        mu_p = mu**p
        a1 = mu_p / alpha
        a2 = mu + a1
        a3 = y + a1
        a4 = p * a1 / mu

        dgpart = digamma(a3) - digamma(a1)
        dgterm = dgpart + np.log(a1 / a2) + 1 - a3 / a2
        # TODO: better name/interpretation for dgterm?

        dparams = (a4 * dgterm -
                   a3 / a2 +
                   y / mu)
        dparams = (self.exog.T * mu * dparams).T
        dalpha = -a1 / alpha * dgterm

        return np.concatenate((dparams, np.atleast_2d(dalpha).T),
                              axis=1)

    def score(self, params):
        """
        Generalized Negative Binomial (NB-P) model score (gradient) vector of the log-likelihood

        Parameters
        ----------
        params : array_like
            The parameters of the model

        Returns
        -------
        score : ndarray, 1-D
            The score vector of the model, i.e. the first derivative of the
            loglikelihood function, evaluated at `params`
        """
        score = np.sum(self.score_obs(params), axis=0)
        if self._transparams:
            score[-1] == score[-1] ** 2
            return score
        else:
            return score

    def score_factor(self, params, endog=None):
        """
        Generalized Negative Binomial (NB-P) model score (gradient) vector of the log-likelihood for each observations.

        Parameters
        ----------
        params : array-like
            The parameters of the model

        Returns
        -------
        score : ndarray, 1-D
            The score vector of the model, i.e. the first derivative of the
            loglikelihood function, evaluated at `params`
        """
        params = np.asarray(params)
        if self._transparams:
            alpha = np.exp(params[-1])
        else:
            alpha = params[-1]

        params = params[:-1]
        p = 2 - self.parameterization
        y = self.endog if endog is None else endog

        mu = self.predict(params)
        mu_p = mu**p
        a1 = mu_p / alpha
        a2 = mu + a1
        a3 = y + a1
        a4 = p * a1 / mu

        dgpart = digamma(a3) - digamma(a1)

        dparams = ((a4 * dgpart -
                   a3 / a2) +
                   y / mu + a4 * (1 - a3 / a2 + np.log(a1 / a2)))
        dparams = (mu * dparams).T
        dalpha = (-a1 / alpha * (dgpart +
                                 np.log(a1 / a2) +
                                 1 - a3 / a2))

        return dparams, dalpha

    def hessian(self, params):
        """
        Generalized Negative Binomial (NB-P) model hessian maxtrix of the log-likelihood

        Parameters
        ----------
        params : array_like
            The parameters of the model

        Returns
        -------
        hessian : ndarray, 2-D
            The hessian matrix of the model.
        """
        if self._transparams:
            alpha = np.exp(params[-1])
        else:
            alpha = params[-1]
        params = params[:-1]

        p = 2 - self.parameterization
        y = self.endog
        exog = self.exog
        mu = self.predict(params)

        mu_p = mu**p
        a1 = mu_p / alpha
        a2 = mu + a1
        a3 = y + a1
        a4 = p * a1 / mu

        prob = a1 / a2
        lprob = np.log(prob)
        dgpart = digamma(a3) - digamma(a1)
        pgpart = polygamma(1, a3) - polygamma(1, a1)

        dim = exog.shape[1]
        hess_arr = np.zeros((dim + 1, dim + 1))

        coeff = mu**2 * (((1 + a4)**2 * a3 / a2**2 -
                          a3 / a2 * (p - 1) * a4 / mu -
                          y / mu**2 -
                          2 * a4 * (1 + a4) / a2 +
                          p * a4 / mu * (lprob + dgpart + 2) -
                          a4 / mu * (lprob + dgpart + 1) +
                          a4**2 * pgpart) +
                         (-(1 + a4) * a3 / a2 +
                          y / mu +
                          a4 * (lprob + dgpart + 1)) / mu)

        for i in range(dim):
            hess_arr[i, :-1] = np.sum(self.exog[:, :].T * self.exog[:, i] * coeff, axis=1)


        hess_arr[-1,:-1] = (self.exog[:, :].T * mu * a1 *
                ((1 + a4) * (1 - a3 / a2) / a2 -
                 p * (lprob + dgpart + 2) / mu +
                 p / mu * (a3 + p * a1) / a2 -
                 a4 * pgpart) / alpha).sum(axis=1)


        da2 = (a1 * (2 * lprob +
                     2 * dgpart + 3 -
                     2 * a3 / a2
                     + a1 * pgpart
                     - 2 * prob +
                     prob * a3 / a2) / alpha**2)

        hess_arr[-1, -1] = da2.sum()

        tri_idx = np.triu_indices(dim + 1, k=1)
        hess_arr[tri_idx] = hess_arr.T[tri_idx]

        return hess_arr

    def hessian_factor(self, params):
        """
        Generalized Negative Binomial (NB-P) model hessian maxtrix of the log-likelihood

        Parameters
        ----------
        params : array-like
            The parameters of the model

        Returns
        -------
        hessian : ndarray, 2-D
            The hessian matrix of the model.
        """
        params = np.asarray(params)
        if self._transparams:
            alpha = np.exp(params[-1])
        else:
            alpha = params[-1]
        params = params[:-1]

        p = 2 - self.parameterization
        y = self.endog
        mu = self.predict(params)

        mu_p = mu**p
        a1 = mu_p / alpha
        a2 = mu + a1
        a3 = y + a1
        a4 = p * a1 / mu
        a5 = a4 * p / mu

        dgpart = digamma(a3) - digamma(a1)

        coeff = mu**2 * (((1 + a4)**2 * a3 / a2**2 -
                          a3 * (a5 - a4 / mu) / a2 -
                          y / mu**2 -
                          2 * a4 * (1 + a4) / a2 +
                          a5 * (np.log(a1) - np.log(a2) + dgpart + 2) -
                          a4 * (np.log(a1) - np.log(a2) + dgpart + 1) / mu -
                          a4**2 * (polygamma(1, a1) - polygamma(1, a3))) +
                         (-(1 + a4) * a3 / a2 +
                          y / mu +
                          a4 * (np.log(a1) - np.log(a2) + dgpart + 1)) / mu)

        hfbb = coeff

        hfba = (mu * a1 *
                ((1 + a4) * (1 - a3 / a2) / a2 -
                 p * (np.log(a1 / a2) + dgpart + 2) / mu +
                 p * (a3 / mu + a4) / a2 +
                 a4 * (polygamma(1, a1) - polygamma(1, a3))) / alpha)

        hfaa = (a1 * (2 * np.log(a1 / a2) +
                     2 * dgpart + 3 -
                     2 * a3 / a2 - a1 * polygamma(1, a1) +
                     a1 * polygamma(1, a3) - 2 * a1 / a2 +
                     a1 * a3 / a2**2) / alpha**2)

        return hfbb, hfba, hfaa

    @Appender(_get_start_params_null_docs)
    def _get_start_params_null(self):
        offset = getattr(self, "offset", 0)
        exposure = getattr(self, "exposure", 0)

        const = (self.endog / np.exp(offset + exposure)).mean()
        params = [np.log(const)]
        mu = const * np.exp(offset + exposure)
        resid = self.endog - mu
        a = self._estimate_dispersion(mu, resid, df_resid=resid.shape[0] - 1)
        params.append(a)

        return np.array(params)

    def _estimate_dispersion(self, mu, resid, df_resid=None):
        q = self.parameterization - 1
        if df_resid is None:
            df_resid = resid.shape[0]
        a = ((resid**2 / mu - 1) * mu**(-q)).sum() / df_resid
        return a

    @Appender(DiscreteModel.fit.__doc__)
    def fit(self, start_params=None, method='bfgs', maxiter=35,
            full_output=1, disp=1, callback=None, use_transparams=False,
            cov_type='nonrobust', cov_kwds=None, use_t=None,
            optim_kwds_prelim=None, **kwargs):
        # TODO: Fix doc string
        """
        use_transparams : bool
            This parameter enable internal transformation to impose
            non-negativity. True to enable. Default is False.
            use_transparams=True imposes the no underdispersion (alpha > 0)
            constraint. In case use_transparams=True and method="newton" or
            "ncg" transformation is ignored.
        """
        if use_transparams and method not in ['newton', 'ncg']:
            self._transparams = True
        else:
            if use_transparams:
                warnings.warn('Parameter "use_transparams" is ignored',
                              RuntimeWarning)
            self._transparams = False
        if start_params is None:
            offset = getattr(self, "offset", 0) + getattr(self, "exposure", 0)
            if np.size(offset) == 1 and offset == 0:
                offset = None
            kwds_prelim = {'disp': 0, 'skip_hessian': True, 'warn_convergence': False}
            if optim_kwds_prelim is not None:
                kwds_prelim.update(optim_kwds_prelim)
            mod_poi = Poisson(self.endog, self.exog, offset=offset)
            with warnings.catch_warnings():
                warnings.simplefilter("always")
                res_poi = mod_poi.fit(**kwds_prelim)
            start_params = res_poi.params
            a = self._estimate_dispersion(res_poi.predict(), res_poi.resid,
                                          df_resid=res_poi.df_resid)
            start_params = np.append(start_params, max(0.05, a))

        if callback is None:
            # work around perfect separation callback #3895
            callback = lambda *x: x

        mlefit = super().fit(start_params=start_params,
                        maxiter=maxiter, method=method, disp=disp,
                        full_output=full_output, callback=callback,
                        **kwargs)
        if optim_kwds_prelim is not None:
            mlefit.mle_settings["optim_kwds_prelim"] = optim_kwds_prelim
        if use_transparams and method not in ["newton", "ncg"]:
            self._transparams = False
            mlefit._results.params[-1] = np.exp(mlefit._results.params[-1])

        nbinfit = NegativeBinomialPResults(self, mlefit._results)
        result = NegativeBinomialPResultsWrapper(nbinfit)

        if cov_kwds is None:
            cov_kwds = {}
        result._get_robustcov_results(cov_type=cov_type,
                                    use_self=True, use_t=use_t, **cov_kwds)
        return result

    @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)

        if np.size(alpha) == 1 and alpha != 0:
            k_params = self.exog.shape[1] + self.k_extra
            alpha = alpha * np.ones(k_params)
            alpha[-1] = 0

        alpha_p = alpha[:-1] if (self.k_extra and np.size(alpha) > 1) else alpha

        self._transparams = False
        if start_params is None:
            offset = getattr(self, "offset", 0) + getattr(self, "exposure", 0)
            if np.size(offset) == 1 and offset == 0:
                offset = None
            mod_poi = Poisson(self.endog, self.exog, offset=offset)
            with warnings.catch_warnings():
                warnings.simplefilter("always")
                start_params = mod_poi.fit_regularized(
                    start_params=start_params, method=method, maxiter=maxiter,
                    full_output=full_output, disp=0, callback=callback,
                    alpha=alpha_p, trim_mode=trim_mode,
                    auto_trim_tol=auto_trim_tol, size_trim_tol=size_trim_tol,
                    qc_tol=qc_tol, **kwargs).params
            start_params = np.append(start_params, 0.1)

        cntfit = super(CountModel, self).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 = L1NegativeBinomialResults(self, cntfit)

        return L1NegativeBinomialResultsWrapper(discretefit)

    @Appender(Poisson.predict.__doc__)
    def predict(self, params, exog=None, exposure=None, offset=None,
                which='mean', y_values=None):

        if exog is None:
            exog = self.exog

        if exposure is None:
            exposure = getattr(self, 'exposure', 0)
        elif exposure != 0:
            exposure = np.log(exposure)

        if offset is None:
            offset = getattr(self, 'offset', 0)

        fitted = np.dot(exog, params[:exog.shape[1]])
        linpred = fitted + exposure + offset

        if which == 'mean':
            return np.exp(linpred)
        elif which == 'linear':
            return linpred
        elif which == 'var':
            mean = np.exp(linpred)
            alpha = params[-1]
            p = self.parameterization  # no `-1` as in GPP
            var_ = mean * (1 + alpha * mean**(p - 1))
            return var_
        elif which == 'prob':
            if y_values is None:
                y_values = np.atleast_2d(np.arange(0, np.max(self.endog)+1))

            mu = self.predict(params, exog, exposure, offset)
            size, prob = self.convert_params(params, mu)
            return nbinom.pmf(y_values, size[:, None], prob[:, None])
        else:
            raise ValueError('keyword "which" = %s not recognized' % which)

    def convert_params(self, params, mu):
        alpha = params[-1]
        p = 2 - self.parameterization

        size = 1. / alpha * mu**p
        prob = size / (size + mu)

        return (size, prob)

    def _deriv_score_obs_dendog(self, params):
        """derivative of score_obs w.r.t. endog

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

        Returns
        -------
        derivative : ndarray_2d
            The derivative of the score_obs with respect to endog.
        """
        from statsmodels.tools.numdiff import _approx_fprime_cs_scalar

        def f(y):
            if y.ndim == 2 and y.shape[1] == 1:
                y = y[:, 0]
            sf = self.score_factor(params, endog=y)
            return np.column_stack(sf)

        dsf = _approx_fprime_cs_scalar(self.endog[:, None], f)
        # deriv is 2d vector
        d1 = dsf[:, :1] * self.exog
        d2 = dsf[:, 1:2]

        return np.column_stack((d1, d2))

    def _var(self, mu, params=None):
        """variance implied by the distribution

        internal use, will be refactored or removed
        """
        alpha = params[-1]
        p = self.parameterization  # no `-1` as in GPP
        var_ = mu * (1 + alpha * mu**(p - 1))
        return var_

    def _prob_nonzero(self, mu, params):
        """Probability that count is not zero

        internal use in Censored model, will be refactored or removed
        """
        alpha = params[-1]
        p = self.parameterization
        prob_nz = 1 - (1 + alpha * mu**(p-1))**(- 1 / alpha)
        return prob_nz

    @Appender(Poisson.get_distribution.__doc__)
    def get_distribution(self, params, exog=None, exposure=None, offset=None):
        """get frozen instance of distribution
        """
        mu = self.predict(params, exog=exog, exposure=exposure, offset=offset)
        size, prob = self.convert_params(params, mu)
        # distr = nbinom(size[:, None], prob[:, None])
        distr = nbinom(size, prob)
        return distr


### Results Class ###

class DiscreteResults(base.LikelihoodModelResults):
    __doc__ = _discrete_results_docs % {"one_line_description" :
        "A results class for the discrete dependent variable models.",
        "extra_attr" : ""}

    def __init__(self, model, mlefit, cov_type='nonrobust', cov_kwds=None,
                 use_t=None):
        #super(DiscreteResults, self).__init__(model, params,
        #        np.linalg.inv(-hessian), scale=1.)
        self.model = model
        self.method = "MLE"
        self.df_model = model.df_model
        self.df_resid = model.df_resid
        self._cache = {}
        self.nobs = model.exog.shape[0]
        self.__dict__.update(mlefit.__dict__)
        self.converged = mlefit.mle_retvals["converged"]

        if not hasattr(self, 'cov_type'):
            # do this only if super, i.e. mlefit did not already add cov_type
            # robust covariance
            if use_t is not None:
                self.use_t = use_t
            if cov_type == 'nonrobust':
                self.cov_type = 'nonrobust'
                self.cov_kwds = {'description' : 'Standard Errors assume that the ' +
                                 'covariance matrix of the errors is correctly ' +
                                 'specified.'}
            else:
                if cov_kwds is None:
                    cov_kwds = {}
                from statsmodels.base.covtype import get_robustcov_results
                get_robustcov_results(self, cov_type=cov_type, use_self=True,
                                           **cov_kwds)


    def __getstate__(self):
        # remove unpicklable methods
        mle_settings = getattr(self, 'mle_settings', None)
        if mle_settings is not None:
            if 'callback' in mle_settings:
                mle_settings['callback'] = None
            if 'cov_params_func' in mle_settings:
                mle_settings['cov_params_func'] = None
        return self.__dict__

    @cache_readonly
    def prsquared(self):
        """
        McFadden's pseudo-R-squared. `1 - (llf / llnull)`
        """
        return 1 - self.llf/self.llnull

    @cache_readonly
    def llr(self):
        """
        Likelihood ratio chi-squared statistic; `-2*(llnull - llf)`
        """
        return -2*(self.llnull - self.llf)

    @cache_readonly
    def llr_pvalue(self):
        """
        The chi-squared probability of getting a log-likelihood ratio
        statistic greater than llr.  llr has a chi-squared distribution
        with degrees of freedom `df_model`.
        """
        return stats.distributions.chi2.sf(self.llr, self.df_model)

    def set_null_options(self, llnull=None, attach_results=True, **kwargs):
        """
        Set the fit options for the Null (constant-only) model.

        This resets the cache for related attributes which is potentially
        fragile. This only sets the option, the null model is estimated
        when llnull is accessed, if llnull is not yet in cache.

        Parameters
        ----------
        llnull : {None, float}
            If llnull is not None, then the value will be directly assigned to
            the cached attribute "llnull".
        attach_results : bool
            Sets an internal flag whether the results instance of the null
            model should be attached. By default without calling this method,
            thenull model results are not attached and only the loglikelihood
            value llnull is stored.
        **kwargs
            Additional keyword arguments used as fit keyword arguments for the
            null model. The override and model default values.

        Notes
        -----
        Modifies attributes of this instance, and so has no return.
        """
        # reset cache, note we need to add here anything that depends on
        # llnullor the null model. If something is missing, then the attribute
        # might be incorrect.
        self._cache.pop('llnull', None)
        self._cache.pop('llr', None)
        self._cache.pop('llr_pvalue', None)
        self._cache.pop('prsquared', None)
        if hasattr(self, 'res_null'):
            del self.res_null

        if llnull is not None:
            self._cache['llnull'] = llnull
        self._attach_nullmodel = attach_results
        self._optim_kwds_null = kwargs

    @cache_readonly
    def llnull(self):
        """
        Value of the constant-only loglikelihood
        """
        model = self.model
        kwds = model._get_init_kwds().copy()
        for key in getattr(model, '_null_drop_keys', []):
            del kwds[key]
        # TODO: what parameters to pass to fit?
        mod_null = model.__class__(model.endog, np.ones(self.nobs), **kwds)
        # TODO: consider catching and warning on convergence failure?
        # in the meantime, try hard to converge. see
        # TestPoissonConstrained1a.test_smoke

        optim_kwds = getattr(self, '_optim_kwds_null', {}).copy()

        if 'start_params' in optim_kwds:
            # user provided
            sp_null = optim_kwds.pop('start_params')
        elif hasattr(model, '_get_start_params_null'):
            # get moment estimates if available
            sp_null = model._get_start_params_null()
        else:
            sp_null = None

        opt_kwds = dict(method='bfgs', warn_convergence=False, maxiter=10000,
                        disp=0)
        opt_kwds.update(optim_kwds)

        if optim_kwds:
            res_null = mod_null.fit(start_params=sp_null, **opt_kwds)
        else:
            # this should be a reasonably method case across versions
            res_null = mod_null.fit(start_params=sp_null, method='nm',
                                    warn_convergence=False,
                                    maxiter=10000, disp=0)
            res_null = mod_null.fit(start_params=res_null.params, method='bfgs',
                                    warn_convergence=False,
                                    maxiter=10000, disp=0)

        if getattr(self, '_attach_nullmodel', False) is not False:
            self.res_null = res_null

        return res_null.llf

    @cache_readonly
    def fittedvalues(self):
        """
        Linear predictor XB.
        """
        return np.dot(self.model.exog, self.params[:self.model.exog.shape[1]])

    @cache_readonly
    def resid_response(self):
        """
        Respnose residuals. The response residuals are defined as
        `endog - fittedvalues`
        """
        return self.model.endog - self.predict()

    @cache_readonly
    def resid_pearson(self):
        """
        Pearson residuals defined as response residuals divided by standard
        deviation implied by the model.
        """
        var_ = self.predict(which="var")
        return self.resid_response / np.sqrt(var_)

    @cache_readonly
    def aic(self):
        """
        Akaike information criterion.  `-2*(llf - p)` where `p` is the number
        of regressors including the intercept.
        """
        k_extra = getattr(self.model, 'k_extra', 0)
        return -2*(self.llf - (self.df_model + 1 + k_extra))

    @cache_readonly
    def bic(self):
        """
        Bayesian information criterion. `-2*llf + ln(nobs)*p` where `p` is the
        number of regressors including the intercept.
        """
        k_extra = getattr(self.model, 'k_extra', 0)
        return -2*self.llf + np.log(self.nobs)*(self.df_model + 1 + k_extra)

    @cache_readonly
    def im_ratio(self):
        return pinfer.im_ratio(self)

    def info_criteria(self, crit, dk_params=0):
        """Return an information criterion for the model.

        Parameters
        ----------
        crit : string
            One of 'aic', 'bic', 'tic' or 'gbic'.
        dk_params : int or float
            Correction to the number of parameters used in the information
            criterion.

        Returns
        -------
        Value of information criterion.

        Notes
        -----
        Tic and gbic

        References
        ----------
        Burnham KP, Anderson KR (2002). Model Selection and Multimodel
        Inference; Springer New York.
        """
        crit = crit.lower()
        k_extra = getattr(self.model, 'k_extra', 0)
        k_params = self.df_model + 1 + k_extra + dk_params

        if crit == "aic":
            return -2 * self.llf + 2 * k_params
        elif crit == "bic":
            nobs = self.df_model + self.df_resid + 1
            bic = -2*self.llf + k_params*np.log(nobs)
            return bic
        elif crit == "tic":
            return pinfer.tic(self)
        elif crit == "gbic":
            return pinfer.gbic(self)
        else:
            raise ValueError("Name of information criterion not recognized.")

    def score_test(self, exog_extra=None, params_constrained=None,
                   hypothesis='joint', cov_type=None, cov_kwds=None,
                   k_constraints=None, observed=True):

        res = pinfer.score_test(self, exog_extra=exog_extra,
                                params_constrained=params_constrained,
                                hypothesis=hypothesis,
                                cov_type=cov_type, cov_kwds=cov_kwds,
                                k_constraints=k_constraints,
                                observed=observed)
        return res

    score_test.__doc__ = pinfer.score_test.__doc__

    def get_prediction(self, exog=None,
                       transform=True, which="mean", linear=None,
                       row_labels=None, average=False,
                       agg_weights=None, y_values=None,
                       **kwargs):
        """
        Compute prediction results when endpoint transformation is valid.

        Parameters
        ----------
        exog : array_like, optional
            The values for which you want to predict.
        transform : bool, optional
            If the model was fit via a formula, do you want to pass
            exog through the formula. Default is True. E.g., if you fit
            a model y ~ log(x1) + log(x2), and transform is True, then
            you can pass a data structure that contains x1 and x2 in
            their original form. Otherwise, you'd need to log the data
            first.
        which : str
            Which statistic is to be predicted. Default is "mean".
            The available statistics and options depend on the model.
            see the model.predict docstring
        linear : bool
            Linear has been replaced by the `which` keyword and will be
            deprecated.
            If linear is True, then `which` is ignored and the linear
            prediction is returned.
        row_labels : list of str or None
            If row_lables are provided, then they will replace the generated
            labels.
        average : bool
            If average is True, then the mean prediction is computed, that is,
            predictions are computed for individual exog and then the average
            over observation is used.
            If average is False, then the results are the predictions for all
            observations, i.e. same length as ``exog``.
        agg_weights : ndarray, optional
            Aggregation weights, only used if average is True.
            The weights are not normalized.
        y_values : None or nd_array
            Some predictive statistics like which="prob" are computed at
            values of the response variable. If y_values is not None, then
            it will be used instead of the default set of y_values.

            **Warning:** ``which="prob"`` for count models currently computes
            the pmf for all y=k up to max(endog). This can be a large array if
            the observed endog values are large.
            This will likely change so that the set of y_values will be chosen
            to limit the array size.
        **kwargs :
            Some models can take additional keyword arguments, such as offset,
            exposure or additional exog in multi-part models like zero inflated
            models.
            See the predict method of the model for the details.

        Returns
        -------
        prediction_results : PredictionResults
            The prediction results instance contains prediction and prediction
            variance and can on demand calculate confidence intervals and
            summary dataframe for the prediction.

        Notes
        -----
        Status: new in 0.14, experimental
        """

        if linear is True:
            # compatibility with old keyword
            which = "linear"

        pred_kwds = kwargs
        # y_values is explicit so we can add it to the docstring
        if y_values is not None:
            pred_kwds["y_values"] = y_values

        res = pred.get_prediction(
            self,
            exog=exog,
            which=which,
            transform=transform,
            row_labels=row_labels,
            average=average,
            agg_weights=agg_weights,
            pred_kwds=pred_kwds
            )
        return res

    def get_distribution(self, exog=None, transform=True, **kwargs):

        exog, _ = self._transform_predict_exog(exog, transform=transform)
        if exog is not None:
            exog = np.asarray(exog)
        distr = self.model.get_distribution(self.params,
                                            exog=exog,
                                            **kwargs
                                            )
        return distr

    def _get_endog_name(self, yname, yname_list):
        if yname is None:
            yname = self.model.endog_names
        if yname_list is None:
            yname_list = self.model.endog_names
        return yname, yname_list

    def get_margeff(self, at='overall', method='dydx', atexog=None,
            dummy=False, count=False):
        """Get marginal effects of the fitted model.

        Parameters
        ----------
        at : str, optional
            Options are:

            - 'overall', The average of the marginal effects at each
              observation.
            - 'mean', The marginal effects at the mean of each regressor.
            - 'median', The marginal effects at the median of each regressor.
            - 'zero', The marginal effects at zero for each regressor.
            - 'all', The marginal effects at each observation. If `at` is all
              only margeff will be available from the returned object.

            Note that if `exog` is specified, then marginal effects for all
            variables not specified by `exog` are calculated using the `at`
            option.
        method : str, optional
            Options are:

            - 'dydx' - dy/dx - No transformation is made and marginal effects
              are returned.  This is the default.
            - 'eyex' - estimate elasticities of variables in `exog` --
              d(lny)/d(lnx)
            - 'dyex' - estimate semi-elasticity -- dy/d(lnx)
            - 'eydx' - estimate semi-elasticity -- d(lny)/dx

            Note that tranformations are done after each observation is
            calculated.  Semi-elasticities for binary variables are computed
            using the midpoint method. 'dyex' and 'eyex' do not make sense
            for discrete variables. For interpretations of these methods
            see notes below.
        atexog : array_like, optional
            Optionally, you can provide the exogenous variables over which to
            get the marginal effects.  This should be a dictionary with the key
            as the zero-indexed column number and the value of the dictionary.
            Default is None for all independent variables less the constant.
        dummy : bool, optional
            If False, treats binary variables (if present) as continuous.  This
            is the default.  Else if True, treats binary variables as
            changing from 0 to 1.  Note that any variable that is either 0 or 1
            is treated as binary.  Each binary variable is treated separately
            for now.
        count : bool, optional
            If False, treats count variables (if present) as continuous.  This
            is the default.  Else if True, the marginal effect is the
            change in probabilities when each observation is increased by one.

        Returns
        -------
        DiscreteMargins : marginal effects instance
            Returns an object that holds the marginal effects, standard
            errors, confidence intervals, etc. See
            `statsmodels.discrete.discrete_margins.DiscreteMargins` for more
            information.

        Notes
        -----
        Interpretations of methods:

        - 'dydx' - change in `endog` for a change in `exog`.
        - 'eyex' - proportional change in `endog` for a proportional change
          in `exog`.
        - 'dyex' - change in `endog` for a proportional change in `exog`.
        - 'eydx' - proportional change in `endog` for a change in `exog`.

        When using after Poisson, returns the expected number of events per
        period, assuming that the model is loglinear.
        """
        if getattr(self.model, "offset", None) is not None:
            raise NotImplementedError("Margins with offset are not available.")
        from statsmodels.discrete.discrete_margins import DiscreteMargins
        return DiscreteMargins(self, (at, method, atexog, dummy, count))

    def get_influence(self):
        """
        Get an instance of MLEInfluence with influence and outlier measures

        Returns
        -------
        infl : MLEInfluence instance
            The instance has methods to calculate the main influence and
            outlier measures as attributes.

        See Also
        --------
        statsmodels.stats.outliers_influence.MLEInfluence
        """
        from statsmodels.stats.outliers_influence import MLEInfluence
        return MLEInfluence(self)

    def summary(self, yname=None, xname=None, title=None, alpha=.05,
                yname_list=None):
        """
        Summarize the Regression Results.

        Parameters
        ----------
        yname : str, optional
            The name of the endog variable in the tables. The default is `y`.
        xname : list[str], optional
            The names for the exogenous variables, default is "var_xx".
            Must match the number of parameters in the model.
        title : str, optional
            Title for the top table. If not None, then this replaces the
            default title.
        alpha : float
            The significance level for the confidence intervals.

        Returns
        -------
        Summary
            Class that holds the summary tables and text, which can be printed
            or converted to various output formats.

        See Also
        --------
        statsmodels.iolib.summary.Summary : Class that hold summary results.
        """

        top_left = [('Dep. Variable:', None),
                     ('Model:', [self.model.__class__.__name__]),
                     ('Method:', [self.method]),
                     ('Date:', None),
                     ('Time:', None),
                     ('converged:', ["%s" % self.mle_retvals['converged']]),
                    ]

        top_right = [('No. Observations:', None),
                     ('Df Residuals:', None),
                     ('Df Model:', None),
                     ('Pseudo R-squ.:', ["%#6.4g" % self.prsquared]),
                     ('Log-Likelihood:', None),
                     ('LL-Null:', ["%#8.5g" % self.llnull]),
                     ('LLR p-value:', ["%#6.4g" % self.llr_pvalue])
                     ]

        if hasattr(self, 'cov_type'):
            top_left.append(('Covariance Type:', [self.cov_type]))

        if title is None:
            title = self.model.__class__.__name__ + ' ' + "Regression Results"

        # boiler plate
        from statsmodels.iolib.summary import Summary
        smry = Summary()
        yname, yname_list = self._get_endog_name(yname, yname_list)

        # for top of table
        smry.add_table_2cols(self, gleft=top_left, gright=top_right,
                             yname=yname, xname=xname, title=title)

        # for parameters, etc
        smry.add_table_params(self, yname=yname_list, xname=xname, alpha=alpha,
                              use_t=self.use_t)

        if hasattr(self, 'constraints'):
            smry.add_extra_txt(['Model has been estimated subject to linear '
                                'equality constraints.'])

        return smry

    def summary2(self, yname=None, xname=None, title=None, alpha=.05,
                 float_format="%.4f"):
        """
        Experimental function to summarize regression results.

        Parameters
        ----------
        yname : str
            Name of the dependent variable (optional).
        xname : list[str], optional
            List of strings of length equal to the number of parameters
            Names of the independent variables (optional).
        title : str, optional
            Title for the top table. If not None, then this replaces the
            default title.
        alpha : float
            The significance level for the confidence intervals.
        float_format : str
            The print format for floats in parameters summary.

        Returns
        -------
        Summary
            Instance that contains the summary tables and text, which can be
            printed or converted to various output formats.

        See Also
        --------
        statsmodels.iolib.summary2.Summary : Class that holds summary results.
        """
        from statsmodels.iolib import summary2
        smry = summary2.Summary()
        smry.add_base(results=self, alpha=alpha, float_format=float_format,
                      xname=xname, yname=yname, title=title)

        if hasattr(self, 'constraints'):
            smry.add_text('Model has been estimated subject to linear '
                          'equality constraints.')

        return smry


class CountResults(DiscreteResults):
    __doc__ = _discrete_results_docs % {
        "one_line_description": "A results class for count data",
        "extra_attr": ""}

    @cache_readonly
    def resid(self):
        """
        Residuals

        Notes
        -----
        The residuals for Count models are defined as

        .. math:: y - p

        where :math:`p = \\exp(X\\beta)`. Any exposure and offset variables
        are also handled.
        """
        return self.model.endog - self.predict()

    def get_diagnostic(self, y_max=None):
        """
        Get instance of class with specification and diagnostic methods.

        experimental, API of Diagnostic classes will change

        Returns
        -------
        CountDiagnostic instance
            The instance has methods to perform specification and diagnostic
            tesst and plots

        See Also
        --------
        statsmodels.statsmodels.discrete.diagnostic.CountDiagnostic
        """
        from statsmodels.discrete.diagnostic import CountDiagnostic
        return CountDiagnostic(self, y_max=y_max)


class NegativeBinomialResults(CountResults):
    __doc__ = _discrete_results_docs % {
        "one_line_description": "A results class for NegativeBinomial 1 and 2",
        "extra_attr": ""}

    @cache_readonly
    def lnalpha(self):
        """Natural log of alpha"""
        return np.log(self.params[-1])

    @cache_readonly
    def lnalpha_std_err(self):
        """Natural log of standardized error"""
        return self.bse[-1] / self.params[-1]

    @cache_readonly
    def aic(self):
        # + 1 because we estimate alpha
        k_extra = getattr(self.model, 'k_extra', 0)
        return -2*(self.llf - (self.df_model + self.k_constant + k_extra))

    @cache_readonly
    def bic(self):
        # + 1 because we estimate alpha
        k_extra = getattr(self.model, 'k_extra', 0)
        return -2*self.llf + np.log(self.nobs)*(self.df_model +
                                                self.k_constant + k_extra)


class NegativeBinomialPResults(NegativeBinomialResults):
    __doc__ = _discrete_results_docs % {
        "one_line_description": "A results class for NegativeBinomialP",
        "extra_attr": ""}


class GeneralizedPoissonResults(NegativeBinomialResults):
    __doc__ = _discrete_results_docs % {
        "one_line_description": "A results class for Generalized Poisson",
        "extra_attr": ""}

    @cache_readonly
    def _dispersion_factor(self):
        p = getattr(self.model, 'parameterization', 0)
        mu = self.predict()
        return (1 + self.params[-1] * mu**p)**2


class L1CountResults(DiscreteResults):
    __doc__ = _discrete_results_docs % {"one_line_description" :
            "A results class for count data fit by l1 regularization",
            "extra_attr" : _l1_results_attr}

    def __init__(self, model, cntfit):
        super().__init__(model, cntfit)
        # self.trimmed is a boolean array with T/F telling whether or not that
        # entry in params has been set zero'd out.
        self.trimmed = cntfit.mle_retvals['trimmed']
        self.nnz_params = (~self.trimmed).sum()

        # Set degrees of freedom.  In doing so,
        # adjust for extra parameter in NegativeBinomial nb1 and nb2
        # extra parameter is not included in df_model
        k_extra = getattr(self.model, 'k_extra', 0)

        self.df_model = self.nnz_params - 1 - k_extra
        self.df_resid = float(self.model.endog.shape[0] - self.nnz_params) + k_extra


class PoissonResults(CountResults):

    def predict_prob(self, n=None, exog=None, exposure=None, offset=None,
                     transform=True):
        """
        Return predicted probability of each count level for each observation

        Parameters
        ----------
        n : array_like or int
            The counts for which you want the probabilities. If n is None
            then the probabilities for each count from 0 to max(y) are
            given.

        Returns
        -------
        ndarray
            A nobs x n array where len(`n`) columns are indexed by the count
            n. If n is None, then column 0 is the probability that each
            observation is 0, column 1 is the probability that each
            observation is 1, etc.
        """
        if n is not None:
            counts = np.atleast_2d(n)
        else:
            counts = np.atleast_2d(np.arange(0, np.max(self.model.endog)+1))
        mu = self.predict(exog=exog, exposure=exposure, offset=offset,
                          transform=transform, which="mean")[:,None]
        # uses broadcasting
        return stats.poisson.pmf(counts, mu)

    @property
    def resid_pearson(self):
        """
        Pearson residuals

        Notes
        -----
        Pearson residuals are defined to be

        .. math:: r_j = \\frac{(y - M_jp_j)}{\\sqrt{M_jp_j(1-p_j)}}

        where :math:`p_j=cdf(X\\beta)` and :math:`M_j` is the total number of
        observations sharing the covariate pattern :math:`j`.

        For now :math:`M_j` is always set to 1.
        """
        # Pearson residuals
        p = self.predict()  # fittedvalues is still linear
        return (self.model.endog - p)/np.sqrt(p)

    def get_influence(self):
        """
        Get an instance of MLEInfluence with influence and outlier measures

        Returns
        -------
        infl : MLEInfluence instance
            The instance has methods to calculate the main influence and
            outlier measures as attributes.

        See Also
        --------
        statsmodels.stats.outliers_influence.MLEInfluence
        """
        from statsmodels.stats.outliers_influence import MLEInfluence
        return MLEInfluence(self)

    def get_diagnostic(self, y_max=None):
        """
        Get instance of class with specification and diagnostic methods

        experimental, API of Diagnostic classes will change

        Returns
        -------
        PoissonDiagnostic instance
            The instance has methods to perform specification and diagnostic
            tesst and plots

        See Also
        --------
        statsmodels.statsmodels.discrete.diagnostic.PoissonDiagnostic
        """
        from statsmodels.discrete.diagnostic import (
            PoissonDiagnostic)
        return PoissonDiagnostic(self, y_max=y_max)


class L1PoissonResults(L1CountResults, PoissonResults):
    pass

class L1NegativeBinomialResults(L1CountResults, NegativeBinomialResults):
    pass

class L1GeneralizedPoissonResults(L1CountResults, GeneralizedPoissonResults):
    pass

class OrderedResults(DiscreteResults):
    __doc__ = _discrete_results_docs % {"one_line_description" : "A results class for ordered discrete data." , "extra_attr" : ""}
    pass

class BinaryResults(DiscreteResults):
    __doc__ = _discrete_results_docs % {"one_line_description" : "A results class for binary data", "extra_attr" : ""}

    def pred_table(self, threshold=.5):
        """
        Prediction table

        Parameters
        ----------
        threshold : scalar
            Number between 0 and 1. Threshold above which a prediction is
            considered 1 and below which a prediction is considered 0.

        Notes
        -----
        pred_table[i,j] refers to the number of times "i" was observed and
        the model predicted "j". Correct predictions are along the diagonal.
        """
        model = self.model
        actual = model.endog
        pred = np.array(self.predict() > threshold, dtype=float)
        bins = np.array([0, 0.5, 1])
        return np.histogram2d(actual, pred, bins=bins)[0]

    @Appender(DiscreteResults.summary.__doc__)
    def summary(self, yname=None, xname=None, title=None, alpha=.05,
                yname_list=None):
        smry = super().summary(yname, xname, title, alpha,
                                                  yname_list)
        fittedvalues = self.model.cdf(self.fittedvalues)
        absprederror = np.abs(self.model.endog - fittedvalues)
        predclose_sum = (absprederror < 1e-4).sum()
        predclose_frac = predclose_sum / len(fittedvalues)

        # add warnings/notes
        etext = []
        if predclose_sum == len(fittedvalues):  # TODO: nobs?
            wstr = "Complete Separation: The results show that there is"
            wstr += "complete separation or perfect prediction.\n"
            wstr += "In this case the Maximum Likelihood Estimator does "
            wstr += "not exist and the parameters\n"
            wstr += "are not identified."
            etext.append(wstr)
        elif predclose_frac > 0.1:  # TODO: get better diagnosis
            wstr = "Possibly complete quasi-separation: A fraction "
            wstr += "%4.2f of observations can be\n" % predclose_frac
            wstr += "perfectly predicted. This might indicate that there "
            wstr += "is complete\nquasi-separation. In this case some "
            wstr += "parameters will not be identified."
            etext.append(wstr)
        if etext:
            smry.add_extra_txt(etext)
        return smry

    @cache_readonly
    def resid_dev(self):
        """
        Deviance residuals

        Notes
        -----
        Deviance residuals are defined

        .. math:: d_j = \\pm\\left(2\\left[Y_j\\ln\\left(\\frac{Y_j}{M_jp_j}\\right) + (M_j - Y_j\\ln\\left(\\frac{M_j-Y_j}{M_j(1-p_j)} \\right) \\right] \\right)^{1/2}

        where

        :math:`p_j = cdf(X\\beta)` and :math:`M_j` is the total number of
        observations sharing the covariate pattern :math:`j`.

        For now :math:`M_j` is always set to 1.
        """
        #These are the deviance residuals
        #model = self.model
        endog = self.model.endog
        #exog = model.exog
        # M = # of individuals that share a covariate pattern
        # so M[i] = 2 for i = two share a covariate pattern
        M = 1
        p = self.predict()
        #Y_0 = np.where(exog == 0)
        #Y_M = np.where(exog == M)
        #NOTE: Common covariate patterns are not yet handled
        res = -(1-endog)*np.sqrt(2*M*np.abs(np.log(1-p))) + \
                endog*np.sqrt(2*M*np.abs(np.log(p)))
        return res

    @cache_readonly
    def resid_pearson(self):
        """
        Pearson residuals

        Notes
        -----
        Pearson residuals are defined to be

        .. math:: r_j = \\frac{(y - M_jp_j)}{\\sqrt{M_jp_j(1-p_j)}}

        where :math:`p_j=cdf(X\\beta)` and :math:`M_j` is the total number of
        observations sharing the covariate pattern :math:`j`.

        For now :math:`M_j` is always set to 1.
        """
        # Pearson residuals
        #model = self.model
        endog = self.model.endog
        #exog = model.exog
        # M = # of individuals that share a covariate pattern
        # so M[i] = 2 for i = two share a covariate pattern
        # use unique row pattern?
        M = 1
        p = self.predict()
        return (endog - M*p)/np.sqrt(M*p*(1-p))

    @cache_readonly
    def resid_response(self):
        """
        The response residuals

        Notes
        -----
        Response residuals are defined to be

        .. math:: y - p

        where :math:`p=cdf(X\\beta)`.
        """
        return self.model.endog - self.predict()


class LogitResults(BinaryResults):
    __doc__ = _discrete_results_docs % {
        "one_line_description": "A results class for Logit Model",
        "extra_attr": ""}

    @cache_readonly
    def resid_generalized(self):
        """
        Generalized residuals

        Notes
        -----
        The generalized residuals for the Logit model are defined

        .. math:: y - p

        where :math:`p=cdf(X\\beta)`. This is the same as the `resid_response`
        for the Logit model.
        """
        # Generalized residuals
        return self.model.endog - self.predict()

    def get_influence(self):
        """
        Get an instance of MLEInfluence with influence and outlier measures

        Returns
        -------
        infl : MLEInfluence instance
            The instance has methods to calculate the main influence and
            outlier measures as attributes.

        See Also
        --------
        statsmodels.stats.outliers_influence.MLEInfluence
        """
        from statsmodels.stats.outliers_influence import MLEInfluence
        return MLEInfluence(self)


class ProbitResults(BinaryResults):
    __doc__ = _discrete_results_docs % {
        "one_line_description": "A results class for Probit Model",
        "extra_attr": ""}

    @cache_readonly
    def resid_generalized(self):
        """
        Generalized residuals

        Notes
        -----
        The generalized residuals for the Probit model are defined

        .. math:: y\\frac{\\phi(X\\beta)}{\\Phi(X\\beta)}-(1-y)\\frac{\\phi(X\\beta)}{1-\\Phi(X\\beta)}
        """
        # generalized residuals
        model = self.model
        endog = model.endog
        XB = self.predict(which="linear")
        pdf = model.pdf(XB)
        cdf = model.cdf(XB)
        return endog * pdf/cdf - (1-endog)*pdf/(1-cdf)

class L1BinaryResults(BinaryResults):
    __doc__ = _discrete_results_docs % {"one_line_description" :
    "Results instance for binary data fit by l1 regularization",
    "extra_attr" : _l1_results_attr}
    def __init__(self, model, bnryfit):
        super().__init__(model, bnryfit)
        # self.trimmed is a boolean array with T/F telling whether or not that
        # entry in params has been set zero'd out.
        self.trimmed = bnryfit.mle_retvals['trimmed']
        self.nnz_params = (~self.trimmed).sum()
        self.df_model = self.nnz_params - 1
        self.df_resid = float(self.model.endog.shape[0] - self.nnz_params)


class MultinomialResults(DiscreteResults):
    __doc__ = _discrete_results_docs % {"one_line_description" :
            "A results class for multinomial data", "extra_attr" : ""}

    def __init__(self, model, mlefit):
        super().__init__(model, mlefit)
        self.J = model.J
        self.K = model.K

    @staticmethod
    def _maybe_convert_ynames_int(ynames):
        # see if they're integers
        issue_warning = False
        msg = ('endog contains values are that not int-like. Uses string '
               'representation of value. Use integer-valued endog to '
               'suppress this warning.')
        for i in ynames:
            try:
                if ynames[i] % 1 == 0:
                    ynames[i] = str(int(ynames[i]))
                else:
                    issue_warning = True
                    ynames[i] = str(ynames[i])
            except TypeError:
                ynames[i] = str(ynames[i])
        if issue_warning:
            warnings.warn(msg, SpecificationWarning)

        return ynames

    def _get_endog_name(self, yname, yname_list, all=False):
        """
        If all is False, the first variable name is dropped
        """
        model = self.model
        if yname is None:
            yname = model.endog_names
        if yname_list is None:
            ynames = model._ynames_map
            ynames = self._maybe_convert_ynames_int(ynames)
            # use range below to ensure sortedness
            ynames = [ynames[key] for key in range(int(model.J))]
            ynames = ['='.join([yname, name]) for name in ynames]
            if not all:
                yname_list = ynames[1:] # assumes first variable is dropped
            else:
                yname_list = ynames
        return yname, yname_list

    def pred_table(self):
        """
        Returns the J x J prediction table.

        Notes
        -----
        pred_table[i,j] refers to the number of times "i" was observed and
        the model predicted "j". Correct predictions are along the diagonal.
        """
        ju = self.model.J - 1  # highest index
        # these are the actual, predicted indices
        #idx = lzip(self.model.endog, self.predict().argmax(1))
        bins = np.concatenate(([0], np.linspace(0.5, ju - 0.5, ju), [ju]))
        return np.histogram2d(self.model.endog, self.predict().argmax(1),
                              bins=bins)[0]

    @cache_readonly
    def bse(self):
        bse = np.sqrt(np.diag(self.cov_params()))
        return bse.reshape(self.params.shape, order='F')

    @cache_readonly
    def aic(self):
        return -2*(self.llf - (self.df_model+self.model.J-1))

    @cache_readonly
    def bic(self):
        return -2*self.llf + np.log(self.nobs)*(self.df_model+self.model.J-1)

    def conf_int(self, alpha=.05, cols=None):
        confint = super(DiscreteResults, self).conf_int(alpha=alpha,
                                                            cols=cols)
        return confint.transpose(2,0,1)

    def get_prediction(self):
        """Not implemented for Multinomial
        """
        raise NotImplementedError

    def margeff(self):
        raise NotImplementedError("Use get_margeff instead")

    @cache_readonly
    def resid_misclassified(self):
        """
        Residuals indicating which observations are misclassified.

        Notes
        -----
        The residuals for the multinomial model are defined as

        .. math:: argmax(y_i) \\neq argmax(p_i)

        where :math:`argmax(y_i)` is the index of the category for the
        endogenous variable and :math:`argmax(p_i)` is the index of the
        predicted probabilities for each category. That is, the residual
        is a binary indicator that is 0 if the category with the highest
        predicted probability is the same as that of the observed variable
        and 1 otherwise.
        """
        # it's 0 or 1 - 0 for correct prediction and 1 for a missed one
        return (self.model.wendog.argmax(1) !=
                self.predict().argmax(1)).astype(float)

    def summary2(self, alpha=0.05, float_format="%.4f"):
        """Experimental function to summarize regression results

        Parameters
        ----------
        alpha : float
            significance level for the confidence intervals
        float_format : str
            print format for floats in parameters summary

        Returns
        -------
        smry : Summary instance
            this holds the summary tables and text, which can be printed or
            converted to various output formats.

        See Also
        --------
        statsmodels.iolib.summary2.Summary : class to hold summary results
        """

        from statsmodels.iolib import summary2
        smry = summary2.Summary()
        smry.add_dict(summary2.summary_model(self))
        # One data frame per value of endog
        eqn = self.params.shape[1]
        confint = self.conf_int(alpha)
        for i in range(eqn):
            coefs = summary2.summary_params((self, self.params[:, i],
                                             self.bse[:, i],
                                             self.tvalues[:, i],
                                             self.pvalues[:, i],
                                             confint[i]),
                                            alpha=alpha)
            # Header must show value of endog
            level_str =  self.model.endog_names + ' = ' + str(i)
            coefs[level_str] = coefs.index
            coefs = coefs.iloc[:, [-1, 0, 1, 2, 3, 4, 5]]
            smry.add_df(coefs, index=False, header=True,
                        float_format=float_format)
            smry.add_title(results=self)
        return smry


class L1MultinomialResults(MultinomialResults):
    __doc__ = _discrete_results_docs % {"one_line_description" :
        "A results class for multinomial data fit by l1 regularization",
        "extra_attr" : _l1_results_attr}
    def __init__(self, model, mlefit):
        super().__init__(model, mlefit)
        # self.trimmed is a boolean array with T/F telling whether or not that
        # entry in params has been set zero'd out.
        self.trimmed = mlefit.mle_retvals['trimmed']
        self.nnz_params = (~self.trimmed).sum()

        # Note: J-1 constants
        self.df_model = self.nnz_params - (self.model.J - 1)
        self.df_resid = float(self.model.endog.shape[0] - self.nnz_params)


#### Results Wrappers ####

class OrderedResultsWrapper(lm.RegressionResultsWrapper):
    pass


wrap.populate_wrapper(OrderedResultsWrapper, OrderedResults)


class CountResultsWrapper(lm.RegressionResultsWrapper):
    pass


wrap.populate_wrapper(CountResultsWrapper, CountResults)


class NegativeBinomialResultsWrapper(lm.RegressionResultsWrapper):
    pass


wrap.populate_wrapper(NegativeBinomialResultsWrapper,
                      NegativeBinomialResults)


class NegativeBinomialPResultsWrapper(lm.RegressionResultsWrapper):
    pass


wrap.populate_wrapper(NegativeBinomialPResultsWrapper,
                      NegativeBinomialPResults)


class GeneralizedPoissonResultsWrapper(lm.RegressionResultsWrapper):
    pass


wrap.populate_wrapper(GeneralizedPoissonResultsWrapper,
                      GeneralizedPoissonResults)


class PoissonResultsWrapper(lm.RegressionResultsWrapper):
    pass


wrap.populate_wrapper(PoissonResultsWrapper, PoissonResults)


class L1CountResultsWrapper(lm.RegressionResultsWrapper):
    pass


class L1PoissonResultsWrapper(lm.RegressionResultsWrapper):
    pass


wrap.populate_wrapper(L1PoissonResultsWrapper, L1PoissonResults)


class L1NegativeBinomialResultsWrapper(lm.RegressionResultsWrapper):
    pass


wrap.populate_wrapper(L1NegativeBinomialResultsWrapper,
                      L1NegativeBinomialResults)


class L1GeneralizedPoissonResultsWrapper(lm.RegressionResultsWrapper):
    pass


wrap.populate_wrapper(L1GeneralizedPoissonResultsWrapper,
                      L1GeneralizedPoissonResults)


class BinaryResultsWrapper(lm.RegressionResultsWrapper):
    _attrs = {"resid_dev": "rows",
              "resid_generalized": "rows",
              "resid_pearson": "rows",
              "resid_response": "rows"
              }
    _wrap_attrs = wrap.union_dicts(lm.RegressionResultsWrapper._wrap_attrs,
                                   _attrs)


wrap.populate_wrapper(BinaryResultsWrapper, BinaryResults)


class L1BinaryResultsWrapper(lm.RegressionResultsWrapper):
    pass


wrap.populate_wrapper(L1BinaryResultsWrapper, L1BinaryResults)


class MultinomialResultsWrapper(lm.RegressionResultsWrapper):
    _attrs = {"resid_misclassified": "rows"}
    _wrap_attrs = wrap.union_dicts(lm.RegressionResultsWrapper._wrap_attrs,
                                   _attrs)
    _methods = {'conf_int': 'multivariate_confint'}
    _wrap_methods = wrap.union_dicts(lm.RegressionResultsWrapper._wrap_methods,
                                     _methods)


wrap.populate_wrapper(MultinomialResultsWrapper, MultinomialResults)


class L1MultinomialResultsWrapper(lm.RegressionResultsWrapper):
    pass


wrap.populate_wrapper(L1MultinomialResultsWrapper, L1MultinomialResults)
