Discrete NegativeBinomialModel regularized_fit ValueError: matrices are not aligned #1453

Closed
KeninSydney opened this Issue Mar 8, 2014 · 5 comments

Projects

None yet

2 participants

@KeninSydney

This code returns:
ValueError: matrices are not aligned
The params array is always one element too short.

import numpy as np
import statsmodels.api as sm
sm.version.full_version
sm.show_versions()
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=False)

mod_nbin = sm.NegativeBinomial(rand_data.endog, rand_exog)
res_nbin = mod_nbin.fit_regularized(disp=False)
print res_nbin.summary()

Output:
INSTALLED VERSIONS
------------------
Python: 2.7.6.final.0
OS: Linux 2.6.32-042stab076.8 #1 SMP Tue May 14 20:38:14 MSK 2013 x86_64
byteorder: little
LC_ALL: None
LANG: None

Statsmodels
===========

Installed: 0.5.0 (/opt/anaconda/envs/np18py27-1.9/lib/python2.7/site-packages/statsmodels)

Required Dependencies
=====================

cython: 0.20 (/opt/anaconda/envs/np18py27-1.9/lib/python2.7/site-packages/Cython)
numpy: 1.8.0 (/opt/anaconda/envs/np18py27-1.9/lib/python2.7/site-packages/numpy)
scipy: 0.13.3 (/opt/anaconda/envs/np18py27-1.9/lib/python2.7/site-packages/scipy)
pandas: 0.13.1 (/opt/anaconda/envs/np18py27-1.9/lib/python2.7/site-packages/pandas)
    dateutil: 1.5 (/opt/anaconda/envs/np18py27-1.9/lib/python2.7/site-packages/dateutil)
patsy: 0.2.1 (/opt/anaconda/envs/np18py27-1.9/lib/python2.7/site-packages/patsy)

Optional Dependencies
=====================

matplotlib: 1.3.1 (/opt/anaconda/envs/np18py27-1.9/lib/python2.7/site-packages/matplotlib)
cvxopt: Not installed

Developer Tools
================

IPython: 2.0.0-dev (/opt/anaconda/envs/np18py27-1.9/lib/python2.7/site-packages/IPython)
    jinja2: 2.7.2 (/opt/anaconda/envs/np18py27-1.9/lib/python2.7/site-packages/jinja2)
sphinx: 1.2.1 (/opt/anaconda/envs/np18py27-1.9/lib/python2.7/site-packages/sphinx)
    pygments: 1.6 (/opt/anaconda/envs/np18py27-1.9/lib/python2.7/site-packages/pygments)
nose: 1.3.0 (/opt/anaconda/envs/np18py27-1.9/lib/python2.7/site-packages/nose)
virtualenv: Not installed


---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-1-5c06d718e7a6> in <module>()
      8 
      9 mod_nbin = sm.NegativeBinomial(rand_data.endog, rand_exog)
---> 10 res_nbin = mod_nbin.fit_regularized(disp=False)
     11 print res_nbin.summary()

/opt/anaconda/envs/np18py27-1.9/lib/python2.7/site-packages/statsmodels/discrete/discrete_model.pyc in fit_regularized(self, start_params, method, maxiter, full_output, disp, callback, alpha, trim_mode, auto_trim_tol, size_trim_tol, qc_tol, **kwargs)
    754                 full_output=full_output, disp=disp, callback=callback,
    755                 alpha=alpha, trim_mode=trim_mode, auto_trim_tol=auto_trim_tol,
--> 756                 size_trim_tol=size_trim_tol, qc_tol=qc_tol, **kwargs)
    757         if method in ['l1', 'l1_cvxopt_cp']:
    758             discretefit = L1CountResults(self, cntfit)

/opt/anaconda/envs/np18py27-1.9/lib/python2.7/site-packages/statsmodels/discrete/discrete_model.pyc in fit_regularized(self, start_params, method, maxiter, full_output, disp, callback, alpha, trim_mode, auto_trim_tol, size_trim_tol, qc_tol, qc_verbose, **kwargs)
    310                 method=method, maxiter=maxiter, full_output=full_output,
    311                 disp=disp, callback=callback, extra_fit_funcs=extra_fit_funcs,
--> 312                 cov_params_func=cov_params_func, **kwargs)
    313 
    314         return mlefit # up to subclasses to wrap results

/opt/anaconda/envs/np18py27-1.9/lib/python2.7/site-packages/statsmodels/base/model.pyc in fit(self, start_params, method, maxiter, full_output, disp, fargs, callback, retall, **kwargs)
    355                              disp=disp, maxiter=maxiter, callback=callback,
    356                              retall=retall, full_output=full_output,
--> 357                              hess=hess)
    358 
    359         if not full_output: # xopt should be None and retvals is argmin

/opt/anaconda/envs/np18py27-1.9/lib/python2.7/site-packages/statsmodels/base/l1_slsqp.pyc in fit_l1_slsqp(f, score, start_params, args, kwargs, disp, maxiter, callback, retall, full_output, hess)
     77         func, x0, f_ieqcons=f_ieqcons_wrap, fprime=fprime_wrap, acc=acc,
     78         iter=maxiter, disp=disp_slsqp, full_output=full_output,
---> 79         fprime_ieqcons=fprime_ieqcons_wrap)
     80     params = np.asarray(results[0][:k_params])
     81 

/opt/anaconda/envs/np18py27-1.9/lib/python2.7/site-packages/scipy/optimize/slsqp.pyc in fmin_slsqp(func, x0, eqcons, f_eqcons, ieqcons, f_ieqcons, bounds, fprime, fprime_eqcons, fprime_ieqcons, args, iter, acc, iprint, disp, full_output, epsilon)
    198 
    199     res = _minimize_slsqp(func, x0, args, jac=fprime, bounds=bounds,
--> 200                           constraints=cons, **opts)
    201     if full_output:
    202         return res['x'], res['fun'], res['nit'], res['status'], res['message']

/opt/anaconda/envs/np18py27-1.9/lib/python2.7/site-packages/scipy/optimize/slsqp.pyc in _minimize_slsqp(func, x0, args, jac, bounds, constraints, maxiter, ftol, iprint, disp, eps, **unknown_options)
    352 
    353             # Compute objective function
--> 354             fx = func(x)
    355             # Compute the constraints
    356             if cons['eq']:

/opt/anaconda/envs/np18py27-1.9/lib/python2.7/site-packages/scipy/optimize/optimize.pyc in function_wrapper(*wrapper_args)
    279     def function_wrapper(*wrapper_args):
    280         ncalls[0] += 1
--> 281         return function(*(wrapper_args + args))
    282 
    283     return ncalls, function_wrapper

/opt/anaconda/envs/np18py27-1.9/lib/python2.7/site-packages/statsmodels/base/l1_slsqp.pyc in <lambda>(x_full)
     68 
     69     ### Wrap up for use in fmin_slsqp
---> 70     func = lambda x_full: _objective_func(f, x_full, k_params, alpha, *args)
     71     f_ieqcons_wrap = lambda x_full: _f_ieqcons(x_full, k_params)
     72     fprime_wrap = lambda x_full: _fprime(score, x_full, k_params, alpha)

/opt/anaconda/envs/np18py27-1.9/lib/python2.7/site-packages/statsmodels/base/l1_slsqp.pyc in _objective_func(f, x_full, k_params, alpha, *args)
    133     x_added = x_full[k_params:]
    134     ## Return
--> 135     return f(x_params, *args) + (alpha * x_added).sum()
    136 
    137 

/opt/anaconda/envs/np18py27-1.9/lib/python2.7/site-packages/statsmodels/base/model.pyc in <lambda>(params, *args)
    327 
    328         nobs = self.endog.shape[0]
--> 329         f = lambda params, *args: -self.loglike(params, *args) / nobs
    330         score = lambda params: -self.score(params) / nobs
    331         try:

/opt/anaconda/envs/np18py27-1.9/lib/python2.7/site-packages/statsmodels/discrete/discrete_model.pyc in loglike(self, params)
   1840 
   1841         """
-> 1842         llf = np.sum(self.loglikeobs(params))
   1843         return llf
   1844 

/opt/anaconda/envs/np18py27-1.9/lib/python2.7/site-packages/statsmodels/discrete/discrete_model.pyc in _ll_nb2(self, params)
   1793         else:
   1794             alpha = params[-1]
-> 1795         return self._ll_nbin(params[:-1], alpha, Q=0)
   1796 
   1797     def _ll_nb1(self, params):

/opt/anaconda/envs/np18py27-1.9/lib/python2.7/site-packages/statsmodels/discrete/discrete_model.pyc in _ll_nbin(self, params, alpha, Q)
   1780     def _ll_nbin(self, params, alpha, Q=0):
   1781         endog = self.endog
-> 1782         mu = np.exp(np.dot(self.exog, params))
   1783         size = 1/alpha * mu**Q
   1784         prob = size/(size+mu)

ValueError: matrices are not aligned
@josef-pkt
Member

Thanks for reporting, I get the same

first guess: start_params have the wrong length

If I use an explicit start_params, as in the following, I don't get any error

mod_nbin = sm.NegativeBinomial(rand_data.endog, rand_exog)
res = mod_nbin.fit(disp=False)
res_nbin = mod_nbin.fit_regularized(start_params=res.params, disp=False)
print res_nbin.summary()

however, the params with fit and fit_regularized lead to the essentially the same parameters

>>> res.params
array([-0.05795413, -0.26781396,  0.04121061, -0.03814017,  0.26895441,
        0.03816294, -0.04411766,  0.01724125,  0.17800361,  0.66356013,
        1.29295351])
>>> res_nbin.params
array([-0.05795374, -0.26781391,  0.04121068, -0.03813948,  0.26895435,
        0.03816315, -0.04411768,  0.01724126,  0.1780036 ,  0.66356014,
        1.29295351])

the unregularized fit says converged=False in the summary
regularized fit has converged=True

Is this correct?
need to try out the fit_regularized options

@josef-pkt
Member

I don't remember what qc check is.

Below I get nans in gopt and hopt. why?

>>> res_nbin2 = mod_nbin.fit_regularized(start_params=res.params,alpha=0.25*amask, qc_verbose=True, retall=True)
  NIT    FC           OBJFUN            GNORM
    1     1     2.148790E+00     5.345830E-05
Optimization terminated successfully.    (Exit mode 0)
            Current function value: 2.1487897904
            Iterations: 1
            Function evaluations: 3
            Gradient evaluations: 1
QC check did not pass for 3 out of 11 parameters
Try increasing solver accuracy or number of iterations, decreasing alpha, or switch solvers
------ Recall the problem was rescaled by 1 / nobs ---
|passed    |alpha     |fprime    |param     |
--------------------------------------------
|True      |1.238e-05 |8.252e-06 |-5.795e-02|
|True      |1.238e-05 |1.305e-07 |-2.678e-01|
|False     |1.238e-05 |4.232e-05 |4.121e-02 |
|False     |1.238e-05 |2.823e-05 |-3.814e-02|
|True      |1.238e-05 |3.439e-06 |2.690e-01 |
|False     |1.238e-05 |8.643e-05 |3.816e-02 |
|True      |1.238e-05 |4.983e-06 |-4.412e-02|
|True      |1.238e-05 |8.597e-08 |1.724e-02 |
|True      |1.238e-05 |6.610e-07 |1.780e-01 |
|True      |1.238e-05 |7.236e-06 |6.636e-01 |
|True      |0.000e+00 |1.472e-08 |1.293e+00 |
Could not trim params automatically due to failed QC check.  Trimming using trim_mode == 'size' will still work.
>>> res_nbin2.mle_retvals
{'fopt': 2.1487897904014321, 'gopt': nan, 'hopt': nan, 'trimmed': array([False, False, False, False, False, False, False, False, False,
       False, False], dtype=bool), 'iterations': 1, 'converged': 'True'}

something looks weird, score doesn't look close to zero for all versions of fit that I tried, e.g. for regularized

>>> res_nbin.model.score(res_nbin.params)
array([-0.08566124, -0.00530527, -0.07763577, -0.10601139, -0.02636937,
       -0.01563745, -0.01031334, -0.01058   , -0.00927419, -0.01385042,
        0.02808423])
@josef-pkt
Member

aside: I think we shouldn't penalize the NegativeBinomial shape parameter in the same way as the coefficients for the linear part, e.g.

>>> amask = np.ones(len(res.params))
>>> amask[-1] = 0
>>> res_nbin2 = mod_nbin.fit_regularized(start_params=res.params, alpha=0.25*amask)
@josef-pkt
Member

using "nb1" works, but it looks difficult to figure out what alpha should be

>>> res_nbin3 = mod_nbin.fit_regularized(start_params=res.params,alpha=100*amask, qc_verbose=True, retall=True)
  NIT    FC           OBJFUN            GNORM
    1     1     2.151538E+00     1.566259E-02
    2     3     2.151532E+00     5.249264E-02
    3     5     2.151495E+00     4.557341E-02
    4     7     2.151429E+00     3.590288E-02
    5     8     2.151119E+00     3.016578E-02
    6     9     2.150752E+00     2.210603E-02
    7    10     2.150661E+00     2.018570E-02
    8    11     2.150639E+00     2.084313E-02
    9    12     2.150569E+00     2.570237E-02
   10    13     2.150544E+00     2.481202E-02
   11    14     2.150518E+00     2.287230E-02
   12    15     2.150501E+00     2.261447E-02
   13    16     2.150468E+00     2.281161E-02
   14    17     2.150457E+00     2.157790E-02
   15    18     2.150453E+00     2.095163E-02
   16    19     2.150452E+00     2.112453E-02
   17    20     2.150452E+00     2.123148E-02
   18    21     2.150452E+00     2.134927E-02
   19    22     2.150451E+00     2.134492E-02
   20    23     2.150451E+00     2.122605E-02
   21    24     2.150451E+00     2.114686E-02
   22    25     2.150451E+00     2.110636E-02
   23    26     2.150451E+00     2.106363E-02
   24    27     2.150451E+00     2.103672E-02
   25    28     2.150451E+00     2.104644E-02
   26    29     2.150451E+00     2.109452E-02
   27    30     2.150451E+00     2.113491E-02
   28    31     2.150451E+00     2.115827E-02
   29    32     2.150451E+00     2.118305E-02
   30    33     2.150450E+00     2.118928E-02
   31    34     2.150450E+00     2.116706E-02
   32    35     2.150450E+00     2.114533E-02
   33    36     2.150450E+00     2.113928E-02
Optimization terminated successfully.    (Exit mode 0)
            Current function value: 2.15045044517
            Iterations: 33
            Function evaluations: 36
            Gradient evaluations: 33
>>> print res_nbin3.summary()
                     NegativeBinomial Regression Results                      
==============================================================================
Dep. Variable:                      y   No. Observations:                20190
Model:               NegativeBinomial   Df Residuals:                    20182
Method:                           MLE   Df Model:                            7
Date:                Sat, 08 Mar 2014   Pseudo R-squ.:                 0.02053
Time:                        16:34:52   Log-Likelihood:                -43292.
converged:                       True   LL-Null:                       -44199.
                                        LLR p-value:                     0.000
==============================================================================
                 coef    std err          z      P>|z|      [95.0% Conf. Int.]
------------------------------------------------------------------------------
x1            -0.0603      0.005    -11.471      0.000        -0.071    -0.050
x2            -0.2536      0.019    -13.073      0.000        -0.292    -0.216
x3             0.0422      0.003     12.466      0.000         0.036     0.049
x4            -0.0325      0.003    -10.945      0.000        -0.038    -0.027
x5             0.1367      0.024      5.811      0.000         0.091     0.183
x6             0.0339      0.001     32.318      0.000         0.032     0.036
x7                  0        nan        nan        nan           nan       nan
x8                  0        nan        nan        nan           nan       nan
x9                  0        nan        nan        nan           nan       nan
const          0.6978      0.021     33.228      0.000         0.657     0.739
alpha          3.7085      0.062     59.518      0.000         3.586     3.831
==============================================================================

which looks pretty good compared to alpha=0, res_nbin,
the least significant params are set to zero, the most of the others don't change much

>>> res_nbin3.params / res_nbin.params
array([ 0.92337179,  0.8565406 ,  1.02416303,  1.01404461,  0.71715496,
        1.06321644, -0.        , -0.        ,  0.        ,  0.92140039,
        0.99392089])

>>> res_nbin.pvalues
array([  3.73543154e-034,   4.95150690e-051,   1.39309038e-034,
         2.04536907e-026,   1.13813373e-015,   1.55134219e-193,
         4.66928683e-002,   1.03639840e-001,   4.08960285e-002,
         2.03136069e-272,   0.00000000e+000])
@josef-pkt
Member

maybe everything is fine except for the starting value, nb2 also looks good if alpha is large enough

mod_nbin = sm.NegativeBinomial(rand_data.endog, rand_exog, loglike_method='nb2')
res = mod_nbin.fit(disp=False, method='bfgs')
res_nbin = mod_nbin.fit_regularized(start_params=res.params, disp=False)
print res_nbin.summary()

res_nbin3 = mod_nbin.fit_regularized(start_params=res.params,alpha=100*amask, qc_verbose=True, retall=True)
print res_nbin3.summary()

score at optimum is large in absolute terms, but not in relative terms

>>> res_nbin.model.score(res_nbin.params)
array([  1.35411669e-01,   2.11502368e-02,  -2.83535542e-01,
         1.39332197e-01,  -5.81716604e-02,  -6.38596202e-01,
        -6.46757111e-02,   5.31774935e-03,  -1.22772331e-02,
        -5.08886427e-02,  -1.03765058e-04])

>>> res_nbin.model.score(res_nbin.params) / res_nbin.llf
array([ -3.12125953e-06,  -4.87516170e-07,   6.53553731e-06,
        -3.21162830e-06,   1.34086561e-06,   1.47197394e-05,
         1.49078496e-06,  -1.22574930e-07,   2.82992087e-07,
         1.17299094e-06,   2.39180035e-09])

IIRC, we used this dataset for the Poisson test case, and it was the reason to normalize loglike by nobs in the optimization, because of the large number of observations and large llf.sum() instead of llf.mean()

@josef-pkt josef-pkt added a commit to josef-pkt/statsmodels that referenced this issue Apr 22, 2014
@josef-pkt josef-pkt BUG NegativeBinomial add fit_regularized closes #1453 closes #1454
    adjust test to handle extra parameter
16dfb7d
@josef-pkt josef-pkt added a commit to josef-pkt/statsmodels that referenced this issue Jul 22, 2014
@josef-pkt josef-pkt BUG NegativeBinomial add fit_regularized closes #1453 closes #1454
    adjust test to handle extra parameter
544400d
@josef-pkt josef-pkt closed this in #1623 Jul 22, 2014
@PierreBdR PierreBdR pushed a commit to PierreBdR/statsmodels that referenced this issue Sep 2, 2014
@josef-pkt josef-pkt BUG NegativeBinomial add fit_regularized closes #1453 closes #1454
    adjust test to handle extra parameter
5be6a86
@yarikoptic yarikoptic added a commit to yarikoptic/statsmodels that referenced this issue Oct 23, 2014
@yarikoptic yarikoptic Merge commit 'v0.5.0-1189-g48f7a21' into releases
* commit 'v0.5.0-1189-g48f7a21': (970 commits)
  REF _fit_newton use np.asarray(hess)
  REF: _fit_newton, add ridge_factor option, default 1e-10
  REF: _fit_newton add Ridge correction to Hessian, see also #953
  Bug comment: commented out code, weigths don't sum to 1 see #1845
  CLN: TST cleanup comment, unused code.
  REF: "is" None (not ==)
  BUG: fix if alpha is scalar, TST: try standardized
  REF: NegativeBinomial fit_regularized, try regularized Poisson for start_params
  DOC: add comment in notes about not penalizing NegBin shape parameter [skip ci]
  TST: TestNegativeBinomialL1Compatability use desired as start_params
  TST: adjust test, precision and start_params (failure on TravisCI)
  BUG: NegativeBinomial fix aic, bic for 'geometric', adjust tests
  REF add k_extra to NegativeBinomial, don't count it in df_model, df_resid
  TST explicitely define k_extra in test class
  CLN: a few cosmetic changes
  TST fix negbin geometric test for fit_regularized
  CLN L1NegativeBinomialResults rename closes #1615,     remove redundant `__init__`
  BUG NegativeBinomial add fit_regularized closes #1453 closes  #1454     adjust test to handle extra parameter
  REF: has_constant remove special code in linear_model, is moved to data
  Fix const_idx with multiple const, more tests
  ...
18a1faa
@yarikoptic yarikoptic added a commit to yarikoptic/statsmodels that referenced this issue Oct 23, 2014
@yarikoptic yarikoptic Merge branch 'releases' into debian-experimental
* releases: (970 commits)
  REF _fit_newton use np.asarray(hess)
  REF: _fit_newton, add ridge_factor option, default 1e-10
  REF: _fit_newton add Ridge correction to Hessian, see also #953
  Bug comment: commented out code, weigths don't sum to 1 see #1845
  CLN: TST cleanup comment, unused code.
  REF: "is" None (not ==)
  BUG: fix if alpha is scalar, TST: try standardized
  REF: NegativeBinomial fit_regularized, try regularized Poisson for start_params
  DOC: add comment in notes about not penalizing NegBin shape parameter [skip ci]
  TST: TestNegativeBinomialL1Compatability use desired as start_params
  TST: adjust test, precision and start_params (failure on TravisCI)
  BUG: NegativeBinomial fix aic, bic for 'geometric', adjust tests
  REF add k_extra to NegativeBinomial, don't count it in df_model, df_resid
  TST explicitely define k_extra in test class
  CLN: a few cosmetic changes
  TST fix negbin geometric test for fit_regularized
  CLN L1NegativeBinomialResults rename closes #1615,     remove redundant `__init__`
  BUG NegativeBinomial add fit_regularized closes #1453 closes  #1454     adjust test to handle extra parameter
  REF: has_constant remove special code in linear_model, is moved to data
  Fix const_idx with multiple const, more tests
  ...
8f2ab7a
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment