我用sklearn.linear_model做了一个多元回归。线性回归并得到这样的回归系数:
import numpy as np
from sklearn import linear_model
clf = linear_model.LinearRegression()
TST = np.vstack([x1,x2,x3,x4])
TST = TST.transpose()
clf.fit (TST,y)
clf.coef_
现在,我需要这些相同系数的标准误差。我该怎么做?多谢。
根据这个统计问题和维基百科,我最好的猜测是:
MSE = np.mean((y - clf.predict(TST).T)**2)
var_est = MSE * np.diag(np.linalg.pinv(np.dot(TST.T,TST)))
SE_est = np.sqrt(var_est)
但是,我的线性代数和统计数据都非常差,所以我可能会错过一些重要的东西。另一种选择可能是引导方差估计值。
MSE = np.mean((y - clf.predict(TST).T)**2)
var_est = MSE * np.diag(np.linalg.pinv(np.dot(TST.T,TST)))
SE_est = np.sqrt(var_est)
我想这个答案并不完全正确。特别是,如果我没记错的话,根据您的代码,sklearn 正在添加常量项以默认计算您的系数。
然后,您需要在矩阵 TST 中包含 1 列。然后,代码是正确的,它将为您提供一个包含所有 SE 的数组
这些代码已经过数据测试。他们是正确的。
找到每个数据集的 X 矩阵,n 是数据集的长度,m 是变量数
X, n, m=arrays(data)
y=***.reshape((n,1))
linear = linear_model.LinearRegression()
linear.fit(X, y , n_jobs=-1) ## delete n_jobs=-1, if it's one variable only.
和平方
s=np.sum((linear.predict(X) - y) ** 2)/(n-(m-1)-1)
标准差,方差-协方差矩阵对角线的平方根(σ向量分解)
sd_alpha=np.sqrt(s*(np.diag(np.linalg.pinv(np.dot(X.T,X)))))
(使用 t 统计量,linear.intercept_一个变量)
t_stat_alpha=linear.intercept_[0]/sd_alpha[0] #( use linear.intercept_ for one variable_
我发现接受的答案有一些数学故障,总共需要超出修改帖子的建议礼仪的编辑。因此,这里有一个解决方案来计算通过线性模型获得的系数的标准误差估计值(使用此处建议的无偏估计值):
# preparation
X = np.concatenate((np.ones(TST.shape[0], 1)), TST), axis=1)
y_hat = clf.predict(TST).T
m, n = X.shape
# computation
MSE = np.sum((y_hat - y)**2)/(m - n)
coef_var_est = MSE * np.diag(np.linalg.pinv(np.dot(X.T,X)))
coef_SE_est = np.sqrt(var_est)
请注意,我们必须添加一列 1 来TST
因为原始帖子以适合截距术语的方式使用了linear_model.LinearRegression
。此外,我们需要像方差分析一样计算均方误差(MSE)。也就是说,我们需要将误差平方和 (SSE) 除以误差的自由度,即 df_error = df_observations - df_features。
生成的数组coef_SE_est
包含截距的标准误差估计值和所有其他系数(以coef_SE_est[0]
和coef_SE_est[1:]
为单位)。要打印出来,您可以使用
print('intercept: coef={:.4f} / std_err={:.4f}'.format(clf.intercept_[0], coef_SE_est[0]))
for i, coef in enumerate(clf.coef_[0,:]):
print('x{}: coef={:.4f} / std_err={:.4f}'.format(i+1, coef, coef_SE_est[i+1]))
文档中的示例显示了如何获取均方误差和解释方差分数:
# Create linear regression object
regr = linear_model.LinearRegression()
# Train the model using the training sets
regr.fit(diabetes_X_train, diabetes_y_train)
# The coefficients
print('Coefficients: n', regr.coef_)
# The mean square error
print("Residual sum of squares: %.2f"
% np.mean((regr.predict(diabetes_X_test) - diabetes_y_test) ** 2))
# Explained variance score: 1 is perfect prediction
print('Variance score: %.2f' % regr.score(diabetes_X_test, diabetes_y_test))
这是否涵盖您的需求?