Skip to content

Commit

Permalink
BUG NegativeBinomial add fit_regularized closes #1453 closes #1454
Browse files Browse the repository at this point in the history 
    adjust test to handle extra parameter
  • Loading branch information
@josef-pkt
josef-pkt committed Jul 22, 2014
1 parent 40bd7e0 commit 544400d
Show file tree
Hide file tree
Showing 2 changed files with 106 additions and 12 deletions.
41 changes: 41 additions & 0 deletions statsmodels/discrete/discrete_model.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2192,6 +2192,38 @@ def fit(self, start_params=None, method='bfgs', maxiter=35,
else:
return mlefit


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

self._transparams = False
if start_params == None:
# Use poisson fit as first guess.
start_params = Poisson(self.endog, self.exog).fit(disp=0).params
if self.loglike_method.startswith('nb'):
start_params = np.append(start_params, 0.1)

if self.loglike_method.startswith('nb') and (np.size(alpha) == 1 and
alpha != 0):
alpha = alpha * np.ones(len(start_params))
alpha[-1] = 0

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)
if method in ['l1', 'l1_cvxopt_cp']:
discretefit = L1NegativeBinomialAncillaryResults(self, cntfit)
else:
raise Exception(
"argument method == %s, which is not handled" % method)
#return discretefit
return L1NegativeBinomialAncillaryResultsWrapper(discretefit)


### Results Class ###

class DiscreteResults(base.LikelihoodModelResults):
Expand Down Expand Up @@ -2547,6 +2579,10 @@ def predict_prob(self, n=None, exog=None, exposure=None, offset=None,
class L1PoissonResults(L1CountResults, PoissonResults):
pass

class L1NegativeBinomialAncillaryResults(L1CountResults,
NegativeBinomialAncillaryResults):
pass

class OrderedResults(DiscreteResults):
__doc__ = _discrete_results_docs % {"one_line_description" : "A results class for ordered discrete data." , "extra_attr" : ""}
pass
Expand Down Expand Up @@ -2934,6 +2970,11 @@ class L1PoissonResultsWrapper(lm.RegressionResultsWrapper):
# _methods)
wrap.populate_wrapper(L1PoissonResultsWrapper, L1PoissonResults)

class L1NegativeBinomialAncillaryResultsWrapper(lm.RegressionResultsWrapper):
pass
wrap.populate_wrapper(L1NegativeBinomialAncillaryResultsWrapper,
L1NegativeBinomialAncillaryResults)

class BinaryResultsWrapper(lm.RegressionResultsWrapper):
_attrs = {"resid_dev" : "rows",
"resid_generalized" : "rows",
Expand Down
77 changes: 65 additions & 12 deletions statsmodels/discrete/tests/test_discrete.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -540,37 +540,45 @@ class CheckL1Compatability(object):
def test_params(self):
m = self.m
assert_almost_equal(
self.res_unreg.params, self.res_reg.params[:m], DECIMAL_4)
self.res_unreg.params[:m], self.res_reg.params[:m], DECIMAL_4)
# The last entry should be close to zero
assert_almost_equal(0, self.res_reg.params[m:], DECIMAL_4)
# handle extra parameter of NegativeBinomial
kvars = self.res_reg.model.exog.shape[1]
assert_almost_equal(0, self.res_reg.params[m:kvars], DECIMAL_4)

def test_cov_params(self):
m = self.m
# The restricted cov_params should be equal
assert_almost_equal(
self.res_unreg.cov_params(), self.res_reg.cov_params()[:m, :m],
self.res_unreg.cov_params()[:m, :m],
self.res_reg.cov_params()[:m, :m],
DECIMAL_1)

def test_df(self):
assert_equal(self.res_unreg.df_model, self.res_reg.df_model)
assert_equal(self.res_unreg.df_resid, self.res_reg.df_resid)
extra = len(self.res_reg.params) == (self.res_reg.model.exog.shape[1] + 1)
assert_equal(self.res_unreg.df_model + extra, self.res_reg.df_model)
assert_equal(self.res_unreg.df_resid - extra, self.res_reg.df_resid)

def test_t_test(self):
m = self.m
kvars = self.kvars
t_unreg = self.res_unreg.t_test(np.eye(m))
t_reg = self.res_reg.t_test(np.eye(kvars))
assert_almost_equal(t_unreg.effect, t_reg.effect[:m], DECIMAL_3)
assert_almost_equal(t_unreg.sd, t_reg.sd[:m], DECIMAL_3)
# handle extra parameter of NegativeBinomial
extra = len(self.res_reg.params) == (self.res_reg.model.exog.shape[1] + 1)
t_unreg = self.res_unreg.t_test(np.eye(len(self.res_unreg.params)))
t_reg = self.res_reg.t_test(np.eye(kvars + extra))
assert_almost_equal(t_unreg.effect[:m], t_reg.effect[:m], DECIMAL_3)
assert_almost_equal(t_unreg.sd[:m], t_reg.sd[:m], DECIMAL_3)
assert_almost_equal(np.nan, t_reg.sd[m])
assert_almost_equal(t_unreg.tvalue, t_reg.tvalue[:m], DECIMAL_3)
assert_almost_equal(t_unreg.tvalue[:m], t_reg.tvalue[:m], DECIMAL_3)
assert_almost_equal(np.nan, t_reg.tvalue[m])

def test_f_test(self):
m = self.m
kvars = self.kvars
f_unreg = self.res_unreg.f_test(np.eye(m))
f_reg = self.res_reg.f_test(np.eye(kvars)[:m])
# handle extra parameter of NegativeBinomial
extra = len(self.res_reg.params) == (self.res_reg.model.exog.shape[1] + 1)
f_unreg = self.res_unreg.f_test(np.eye(len(self.res_unreg.params))[:m])
f_reg = self.res_reg.f_test(np.eye(kvars + extra)[:m])
assert_almost_equal(f_unreg.fvalue, f_reg.fvalue, DECIMAL_2)
assert_almost_equal(f_unreg.pvalue, f_reg.pvalue, DECIMAL_3)

Expand Down Expand Up @@ -599,6 +607,51 @@ def setupClass(cls):
trim_mode='auto')


class TestNegativeBinomialL1Compatability(CheckL1Compatability):
@classmethod
def setupClass(cls):
cls.kvars = 10 # Number of variables
cls.m = 7 # Number of unregularized parameters
rand_data = sm.datasets.randhie.load()
rand_exog = rand_data.exog.view(float).reshape(len(rand_data.exog), -1)
rand_exog = sm.add_constant(rand_exog, prepend=True)
# Drop some columns and do an unregularized fit
exog_no_PSI = rand_exog[:, :cls.m]
cls.res_unreg = sm.NegativeBinomial(
rand_data.endog, exog_no_PSI).fit(method="newton", disp=False)
# Do a regularized fit with alpha, effectively dropping the last column
alpha = 10 * len(rand_data.endog) * np.ones(cls.kvars + 1)
alpha[:cls.m] = 0
alpha[-1] = 0 # don't penalize alpha
mod_reg = sm.NegativeBinomial(rand_data.endog, rand_exog)
cls.res_reg = mod_reg.fit_regularized(
method='l1', alpha=alpha, disp=False, acc=1e-10, maxiter=2000,
trim_mode='auto')


class TestNegativeBinomialGeoL1Compatability(CheckL1Compatability):
@classmethod
def setupClass(cls):
cls.kvars = 10 # Number of variables
cls.m = 7 # Number of unregularized parameters
rand_data = sm.datasets.randhie.load()
rand_exog = rand_data.exog.view(float).reshape(len(rand_data.exog), -1)
rand_exog = sm.add_constant(rand_exog, prepend=True)
# Drop some columns and do an unregularized fit
exog_no_PSI = rand_exog[:, :cls.m]
cls.res_unreg = sm.NegativeBinomial(
rand_data.endog, exog_no_PSI).fit(loglike_method='geometric',
method="newton", disp=False)
# Do a regularized fit with alpha, effectively dropping the last column
alpha = 10 * len(rand_data.endog) * np.ones(cls.kvars + 1)
alpha[:cls.m] = 0
alpha[-1] = 0 # don't penalize alpha
mod_reg = sm.NegativeBinomial(rand_data.endog, rand_exog)
cls.res_reg = mod_reg.fit_regularized(loglike_method='geometric',
method='l1', alpha=alpha, disp=False, acc=1e-10, maxiter=2000,
trim_mode='auto')


class TestLogitL1Compatability(CheckL1Compatability):
@classmethod
def setupClass(cls):
Expand Down

0 comments on commit 544400d

Please sign in to comment.