Погрешности данных и весов

Вероятностная формулировка ошибки

Мы уже узнали, что данные практически никогда не бывают «чистыми» - всегда будет какая-то доля ошибки (погрешность) в исходной выборке для регрессии \(\left\{\left(x_i, y_i\right)\right\}_{i=1}^N\). Хотим как-то описать эту ситуацию на языке математики.

Самым естественным способом видется описание элемента выборки как случайной величины с какой-то плотностью распределения \((x_i, y_i) \sim P(x, y)\). Т.е. для всей выборки с учётом предположения о независимости отдельных точек:

\[ P\left(x_1, \ldots, x_N, y_1, \ldots, y_N\right)=\prod_{i=1}^N P\left(x_i, y_i\right) \]

Таким образом мы можем учесть случайную погрешность, объясняющую несовпадение нашей модели и исходных данных. В частности, это может выглядеть следующим образом в предположении (гипотезе) о том, что наша модель верно описывает исходную зависимость:

\[ y_i=X_{i \alpha} \mathrm{w}_\alpha+\epsilon_i, \quad \epsilon_i \sim \mathcal{N}\left(0, \sigma^2\right) \]

Здесь \(\epsilon_i\) - это случайный гауссов шум, объясняющий, почему исходные точки не лежат точно на нашей модели. Отметим, что \(\sigma\) нам неизвестна, но одинакова для всех точек выборки (хотя это может быть и не так).

Также ещё есть систематическая погрешность данных (она, как правило, описывается крестами ошибок). Даже если все ваши точки идеально легли на модель (нулевая случайная погрешность), сами точки были измерены с какой-то точностью, что также внесёт свой вклад в погрешность весов модели.

\[ \widetilde{y}_i \sim \mathcal{N}\left(y_i, s^2\right) \]

Наша цель - понять вклад обоих видов погрешностей данных в итоговую погрешность весов модели.

Учёт систематической погрешности

Воспользуемся процедурой усреднения по ансамблю приборов. Наши данные \((x_i, y_i)\) знаем с какой-то точностью, обусловленной «прибором», который их измеряет/хранит/воводит. Если бы мы взяли эти данные из другого «прибора», значения бы поменялись на какую-то одинаковую величину (систематическую погрешность данных). Для простоты, рассмотрим только \(y_i\), т.е. систематическую погрешность \(x_i\) примем за ноль. Максимально абстрагируясь, истинные данные для описания - случайные величины на вероятностном пространстве приборов:

\[ \widetilde{y}_i \sim \mathcal{N}\left(y_i, s^2\right) \]

где мы приняли наши имеющиеся данные за среднее, а саму зависимость положили гауссовой со дисперсией \(s^2\).

Иными словами, плотность распределения истинных игреков равна:

\[ p\left(\tilde{y}_i\right)=\frac{1}{\sqrt{2 \pi s^2}} \exp \left(-\frac{\left(\tilde{y}_i-y_i\right)^2}{2 s^2}\right) \]

А тогда плотность распределения всей выборки (то есть «всех N игреков») в случае их независимости (т.е. случайные величины \(\tilde{y}_i\) независимы между собой):

\[ P(\tilde{\mathbf{y}})=\prod_{i=1}^l p\left[\tilde{y}_i\right] \]

где мы просто перемножили все вероятности.

А в случае их зависимости (скоррелированности разных точек из выборки между собой = недиагональная ковариационная матрица):

\[ P(\tilde{\mathbf{y}})= \frac{1}{(2 \pi)^{N / 2} \operatorname{det} A} \exp \left(-\frac{(\tilde{\mathbf{y}}-\mathbf{y})^T A^{-1} (\tilde{\mathbf{y}}-\mathbf{y})}{2}\right) \]

\[ \text{cov} \left( \tilde{y}_i \tilde{y}_j \right) \stackrel{\text{def}}{=} \left\langle\tilde{y}_i \tilde{y}_j\right\rangle = A_{ij} \]

где \(A\) - матрица ковариации точек выборки. Напоминаем, что мы рассматриваем в этом пункте только систематическую погрешность.

Хотим получить погрешность весов. Для этого нужно взять корень из дисперсии \(\tilde{\mathbf{w}}_i\) - случайных величин, получающихся по уже знакомой формуле:

\[ \tilde{\mathbf{w}}=\underbrace{\left(X^{\top} X\right)^{-1} X^{\top}}_Q \cdot \mathbf{y} \]

\[ \tilde{\mathrm{w}}_\alpha=Q_{\alpha k} \cdot \mathbf{y}_k \]

А тогда их дисперсия буквально зависит от корреляторов игреков, которые мы уже знаем (нет суммирования по \(\alpha\)):

\[ \sigma_{\mathrm{w}_\alpha}^2\stackrel{\text{def}}{=}\left\langle\tilde{\mathrm{w}}_\alpha \tilde{\mathrm{w}}_\alpha\right\rangle=Q_{\alpha i} \cdot Q_{\alpha j} \cdot \left\langle\tilde{y}_i \tilde{y}_j\right\rangle = \left(QA Q^{\top}\right)_{\alpha \alpha} \]

В случае диагональной A (независимые признаки с разными дисперсиями \(s_i^2\)) с учётом \(QQ^\top = \left(X^{\top} X\right)^{-1}\):

\[ \sigma_{\mathrm{w}_\alpha}^2 = \left(X^{\top} X\right)^{-1}_{\alpha \alpha} \cdot s_i^2 \]

Ну и в случае скалярной A (независимые признаки с одинаковой дисперсией \(s^2\)):

\[ \sigma_{\mathrm{w}_\alpha}^2 = \left(X^{\top} X\right)^{-1}_{\alpha \alpha} \cdot s^2 \]

Итого, систематическая погрешность весов при одинаковой систематической ошибке:

\[ \Delta_{\text{сист}} \mathrm{w}_\alpha=\sqrt{\left(X^{T} X\right)^{-1}_{ \alpha \alpha }\cdot s^2} \]

Примечание. Выше нигде нет суммирования по \(\alpha\)! Просто берём диагональный элемент.

Учёт случайной погрешности

Гауссова вероятностная модель

Пусть нам дан набор точек \(\left\{\left(x_i, y_i\right)\right\}_{i=1}^N\) и пусть их систематическая (приборная) погрешность равна нулю. Опять-таки надо найти наилучшую аппроксимирующую модель.

Как и раньше, будем решать эту задачу линейной регрессией. Объясним несовпадение истинных данных \(y_i\) и наших предсказаний \(\check{y}_i\) некоторым случайным шумом

\[ y_i=X_{i \alpha} \mathrm{w}_\alpha+\epsilon_i, \quad \epsilon_i \sim \mathcal{N}\left(0, \sigma^2\right) \]

Здесь \(\epsilon_i\) - это случайный гауссов шум, объясняющий, почему исходные точки не лежат точно на нашей модели. Отметим, что \(\sigma\) нам неизвестна, но предполагаем её одинаковость для всех точек выборки.

Принцип максимизации правдоподобия (likelyhood)

В отличие от эмпирически выбранной функции потерь, будем пользоваться принципом максимизации правдоподобия (likelyhood). Он заключается в максимизации вероятности пронаблюдать истинные \(y_i\) в нашей нашей вероятностной модели (при условии, что выполнена гипотеза истинности модели для фиттируемых данных).

Иными словами, нам надо подобрать такие \(\mathbf{w}\), чтобы вероятность \(P_\mathbf{w}(\mathbf{y})\) нашей модели \(\check{y}_i=X_{i \alpha} \mathrm{w}_\alpha\) с учётом шума \(\epsilon_i\) выдать истинные \(y_i\) была максимальной. При этом выборка фиксирована. Эта вероятность называется правдоподобием данных.

\[P_\mathbf{w}(\mathbf{y}) \rightarrow \max _\mathbf{w}\]

Т.к. вероятность независимость событий есть произведение их вероятностей, суммарная верность будет большим произведением, что не очень удобно. Поэтому переходят к эквивалентной задаче, беря логарифм от этой вероятности (это эквивалентный шаг, т.к. максимизация \(\log P_\mathbf{w}(\mathbf{y})\) равносильна максимизации \(P_\mathbf{w}(\mathbf{y})\)). Для удобства, эту штуку также называют правдоподобием.

\[\log P_\mathbf{w}(y) \]

Распишем подробнее,

\[ P_\mathbf{w}(\mathbf{y})=\prod_{i=1}^N \mathcal{N}\left(0, \sigma^2\right)\left[\epsilon_i\right]=\prod_{i=1}^N \frac{1}{\sqrt{2 \pi \sigma^2}} \exp \left(-\frac{\left(y_i-X_{i \alpha} \mathrm{w}_\alpha\right)^2}{2 \sigma^2}\right) \]

Логарифмируем,

\[ \log P_\mathbf{w}(\mathbf{y})=\sum_{i=1}^N\left\{-\frac{1}{2} \log 2 \pi \sigma^2-\frac{\left(y_i-X_{i\alpha} \mathrm{w}_\alpha\right)^2}{2 \sigma^2}\right\} \rightarrow \max _\mathbf{w} \]

Заметим, что первые слагаемые от \(\mathbf{w}\) не зависят, а вторые слагаемые в этой сумме образуют сумму квадратов отклонений, т.е. они образуют функцию потерь в методе наименьших квадратов.

Получается, что МНК эквивалентен принципу максимизации правдоподобия с гауссовой вероятностной моделью. Это и есть вероятностный смысл МНК.

Отметим, что от \(\sigma\) решение \(\mathbf{w}\) не зависит, как и в МНК - главное, чтобы она была одинаковой. Тем не менее, эту \(\sigma\) можно оценить из

\[ \frac{\partial \log P_\mathbf{w}(g)}{\partial \sigma^2}=-\frac{N}{\sigma^2}+\frac{1}{2 \sigma^4} \cdot R S S=0 \Rightarrow {\sigma}^2 =\frac{R S S}{N} = MSE \]

где введена полная сумма квадратов отклонений \(RSS \stackrel{\text{def}}{=} MSE \cdot N\).

Отметим, что есть более точные оценки для общей задачи с количеством параметров \(m\) (т.е. искомых весов \(w_i\) коих также \(m\) штук) оценка на случайную погрешность данных будет

\[ \sigma = \frac{RSS}{N-m} \]

Остюда, при \(N \gg m\) можно спокойно оценивать как \(\sigma^2 \approx \frac{R S S}{N} = MSE\).

Повторяя шаги из прошлого пункта, можно написать случайную погрешность весов для независимого шума в гауссовой вероятностной модели:

\[ \Delta_{\text{случ}} \mathrm{w}_\alpha=\sqrt{\left(X^{T} X\right)^{-1}_{ \alpha \alpha }\cdot MSE} \]

Итого суммарная погрешность весов в самом простом случае («когда всё хорошо» на жаргоне) берётся из корня суммы квадратов, т.к. работаем с вероятностями, а не с максимальными оценками:

\[ \Delta \mathrm{w}_\alpha=\sqrt{\left(X^{T} X\right)_{\alpha \alpha}^{-1} \cdot\left(s^2+MSE\right)} \]

где \(s\) - систематическая погрешность данных, MSE - ошибка модели.

# Визуализация доверительного интервала w_i +- \Delta w_i
# Для зашумлённых данных построим предсказания:
# 1) Модели со средними весами жирным
# 2) Нескольких моделей с весами из доверительного интервала прозрачным

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error

def create_X(x, m):
    X = np.ones(x.shape) # Сразу же добавим intercept
    for i in range(1, m + 1):
        X = np.concatenate([X, x**i], axis = 1)
    return X

def predict(x, w):
    X_predict = create_X(x, len(w) - 1)
    return X_predict @ w

m = 10

xs = np.arange(0, 1, 0.01).reshape((-1, 1))
noise = np.random.normal(0., 0.05, size = xs.shape[0]).reshape((-1, 1))
ys = np.exp(-5*xs**2)*np.sin(20*xs) + noise

X_train = create_X(xs, m)
w_middle = np.linalg.inv(X_train.T @ X_train) @ X_train.T @ ys

delta_w = np.diagonal(X_train.T @ X_train) * mean_squared_error(ys, predict(xs, w_middle))

plt.figure(figsize=(10, 5))

plt.ylim((-1, 2))


x_vis = np.arange(0, 1, 0.001).reshape((-1, 1))

for i in range(1000):
    w = np.random.normal(loc = w_middle ,
                          scale = delta_w)
    plt.plot(x_vis, predict(x_vis, w), lw = 0.75, c = "b", alpha = 0.01)
 

plt.scatter(xs, ys, c="r", zorder = 3)
plt.plot(x_vis, predict(x_vis, w_middle), lw = 3, c = "g", zorder = 4)

[<matplotlib.lines.Line2D at 0x7f3105218e80>]