§2. Линейная регрессионная модель

Предположим, что (с точностью до случайных ошибок) целевая переменная $Y$ есть линейная комбинация $\theta_1X_1 + \dots + \theta_mX_m$ предикторных переменных $X_1, \dots, X_m$ с неизвестными коэффициентами $\theta_1, \dots, \theta_m$.

Измерения $\eta_i$ переменной $Y$ ($i = 1, \dots, n$, где $n \ge m$), отвечающие заданным (не обязательно различным) значениям $x_{i1}, \dots, x_{im}$ предикторных переменных, имеют вид

где $\varepsilon_i$ — случайные ошибки. Вводя для векторов и матриц обозначения $\boldsymbol\eta = (\eta_1, \dots, \eta_n)^T$, $\boldsymbol X = \|x_{il}\|_{n\times m}$, $\boldsymbol\theta = (\theta_1, \dots, \theta_m)^T$, $\boldsymbol\varepsilon = (\varepsilon_1, \dots, \varepsilon_n)^T$, можно записать модель в матричной форме:

Матрицу $\boldsymbol X$ называют матрицей плана эксперимента.

Будем предполагать, что для модели (5) выполняется допущение

Д1. Столбцы $\boldsymbol x_l = (x_{1l}, \dots, x_{nl})^T$, $l = 1, \dots, m$, матрицы $\boldsymbol X$ линейно независимы. Иными словами, ввиду выполнения неравенства $n \ge m$ матрица $\boldsymbol X$ имеет ранг $m$ (см. П10).

Смысл по-пацански: «линейная» — это про параметры $\theta$, а не про $X$. Хочешь параболу — бери столбцы $1, u, u^2$; хочешь синусоиду — бери $\sin\pi u$. Главное, чтобы неизвестные коэффициенты входили «множителями», а не сидели внутри экспоненты.

Пример 2. Подгонка полинома. Рассмотрим $\boldsymbol x_1 = (1, \dots, 1)^T$ и $\boldsymbol x_{l+1} = (u_1^l, \dots, u_n^l)^T$, $l = 1, \dots, m - 1$. Здесь $u_1 < \dots < u_n$ — так называемые «узлы», в которых вычисляются значения многочлена

и «зашумляются» случайными ошибками $\varepsilon_i$. Матрица $\boldsymbol X$ имеет вид

Определитель подматрицы $\tilde{\boldsymbol X}$, образованной первыми $m$ строками $\boldsymbol X$ ($m \le n$), является известным из линейной алгебры определителем Вандермонда: $\det\tilde{\boldsymbol X} = \prod_{1\le i<j\le m}(u_j - u_i)$. Он отличен от нуля, поскольку «узлы» различны. Поэтому ранг подматрицы $\tilde{\boldsymbol X}$ (и матрицы $\boldsymbol X$) равен $m$, а столбцы $\boldsymbol x_1, \dots, \boldsymbol x_m$ линейно независимы.

В дальнейшем при изучении статистических свойств оценок параметров $\theta_1, \dots, \theta_m$ (см. §3) будем считать выполняющимся также допущение

Д2. Случайные величины $\varepsilon_1, \dots, \varepsilon_n$ одинаково распределены с $\mathbf M\varepsilon_i = 0$, $\mathbf D\varepsilon_i = \sigma^2$ (параметр $0 < \sigma < \infty$ также неизвестен) и некоррелированы: $\mathbf M\varepsilon_i\varepsilon_j = 0$ при $i \ne j$.

МНК как проекция

Оценим параметры $\theta_1, \dots, \theta_m$ методом наименьших квадратов, минимизируя по $\boldsymbol\theta$ функцию $F(\boldsymbol\theta) = \sum_{i=1}^n(\eta_i - \theta_1x_{i1} - \dots - \theta_mx_{im})^2$. Точка её минимума $\hat{\boldsymbol\theta}$ называется МНК-оценкой, вектор $\hat{\boldsymbol\delta}$ с компонентами $\hat\delta_i = \eta_i - \hat\theta_1x_{i1} - \dots - \hat\theta_mx_{im}$, $i = 1, \dots, n$, — вектором остатков, а значение в точке минимума $RSS = F(\hat{\boldsymbol\theta})$ — остаточной суммой квадратов (в обозначении использованы первые буквы английского термина Residual Sum of Squares).

Для нахождения оценки $\hat{\boldsymbol\theta}$ интерпретируем задачу минимизации функции $F(\boldsymbol\theta)$ в терминах пространства $\mathbb R^n$. Рассмотрим в $\mathbb R^n$ подпространство $L(\boldsymbol X)$, порождаемое столбцами $\boldsymbol x_1, \dots, \boldsymbol x_m$ матрицы $\boldsymbol X$. Очевидно, что вектор $\boldsymbol X\boldsymbol\theta = \theta_1\boldsymbol x_1 + \dots + \theta_m\boldsymbol x_m$ пробегает $L(\boldsymbol X)$, когда $\boldsymbol\theta$ пробегает $\mathbb R^m$. Поскольку

видим, что минимизация по $\boldsymbol\theta$ равносильна нахождению в подпространстве $L(\boldsymbol X)$ вектора $\boldsymbol X\hat{\boldsymbol\theta}$, наименее удалённого от $\boldsymbol\eta$. Как известно, таким вектором служит ортогональная проекция $\boldsymbol\eta$ на $L(\boldsymbol X)$. Следовательно, вектор остатков $\hat{\boldsymbol\delta} = \boldsymbol\eta - \boldsymbol X\hat{\boldsymbol\theta}$ должен быть ортогонален векторам $\boldsymbol x_1, \dots, \boldsymbol x_m$, порождающим подпространство $L(\boldsymbol X)$, то есть

Соотношение (6) представляет собой систему линейных уравнений относительно $\hat{\boldsymbol\theta}$ с положительно определённой (задача 3) матрицей $\boldsymbol B = \boldsymbol X^T\boldsymbol X$, которую называют информационной. Решить систему (6) можно, например, методом Холецкого (П10). Так как ввиду положительной определённости матрица $\boldsymbol B$ является невырожденной, для МНК-оценки $\hat{\boldsymbol\theta}$ получаем представление

Смысл по-пацански: вектор наблюдений $\boldsymbol\eta$ живёт в огромном $\mathbb R^n$, а модель умеет порождать только векторы из маленькой «плоскости» $L(\boldsymbol X)$. МНК просто роняет перпендикуляр с $\boldsymbol\eta$ на эту плоскость. Остаток — это перпендикуляр, поэтому он ортогонален всем столбцам. Вся линейная регрессия — это одна картинка из школьной стереометрии.

Пример 3. На рис. 7 у Лагутина точками изображены результаты эксперимента по изучению зависимости между скоростью автомобиля $V$ (в милях/час) и расстоянием $Y$ (в футах), пройденным им после сигнала об остановке. Для каждого отдельного случая результат определяется в основном тремя факторами: скоростью $V$ в момент подачи сигнала, временем реакции $\theta_1$ водителя на этот сигнал и тормозами автомобиля. Автомобиль успеет проехать путь $\theta_1V$ до момента включения водителем тормозов и ещё $\theta_2V^2$ после этого момента, поскольку согласно элементарным физическим законам теоретическое расстояние, пройденное до остановки с момента торможения, пропорционально квадрату скорости (убедитесь!).

По кирпичикам: при постоянном замедлении $w$ скорость падает до нуля за время $V/w$, а путь равен $V^2/(2w)$ — квадрат скорости.

Таким образом, в качестве модели годится $Y = \theta_1V + \theta_2V^2$. Для экспериментальных данных по формуле (7) были подсчитаны значения $\hat\theta_1 = 0{,}76$ и $\hat\theta_2 = 0{,}056$.

Самих данных в книге нет, поэтому сымитируем эксперимент и проверим, что формула (7) восстанавливает коэффициенты:

(Последняя строка: 1 миля/час ≈ 1,467 фута/с, так что $\theta_1 \approx 0{,}76$ фута на милю/час — это примерно полсекунды реакции.)

Ортогональные столбцы

Замечание 3. Удобно иметь дело с матрицей $\boldsymbol X$, столбцы которой ортогональны друг другу. В этом случае формула (7) упрощается: так как матрица $\boldsymbol B$ оказывается диагональной с элементами $|\boldsymbol x_l|^2$ ($l = 1, \dots, m$) на главной диагонали, из представления (7) получаем выражения

$$\hat\theta_l = \boldsymbol x_l^T\boldsymbol\eta/|\boldsymbol x_l|^2, \qquad l = 1, \dots, m. \tag{8}$$

Если исходный набор линейно независимых векторов $\boldsymbol x_1, \dots, \boldsymbol x_m$ таким свойством не обладает, то его можно ортогонализировать с помощью стандартной процедуры Грама — Шмидта, вводя новые векторы

$$\boldsymbol x_l' = \boldsymbol x_l - a_{(l-1)l}\boldsymbol x_{l-1}' - \dots - a_{1l}\boldsymbol x_1', \qquad l = 1, \dots, m, \tag{9}$$

где $a_{il} = \boldsymbol x_l^T\boldsymbol x_i'/|\boldsymbol x_i'|^2$. Выражая $\boldsymbol x_l$ через $\boldsymbol x_1', \dots, \boldsymbol x_l'$, получаем

$$\boldsymbol x_l = \boldsymbol x_l' + a_{(l-1)l}\boldsymbol x_{l-1}' + \dots + a_{1l}\boldsymbol x_1', \qquad l = 1, \dots, m. \tag{10}$$

Соотношения (10) можно записать в матричной форме: $\boldsymbol X = \boldsymbol X'\boldsymbol A$, где $\boldsymbol A$ — верхнетреугольная матрица с единицами на главной диагонали (она, очевидно, обратима). Подставляя это выражение в формулу (5), приходим к модели $\boldsymbol\eta = \boldsymbol X'\boldsymbol\theta' + \boldsymbol\varepsilon$, где $\boldsymbol\theta' = \boldsymbol A\boldsymbol\theta$. Исходный вектор $\boldsymbol\theta$ восстанавливается по $\boldsymbol\theta'$ с помощью формулы

$$\boldsymbol\theta = \boldsymbol A^{-1}\boldsymbol\theta'. \tag{11}$$

Внимание: в книге в формулах после (10) стоит $\boldsymbol\theta' = \boldsymbol A^{-1}\boldsymbol\theta$ и $\boldsymbol\theta = \boldsymbol A\boldsymbol\theta'$. Проверим: $\boldsymbol X\boldsymbol\theta = \boldsymbol X'(\boldsymbol A\boldsymbol\theta)$, значит, новый вектор коэффициентов — $\boldsymbol\theta' = \boldsymbol A\boldsymbol\theta$, а обратно $\boldsymbol\theta = \boldsymbol A^{-1}\boldsymbol\theta'$. Код ниже проверяет это численно.

Смысл по-пацански: в ортогональном базисе каждый коэффициент считается отдельно, без решения системы: «спроецировал — поделил». Грам — Шмидт превращает любые столбцы в ортогональные, а треугольная матрица $\boldsymbol A$ позволяет потом вернуться к исходным коэффициентам.

Проверь себя