在Python中建模时检测mulicollinear或具有线性组合的列：LinAlgError

debugcn 发表于 Dev

埃克塔

我正在为具有34个因变量的logit模型的数据建模，并且不断抛出奇异矩阵错误，如下所示：

Traceback (most recent call last):
  File "<pyshell#1116>", line 1, in <module>
    test_scores  = smf.Logit(m['event'], train_cols,missing='drop').fit()
  File "/usr/local/lib/python2.7/site-packages/statsmodels-0.5.0-py2.7-linux-i686.egg/statsmodels/discrete/discrete_model.py", line 1186, in fit
    disp=disp, callback=callback, **kwargs)
  File "/usr/local/lib/python2.7/site-packages/statsmodels-0.5.0-py2.7-linux-i686.egg/statsmodels/discrete/discrete_model.py", line 164, in fit
    disp=disp, callback=callback, **kwargs)
  File "/usr/local/lib/python2.7/site-packages/statsmodels-0.5.0-py2.7-linux-i686.egg/statsmodels/base/model.py", line 357, in fit
    hess=hess)
  File "/usr/local/lib/python2.7/site-packages/statsmodels-0.5.0-py2.7-linux-i686.egg/statsmodels/base/model.py", line 405, in _fit_mle_newton
    newparams = oldparams - np.dot(np.linalg.inv(H),
  File "/usr/local/lib/python2.7/site-packages/numpy/linalg/linalg.py", line 445, in inv
    return wrap(solve(a, identity(a.shape[0], dtype=a.dtype)))
  File "/usr/local/lib/python2.7/site-packages/numpy/linalg/linalg.py", line 328, in solve
    raise LinAlgError, 'Singular matrix'
LinAlgError: Singular matrix

那是我在用这种方法将矩阵简化为独立列时的情况

def independent_columns(A, tol = 0):#1e-05):
    """
    Return an array composed of independent columns of A.

    Note the answer may not be unique; this function returns one of many
    possible answers.

    https://stackoverflow.com/q/13312498/190597 (user1812712)
    http://math.stackexchange.com/a/199132/1140 (Gerry Myerson)
    http://mail.scipy.org/pipermail/numpy-discussion/2008-November/038705.html
        (Anne Archibald)

    >>> A = np.array([(2,4,1,3),(-1,-2,1,0),(0,0,2,2),(3,6,2,5)])
    2 4 1 3
    -1 -2 1 0
    0 0 2 2
    3 6 2 5
    # try with checking the rank of matrixs 
    >>> independent_columns(A)
    np.array([[1, 4],
              [2, 5],
              [3, 6]])
    """
    Q, R = linalg.qr(A)
    independent = np.where(np.abs(R.diagonal()) > tol)[0]
    #print independent
    return A[:, independent], independent


A,independent_col_indexes=independent_columns(train_cols.as_matrix(columns=None)) 
#train_cols will not be converted back from a df to a  matrix object,so doing this explicitly
A2=pd.DataFrame(A, columns=train_cols.columns[independent_col_indexes])

test_scores = smf.Logit(m['event'],A2,missing='drop').fit()

我仍然收到LinAlgError，尽管我希望现在可以降低矩阵等级。

另外，我看到了np.linalg.matrix_rank(train_cols)返回33（即，在调用Independent_columns函数之前，“ x”列总数为34（即len(train_cols.ix[0])=34），这意味着我没有完整的排名矩阵），而np.linalg.matrix_rank(A2)返回33（即，我删除了一个列，但是我仍然看到LinAlgError，当我运行时test_scores = smf.Logit(m['event'],A2,missing='drop').fit()，我想念的是什么？

参考上面的代码-如何在协方差矩阵中查找退化的行/列

我试图通过一次引入每个变量来开始建立模型，这不会给我单一的矩阵误差，但是我宁愿有一个确定性的方法，并让我知道我在做什么错＆如何消除这些列。

编辑（通过以下@ user333700发布建议更新）

1.是的，“ A2”的排名降低为33。IE。len(A2.ix[0]) =34->表示可能的共线列未删除-我应该增加“ tol”的公差以获取A2的排名（及其列数）为33。如果我将tol更改为上面的“ 1e-05”，然后我得到了len(A2.ix[0]) =33，这对我来说建议tol> 0（严格来说）是一个指标。在此之后，我只是做了相同的操作test_scores = smf.Logit(m['event'],A2,missing='drop').fit()，没有nm获得收敛。

2.尝试“ nm”方法后出错。不过，奇怪的是，如果我只进行20,000行，我确实会得到结果。由于它没有显示出内存错误，而是“ Inverting hessian failed, no bse or cov_params available”-我假设是有多个几乎相似的记录-您会说什么？

m  = smf.Logit(data['event_custom'].ix[0:1000000] , train_cols.ix[0:1000000],missing='drop')
test_scores=m.fit(start_params=None,method='nm',maxiter=200,full_output=1)
Warning: Maximum number of iterations has been exceeded

Warning (from warnings module):
  File "/usr/local/lib/python2.7/site-packages/statsmodels-0.5.0-py2.7-linux-i686.egg/statsmodels/base/model.py", line 374
    warn(warndoc, Warning)
Warning: Inverting hessian failed, no bse or cov_params available


test_scores.summary()

Traceback (most recent call last):
  File "<pyshell#17>", line 1, in <module>
    test_scores.summary()
  File "/usr/local/lib/python2.7/site-packages/statsmodels-0.5.0-py2.7-linux-i686.egg/statsmodels/discrete/discrete_model.py", line 2396, in summary
    yname_list)
  File "/usr/local/lib/python2.7/site-packages/statsmodels-0.5.0-py2.7-linux-i686.egg/statsmodels/discrete/discrete_model.py", line 2253, in summary
    use_t=False)
  File "/usr/local/lib/python2.7/site-packages/statsmodels-0.5.0-py2.7-linux-i686.egg/statsmodels/iolib/summary.py", line 826, in add_table_params
    use_t=use_t)
  File "/usr/local/lib/python2.7/site-packages/statsmodels-0.5.0-py2.7-linux-i686.egg/statsmodels/iolib/summary.py", line 447, in summary_params
    std_err = results.bse
  File "/usr/local/lib/python2.7/site-packages/statsmodels-0.5.0-py2.7-linux-i686.egg/statsmodels/tools/decorators.py", line 95, in __get__
    _cachedval = self.fget(obj)
  File "/usr/local/lib/python2.7/site-packages/statsmodels-0.5.0-py2.7-linux-i686.egg/statsmodels/base/model.py", line 1037, in bse
    return np.sqrt(np.diag(self.cov_params()))
  File "/usr/local/lib/python2.7/site-packages/statsmodels-0.5.0-py2.7-linux-i686.egg/statsmodels/base/model.py", line 1102, in cov_params
    raise ValueError('need covariance of parameters for computing '
ValueError: need covariance of parameters for computing (unnormalized) covariances

编辑2 ：（更新后由@ user333700发布建议）

重申我要建模的内容-不到约1％的总用户“转化”（成功的结果）-因此，我采取了35（+ ve）/ 65（-ve）的均衡样本

我怀疑该模型虽然收敛，但仍不够鲁棒。因此，将使用“ start_params”作为来自不同数据集的早期迭代的参数。此编辑是为了确认是否可以将“ start_params”输入到结果中，如下所示：

A,independent_col_indexes=independent_columns(train_cols.as_matrix(columns=None))
A2=pd.DataFrame(A, columns=train_cols.columns[independent_col_indexes])
m  = smf.Logit(data['event_custom'], A2,missing='drop')
#m  = smf.Logit(data['event_custom'], train_cols,missing='drop')#,method='nm').fit()#This doesnt work, so tried 'nm' which work, but used lasso, as nm did not converge.
test_scores=m.fit_regularized(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=0.0001, qc_tol=0.03)

a_good_looking_previous_result.params=test_scores.params #storing the parameters of pass1 to feed into pass2

test_scores.params
bidfloor_Quartile_modified_binned_0               0.305765
connectiontype_binned_0                          -0.436798
day_custom_binned_Fri                            -0.040269
day_custom_binned_Mon                             0.138599
day_custom_binned_Sat                            -0.319997
day_custom_binned_Sun                            -0.236507
day_custom_binned_Thu                            -0.058922
user_agent_device_family_binned_iPad            -10.793270
user_agent_device_family_binned_iPhone           -8.483099
user_agent_masterclass_binned_apple               9.038889
user_agent_masterclass_binned_generic            -0.760297
user_agent_masterclass_binned_samsung            -0.063522
log_height_width                                  0.593199
log_height_width_ScreenResolution                -0.520836
productivity                                     -1.495373
games                                             0.706340
entertainment                                    -1.806886
IAB24                                             2.531467
IAB17                                             0.650327
IAB14                                             0.414031
utilities                                         9.968253
IAB1                                              1.850786
social_networking                                -2.814148
IAB3                                             -9.230780
music                                             0.019584
IAB9                                             -0.415559
C(time_day_modified)[(6, 12]]:C(country)[AUS]    -0.103003
C(time_day_modified)[(0, 6]]:C(country)[HKG]      0.769272
C(time_day_modified)[(6, 12]]:C(country)[HKG]     0.406882
C(time_day_modified)[(0, 6]]:C(country)[IDN]      0.073306
C(time_day_modified)[(6, 12]]:C(country)[IDN]    -0.207568
C(time_day_modified)[(0, 6]]:C(country)[IND]      0.033370
... more params here

现在，在另一个数据集（pass2，用于索引）上，我与下面的模型相同-：即。我读取了一个新的数据框，进行了所有变量转换，然后通过Logit进行建模，如之前所述。

m_pass2  = smf.Logit(data['event_custom'], A2_pass2,missing='drop')
test_scores_pass2=m_pass2.fit_regularized(start_params=a_good_looking_previous_result.params, 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=0.0001, qc_tol=0.03)

并且可以通过从较早的过程中获取“ start_params”来保持迭代。

约瑟夫

几点要点：

您需要tol> 0才能检测到接近完美的共线性，这也可能在以后的计算中引起数值问题。检查的列数，A2以查看是否确实删除了列。

Logit需要使用exog进行一些非线性计算，因此，即使设计矩阵与理想共线性不是十分接近，对数似然，导数或Hessian计算的转换变量仍可能最终会遇到数值问题，例如奇异的黑森州。

（当我们接近浮点精度1e-15、1e-16时，所有这些都是浮点问题。matrix_rank和类似的linalg函数的默认阈值有时会有所不同，这可能意味着在某些情况下，一个函数会将其标识为单数，而另一种则没有。）

离散模型（包括Logit）的默认优化方法是简单的牛顿法，该方法在相当好的情况下速度很快，但在条件恶劣的情况下可能会失败。您可以尝试使用其他优化器之一，该优化器将是scipy.optimize中的优化器之一。method='nm'通常它非常健壮但运行缓慢，method='bfgs'在许多情况下都可以正常工作，但也可能会遇到收敛问题。

尽管如此，即使其他优化方法之一成功了，仍然有必要检查结果。通常，一种方法的失败意味着可能无法很好地定义模型或估计问题。

检查是仅仅是起始值不好的问题还是规范问题的一种好方法是先运行method='nm'，然后运行更准确的方法之一，例如将估算值用作起始值，newton或者bfgs使用nm估算值作为起始值，然后查看它是否从良好的起始值成功。

本文收集自互联网，转载请注明来源。

如有侵权，请联系[email protected] 删除。

编辑于2021-06-7

我来说两句

0条评论

登录后参与评论

来自分类Dev

Related 相关文章

文章