我想使用Statsmodels在Python中进行逻辑回归。
X和y各有750行,y是二进制结果,X是10个特征(包括intecept)。
下面是X的前12行(最后一列是截距):

      lngdp_      lnpop    sxp      sxp2    gy1    frac  etdo4590  geogia  \
0   7.367709  16.293980  0.190  0.036100 -1.682   132.0         1   0.916
1   7.509883  16.436258  0.193  0.037249  2.843   132.0         1   0.916
2   7.759187  16.589224  0.269  0.072361  4.986   132.0         1   0.916
3   7.922261  16.742384  0.368  0.135424  3.261   132.0         1   0.916
4   8.002359  16.901037  0.170  0.028900  1.602   132.0         1   0.916
5   7.929126  17.034786  0.179  0.032041 -1.465   132.0         1   0.916
6   6.594413  15.627563  0.360  0.129600 -9.321  4134.0         0   0.648
7   6.448889  16.037861  0.476  0.226576 -2.356  3822.0         0   0.648
8   8.520786  16.919334  0.048  0.002304  2.349   434.0         1   0.858
9   8.637107  16.991980  0.050  0.002500  2.326   434.0         1   0.858
10  8.708144  17.075489  0.042  0.001764  1.421   465.0         1   0.858
11  8.780480  17.151779  0.080  0.006400  1.447   496.0         1   0.858

    peace  intercept
0    24.0        1.0
1    84.0        1.0
2   144.0        1.0
3   204.0        1.0
4   264.0        1.0
5   324.0        1.0
6     1.0        1.0
7    16.0        1.0
8   112.0        1.0
9   172.0        1.0
10  232.0        1.0
11  292.0        1.0

这是我的代码:
import statsmodels.api as sm

logit = sm.Logit(y, X, missing='drop')
result = logit.fit()
print(result.summary())

这是输出:
     Optimization terminated successfully.

     Current function value: inf

     Iterations 9

/家庭/ipattern/anaconda3/lib/python3.6/site packages/statsmodels/discrete/discrete_model.py:1214:
运行时警告:exp中遇到溢出
返回1/(1+np.exp(-X))
/家庭/ipattern/anaconda3/lib/python3.6/site packages/statsmodels/discrete/discrete_model.py:1264:
RuntimeWarning:除以日志中遇到的零
返回
np.sum(np.log(self.cdf(q*np.dot(X,params)))
                           Logit Regression Results
==============================================================================
Dep. Variable:                  warsa   No. Observations:                  750
Model:                          Logit   Df Residuals:                      740
Method:                           MLE   Df Model:                            9
Date:                Tue, 12 Sep 2017   Pseudo R-squ.:                    -inf
Time:                        11:16:58   Log-Likelihood:                   -inf
converged:                       True   LL-Null:                   -4.6237e+05
                                        LLR p-value:                     1.000
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
lngdp_        -0.9504      0.245     -3.872      0.000      -1.431      -0.469
lnpop          0.5105      0.128      3.975      0.000       0.259       0.762
sxp           16.7734      5.206      3.222      0.001       6.569      26.978
sxp2         -23.8004     10.040     -2.371      0.018     -43.478      -4.123
gy1           -0.0980      0.041     -2.362      0.018      -0.179      -0.017
frac          -0.0002    9.2e-05     -2.695      0.007      -0.000   -6.76e-05
etdo4590       0.4801      0.328      1.463      0.144      -0.163       1.124
geogia        -0.9919      0.909     -1.091      0.275      -2.774       0.790
peace         -0.0038      0.001     -3.808      0.000      -0.006      -0.002
intercept     -3.4375      2.486     -1.383      0.167      -8.310       1.435
==============================================================================

底部的系数、标准差、p值等是正确的(我知道这是因为我有“解”)。
但正如你所看到的,我认为这是错误的。
我收到两个警告。显然statsmodels在某些地方做np.exp(BIGNUMBER),例如np.exp(999)和np.log(0)。
还有Current function value is infPseudo R-squ. is -inf,我认为不应该是Log-Likelihood is -inf
我做错什么了?
编辑:
X.描述():
           lngdp_       lnpop         sxp        sxp2         gy1  \
count  750.000000  750.000000  750.000000  750.000000  750.000000
mean     7.766948   15.702191    0.155329    0.043837    1.529772
std      1.045121    1.645154    0.140486    0.082838    3.546621
min      5.402678   11.900227    0.002000    0.000004  -13.088000
25%      6.882694   14.723123    0.056000    0.003136   -0.411250
50%      7.696212   15.680984    0.111000    0.012321    1.801000
75%      8.669355   16.652981    0.203000    0.041209    3.625750
max      9.851826   20.908354    0.935000    0.874225   14.409000

              frac    etdo4590      geogia       peace  intercept
count   750.000000  750.000000  750.000000  750.000000      750.0
mean   1812.777333    0.437333    0.600263  348.209333        1.0
std    1982.106029    0.496388    0.209362  160.941996        0.0
min      12.000000    0.000000    0.000000    1.000000        1.0
25%     176.000000    0.000000    0.489250  232.000000        1.0
50%     864.000000    0.000000    0.608000  352.000000        1.0
75%    3375.000000    1.000000    0.763000  472.000000        1.0
max    6975.000000    1.000000    0.971000  592.000000        1.0

logit.loglikeobs(结果.params):
array([ -4.61803704e+01,  -2.26983454e+02,  -2.66741244e+02,
        -2.60206733e+02,  -4.75585266e+02,  -1.76454554e+00,
        -4.86048292e-01,  -8.02300533e-01,             -inf,
                   -inf,             -inf,             -inf,
                   -inf,             -inf,             -inf,
                   -inf,             -inf,             -inf,
                   -inf,             -inf,             -inf,
                   -inf,             -inf,             -inf,
                   -inf,             -inf,             -inf,
                   -inf,             -inf,             -inf,
                   -inf,             -inf,  -6.02780923e+02,
        -4.12209348e+02,  -6.42901288e+02,  -6.94331125e+02,
                   -inf,             -inf,             -inf,
                   -inf,             -inf,             -inf,
                   -inf,             -inf,             -inf, ...

(logit.exog*np.array(result.params)).min(0):
array([ -9.36347474,   6.07506083,   0.03354677, -20.80694575,
        -1.41162588,  -1.72895247,   0.        ,  -0.9631801 ,
        -2.23188846,  -3.4374963 ])

数据集:
X:https://pastebin.com/VRNSepBg
是:https://pastebin.com/j2Udyc7m

最佳答案

我很惊讶在这种情况下它仍然收敛。
当x值较大时,Logit或Poisson中使用的exp函数可能存在溢出的收敛问题。这通常可以通过重新缩放回归器来避免。
但是,在这种情况下,我的猜测是x中的异常值。第6列的值类似于4134.0,而其他列的值要小得多。
您可以检查每个观察值logit.loglikeobs(result.params)的对数可能性,以查看哪些观察值可能导致问题,其中logit是引用模型的名称
另外,每个预测因子的贡献也可能有帮助,例如
np.argmax(np.abs(logit.exog * result.params), 0)

(logit.exog * result.params).min(0)
如果只是一个或几个观察结果,那么放弃它们可能会有帮助。重新缩放exog很可能没有帮助,因为在收敛时,它只会被估计系数的重新缩放所补偿。
还要检查是否有编码错误或大值作为丢失值的占位符。
编辑
鉴于logLoops中的-inf的数量似乎很大,我认为可能存在比离群值更为根本的问题,因为logit模型不是该数据集的正确指定的最大似然模型。
一般来说有两种可能性(因为我没有看到数据集):
完全分离:Logit假设预测的概率远离0和1。在某些情况下,解释变量或它们的组合允许完美地预测因变量。在这种情况下,参数要么没有被识别,要么转到正负无穷大。实际参数估计取决于优化的收敛准则。Statsmodels Logit检测到一些这样的情况,然后引发并完善分离异常,但它并没有检测到所有部分分离的情况。
Logit或GLM二项式属于单参数线性指数族。在这种情况下,参数估计仅依赖于指定的均值函数和隐含方差。它不要求正确指定似然函数。因此,即使似然函数对给定数据集不正确,也有可能得到良好(一致)的估计。在这种情况下,该解是一个准极大似然估计,但对数似然值无效。
这可能会导致收敛性和数值稳定性方面的结果取决于如何处理边缘或极端情况的计算细节。Statsmodels正在剪裁这些值,以使它们在某些情况下远离边界,但还不是在任何地方。
困难在于找出如何处理数值问题,避免在底层模型不适合或与数据不兼容时返回“一些”数字而不警告用户。
在这种情况下,也许llf = -inf是“正确”的答案,而任何有限数都只是对-in的近似。也许这只是一个数值问题,因为函数是以双精度的方式实现的。

09-27 03:04