Python scikit learn线性模型参数标准误差

Question

我正在使用sklearn库，具体来说是linear_model模块。在进行简单线性拟合后，如下所示：

import pandas as pd
import numpy as np
from sklearn import linear_model
randn = np.random.randn

X = pd.DataFrame(randn(10,3), columns=['X1','X2','X3'])
y = pd.DataFrame(randn(10,1), columns=['Y'])        

model = linear_model.LinearRegression()
model.fit(X=X, y=y)

我知道可以通过coef_和intercept_访问系数和截距，并且预测也很简单。我想访问这个简单模型参数的方差协方差矩阵和这些参数的标准误差。我熟悉R和vcov()函数，看起来scipy.optimize也有一些功能(在Python中使用optimize.leastsq方法获取拟合参数的标准误差)——sklearn是否具有访问这些统计数据的功能？感谢任何关于此的帮助。

-Ryan

- Ryan

Answer 1

tl;dr

虽然scikit-learn不能计算，但是您可以使用一些线性代数手动计算。下面是您的示例代码。

此外，这里有一个带有此代码的jupyter笔记本：https://gist.github.com/grisaitis/cf481034bb413a14d3ea851dab201d31

what and why

您的估计标准误差只是您估计方差的平方根。您的估计方差是多少？如果假设您的模型具有高斯误差，则为：

Var(beta_hat) = inverse(X.T @ X) * sigma_squared_hat

然后，beta_hat[i]的标准误差为Var(beta_hat)[i, i] ** 0.5。

您只需要计算sigma_squared_hat。这是您模型的高斯误差的估计值。这不是先验已知的，但可以用残差的样本方差进行估计。

此外，您需要在数据矩阵中添加拦截项。Scikit-learn使用LinearRegression类自动执行此操作。因此，如果要手动计算，则需要将其添加到X矩阵或数据帧中。

how

从您的代码开始，

show your scikit-learn results

print(model.intercept_)
print(model.coef_)

[-0.28671532]
[[ 0.17501115 -0.6928708   0.22336584]]

用线性代数实现这个问题

N = len(X)
p = len(X.columns) + 1  # plus one because LinearRegression adds an intercept term

X_with_intercept = np.empty(shape=(N, p), dtype=np.float)
X_with_intercept[:, 0] = 1
X_with_intercept[:, 1:p] = X.values

beta_hat = np.linalg.inv(X_with_intercept.T @ X_with_intercept) @ X_with_intercept.T @ y.values
print(beta_hat)

[[-0.28671532]
 [ 0.17501115]
 [-0.6928708 ]
 [ 0.22336584]]

计算参数估计的标准误差

y_hat = model.predict(X)
residuals = y.values - y_hat
residual_sum_of_squares = residuals.T @ residuals
sigma_squared_hat = residual_sum_of_squares[0, 0] / (N - p)
var_beta_hat = np.linalg.inv(X_with_intercept.T @ X_with_intercept) * sigma_squared_hat
for p_ in range(p):
    standard_error = var_beta_hat[p_, p_] ** 0.5
    print(f"SE(beta_hat[{p_}]): {standard_error}")

SE(beta_hat[0]): 0.2468580488280805
SE(beta_hat[1]): 0.2965501221823944
SE(beta_hat[2]): 0.3518847753610169
SE(beta_hat[3]): 0.3250760291745124

使用statsmodels进行确认

import statsmodels.api as sm
ols = sm.OLS(y.values, X_with_intercept)
ols_result = ols.fit()
ols_result.summary()

...
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const         -0.2867      0.247     -1.161      0.290      -0.891       0.317
x1             0.1750      0.297      0.590      0.577      -0.551       0.901
x2            -0.6929      0.352     -1.969      0.096      -1.554       0.168
x3             0.2234      0.325      0.687      0.518      -0.572       1.019
==============================================================================

太好了，完成了！

- william_grisaitis

---

太棒了。非常感谢！ - TravelTrader

我的数据集在 sigma_squared_hat = residual_sum_of_squares [0，0] /（N-p）处报错无效的标量变量索引。residual_sum_of_squares计算结果为numpy.float64。我漏掉了什么？ - Bharat

@Bharat，你生成“残差平方和”代码是什么？ - william_grisaitis

那么，当您使用弹性网来收缩系数时呢... - thistleknot

@thistleknot 对于弹性网络等其他非MLE的事情，我认为你唯一的选择是引导系数的标准误差，但我可能错了。 - william_grisaitis

Answer 2

不，scikit-learn没有内置的误差估计用于推理。但Statsmodels有。

import statsmodels.api as sm
ols = sm.OLS(y, X)
ols_result = ols.fit()
# Now you have at your disposition several error estimates, e.g.
ols_result.HC0_se
# and covariance estimates
ols_result.cov_HC0

查看文档

- eickenberg

---

有没有办法使用scikit回归模型中的任何数据计算标准误差？我知道statsmodels提供这些数据，但我需要l2惩罚项，而statsmodels没有。 - TheDude

据我所知没有。对于L2惩罚项和n > p，我猜你可以写出公式。对于n < p，这其实是非常不平凡的，只有最近才开始有人着手解决这个问题。 - eickenberg

这并没有直接回答问题，但是对于预测误差，你可以像这里所述获取均方误差，这是预测标准误差的一步。 - ryanwc

有关 @eickenberg 的答案的更详细版本，请参见：https://dev59.com/MlwZ5IYBdhLWcg3wYfdu - Ambareesh

Answer 3

每个预测器列都是随机的相同格式。因此，这就像运行三个模拟：

import pandas as pd
import numpy as np
from sklearn import linear_model
randn = np.random.randn

X = pd.DataFrame(randn(10,1))
y = pd.DataFrame(randn(10,1)) 
model = linear_model.LinearRegression()
model.fit(X=X, y=y)
y_pred = model.predict(X)
print(y)
print(y_pred)
residuals = y - y_pred
residuals['c'] = residuals.iloc[:, 0]**2
sq = residuals['c']
print(sq)
standard_error = (sum(sq)/(10-2))**0.5
print(standard_error)

- Candace Agonafir