为何用sklearn与statsmodels做Bootstrap得到的系数不同？

使用ISLR的Default数据集，分别用sklearn和statsmodels执行逻辑回归，非Bootstrap时两者系数一致，但通过Bootstrap估计系数和标准误时，sklearn Bootstrap得到的系数与其他结果（非Bootstrap sklearn、非Bootstrap statsmodels、Bootstrap statsmodels）差异显著。

1. Statsmodels非Bootstrap系数结果 #

代码：

endog = (Default["default"] == "Yes").astype(int)
exog = sm.add_constant(Default[["income", "balance"]])
mod = sm.Logit(endog, exog)
res = mod.fit()
res.summary()

输出摘要中的系数部分：

==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const        -11.5405      0.435    -26.544      0.000     -12.393     -10.688
income      2.081e-05   4.99e-06      4.174      0.000     1.1e-05    3.06e-05
balance        0.0056    0.00022     24.835      0.000       0.005       0.006
==============================================================================

核心系数：

• income系数：2.081e-05，标准误：4.99e-06
• balance系数：0.0056，标准误：0.00022

2. Bootstrap实现代码 #

2.1 sklearn的Bootstrap代码 #

首先定义boot_fn：

def boot_fn(data, clf):
    clf.fit(data[["income", "balance"]], data["default"])
    return clf.coef_

Bootstrap循环：

np.random.seed(17)
num_estimates = 1000
boot_estimates = np.empty((num_estimates, 2))

for i in range(num_estimates):
    clf = LogisticRegression(penalty = "none", solver = "lbfgs")
    sample = resample(Default)
    coefs = boot_fn(sample, clf)
    boot_estimates[i, 0] = coefs[0, 0]
    boot_estimates[i, 1] = coefs[0, 1]
    print(f"{i+1} out of {num_estimates}", end="\r")

2.2 statsmodels的Bootstrap代码 #

np.random.seed(17)
num_estimates = 1000
boot_estimates = np.empty((num_estimates, 2))
for i in range(num_estimates):
    sample = resample(Default)
    endog = (sample["default"] == "Yes").astype(int)
    exog = sm.add_constant(sample[["income", "balance"]])
    mod = sm.Logit(endog, exog)
    res = mod.fit(disp = False)
    boot_estimates[i, 0] = res.params["income"]
    boot_estimates[i, 1] = res.params["balance"]
    print(f'{i+1} of {num_estimates}', end='\r')

3. 所有结果对比 #

通过.mean()和.sem()计算系数和标准误，整理如下：

4. 补充：非Bootstrap的sklearn系数结果 #

代码：

X = Default[["balance", "income"]]
y = Default["default"]
clf = LogisticRegression(penalty = "none", solver = "lbfgs")
clf.fit(X, y)
clf.coef_

输出：

array([[5.64710291e-03, 2.08089921e-05]])

可以看到，statsmodels非Bootstrap、sklearn非Bootstrap、statsmodels Bootstrap的系数结果均一致，但sklearn Bootstrap的系数差异明显，请问这是什么原因？

内容的提问来源于stack exchange，提问作者guin0x