Source code for statsmodels.discrete.discrete_model

"""Limited dependent variable and qualitative variables.Includes binary outcomes, count data, (ordered) ordinal data and limiteddependent variables.General References--------------------A.C. Cameron and P.K. Trivedi. `Regression Analysis of Count Data`. Cambridge, 1998G.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"]importnumpyasnpfrompandasimportget_dummies,MultiIndexfromscipy.specialimportgammaln,digamma,polygamma,loggammafromscipyimportstats,specialfromscipy.statsimportnbinomfromstatsmodels.compat.pandasimportAppenderimportstatsmodels.tools.toolsastoolsfromstatsmodels.toolsimportdataasdata_toolsfromstatsmodels.tools.decoratorsimportcache_readonlyfromstatsmodels.tools.sm_exceptionsimport(PerfectSeparationError,SpecificationWarning)fromstatsmodels.tools.numdiffimportapprox_fprime_csimportstatsmodels.base.modelasbasefromstatsmodels.base.dataimporthandle_data# for mnlogitimportstatsmodels.regression.linear_modelaslmimportstatsmodels.base.wrapperaswrapfromstatsmodels.base.l1_slsqpimportfit_l1_slsqpfromstatsmodels.distributionsimportgenpoisson_ptry:importcvxopt# noqa:F401have_cvxopt=TrueexceptImportError:have_cvxopt=Falseimportwarnings#TODO: When we eventually get user-settable precision, we need to change# thisFLOAT_EPS=np.finfo(float).eps#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) modelThis is a preliminary estimator used as start_params.Returns-------params : ndarray parameter estimate based one one-step moment matching"""# helper for MNLogit (will be generally useful later)def_numpy_to_dummies(endog):ifendog.dtype.kindin['S','O']:endog_dummies,ynames=tools.categorical(endog,drop=True,dictnames=True)elifendog.ndim==2:endog_dummies=endogynames=range(endog.shape[1])else:endog_dummies,ynames=tools.categorical(endog,drop=True,dictnames=True)returnendog_dummies,ynamesdef_pandas_to_dummies(endog):ifendog.ndim==2:ifendog.shape[1]==1:yname=endog.columns[0]endog_dummies=get_dummies(endog.iloc[:,0])else:# seriesyname='y'endog_dummies=endogelse:yname=endog.nameendog_dummies=get_dummies(endog)ynames=endog_dummies.columns.tolist()returnendog_dummies,ynames,ynamedef_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 """ifmethodnotin['l1','l1_cvxopt_cp']:raiseValueError('`method` = {method} is not supported, use either ''"l1" or "l1_cvxopt_cp"'.format(method=method))#### Private Model Classes ####classDiscreteModel(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,**kwargs):super(DiscreteModel,self).__init__(endog,exog,**kwargs)self.raise_on_perfect_prediction=Truedefinitialize(self):""" Initialize is called by statsmodels.model.LikelihoodModel.__init__ and should contain any preprocessing that needs to be done for a model. """# assumes constantrank=np.linalg.matrix_rank(self.exog)self.df_model=float(rank-1)self.df_resid=float(self.exog.shape[0]-rank)defcdf(self,X):""" The cumulative distribution function of the model. """raiseNotImplementedErrordefpdf(self,X):""" The probability density (mass) function of the model. """raiseNotImplementedErrordef_check_perfect_pred(self,params,*args):endog=self.endogfittedvalues=self.cdf(np.dot(self.exog,params[:self.exog.shape[1]]))if(self.raise_on_perfect_predictionandnp.allclose(fittedvalues-endog,0)):msg="Perfect separation detected, results not available"raisePerfectSeparationError(msg)@Appender(base.LikelihoodModel.fit.__doc__)deffit(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 """ifcallbackisNone:callback=self._check_perfect_predelse:pass# TODO: make a function factory to have multiple call-backsmlefit=super(DiscreteModel,self).fit(start_params=start_params,method=method,maxiter=maxiter,full_output=full_output,disp=disp,callback=callback,**kwargs)returnmlefit# It is up to subclasses to wrap resultsdeffit_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 methodcov_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 timesalpha=np.array(alpha)assertalpha.min()>=0try:kwargs['alpha']=alphaexceptTypeError:kwargs=dict(alpha=alpha)kwargs['alpha_rescaled']=kwargs['alpha']/float(self.endog.shape[0])kwargs['trim_mode']=trim_modekwargs['size_trim_tol']=size_trim_tolkwargs['auto_trim_tol']=auto_trim_tolkwargs['qc_tol']=qc_tolkwargs['qc_verbose']=qc_verbose### Define default keyword arguments to be passed to super(...).fit()ifmaxiter=='defined_by_method':ifmethod=='l1':maxiter=1000elifmethod=='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}ifhave_cvxoptandmethod=='l1_cvxopt_cp':fromstatsmodels.base.l1_cvxoptimportfit_l1_cvxopt_cpextra_fit_funcs['l1_cvxopt_cp']=fit_l1_cvxopt_cpelifmethod.lower()=='l1_cvxopt_cp':raiseValueError("Cannot use l1_cvxopt_cp as cvxopt ""was not found (install it, or use method='l1' instead)")ifcallbackisNone:callback=self._check_perfect_predelse:pass# make a function factory to have multiple call-backsmlefit=super(DiscreteModel,self).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)returnmlefit# up to subclasses to wrap resultsdefcov_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()ifnnz_params>0:H_restricted=H[nz_idx[:,None],nz_idx]# Covariance estimate for the nonzero paramsH_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_invreturncov_paramsdefpredict(self,params,exog=None,linear=False):""" Predict response variable of a model given exogenous variables. """raiseNotImplementedErrordef_derivative_exog(self,params,exog=None,dummy_idx=None,count_idx=None):""" This should implement the derivative of the non-linear function """raiseNotImplementedErrordef_derivative_exog_helper(self,margeff,params,exog,dummy_idx,count_idx,transform):""" Helper for _derivative_exog to wrap results appropriately """from.discrete_marginsimport_get_count_effects,_get_dummy_effectsifcount_idxisnotNone:margeff=_get_count_effects(margeff,exog,count_idx,transform,self,params)ifdummy_idxisnotNone:margeff=_get_dummy_effects(margeff,exog,dummy_idx,transform,self,params)returnmargeff

[docs]classBinaryModel(DiscreteModel):def__init__(self,endog,exog,**kwargs):super(BinaryModel,self).__init__(endog,exog,**kwargs)if(notissubclass(self.__class__,MultinomialModel)andnotnp.all((self.endog>=0)&(self.endog<=1))):raiseValueError("endog must be in the unit interval.")

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. """ifexogisNone:exog=self.exogdF=self.pdf(np.dot(exog,params))[:,None]*exogif'ey'intransform:dF/=self.predict(params,exog)[:,None]returndFdef_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. """# Note: this form should be appropriate for# group 1 probit, logit, logistic, cloglog, heckprob, xtprobitifexogisNone:exog=self.exogmargeff=np.dot(self.pdf(np.dot(exog,params))[:,None],params[None,:])if'ex'intransform:margeff*=exogif'ey'intransform:margeff/=self.predict(params,exog)[:,None]returnself._derivative_exog_helper(margeff,params,exog,dummy_idx,count_idx,transform)

[docs]classMultinomialModel(BinaryModel):def_handle_data(self,endog,exog,missing,hasconst,**kwargs):ifdata_tools._is_using_ndarray_type(endog,None):endog_dummies,ynames=_numpy_to_dummies(endog)yname='y'elifdata_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'ifnotisinstance(ynames,dict):ynames=dict(zip(range(endog_dummies.shape[1]),ynames))self._ynames_map=ynamesdata=handle_data(endog_dummies,exog,missing,hasconst,**kwargs)data.ynames=yname# overwrite this to single endog namedata.orig_endog=endogself.wendog=data.endog# repeating from upstream...forkeyinkwargs:ifkeyin['design_info','formula']:# leave attached to datacontinuetry:setattr(self,key,data.__dict__.pop(key))exceptKeyError:passreturndata

[docs]definitialize(self):""" Preprocesses the data for MNLogit. """super(MultinomialModel,self).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 idxself.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)

[docs]defpredict(self,params,exog=None,linear=False):""" 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. linear : bool, optional If True, returns the linear predictor dot(exog,params). Else, returns the value of the cdf at the linear predictor. Notes ----- Column 0 is the base case, the rest conform to the rows of params shifted up one for the base case. """ifexogisNone:# do here to accommodate user-given exogexog=self.exogifexog.ndim==1:exog=exog[None]pred=super(MultinomialModel,self).predict(params,exog,linear)iflinear:pred=np.column_stack((np.zeros(len(exog)),pred))returnpred

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. """ifexogisNone:exog=self.exogifparams.ndim==1:# will get flatted from approx_fprimeparams=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 levelF0=-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 indexother_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'intransform:dFdX/=self.predict(params,exog)[:,:,None]returndFdXdef_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 choicesK=int(self.K)# number of variables# Note: this form should be appropriate for# group 1 probit, logit, logistic, cloglog, heckprob, xtprobitifexogisNone:exog=self.exogifparams.ndim==1:# will get flatted from approx_fprimeparams=params.reshape(K,J-1,order='F')zeroparams=np.c_[np.zeros(K),params]# add base incdf=self.cdf(np.dot(exog,params))# TODO: meaningful interpretation for `iterm`?iterm=np.array([cdf[:,[i]]*zeroparams[:,i]foriinrange(int(J))]).sum(0)margeff=np.array([cdf[:,[j]]*(zeroparams[:,j]-iterm)forjinrange(J)])# swap the axes to make sure margeff are in order nobs, K, Jmargeff=np.transpose(margeff,(1,2,0))if'ex'intransform:margeff*=exogif'ey'intransform:margeff/=self.predict(params,exog)[:,None,:]margeff=self._derivative_exog_helper(margeff,params,exog,dummy_idx,count_idx,transform)returnmargeff.reshape(len(exog),-1,order='F')

classCountModel(DiscreteModel):def__init__(self,endog,exog,offset=None,exposure=None,missing='none',**kwargs):super(CountModel,self).__init__(endog,exog,missing=missing,offset=offset,exposure=exposure,**kwargs)ifexposureisnotNone:self.exposure=np.log(self.exposure)self._check_inputs(self.offset,self.exposure,self.endog)ifoffsetisNone:delattr(self,'offset')ifexposureisNone:delattr(self,'exposure')# promote dtype to float64 if neededdt=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):ifoffsetisnotNoneandoffset.shape[0]!=endog.shape[0]:raiseValueError("offset is not the same length as endog")ifexposureisnotNoneandexposure.shape[0]!=endog.shape[0]:raiseValueError("exposure is not the same length as endog")def_get_init_kwds(self):# this is a temporary fixup because exposure has been transformed# see #1609kwds=super(CountModel,self)._get_init_kwds()if'exposure'inkwdsandkwds['exposure']isnotNone:kwds['exposure']=np.exp(kwds['exposure'])returnkwdsdefpredict(self,params,exog=None,exposure=None,offset=None,linear=False):""" 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. linear : bool If True, returns the linear predicted values. If False, returns the exponential of the linear predicted value. Notes ----- If exposure is specified, then it will be logged by the method. The user does not need to log it first. """# the following is copied from GLM predict (without family/link check)# Use fit offset if appropriateifoffsetisNoneandexogisNoneandhasattr(self,'offset'):offset=self.offsetelifoffsetisNone:offset=0.# Use fit exposure if appropriateifexposureisNoneandexogisNoneandhasattr(self,'exposure'):# Already loggedexposure=self.exposureelifexposureisNone:exposure=0.else:exposure=np.log(exposure)ifexogisNone:exog=self.exogfitted=np.dot(exog,params[:exog.shape[1]])linpred=fitted+exposure+offsetifnotlinear:returnnp.exp(linpred)# not cdfelse:returnlinpreddef_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. """ifexogisNone:exog=self.exog#NOTE: this handles offset and exposuredF=self.predict(params,exog)[:,None]*exogif'ey'intransform:dF/=self.predict(params,exog)[:,None]returndFdef_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, zinbifexogisNone:exog=self.exogk_extra=getattr(self,'k_extra',0)params_exog=paramsifk_extra==0elseparams[:-k_extra]margeff=self.predict(params,exog)[:,None]*params_exog[None,:]if'ex'intransform:margeff*=exogif'ey'intransform:margeff/=self.predict(params,exog)[:,None]returnself._derivative_exog_helper(margeff,params,exog,dummy_idx,count_idx,transform)@Appender(DiscreteModel.fit.__doc__)deffit(self,start_params=None,method='newton',maxiter=35,full_output=1,disp=1,callback=None,**kwargs):cntfit=super(CountModel,self).fit(start_params=start_params,method=method,maxiter=maxiter,full_output=full_output,disp=disp,callback=callback,**kwargs)discretefit=CountResults(self,cntfit)returnCountResultsWrapper(discretefit)@Appender(DiscreteModel.fit_regularized.__doc__)deffit_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=L1CountResults(self,cntfit)returnL1CountResultsWrapper(discretefit)classOrderedModel(DiscreteModel):pass# Public Model ClassesclassPoisson(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}@propertydeffamily(self):fromstatsmodels.genmodimportfamiliesreturnfamilies.Poisson()defcdf(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.endogreturnstats.poisson.cdf(y,np.exp(X))defpdf(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.endogreturnnp.exp(stats.poisson.logpmf(y,np.exp(X)))defloglike(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+exposureendog=self.endogreturnnp.sum(-np.exp(XB)+endog*XB-gammaln(endog+1))defloglikeobs(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+exposureendog=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)]returnparams@Appender(DiscreteModel.fit.__doc__)deffit(self,start_params=None,method='newton',maxiter=35,full_output=1,disp=1,callback=None,**kwargs):ifstart_paramsisNoneandself.data.const_idxisnotNone:# 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]cntfit=super(CountModel,self).fit(start_params=start_params,method=method,maxiter=maxiter,full_output=full_output,disp=disp,callback=callback,**kwargs)if'cov_type'inkwargs:cov_kwds=kwargs.get('cov_kwds',{})kwds={'cov_type':kwargs['cov_type'],'cov_kwds':cov_kwds}else:kwds={}discretefit=PoissonResults(self,cntfit,**kwds)returnPoissonResultsWrapper(discretefit)@Appender(DiscreteModel.fit_regularized.__doc__)deffit_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)returnL1PoissonResultsWrapper(discretefit)deffit_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 importsfrompatsyimportDesignInfofromstatsmodels.base._constraintsimport(fit_constrained,LinearConstraints)# same pattern as in base.LikelihoodModel.t_testlc=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 transformationparams,cov,res_constr=fit_constrained(self,R,q,start_params=start_params,fit_kwds=fit_kwds)#create dummy results Instance, TODO: wire up properlyres=self.fit(maxiter=0,method='nm',disp=0,warn_convergence=False)# we get a wrapper backres.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=paramsres._results.cov_params_default=covcov_type=fit_kwds.get('cov_type','nonrobust')ifcov_type!='nonrobust':res._results.normalized_cov_params=cov# assume scale=1else:res._results.normalized_cov_params=Nonek_constr=len(q)res._results.df_resid+=k_constrres._results.df_model-=k_constrres._results.constraints=LinearConstraints.from_patsy(lc)res._results.k_constr=k_constrres._results.results_constrained=res_constrreturnresdefscore(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.exogL=np.exp(np.dot(X,params)+offset+exposure)returnnp.dot(self.endog-L,X)defscore_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.exogL=np.exp(np.dot(X,params)+offset+exposure)return(self.endog-L)[:,None]*Xdefscore_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.exogL=np.exp(np.dot(X,params)+offset+exposure)return(self.endog-L)defhessian(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.exogL=np.exp(np.dot(X,params)+exposure+offset)return-np.dot(L*X.T,X)defhessian_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.exogL=np.exp(np.dot(X,params)+exposure+offset)returnLclassGeneralizedPoisson(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}def__init__(self,endog,exog,p=1,offset=None,exposure=None,missing='none',**kwargs):super(GeneralizedPoisson,self).__init__(endog,exog,offset=offset,exposure=exposure,missing=missing,**kwargs)self.parameterization=p-1self.exog_names.append('alpha')self.k_extra=1self._transparams=Falsedef_get_init_kwds(self):kwds=super(GeneralizedPoisson,self)._get_init_kwds()kwds['p']=self.parameterization+1returnkwdsdefloglike(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] """returnnp.sum(self.loglikeobs(params))defloglikeobs(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` """ifself._transparams:alpha=np.exp(params[-1])else:alpha=params[-1]params=params[:-1]p=self.parameterizationendog=self.endogmu=self.predict(params)mu_p=np.power(mu,p)a1=1+alpha*mu_pa2=mu+(a1-1)*endogreturn(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-mua=self._estimate_dispersion(mu,resid,df_resid=resid.shape[0]-1)params.append(a)returnnp.array(params)def_estimate_dispersion(self,mu,resid,df_resid=None):q=self.parameterizationifdf_residisNone:df_resid=resid.shape[0]a=((np.abs(resid)/np.sqrt(mu)-1)*mu**(-q)).sum()/df_residreturna@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__)deffit(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,**kwargs):ifuse_transparamsandmethodnotin['newton','ncg']:self._transparams=Trueelse:ifuse_transparams:warnings.warn('Parameter "use_transparams" is ignored',RuntimeWarning)self._transparams=Falseifstart_paramsisNone:offset=getattr(self,"offset",0)+getattr(self,"exposure",0)ifnp.size(offset)==1andoffset==0:offset=Noneoptim_kwds_prelim={'disp':0,'skip_hessian':True,'warn_convergence':False}optim_kwds_prelim.update(kwargs.get('optim_kwds_prelim',{}))mod_poi=Poisson(self.endog,self.exog,offset=offset)res_poi=mod_poi.fit(**optim_kwds_prelim)start_params=res_poi.paramsa=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))ifcallbackisNone:# work around perfect separation callback #3895callback=lambda*x:xmlefit=super(GeneralizedPoisson,self).fit(start_params=start_params,maxiter=maxiter,method=method,disp=disp,full_output=full_output,callback=callback,**kwargs)ifuse_transparamsandmethodnotin["newton","ncg"]:self._transparams=Falsemlefit._results.params[-1]=np.exp(mlefit._results.params[-1])gpfit=GeneralizedPoissonResults(self,mlefit._results)result=GeneralizedPoissonResultsWrapper(gpfit)ifcov_kwdsisNone:cov_kwds={}result._get_robustcov_results(cov_type=cov_type,use_self=True,use_t=use_t,**cov_kwds)returnresult@Appender(DiscreteModel.fit_regularized.__doc__)deffit_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)ifnp.size(alpha)==1andalpha!=0:k_params=self.exog.shape[1]+self.k_extraalpha=alpha*np.ones(k_params)alpha[-1]=0alpha_p=alpha[:-1]if(self.k_extraandnp.size(alpha)>1)elsealphaself._transparams=Falseifstart_paramsisNone:offset=getattr(self,"offset",0)+getattr(self,"exposure",0)ifnp.size(offset)==1andoffset==0:offset=Nonemod_poi=Poisson(self.endog,self.exog,offset=offset)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).paramsstart_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)returnL1GeneralizedPoissonResultsWrapper(discretefit)defscore_obs(self,params):ifself._transparams:alpha=np.exp(params[-1])else:alpha=params[-1]params=params[:-1]p=self.parameterizationexog=self.exogy=self.endog[:,None]mu=self.predict(params)[:,None]mu_p=np.power(mu,p)a1=1+alpha*mu_pa2=mu+alpha*mu_p*ya3=alpha*p*mu**(p-1)a4=a3*ydmudb=mu*exogdalpha=(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)returnnp.concatenate((dparams,np.atleast_2d(dalpha)),axis=1)defscore(self,params):score=np.sum(self.score_obs(params),axis=0)ifself._transparams:score[-1]==score[-1]**2returnscoreelse:returnscoredef_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`. """ifself._transparams:alpha=np.exp(params[-1])else:alpha=params[-1]params=params[:-1]p=self.parameterizationexog=self.exogy=self.endog[:,None]mu=self.predict(params)[:,None]mu_p=np.power(mu,p)a1=1+alpha*mu_pa2=mu+alpha*mu_p*ydp=np.sum((np.log(mu)*((a2-mu)*((y-1)/a2-2/a1)+(a1-1)*a2/a1**2)))returndpdefhessian(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` """ifself._transparams:alpha=np.exp(params[-1])else:alpha=params[-1]params=params[:-1]p=self.parameterizationexog=self.exogy=self.endog[:,None]mu=self.predict(params)[:,None]mu_p=np.power(mu,p)a1=1+alpha*mu_pa2=mu+alpha*mu_p*ya3=alpha*p*mu**(p-1)a4=a3*ya5=p*mu**(p-1)dmudb=mu*exog# for dl/dparams dparamsdim=exog.shape[1]hess_arr=np.empty((dim+1,dim+1))foriinrange(dim):forjinrange(i+1):hess_arr[i,j]=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)tri_idx=np.triu_indices(dim,k=1)hess_arr[tri_idx]=hess_arr.T[tri_idx]# for dl/dparams dalphadldpda=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]=dldpdahess_arr[:-1,-1]=dldpda# for dl/dalpha dalphadldada=mu_p**2*(3*y/a1**2-(y/a2)**2.*(y-1)-2*a2/a1**3)hess_arr[-1,-1]=dldada.sum()returnhess_arrdefpredict(self,params,exog=None,exposure=None,offset=None,which='mean'):""" Predict response variable of a count model given exogenous variables. Notes ----- If exposure is specified, then it will be logged by the method. The user does not need to log it first. """ifexogisNone:exog=self.exogifexposureisNone:exposure=getattr(self,'exposure',0)elifexposure!=0:exposure=np.log(exposure)ifoffsetisNone:offset=getattr(self,'offset',0)fitted=np.dot(exog,params[:exog.shape[1]])linpred=fitted+exposure+offsetifwhich=='mean':returnnp.exp(linpred)elifwhich=='linear':returnlinpredelifwhich=='prob':counts=np.atleast_2d(np.arange(0,np.max(self.endog)+1))mu=self.predict(params,exog=exog,exposure=exposure,offset=offset)[:,None]returngenpoisson_p.pmf(counts,mu,params[-1],self.parameterization+1)else:raiseValueError('keyword \'which\' not recognized')

[docs]defloglikeobs(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-1X=self.exogreturnnp.log(self.cdf(q*np.dot(X,params)))

classProbit(BinaryModel):__doc__=""" Probit 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':base._missing_param_doc}defcdf(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 """returnstats.norm._cdf(X)defpdf(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)returnstats.norm._pdf(X)defloglike(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-1X=self.exogreturnnp.sum(np.log(np.clip(self.cdf(q*np.dot(X,params)),FLOAT_EPS,1)))defloglikeobs(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-1X=self.exogreturnnp.log(np.clip(self.cdf(q*np.dot(X,params)),FLOAT_EPS,1))defscore(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.endogX=self.exogXB=np.dot(X,params)q=2*y-1# clip to get rid of invalid divide complaintL=q*self.pdf(q*XB)/np.clip(self.cdf(q*XB),FLOAT_EPS,1-FLOAT_EPS)returnnp.dot(L,X)defscore_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.endogX=self.exogXB=np.dot(X,params)q=2*y-1# clip to get rid of invalid divide complaintL=q*self.pdf(q*XB)/np.clip(self.cdf(q*XB),FLOAT_EPS,1-FLOAT_EPS)returnL[:,None]*Xdefhessian(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.exogXB=np.dot(X,params)q=2*self.endog-1L=q*self.pdf(q*XB)/self.cdf(q*XB)returnnp.dot(-L*(L+XB)*X.T,X)@Appender(DiscreteModel.fit.__doc__)deffit(self,start_params=None,method='newton',maxiter=35,full_output=1,disp=1,callback=None,**kwargs):bnryfit=super(Probit,self).fit(start_params=start_params,method=method,maxiter=maxiter,full_output=full_output,disp=disp,callback=callback,**kwargs)discretefit=ProbitResults(self,bnryfit)returnBinaryResultsWrapper(discretefit)

[docs]classMNLogit(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)s 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. """%{'extra_params':base._missing_param_doc}def__init__(self,endog,exog,**kwargs):super(MNLogit,self).__init__(endog,exog,**kwargs)# Override cov_names since multivariate modelyname=self.endog_namesynames=self._ynames_mapynames=MultinomialResults._maybe_convert_ynames_int(ynames)# use range below to ensure sortednessynames=[ynames[key]forkeyinrange(int(self.J))]idx=MultiIndex.from_product((ynames[1:],self.data.xnames),names=(yname,None))self.data.cov_names=idx

[docs]defscore(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 reshapedreturnnp.dot(firstterm.T,self.exog).flatten()

[docs]defloglike_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()returnloglike_value,score_array

[docs]defscore_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 reshapedreturn(firstterm[:,:,None]*self.exog[:,None,:]).reshape(self.exog.shape[0],-1)

[docs]defhessian(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.exogpr=self.cdf(np.dot(X,params))partials=[]J=self.JK=self.Kforiinrange(J-1):forjinrange(J-1):# this loop assumes we drop the first col.ifi==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 upH=np.transpose(H.reshape(J-1,J-1,K,K),(0,2,1,3)).reshape((J-1)*K,(J-1)*K)returnH

[docs]deffit(self,start_params=None,method='bfgs',maxiter=35,full_output=1,disp=1,callback=None,cov_type='nonrobust',cov_kwds=None,use_t=None,**kwargs):# Note: do not let super handle robust covariance because it has# transformed paramsself._transparams=False# always define attributeifself.loglike_method.startswith('nb')andmethodnotin['newton','ncg']:self._transparams=True# in case same Model instance is refitelifself.loglike_method.startswith('nb'):# method is newton/ncgself._transparams=False# because we need to step in alpha spaceifstart_paramsisNone:# Use poisson fit as first guess.#TODO, Warning: this assumes exposure is loggedoffset=getattr(self,"offset",0)+getattr(self,"exposure",0)ifnp.size(offset)==1andoffset==0:offset=Noneoptim_kwds_prelim={'disp':0,'skip_hessian':True,'warn_convergence':False}optim_kwds_prelim.update(kwargs.get('optim_kwds_prelim',{}))mod_poi=Poisson(self.endog,self.exog,offset=offset)res_poi=mod_poi.fit(**optim_kwds_prelim)start_params=res_poi.paramsifself.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:ifself._transparamsisTrue:# transform user provided start_params dispersion, see #3918start_params=np.array(start_params,copy=True)start_params[-1]=np.log(start_params[-1])ifcallbackisNone:# work around perfect separation callback #3895callback=lambda*x:xmlefit=super(NegativeBinomial,self).fit(start_params=start_params,maxiter=maxiter,method=method,disp=disp,full_output=full_output,callback=callback,**kwargs)# TODO: Fix NBin _check_perfect_predifself.loglike_method.startswith('nb'):# mlefit is a wrapped counts resultsself._transparams=False# do not need to transform anymore now# change from lnalpha to alphaifmethodnotin["newton","ncg"]:mlefit._results.params[-1]=np.exp(mlefit._results.params[-1])nbinfit=NegativeBinomialResults(self,mlefit._results)result=NegativeBinomialResultsWrapper(nbinfit)else:result=mlefitifcov_kwdsisNone:cov_kwds={}#TODO: make this unnecessary ?result._get_robustcov_results(cov_type=cov_type,use_self=True,use_t=use_t,**cov_kwds)returnresult

[docs]deffit_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)ifself.loglike_method.startswith('nb')and(np.size(alpha)==1andalpha!=0):# do not penalize alpha if alpha is scalark_params=self.exog.shape[1]+self.k_extraalpha=alpha*np.ones(k_params)alpha[-1]=0# alpha for regularized poisson to get starting valuesalpha_p=alpha[:-1]if(self.k_extraandnp.size(alpha)>1)elsealphaself._transparams=Falseifstart_paramsisNone:# Use poisson fit as first guess.#TODO, Warning: this assumes exposure is loggedoffset=getattr(self,"offset",0)+getattr(self,"exposure",0)ifnp.size(offset)==1andoffset==0:offset=Nonemod_poi=Poisson(self.endog,self.exog,offset=offset)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).paramsifself.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)returnL1NegativeBinomialResultsWrapper(discretefit)

[docs]defpredict(self,params,exog=None,exposure=None,offset=None,which='mean'):""" 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, 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. linear : bool, optional If True, returns the linear predictor dot(exog,params). Else, returns the value of the cdf at the linear predictor. offset : array_like, optional Offset is added to the linear prediction with coefficient equal to 1. exposure : array_like, optional Log(exposure) is added to the linear prediction with coefficient equal to 1. which : 'mean', 'linear', 'prob', optional. 'mean' returns the exp of linear predictor exp(dot(exog,params)). 'linear' returns the linear predictor dot(exog,params). 'prob' return probabilities for counts from 0 to max(endog). Default is 'mean'. Notes ----- """ifexogisNone:exog=self.exogifexposureisNone:exposure=getattr(self,'exposure',0)elifexposure!=0:exposure=np.log(exposure)ifoffsetisNone:offset=getattr(self,'offset',0)fitted=np.dot(exog,params[:exog.shape[1]])linpred=fitted+exposure+offsetifwhich=='mean':returnnp.exp(linpred)elifwhich=='linear':returnlinpredelifwhich=='prob':counts=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)returnnbinom.pmf(counts,size[:,None],prob[:,None])else:raiseValueError('keyword "which" = %s not recognized'%which)

### Results Class ###classDiscreteResults(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=modelself.df_model=model.df_modelself.df_resid=model.df_residself._cache={}self.nobs=model.exog.shape[0]self.__dict__.update(mlefit.__dict__)ifnothasattr(self,'cov_type'):# do this only if super, i.e. mlefit did not already add cov_type# robust covarianceifuse_tisnotNone:self.use_t=use_tifcov_type=='nonrobust':self.cov_type='nonrobust'self.cov_kwds={'description':'Standard Errors assume that the '+'covariance matrix of the errors is correctly '+'specified.'}else:ifcov_kwdsisNone:cov_kwds={}fromstatsmodels.base.covtypeimportget_robustcov_resultsget_robustcov_results(self,cov_type=cov_type,use_self=True,**cov_kwds)def__getstate__(self):# remove unpicklable methodsmle_settings=getattr(self,'mle_settings',None)ifmle_settingsisnotNone:if'callback'inmle_settings:mle_settings['callback']=Noneif'cov_params_func'inmle_settings:mle_settings['cov_params_func']=Nonereturnself.__dict__@cache_readonlydefprsquared(self):""" McFadden's pseudo-R-squared. `1 - (llf / llnull)` """return1-self.llf/self.llnull@cache_readonlydefllr(self):""" Likelihood ratio chi-squared statistic; `-2*(llnull - llf)` """return-2*(self.llnull-self.llf)@cache_readonlydefllr_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`. """returnstats.distributions.chi2.sf(self.llr,self.df_model)defset_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)ifhasattr(self,'res_null'):delself.res_nullifllnullisnotNone:self._cache['llnull']=llnullself._attach_nullmodel=attach_resultsself._optim_kwds_null=kwargs@cache_readonlydefllnull(self):""" Value of the constant-only loglikelihood """model=self.modelkwds=model._get_init_kwds().copy()forkeyingetattr(model,'_null_drop_keys',[]):delkwds[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_smokeoptim_kwds=getattr(self,'_optim_kwds_null',{}).copy()if'start_params'inoptim_kwds:# user providedsp_null=optim_kwds.pop('start_params')elifhasattr(model,'_get_start_params_null'):# get moment estimates if availablesp_null=model._get_start_params_null()else:sp_null=Noneopt_kwds=dict(method='bfgs',warn_convergence=False,maxiter=10000,disp=0)opt_kwds.update(optim_kwds)ifoptim_kwds:res_null=mod_null.fit(start_params=sp_null,**opt_kwds)else:# this should be a reasonably method case across versionsres_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)ifgetattr(self,'_attach_nullmodel',False)isnotFalse:self.res_null=res_nullreturnres_null.llf@cache_readonlydeffittedvalues(self):""" Linear predictor XB. """returnnp.dot(self.model.exog,self.params[:self.model.exog.shape[1]])@cache_readonlydefresid_response(self):""" Respnose residuals. The response residuals are defined as `endog - fittedvalues` """returnself.model.endog-self.predict()@cache_readonlydefaic(self):""" Akaike information criterion. `-2*(llf - p)` where `p` is the number of regressors including the intercept. """return-2*(self.llf-(self.df_model+1))@cache_readonlydefbic(self):""" Bayesian information criterion. `-2*llf + ln(nobs)*p` where `p` is the number of regressors including the intercept. """return-2*self.llf+np.log(self.nobs)*(self.df_model+1)def_get_endog_name(self,yname,yname_list):ifynameisNone:yname=self.model.endog_namesifyname_listisNone:yname_list=self.model.endog_namesreturnyname,yname_listdefget_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. """fromstatsmodels.discrete.discrete_marginsimportDiscreteMarginsreturnDiscreteMargins(self,(at,method,atexog,dummy,count))defsummary(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:',['MLE']),('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])]ifhasattr(self,'cov_type'):top_left.append(('Covariance Type:',[self.cov_type]))iftitleisNone:title=self.model.__class__.__name__+' '+"Regression Results"# boiler platefromstatsmodels.iolib.summaryimportSummarysmry=Summary()yname,yname_list=self._get_endog_name(yname,yname_list)# for top of tablesmry.add_table_2cols(self,gleft=top_left,gright=top_right,yname=yname,xname=xname,title=title)# for parameters, etcsmry.add_table_params(self,yname=yname_list,xname=xname,alpha=alpha,use_t=self.use_t)ifhasattr(self,'constraints'):smry.add_extra_txt(['Model has been estimated subject to linear ''equality constraints.'])returnsmrydefsummary2(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. """fromstatsmodels.iolibimportsummary2smry=summary2.Summary()smry.add_base(results=self,alpha=alpha,float_format=float_format,xname=xname,yname=yname,title=title)ifhasattr(self,'constraints'):smry.add_text('Model has been estimated subject to linear ''equality constraints.')returnsmryclassCountResults(DiscreteResults):__doc__=_discrete_results_docs%{"one_line_description":"A results class for count data","extra_attr":""}@cache_readonlydefresid(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. """returnself.model.endog-self.predict()

classGeneralizedPoissonResults(NegativeBinomialResults):__doc__=_discrete_results_docs%{"one_line_description":"A results class for Generalized Poisson","extra_attr":""}@cache_readonlydef_dispersion_factor(self):p=getattr(self.model,'parameterization',0)mu=self.predict()return(1+self.params[-1]*mu**p)**2classL1CountResults(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(L1CountResults,self).__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_modelk_extra=getattr(self.model,'k_extra',0)self.df_model=self.nnz_params-1-k_extraself.df_resid=float(self.model.endog.shape[0]-self.nnz_params)+k_extraclassPoissonResults(CountResults):defpredict_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. """ifnisnotNone: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,linear=False)[:,None]# uses broadcastingreturnstats.poisson.pmf(counts,mu)@propertydefresid_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 residualsp=self.predict()# fittedvalues is still linearreturn(self.model.endog-p)/np.sqrt(p)classL1PoissonResults(L1CountResults,PoissonResults):passclassL1NegativeBinomialResults(L1CountResults,NegativeBinomialResults):passclassL1GeneralizedPoissonResults(L1CountResults,GeneralizedPoissonResults):passclassOrderedResults(DiscreteResults):__doc__=_discrete_results_docs%{"one_line_description":"A results class for ordered discrete data.","extra_attr":""}pass

[docs]classBinaryResults(DiscreteResults):__doc__=_discrete_results_docs%{"one_line_description":"A results class for binary data","extra_attr":""}defpred_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.modelactual=model.endogpred=np.array(self.predict()>threshold,dtype=float)bins=np.array([0,0.5,1])returnnp.histogram2d(actual,pred,bins=bins)[0]@Appender(DiscreteResults.summary.__doc__)defsummary(self,yname=None,xname=None,title=None,alpha=.05,yname_list=None):smry=super(BinaryResults,self).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/notesetext=[]ifpredclose_sum==len(fittedvalues):# TODO: nobs?wstr="Complete Separation: The results show that there is"wstr+="complete separation.\n"wstr+="In this case the Maximum Likelihood Estimator does "wstr+="not exist and the parameters\n"wstr+="are not identified."etext.append(wstr)elifpredclose_frac>0.1:# TODO: get better diagnosiswstr="Possibly complete quasi-separation: A fraction "wstr+="%4.2f of observations can be\n"%predclose_fracwstr+="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)ifetext:smry.add_extra_txt(etext)returnsmry@cache_readonlydefresid_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.modelendog=self.model.endog#exog = model.exog# M = # of individuals that share a covariate pattern# so M[i] = 2 for i = two share a covariate patternM=1p=self.predict()#Y_0 = np.where(exog == 0)#Y_M = np.where(exog == M)#NOTE: Common covariate patterns are not yet handledres=-(1-endog)*np.sqrt(2*M*np.abs(np.log(1-p)))+ \
endog*np.sqrt(2*M*np.abs(np.log(p)))returnres@cache_readonlydefresid_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.modelendog=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=1p=self.predict()return(endog-M*p)/np.sqrt(M*p*(1-p))@cache_readonlydefresid_response(self):""" The response residuals Notes ----- Response residuals are defined to be .. math:: y - p where :math:`p=cdf(X\\beta)`. """returnself.model.endog-self.predict()

[docs]classLogitResults(BinaryResults):__doc__=_discrete_results_docs%{"one_line_description":"A results class for Logit Model","extra_attr":""}@cache_readonlydefresid_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 residualsreturnself.model.endog-self.predict()

classProbitResults(BinaryResults):__doc__=_discrete_results_docs%{"one_line_description":"A results class for Probit Model","extra_attr":""}@cache_readonlydefresid_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 residualsmodel=self.modelendog=model.endogXB=self.predict(linear=True)pdf=model.pdf(XB)cdf=model.cdf(XB)returnendog*pdf/cdf-(1-endog)*pdf/(1-cdf)classL1BinaryResults(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(L1BinaryResults,self).__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-1self.df_resid=float(self.model.endog.shape[0]-self.nnz_params)

[docs]classMultinomialResults(DiscreteResults):__doc__=_discrete_results_docs%{"one_line_description":"A results class for multinomial data","extra_attr":""}def__init__(self,model,mlefit):super(MultinomialResults,self).__init__(model,mlefit)self.J=model.Jself.K=model.K@staticmethoddef_maybe_convert_ynames_int(ynames):# see if they're integersissue_warning=Falsemsg=('endog contains values are that not int-like. Uses string ''representation of value. Use integer-valued endog to ''suppress this warning.')foriinynames:try:ifynames[i]%1==0:ynames[i]=str(int(ynames[i]))else:issue_warning=Trueynames[i]=str(ynames[i])exceptTypeError:ynames[i]=str(ynames[i])ifissue_warning:importwarningswarnings.warn(msg,SpecificationWarning)returnynamesdef_get_endog_name(self,yname,yname_list,all=False):""" If all is False, the first variable name is dropped """model=self.modelifynameisNone:yname=model.endog_namesifyname_listisNone:ynames=model._ynames_mapynames=self._maybe_convert_ynames_int(ynames)# use range below to ensure sortednessynames=[ynames[key]forkeyinrange(int(model.J))]ynames=['='.join([yname,name])fornameinynames]ifnotall:yname_list=ynames[1:]# assumes first variable is droppedelse:yname_list=ynamesreturnyname,yname_list

[docs]defpred_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]))returnnp.histogram2d(self.model.endog,self.predict().argmax(1),bins=bins)[0]

@cache_readonlydefresid_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 onereturn(self.model.wendog.argmax(1)!=self.predict().argmax(1)).astype(float)

[docs]defsummary2(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 """fromstatsmodels.iolibimportsummary2smry=summary2.Summary()smry.add_dict(summary2.summary_model(self))# One data frame per value of endogeqn=self.params.shape[1]confint=self.conf_int(alpha)foriinrange(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 endoglevel_str=self.model.endog_names+' = '+str(i)coefs[level_str]=coefs.indexcoefs=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)returnsmry

classL1MultinomialResults(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(L1MultinomialResults,self).__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 constantsself.df_model=self.nnz_params-(self.model.J-1)self.df_resid=float(self.model.endog.shape[0]-self.nnz_params)#### Results Wrappers ####classOrderedResultsWrapper(lm.RegressionResultsWrapper):passwrap.populate_wrapper(OrderedResultsWrapper,OrderedResults)classCountResultsWrapper(lm.RegressionResultsWrapper):passwrap.populate_wrapper(CountResultsWrapper,CountResults)classNegativeBinomialResultsWrapper(lm.RegressionResultsWrapper):passwrap.populate_wrapper(NegativeBinomialResultsWrapper,NegativeBinomialResults)classGeneralizedPoissonResultsWrapper(lm.RegressionResultsWrapper):passwrap.populate_wrapper(GeneralizedPoissonResultsWrapper,GeneralizedPoissonResults)classPoissonResultsWrapper(lm.RegressionResultsWrapper):passwrap.populate_wrapper(PoissonResultsWrapper,PoissonResults)classL1CountResultsWrapper(lm.RegressionResultsWrapper):passclassL1PoissonResultsWrapper(lm.RegressionResultsWrapper):passwrap.populate_wrapper(L1PoissonResultsWrapper,L1PoissonResults)classL1NegativeBinomialResultsWrapper(lm.RegressionResultsWrapper):passwrap.populate_wrapper(L1NegativeBinomialResultsWrapper,L1NegativeBinomialResults)classL1GeneralizedPoissonResultsWrapper(lm.RegressionResultsWrapper):passwrap.populate_wrapper(L1GeneralizedPoissonResultsWrapper,L1GeneralizedPoissonResults)classBinaryResultsWrapper(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)classL1BinaryResultsWrapper(lm.RegressionResultsWrapper):passwrap.populate_wrapper(L1BinaryResultsWrapper,L1BinaryResults)classMultinomialResultsWrapper(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)classL1MultinomialResultsWrapper(lm.RegressionResultsWrapper):passwrap.populate_wrapper(L1MultinomialResultsWrapper,L1MultinomialResults)