Глава 5. РЕГРЕССИЯ
5.1. Гипотеза линейной регрессии
\[ \hat{y} = X\boldsymbol{\beta} + \epsilon \quad \text{(или } \mathbf{y} = X\boldsymbol{\beta} + \boldsymbol{\epsilon}) \] (где \(\hat{y}\) или \(\mathbf{y}\) - предсказанные или истинные значения, X - матрица признаков (с добавленным столбцом единиц для свободного члена), \(\boldsymbol{\beta}\) - вектор коэффициентов, \(\epsilon\) или \(\boldsymbol{\epsilon}\) - ошибка/шум)Объяснение: Гипотеза линейной регрессии предполагает, что целевая переменная \(y\) является линейной комбинацией признаков X, коэффициентов \(\boldsymbol{\beta}\) и члена ошибки \(\epsilon\).
Пример: Для \(y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \epsilon\), предсказать \(y\) как линейную функцию \(x_1\) и \(x_2\). Если \(\boldsymbol{\beta} = \begin{bmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \end{bmatrix} = \begin{bmatrix} 1 \\ 2 \\ 3 \end{bmatrix}\) и \(X = \begin{bmatrix} 1 & x_{11} & x_{12} \\ 1 & x_{21} & x_{22} \end{bmatrix}\), то \(\hat{\mathbf{y}} = X\boldsymbol{\beta}\).
Реализация (предсказание):
import numpy as np
# X - матрица признаков (n_samples, n_features)
# beta - вектор коэффициентов (n_features,) или (n_features, 1)
# Предположим, X уже включает столбец единиц, если нужен свободный член
# Пример данных
X = np.array([[1, 1, 2], [1, 3, 4]]) # Добавлен столбец единиц
beta = np.array([1, 2, 3]) # beta_0=1, beta_1=2, beta_2=3
y_pred = X @ beta # Матричное умножение
print(f"Predictions: {y_pred}") # [1*1+1*2+2*3, 1*1+3*2+4*3] = [1+2+6, 1+6+12] = [9, 19]
5.2. Метод наименьших квадратов (МНК / OLS - Ordinary Least Squares)
\[ \boldsymbol{\hat{\beta}} = (X^T X)^{-1} X^T \mathbf{y} \] (Решение нормальных уравнений)Объяснение: МНК находит вектор коэффициентов \(\boldsymbol{\beta}\), который минимизирует сумму квадратов остатков (разниц) между предсказанными и фактическими значениями.
Пример: Для \(X = \begin{bmatrix} 1 & 2 \\ 3 & 4 \end{bmatrix}\) (предположим, это признаки без столбца единиц) и \(\mathbf{y} = \begin{bmatrix} 5 \\ 11 \end{bmatrix}\), вычислить \(\boldsymbol{\beta}\).
Реализация:
import numpy as np
# X - матрица признаков (n_samples, n_features)
# y - вектор целевой переменной (n_samples,)
X = np.array([[1, 2], [3, 4]])
y = np.array([5, 11])
try:
# Вычисляем (X^T X)^(-1)
XTX_inv = np.linalg.inv(X.T @ X)
# Вычисляем X^T y
XTy = X.T @ y
# Находим beta
beta = XTX_inv @ XTy
print(f"OLS coefficients beta: {beta}")
# X^T X = [[1,3],[2,4]] @ [[1,2],[3,4]] = [[10, 14],[14, 20]]
# det(X^T X) = 200 - 196 = 4
# (X^T X)^-1 = 1/4 * [[20, -14],[-14, 10]] = [[5, -3.5],[-3.5, 2.5]]
# X^T y = [[1,3],[2,4]] @ [5, 11] = [5+33, 10+44] = [38, 54]
# beta = [[5, -3.5],[-3.5, 2.5]] @ [38, 54] = [5*38-3.5*54, -3.5*38+2.5*54]
# = [190 - 189, -133 + 135] = [1, 2]
# Значит beta = [1, 2]
except np.linalg.LinAlgError:
print("Матрица X^T X вырождена. Используйте псевдообратную или регуляризацию.")
# Решение с использованием псевдообратной матрицы:
beta_pinv = np.linalg.pinv(X.T @ X) @ X.T @ y
print(f"OLS coefficients beta (using pseudo-inverse): {beta_pinv}")
# Или через np.linalg.lstsq, что предпочтительнее численно:
beta_lstsq, residuals, rank, s = np.linalg.lstsq(X, y, rcond=None)
print(f"OLS coefficients beta (using lstsq): {beta_lstsq}")
5.3. Среднеквадратичная ошибка (MSE - Mean Squared Error)
\[ \text{MSE} = \frac{1}{n} \sum_{i=1}^{n} (y_i - \hat{y}_i)^2 \] (где n - количество наблюдений)Объяснение: MSE количественно определяет среднюю квадратичную разницу между фактическими и предсказанными значениями. Это стандартная функция потерь в регрессии.
Пример: Для \(y = [1, 2, 3]\) и \(\hat{y} = [1.1, 1.9, 3.2]\), вычислить MSE. MSE \( = \frac{1}{3} [(1-1.1)^2 + (2-1.9)^2 + (3-3.2)^2] = \frac{1}{3} [(-0.1)^2 + (0.1)^2 + (-0.2)^2] = \frac{1}{3} [0.01 + 0.01 + 0.04] = \frac{0.06}{3} = 0.02\).
Реализация:
import numpy as np
y_true = np.array([1, 2, 3])
y_pred = np.array([1.1, 1.9, 3.2])
mse = np.mean((y_true - y_pred)**2)
print(f"MSE: {mse}") # Должно быть 0.02
5.4. Градиент потерь MSE
\[ \nabla_{\boldsymbol{\beta}} \text{MSE} = \frac{1}{n} \nabla_{\boldsymbol{\beta}} \|\mathbf{y} - X\boldsymbol{\beta}\|_2^2 = \frac{2}{n} X^T (X\boldsymbol{\beta} - \mathbf{y}) = -\frac{2}{n} X^T (\mathbf{y} - X\boldsymbol{\beta}) \]Объяснение: Градиент MSE по отношению к \(\boldsymbol{\beta}\) используется в градиентных алгоритмах оптимизации, таких как градиентный спуск, для подбора коэффициентов регрессии.
Пример: Вычислить градиент для \(X = \begin{bmatrix} 1 & 2 \\ 3 & 4 \end{bmatrix}\), \(\mathbf{y} = [5, 11]\) и \(\boldsymbol{\beta} = [1, 1]\).
Реализация:
import numpy as np
X = np.array([[1, 2], [3, 4]])
y = np.array([5, 11])
beta = np.array([1, 1]) # Текущие коэффициенты
n = len(y)
# Вычисляем предсказания
y_pred = X @ beta # [1*1+2*1, 3*1+4*1] = [3, 7]
# Вычисляем вектор ошибки
error = y_pred - y # [3-5, 7-11] = [-2, -4]
# Вычисляем градиент
grad = (2/n) * X.T @ error
# grad = (2/2) * [[1, 3], [2, 4]] @ [-2, -4]
# = [[1, 3], [2, 4]] @ [-2, -4]
# = [1*(-2)+3*(-4), 2*(-2)+4*(-4)] = [-2-12, -4-16] = [-14, -20]
print(f"Gradient of MSE loss: {grad}") # Должно быть [-14, -20]
# Вариант с -2/n * X^T(y - X*beta)
error_alt = y - y_pred # [5-3, 11-7] = [2, 4]
grad_alt = (-2/n) * X.T @ error_alt
# grad_alt = (-2/2) * [[1, 3], [2, 4]] @ [2, 4]
# = -1 * [1*2+3*4, 2*2+4*4] = -1 * [14, 20] = [-14, -20]
print(f"Gradient (alternative formula): {grad_alt}")
5.5. Коэффициент детерминации (R²)
\[ R^2 = 1 - \frac{\sum_{i=1}^{n} (y_i - \hat{y}_i)^2}{\sum_{i=1}^{n} (y_i - \bar{y})^2} = 1 - \frac{\text{SSE}}{\text{SST}} \] (SSE - Sum of Squared Errors/Residuals, SST - Total Sum of Squares)Объяснение: R² измеряет долю дисперсии целевой переменной, объясненную моделью. Значение, близкое к 1, указывает на хорошее соответствие модели данным.
Пример: Для \(y = [1, 2, 3]\) и \(\hat{y} = [1.1, 1.9, 3.2]\), вычислить R². \(\bar{y} = (1+2+3)/3 = 2\). SSE \( = (1-1.1)^2 + (2-1.9)^2 + (3-3.2)^2 = 0.01 + 0.01 + 0.04 = 0.06\). SST \( = (1-2)^2 + (2-2)^2 + (3-2)^2 = (-1)^2 + 0^2 + 1^2 = 1 + 0 + 1 = 2\). \(R^2 = 1 - \frac{0.06}{2} = 1 - 0.03 = 0.97\).
Реализация:
import numpy as np
y_true = np.array([1, 2, 3])
y_pred = np.array([1.1, 1.9, 3.2])
# Сумма квадратов ошибок (остатков)
sse = np.sum((y_true - y_pred)**2)
# Общая сумма квадратов
sst = np.sum((y_true - np.mean(y_true))**2)
# Избегаем деления на ноль, если все y_true одинаковы
if sst == 0:
r2 = 1.0 if sse == 0 else 0.0 # Или неопределен
else:
r2 = 1 - (sse / sst)
print(f"R-squared: {r2}") # Должно быть 0.97
5.6. Скорректированный R² (Adjusted R²)
\[ R^2_{\text{adj}} = 1 - \frac{(1 - R^2)(n - 1)}{n - p - 1} \] (где n - количество наблюдений, p - количество предикторов/признаков в модели, не считая свободного члена)Объяснение: Скорректированный R² учитывает количество предикторов \(p\) в модели, штрафуя за переобучение (добавление бесполезных признаков).
Пример: Для \(R^2 = 0.97\), \(n = 50\) наблюдений и \(p = 5\) предикторов, вычислить \(R^2_{\text{adj}}\). \(R^2_{\text{adj}} = 1 - \frac{(1 - 0.97)(50 - 1)}{50 - 5 - 1} = 1 - \frac{0.03 \times 49}{44} = 1 - \frac{1.47}{44} \approx 1 - 0.0334 = 0.9666\).
Реализация:
import numpy as np
r2 = 0.97
n = 50 # Количество наблюдений
p = 5 # Количество предикторов
# Убедимся, что знаменатель не ноль или отрицательный
denominator = n - p - 1
if denominator <= 0:
# Adjusted R^2 не определен или бессмысленен
adjusted_r2 = np.nan
else:
adjusted_r2 = 1 - (1 - r2) * (n - 1) / denominator
print(f"Adjusted R-squared: {adjusted_r2}") # ~ 0.96659
5.7. Средняя абсолютная ошибка (MAE - Mean Absolute Error)
\[ \text{MAE} = \frac{1}{n} \sum_{i=1}^{n} |y_i - \hat{y}_i| \]Объяснение: MAE измеряет среднюю абсолютную величину ошибок предсказания. Она менее чувствительна к выбросам по сравнению с MSE.
Пример: Для \(y = [1, 2, 3]\) и \(\hat{y} = [1.1, 1.9, 3.2]\), вычислить MAE. MAE \( = \frac{1}{3} [|1-1.1| + |2-1.9| + |3-3.2|] = \frac{1}{3} [|-0.1| + |0.1| + |-0.2|] = \frac{1}{3} [0.1 + 0.1 + 0.2] = \frac{0.4}{3} \approx 0.133\).
Реализация:
import numpy as np
y_true = np.array([1, 2, 3])
y_pred = np.array([1.1, 1.9, 3.2])
mae = np.mean(np.abs(y_true - y_pred))
print(f"MAE: {mae}") # ~ 0.1333...
5.8. Взвешенный метод наименьших квадратов (WLS - Weighted Least Squares)
\[ \boldsymbol{\hat{\beta}} = (X^T W X)^{-1} X^T W \mathbf{y} \] (где W - диагональная матрица весов)Объяснение: WLS минимизирует сумму взвешенных квадратов остатков, что позволяет учитывать гетероскедастичность (неодинаковую дисперсию ошибок) в данных.
Пример: Для \(W = \text{diag}([w_1, w_2, \dots, w_n])\), где веса \(w_i\) обратно пропорциональны дисперсии ошибки для i-го наблюдения, вычислить \(\boldsymbol{\beta}\).
Реализация:
import numpy as np
# X - матрица признаков (n_samples, n_features)
# y - вектор целевой переменной (n_samples,)
# weights - вектор весов (n_samples,)
X = np.array([[1, 2], [3, 4]])
y = np.array([5, 11])
weights = np.array([1, 2]) # Пример весов
# Создаем диагональную матрицу W
W = np.diag(weights)
try:
# Вычисляем (X^T W X)^(-1)
XTWX_inv = np.linalg.inv(X.T @ W @ X)
# Вычисляем X^T W y
XTWy = X.T @ W @ y
# Находим beta
beta_wls = XTWX_inv @ XTWy
print(f"WLS coefficients beta: {beta_wls}")
except np.linalg.LinAlgError:
print("Матрица X^T W X вырождена.")
# Можно использовать псевдообратную
beta_wls_pinv = np.linalg.pinv(X.T @ W @ X) @ X.T @ W @ y
print(f"WLS coefficients beta (using pseudo-inverse): {beta_wls_pinv}")
5.9. Гипотеза полиномиальной регрессии
\[ \hat{y} = \beta_0 + \beta_1 x + \beta_2 x^2 + \dots + \beta_n x^n \]Объяснение: Полиномиальная регрессия моделирует зависимость между \(x\) и \(y\) как полином степени \(n\). Она обобщает линейную регрессию для нелинейных зависимостей.
Пример: Подогнать модель \(y = 2x + x^2\).
Реализация (используя scikit-learn):
import numpy as np
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression
# Пример данных
X = np.array([1, 2, 3, 4, 5]).reshape(-1, 1)
y = 2 * X.flatten() + X.flatten()**2 + np.random.randn(5) # y = 2x + x^2 + шум
# Создаем полиномиальные признаки (степень 2)
poly = PolynomialFeatures(degree=2, include_bias=True) # include_bias=True добавит столбец для beta_0
X_poly = poly.fit_transform(X)
# Обучаем линейную регрессию на полиномиальных признаках
model = LinearRegression(fit_intercept=False) # Перехват уже включен в X_poly, если include_bias=True
model.fit(X_poly, y)
# Предсказание
X_new = np.array([2.5]).reshape(-1, 1)
X_new_poly = poly.transform(X_new)
y_pred = model.predict(X_new_poly)
print(f"Coefficients (beta_0, beta_1, beta_2): {model.coef_}") # Должны быть близки к [?, 2, 1]
print(f"Prediction for x=2.5: {y_pred}")
5.10. Нелинейная регрессия
\[ \hat{y} = f(X, \boldsymbol{\beta}) + \epsilon \] (где \(f\) - нелинейная функция параметров \(\boldsymbol{\beta}\))Объяснение: Нелинейная регрессия моделирует зависимости, где целевая переменная является нелинейной функцией параметров.
Пример: Подогнать модель \(y = a e^{bx}\) с помощью оптимизации для нахождения параметров \(a\) и \(b\).
Реализация (используя scipy.optimize.curve_fit):
import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
# Определяем нелинейную модель
def model_func(x, a, b):
"""Пример нелинейной модели: a * exp(b*x)."""
return a * np.exp(b * x)
# Генерируем пример данных
x_data = np.linspace(0, 4, 50)
a_true, b_true = 2.5, 1.3
y_data = model_func(x_data, a_true, b_true) + np.random.normal(size=x_data.shape)*1.5 # Добавляем шум
# Подгоняем модель к данным
try:
params, params_covariance = curve_fit(model_func, x_data, y_data, p0=[2, 1]) # p0 - начальное приближение
a_fit, b_fit = params
print(f"Fitted parameters: a={a_fit:.3f}, b={b_fit:.3f}")
# Визуализация (опционально)
plt.figure(figsize=(6, 4))
plt.scatter(x_data, y_data, label='Data')
plt.plot(x_data, model_func(x_data, a_fit, b_fit), color='red', label=f'Fit: a={a_fit:.2f}, b={b_fit:.2f}')
plt.legend()
plt.xlabel("x")
plt.ylabel("y")
plt.title("Non-linear Regression Fit")
plt.show()
except RuntimeError:
print("Could not find optimal parameters.")
5.11. Оценка максимального правдоподобия (MLE) для регрессии
\[ \boldsymbol{\hat{\beta}}_{\text{MLE}} = \arg \max_{\boldsymbol{\beta}} \prod_{i=1}^{n} p(y_i | X_i, \boldsymbol{\beta}) \quad \Leftrightarrow \quad \arg \min_{\boldsymbol{\beta}} \sum_{i=1}^{n} [-\log p(y_i | X_i, \boldsymbol{\beta})] \] (где \(p(y|X, \beta)\) - условная плотность вероятности y при заданных X и \(\beta\))Объяснение: MLE оценивает параметры, которые максимизируют правдоподобие наблюдения данных при заданной вероятностной модели.
Пример: Оценить \(\boldsymbol{\beta}\), предполагая Гауссов шум: \(y_i \sim \mathcal{N}(X_i \boldsymbol{\beta}, \sigma^2)\). Максимизация правдоподобия (или минимизация отрицательного логарифма правдоподобия) в этом случае эквивалентна минимизации MSE (методу наименьших квадратов).
Реализация (минимизация отрицательного лог. правдоподобия для Гауссова шума):
import numpy as np
from scipy.optimize import minimize
# X - матрица признаков
# y - вектор целевой переменной
X = np.array([[1, 2], [3, 4]])
y = np.array([5, 11])
n_features = X.shape[1]
# Определяем функцию отрицательного логарифма правдоподобия
# Для Гауссова шума с постоянной дисперсией, это пропорционально сумме квадратов ошибок
def neg_log_likelihood_gaussian(beta, X_data, y_data):
"""Отриц. лог. правдоподобия для лин. регрессии с Гауссовым шумом."""
residuals = y_data - X_data @ beta
# Константы (связанные с sigma^2 и 2pi) опускаются, т.к. они не влияют на argmin
return np.sum(residuals**2) # / (2 * sigma**2) + n/2 * log(2 * pi * sigma**2)
# Начальное приближение для beta
initial_beta = np.zeros(n_features)
# Минимизируем функцию
# minimize использует градиентные методы, градиент -negLL пропорционален градиенту MSE
result = minimize(neg_log_likelihood_gaussian, initial_beta, args=(X, y), method='L-BFGS-B')
if result.success:
beta_mle = result.x
print(f"MLE coefficients beta (Gaussian noise): {beta_mle}") # Должно быть [1, 2]
else:
print(f"Optimization failed: {result.message}")
5.12. Минимизация эмпирического риска (ERM - Empirical Risk Minimization)
\[ \hat{\boldsymbol{\theta}} = \arg \min_{\boldsymbol{\theta}} \frac{1}{n} \sum_{i=1}^{n} \ell(y_i, f(X_i, \boldsymbol{\theta})) \] (где \(\ell\) - функция потерь, \(f\) - модель)Объяснение: ERM минимизирует средние потери по обучающим данным для оценки параметров модели.
Пример: Минимизировать потери MSE для линейной регрессии: \(\hat{\boldsymbol{\beta}} = \arg \min_{\boldsymbol{\beta}} \frac{1}{n} \sum (y_i - X_i \boldsymbol{\beta})^2\).
Реализация (вычисление эмпирического риска):
import numpy as np
# theta - параметры модели (например, beta)
# X - матрица признаков
# y - вектор истинных значений
# loss_function - функция потерь l(y_true, y_pred)
# model_predict - функция модели f(X_i, theta), возвращающая предсказание
def empirical_risk(theta, X_data, y_data, loss_function, model_predict):
"""Вычисляет эмпирический риск."""
n = len(y_data)
total_loss = 0
for i in range(n):
# Для некоторых моделей предсказание делается сразу для всех X
# Здесь пример с предсказанием для одного X_i
y_pred_i = model_predict(X_data[i], theta)
total_loss += loss_function(y_data[i], y_pred_i)
return total_loss / n
# Пример использования: Линейная регрессия с MSE
def mse_loss_single(y_true, y_pred):
return (y_true - y_pred)**2
def linear_model_predict(X_i, beta):
return np.dot(X_i, beta) # Предсказание для одного примера X_i
X = np.array([[1, 2], [3, 4]])
y = np.array([5, 11])
beta_test = np.array([1.0, 1.5]) # Пробные параметры
risk = empirical_risk(beta_test, X, y, mse_loss_single, linear_model_predict)
# Предсказания: [1*1+2*1.5, 3*1+4*1.5] = [4, 9]
# Потери: (5-4)^2=1, (11-9)^2=4
# Риск: (1+4)/2 = 2.5
print(f"Empirical risk (MSE) for beta={beta_test}: {risk}")