`statsmodels`与`sklearn`中的Logit估计器

Question

我觉得这应该是一个功能，而不是一个错误，但我想知道有没有办法让 sklearn 和 statsmodels 的逻辑回归估计结果一致。这里有一个很简单的例子：

import numpy as np
import statsmodels.formula.api as sm
from sklearn.linear_model import LogisticRegression

np.random.seed(123)

n = 100
y = np.random.random_integers(0, 1, n)
x = np.random.random((n, 2))
# Constant term
x[:, 0] = 1.

使用 statsmodels 得到的估计结果：

sm_lgt = sm.Logit(y, x).fit()
    Optimization terminated successfully.
             Current function value: 0.675320
             Iterations 4
print sm_lgt.params
    [ 0.38442   -1.1429183]

而使用 sklearn 得到的估计结果：

sk_lgt = LogisticRegression(fit_intercept=False).fit(x, y)
print sk_lgt.coef_
    [[ 0.16546794 -0.72637982]]

我认为这和 sklearn 的实现有关，因为它使用了一种正则化的方法。有没有选项可以像 statsmodels 那样估计一个简单的逻辑回归（这样速度会快很多，而且扩展性也更好）？另外，sklearn 是否提供推断（标准误差）或边际效应的计算？

Answer 1

另外补充一下，我之前在处理矩阵时遇到了一些结果差异，特别是当我的矩阵是共线的时候。显然，这意味着我们需要一些额外的预处理步骤来获得可靠的结果，但我还是想弄清楚为什么在使用sklearn时得到了结果，而在statsmodels中却出错了。

简单来说：在statsmodels中调用fit时，如果设置solver='bfgs'，那么即使在变量共线的情况下，得到的结果也几乎和sklearn模型一样（当然前提是要注意，statsmodels的默认设置是不包含截距，而sklearn的默认设置是包含截距的）。

下面是一个例子（改编自一个关于OLS的类似问题）：

import numpy as np
import statsmodels.api as sm
from sklearn.linear_model import LogisticRegression

np.random.seed = 237
num_samples=1000
X=np.random.random((num_samples, 2))
X[:, 1] = 2*X[:, 0]
X_sm = sm.add_constant(X)

beta = [1, -2, .5]
error = np.random.random(num_samples)
y = np.round(1/(1+np.exp( -(np.dot(X_sm, beta)) + error   )))  # y = 1/(1+exp(-beta*x))

lr = LogisticRegression(C=1e9).fit(X, y)

print "sklearn:"
print lr.intercept_
print lr.coef_

print "statsmodels:"
print sm.Logit(y, X_sm).fit(method='bfgs').params  # method='nm' or default method errors out

（顺便说一下，如果有人对这两种求解器的数学原理和结果的可靠性有意见，我很想听听！我觉得有趣的是，sklearn对此甚至没有发出警告……）

Answer 2

有没有办法像在 statsmodels 中那样估计一个简单的logit模型？

你可以把 C（反向正则化强度）这个参数设置成一个非常大的数，只要这个数是有限的：

>>> sk_lgt = LogisticRegression(fit_intercept=False, C=1e9).fit(x, y)
>>> print(sk_lgt.coef_)
[[ 0.38440594 -1.14287175]]

不过，关闭正则化是不可能的，因为底层的求解器Liblinear不支持这个功能。

另外，sklearn 有没有提供推断（标准误差）或者边际效应的功能？

没有。目前有一个提议要添加这个功能，但还没有在主代码库中实现。