Линейная регрессия с использованием линейной алгебры в Python - PullRequest
Новая
       44

Линейная регрессия с использованием линейной алгебры в Python

0 голосов
/ 26 апреля 2019

Я интерпретирую эти формулы в Википедии (https://en.wikipedia.org/wiki/Coefficient_of_determination) неправильно в Python? Ниже я попробовал.

ssres

def ss_res(X, y, theta):

    y_diff=[]
    y_pred = X.dot(theta)

    for i in range(0, len(y)):
        y_diff.append((y[i]-y_pred[i])**2)

    return np.sum(y_diff)

вывод выглядит корректно, но цифры немного не совпадают ... как несколько десятичных знаков.

* +1012 *stderror 
def std_error(X, y, theta):


    delta = (1/(len(y)-X.shape[1]+1))*(ss_res(X,y,theta))
    matrix1=matrix_power((X.T.dot(X)),-1)
    thing2=delta*matrix1
    thing3=scipy.linalg.sqrtm(thing2)

    res=np.diag(thing3)
    serr=np.reshape(res, (6, 1))
    return serr

std_error_array=std_error(X,y,theta)

1 Ответ

2 голосов
/ 26 апреля 2019

Вы можете или не можете хотеть +1 в том, что вы называете delta, зависит от того, включает ли ваш X столбец «константы» (т.е. все значения = 1)

В противном случае это выглядит нормально, если немного не пифоново. Я хотел бы написать их больше как: 

import numpy as np
from numpy.linalg import inv
from scipy.linalg import sqrtm

def solve_theta(X, Y):
    return np.linalg.solve(X.T @ X, X.T @ Y)

def ss_res(X, Y, theta):
    res = Y - (X @ theta)
    return np.sum(res ** 2)

def std_error(X, Y, theta):
    nr, rank = X.shape
    resid_df = nr - rank
    residvar = ss_res(X, Y, theta) / resid_df
    var_theta = residvar * inv(X.T @ X)
    return np.diag(sqrtm(var_theta))[:,None]

Примечание: здесь используется оператор умножения матрицы в стиле Python 3.5 @ вместо записи .dot()

Числовая стабильность такого алгоритма не удивительна, вы можете посмотреть на использование SVD или QR-разложений. есть доступное описание того, как бы вы сделали это с SVD:

Джон Мандель (1982) «Использование разложения по сингулярным числам в регрессионном анализе» 10.1080 / 00031305.1982.10482771

Мы можем проверить это, создав несколько фиктивных данных: 

np.random.seed(42)

N = 20
K = 3

true_theta = np.random.randn(K, 1) * 5
X = np.random.randn(N, K)
Y = np.random.randn(N, 1) + X @ true_theta

и запустите указанный выше код: 

theta = solve_theta(X, Y)
sse = std_error(X, Y, theta)

print(np.column_stack((theta, sse)))

, что дает: 

[[ 2.23556391  0.35678574]
 [-0.40643163  0.24751913]
 [ 3.14687637  0.26461827]]

Мы можем проверить это с помощью statsmodels: 

import statsmodels.api as sm

sm.OLS(Y, X).fit().summary()

, что дает: 

                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
x1             2.2356      0.358      6.243      0.000       1.480       2.991
x2            -0.4064      0.248     -1.641      0.119      -0.929       0.116
x3             3.1469      0.266     11.812      0.000       2.585       3.709

что довольно близко.

...