BUG Negbin fit regularized #1623

Merged
merged 14 commits into from Jul 22, 2014

Projects

None yet

3 participants

@josef-pkt
Member

adds fit_regularized method to NegativeBinomial
fixes #1453 and #1454

note:
code is mostly copy paste see #1615
took me a while to figure out copy/paste/adjust errors for example #1617

one issue not clear:
NegativeBinomial transforms parameters in some cases, I chose not to: self._transparams = False

TODO: how can I adjust the docstring?
by default if alpha is scalar, then the penalization is not applied to the extra, shape parameter alpha of nb1 and nb2

unittests follow the same pattern as the other L1 test, compare reduced unregularized estimation with regularized estimation that doesn't penalize the original parameters

@jseabold jseabold commented on the diff Apr 22, 2014
statsmodels/discrete/tests/test_discrete.py
@@ -524,37 +524,45 @@ class CheckL1Compatability(object):
def test_params(self):
m = self.m
@jseabold
jseabold Apr 22, 2014 Member

Can use use self.res_reg.model.k_exog here to make sure it's being used systematically?

@josef-pkt
josef-pkt Apr 22, 2014 Member

m is a setting for the test, it is here the number unpenalized beta params,
m=k_exog in the reduced version but m < k_exog in the regularized version.

in my interpreter session I don't have a model.k_exog. where is k_exog set?

@jseabold
jseabold Apr 22, 2014 Member

Right. From the description I thought that the tests were testing the unregularized vs. regularized without any actual regularization.

It looks like k_exog is only systematically set in ARIMA so far. Should be moved up to base. Still TODO in #487.

@josef-pkt josef-pkt commented on an outdated diff Apr 22, 2014
statsmodels/discrete/tests/test_discrete.py
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)
@josef-pkt
josef-pkt Apr 22, 2014 Member

this looks wrong both, reg and unreg, have the extra parameter
L1 counts shape parameter in df_model and df_resid, NegativeBinomial does't
?

@josef-pkt
Member

I'm giving up on #1624 for this PR

temporary solution:
add k_extra as attribute to NegativeBinomialModel, and adjust df_model and df_resid in L1CountResults so it isn't included.
that leaves df_model and df_resid at it's current non-L1 definition for nb1 and nb2

@josef-pkt
Member

TravisCI fails in all versions
looks like a optimization problem in unregularized fit with nb2, gets a nan in params
second failure is small precision issue

two possible problems in this PR for optimization:

  • we use the full Poisson fit as starting value, which might not be very good if there are many extra variables.
  • we don't transform alpha during optimization: I haven't checked if we avoid negative alpha.
@josef-pkt
Member

still red
precision failure is 2405.88524487 / 2405.85274708 - 1 = 1.3507805097123793e-05 ("amended")
why do we check f_test anyway? just to make sure some results are correct (we use normal distribution)

no idea about optimization problem, which param is nan ?

@josef-pkt
Member

all green when I'm cheating with the start_params

@josef-pkt josef-pkt referenced this pull request Apr 22, 2014
Open

Poisson start_params #1511

@josef-pkt
Member

ready to merge, I leave the start_params to another PR.

according to the network graph, this didn't branch off from master, might need rebase before merging

@josef-pkt josef-pkt added this to the 0.6 milestone May 22, 2014
@coveralls

Coverage Status

Coverage increased (+1.63%) when pulling 0739786 on josef-pkt:negbin_fit_regularized into b52bc09 on statsmodels:master.

@josef-pkt
Member

That's better, after standardizing exog (mean=0, std=1) TravisCI manages to estimate correctly. No more exceptions.
Note: data transformation changes the penalization here because I didn't correct for it by adjusting alpha.

#1649 example with convergence or optimization failure

another review of this PR and it is ready to merge

@josef-pkt josef-pkt commented on the diff Jul 21, 2014
statsmodels/discrete/discrete_model.py
+ if start_params == None:
+ # Use poisson fit as first guess.
+ start_params = Poisson(self.endog, self.exog).fit_regularized(
+ start_params=start_params, method=method, maxiter=maxiter,
+ full_output=full_output, disp=0, callback=callback,
+ alpha=alpha_p, trim_mode=trim_mode, auto_trim_tol=auto_trim_tol,
+ size_trim_tol=size_trim_tol, qc_tol=qc_tol, **kwargs).params
+ if self.loglike_method.startswith('nb'):
+ start_params = np.append(start_params, 0.1)
+
+ cntfit = super(CountModel, self).fit_regularized(
+ start_params=start_params, method=method, maxiter=maxiter,
+ full_output=full_output, disp=disp, callback=callback,
+ alpha=alpha, trim_mode=trim_mode, auto_trim_tol=auto_trim_tol,
+ size_trim_tol=size_trim_tol, qc_tol=qc_tol, **kwargs)
+ if method in ['l1', 'l1_cvxopt_cp']:
@josef-pkt
josef-pkt Jul 21, 2014 Member

input/argument checking should be at the beginning before we calculate anything.

@josef-pkt josef-pkt commented on the diff Jul 21, 2014
statsmodels/discrete/discrete_model.py
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):
+
+ if self.loglike_method.startswith('nb') and (np.size(alpha) == 1 and
+ alpha != 0):
+ # don't penalize alpha if alpha is scalar
@josef-pkt
josef-pkt Jul 21, 2014 Member

2 different alpha !

@josef-pkt josef-pkt commented on an outdated diff Jul 21, 2014
statsmodels/discrete/discrete_model.py
+ def fit_regularized(self, start_params=None, method='l1',
+ maxiter='defined_by_method', full_output=1, disp=1, callback=None,
+ alpha=0, trim_mode='auto', auto_trim_tol=0.01, size_trim_tol=1e-4,
+ qc_tol=0.03, **kwargs):
+
+ if self.loglike_method.startswith('nb') and (np.size(alpha) == 1 and
+ alpha != 0):
+ # don't penalize alpha if alpha is scalar
+ alpha = alpha * np.ones(len(start_params))
+ alpha[-1] = 0
+
+ # alpha for regularized poisson to get starting values
+ alpha_p = alpha[:-1] if (self.k_extra and np.size(alpha) > 1) else alpha
+
+ self._transparams = False
+ if start_params == None:
@josef-pkt josef-pkt commented on the diff Jul 21, 2014
statsmodels/discrete/discrete_model.py
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):
+
+ if self.loglike_method.startswith('nb') and (np.size(alpha) == 1 and
@josef-pkt
josef-pkt Jul 21, 2014 Member

we have now self.k_extra as counter for extra params, instead of self.loglike_method.startswith('nb')

@josef-pkt josef-pkt commented on an outdated diff Jul 21, 2014
statsmodels/discrete/tests/test_discrete.py
+ cls.res_unreg = mod_unreg.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
+ # try constant only
+ start_params = np.zeros(len(alpha))
+ start_params[0] = rand_data.endog.mean()
+ start_params[-1] = 0.5
+ # cheating: use params as starting_values
+ start_params = np.zeros(len(alpha))
+ start_params[:cls.m] = cls.res_unreg.params[:-1]
+ start_params[-1] = cls.res_unreg.params[-1]
+
+ mod_reg = sm.NegativeBinomial(rand_data.endog, rand_exog)
+ cls.res_reg = mod_reg.fit_regularized(#start_params=start_params,
@josef-pkt
josef-pkt Jul 21, 2014 Member

commented out start_params - clean up
start_params are not used anymore in here.

@coveralls

Coverage Status

Coverage increased (+0.02%) when pulling 6544287 on josef-pkt:negbin_fit_regularized into 40bd7e0 on statsmodels:master.

@josef-pkt josef-pkt merged commit a40711c into statsmodels:master Jul 22, 2014

2 checks passed

continuous-integration/appveyor AppVeyor build succeeded
Details
continuous-integration/travis-ci The Travis CI build passed
Details
@josef-pkt josef-pkt deleted the josef-pkt:negbin_fit_regularized branch Jul 22, 2014
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment