Python scikit изучает линейную модель Параметр Стандартная ошибка

Я работаю с 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 имеет некоторые функциональные возможности для этого (получение стандартных ошибок по установленным параметрам с использованием метода optimize.leastsq в python) - имеет ли sklearn какие-либо функции для доступа к этим статистическим данным??

Цените любую помощь по этому поводу.

-Ryan

Ответ 2

ТЛ; др

не с scikit-learn, но вы можете вычислить это вручную с помощью некоторой линейной алгебры. Я делаю это для вашего примера ниже.

также здесь записная книжка Jupyter с этим кодом: https://gist.github.com/grisaitis/cf481034bb413a14d3ea851dab201d31

что и почему

стандартные ошибки ваших оценок - это просто квадратный корень из отклонений ваших оценок. какая разница по твоей оценке? если вы предполагаете, что ваша модель имеет ошибку Гаусса, это:

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 или фрейм данных.

как

Начиная с вашего кода,

покажите свои результаты на уроках

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
==============================================================================

да, готово!

Ответ 1