§5. Взвешенный МНК

Почти все величайшие открытия в астрономии вытекают из рассмотрения того, что мы уже раньше назвали качественными или численными остаточными феноменами, иначе говоря, они вытекают из анализа той части численных или качественных результатов наблюдения, которая «торчит» и остаётся необъяснённой после выделения и учёта всех известных причин, что согласуется со строгим применением известных методов.

— Дж. Гершель, «Основы астрономии», 1849 г.

После оценивания параметров регрессионной модели (5) полезно изучить поведение реализации вектора остатков $\hat{\boldsymbol\delta} = \boldsymbol\eta - \boldsymbol X\hat{\boldsymbol\theta}$. В случае хорошего соответствия данных и модели случайные величины $\hat\delta_1, \dots, \hat\delta_n$ должны быть похожими на выборку из закона распределения ошибок $\varepsilon_i$. Однако нередко наблюдается увеличение разброса значений $\hat\delta_i$ с ростом $i$ («воронка») или чередование серий из положительных и отрицательных значений, возникающее из-за сильной коррелированности соседних остатков (см. §5 гл. 15).

В связи с этим рассмотрим два обобщения модели (5), возникающие в результате отказа от допущений одинаковой распределённости и некоррелированности ошибок $\varepsilon_i$.

Смысл по-пацански: остатки — это «что не объяснила модель». Если в остатках виден узор — воронка, волны, серии — модель что-то упустила. Чистые остатки должны выглядеть как бессмысленный шум.

1. Неравноточные наблюдения

Допустим, что в модели (5)

где $\boldsymbol V$ — известная диагональная матрица с элементами $v_i$ на главной диагонали, $\sigma$ — неизвестный параметр. Другими словами, ошибки некоррелированы и $\mathbf D\varepsilon_i = \sigma^2v_i$.

Положим $w_i = 1/v_i$ и, меняя масштаб с помощью преобразования $\tilde\eta_i = \sqrt{w_i}\,\eta_i$, получим новую модель:

где $\tilde x_{ij} = \sqrt{w_i}\,x_{ij}$, $\tilde\varepsilon_i = \sqrt{w_i}\,\varepsilon_i$. Из формул (36) и (37) следует, что $\mathbf M\tilde\varepsilon_i = 0$, $\mathbf D\tilde\varepsilon_i = \sigma^2$, $\operatorname{cov}(\tilde\varepsilon_i, \tilde\varepsilon_j) = 0$ при $i \ne j$, то есть для новой модели выполняется допущение Д2 из §2.

МНК-оценки параметров $\theta_1, \dots, \theta_m$ находятся (согласно §2) путём минимизации функции $F(\boldsymbol\theta) = \sum_{i=1}^n(\tilde\eta_i - \theta_1\tilde x_{i1} - \dots - \theta_m\tilde x_{im})^2$. А для исходной модели имеем

Выходит, что теперь мы будем «взвешивать» каждый член суммы, умножая его на $w_i = 1/v_i$. Это придаёт больший вес наблюдениям, дисперсии ошибок которых меньше, что интуитивно представляется вполне разумным. Такая процедура называется взвешенным методом наименьших квадратов.

Нетрудно видеть, что матрица $\tilde{\boldsymbol X}$ из формулы (37) представляется в виде $\boldsymbol V^{-1/2}\boldsymbol X$ (здесь $\boldsymbol V^{-1/2}$ — диагональная матрица с элементами $1/\sqrt{v_i}$ на диагонали). Аналогично, $\tilde{\boldsymbol\eta} = \boldsymbol V^{-1/2}\boldsymbol\eta$. Подставляя эти выражения в формулу (7), находим МНК-оценку

Это и есть значение $\boldsymbol\theta$, минимизирующее сумму квадратов (38), которая в матричной записи имеет вид

Согласно утверждению 1 из §3, ковариационной матрицей случайного вектора $\hat{\boldsymbol\theta}$ служит

Смысл по-пацански: свидетель, который стоял рядом, — его слово весит больше; тот, что видел издалека, — меньше. Вес — обратно пропорционален «мутности» показаний. Делим каждое уравнение на его $\sqrt{v_i}$ — и снова попадаем в старый добрый МНК.

Разберём несколько приложений рассмотренной модели.

Пример 5. Взвешенное среднее. Если модель (5) содержит единственный неизвестный параметр $\theta$ ($m = 1$), то из формулы (39) несложно получить, что его МНК-оценка с учётом допущений (36) выглядит так:

В случае, когда все $x_i = 1$, она представляет собой просто взвешенное среднее наблюдений $\eta_i$ с весами $w_i$.

Пример 6. Линейный закон — прямая, проходящая через начало координат. Пусть $Y$ — потери тепла в (стандартном) коттедже, а $X$ обозначает разность между внутренней и наружной температурой. Рассмотрим линейную модель $Y = \theta X$, при этом дисперсия наблюдений $\eta_i = \theta x_i + \varepsilon_i$ ($i = 1, \dots, n$) будет, скорее всего, возрастать с увеличением $x_i$. Простейший вид такой зависимости — это $\mathbf D\varepsilon_i = cx_i$, то есть дисперсия пропорциональна разности температур. В этом случае из формулы (42) при $w_i = \sigma^2/(cx_i)$ находим, что

Ввиду соотношения (41) $\mathbf D\hat\theta = \sigma^2/\sum w_ix_i^2 = c/\sum x_i$. Если же допустить, что $\mathbf D\varepsilon_i = cx_i^2$, то приходим к оценке $\hat\theta = \frac1n\sum(\eta_i/x_i)$. Чему равна дисперсия $\mathbf D\hat\theta$ в этом случае — вопрос 3.

При $\mathbf D\varepsilon_i \propto x_i^2$ и $x$ от 1 до 10 взвешенный МНК даёт дисперсию оценки примерно в 1,65 раза меньше, чем обычный. При $\mathbf D\varepsilon_i \propto x_i$ выигрыш скромнее — около 1,12 раза: веса различаются меньше.

Пример 7. Повторные наблюдения. Пусть $\eta_i$ — среднее из $k_i$ наблюдений, причём все они имеют одинаковое математическое ожидание $\sum_{j=1}^mx_{ij}\theta_j$ и общую дисперсию $\sigma^2$, $i = 1, \dots, n$. Если все наблюдения независимы, то $\mathbf D\eta_i = \sigma^2/k_i$, так что в формуле (36) надо положить $v_i = 1/k_i$.

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

2. Коррелированные наблюдения

Рассмотрим обобщение модели (36), допуская, что наблюдения могут не только иметь разную точность, но и быть зависимыми между собой. Для этого матрицу $\boldsymbol V$ в формуле (36) будем считать известной положительно определённой (см. П10) матрицей общего вида. Ковариации ошибок связаны с её элементами $v_{ij}$ равенствами $\operatorname{cov}(\varepsilon_i, \varepsilon_j) = \sigma^2v_{ij}$.

Для такой матрицы существует единственная положительно определённая матрица $\boldsymbol V^{1/2}$, называемая квадратным корнем из матрицы $\boldsymbol V$, такая, что $\boldsymbol V^{1/2}\boldsymbol V^{1/2} = \boldsymbol V$ (см. П10).

Умножая обе части модели (5) на обратную к $\boldsymbol V^{1/2}$ матрицу $\boldsymbol V^{-1/2}$, приходим к новой модели $\tilde{\boldsymbol\eta} = \tilde{\boldsymbol X}\boldsymbol\theta + \tilde{\boldsymbol\varepsilon}$, где $\tilde{\boldsymbol\eta} = \boldsymbol V^{-1/2}\boldsymbol\eta$, $\tilde{\boldsymbol X} = \boldsymbol V^{-1/2}\boldsymbol X$, $\tilde{\boldsymbol\varepsilon} = \boldsymbol V^{-1/2}\boldsymbol\varepsilon$. В силу формул (12) и симметричности матрицы $\boldsymbol V^{-1/2}$ имеем:

то есть видим, что для новой модели опять выполняется допущение Д2 из §2. Нетрудно убедиться, что формулы (39)–(41) остаются без изменений (благодаря симметричности $\boldsymbol V^{-1/2}$) и в случае любой положительно определённой матрицы $\boldsymbol V$.

Смысл по-пацански: матрица $\boldsymbol V^{-1/2}$ «отбеливает» шум: перемешивает наблюдения так, чтобы ошибки стали независимыми и одинаковыми. После этого — обычный МНК.

Оба метода несмещены, но обобщённый МНК (39) точнее — и его разброс совпадает с формулой (41).

Пример 8. Сдвиг и масштаб по порядковым статистикам

Пусть элементы выборки $X_1, \dots, X_n$ имеют функцию распределения $F((x - \mu)/\sigma)$, где $F(u)$ — известная функция распределения, $\mu \in \mathbb R$ и $\sigma > 0$ — неизвестные параметры сдвига и масштаба соответственно. Рассмотрим стандартизованные случайные величины $U_i = (X_i - \mu)/\sigma$, $i = 1, \dots, n$. Тогда $U_1, \dots, U_n$ — выборка из закона с функцией распределения $F(u)$. При этом порядковые статистики (см. §2 гл. 5) двух выборок связаны формулами

и так как функция распределения $F(u)$ случайной величины $U_i$ известна, то в принципе любые вероятностные характеристики порядковых статистик $U_{(i)}$ могут быть рассчитаны. Положим

Учитывая формулы (43), запишем соотношения

Вводя обозначения $\eta_i = X_{(i)}$ и $\varepsilon_i = \sigma(U_{(i)} - \alpha_i)$, представим (44) в виде обобщённой линейной регрессионной модели: $\eta_i = \mu + \sigma\alpha_i + \varepsilon_i$, где $\mathbf M\varepsilon_i = 0$, $\operatorname{cov}(\varepsilon_i, \varepsilon_j) = \sigma^2v_{ij}$, или в матричной форме:

Согласно замечанию 5, МНК-оценки параметров $\mu$ и $\sigma$ в равенстве (45), вычисляемые по формуле (39), имеют минимальную дисперсию в классе линейных несмещённых оценок, построенных по порядковым статистикам $X_{(1)}, \dots, X_{(n)}$. Приведём некоторые результаты вычислений для равномерного и показательного законов.

Равномерный закон. Для $U_i$, равномерно распределённых на отрезке $[0, 1]$, в задаче 1 гл. 5 были подсчитаны $\alpha_i = \mathbf MU_{(i)} = i/(n + 1)$. Известно, что

Матрица $\boldsymbol V^{-1}$ имеет элементы $(n + 1)(n + 2)p_{ij}$, где $p_{ii} = 2$, $p_{i,i-1} = p_{i-1,i} = -1$, $p_{ij} = 0$ при $|i - j| > 1$, а МНК-оценки параметров масштаба и сдвига выглядят так:

Они (с точностью до репараметризации модели) совпадают с оценками метода спейсингов из задачи 4 гл. 9.

Показательный закон. Для показательно распределённых $U_i$ ($F(u) = (1 - e^{-u})I_{\{u\ge0\}}$) при решении задачи 5 гл. 4 было установлено, что

независимы и распределены так же, как и $U_i$. Поскольку $\mathbf MZ_i = 1$,

Матрица $\boldsymbol V^{-1}$ имеет элементы $q_{ij}$, где $q_{ii} = (n - i + 1)^2 + (n - i)^2$, $q_{i,i-1} = q_{i-1,i} = -(n - i + 1)^2$ ($q_{ij} = 0$ при $|i - j| > 1$), а МНК-оценки параметров $\sigma$ и $\mu$ таковы:

В случае нормального закона для моментов порядковых статистик $U_{(i)}$ нет простых формул. Известно, что $\hat\mu = \bar X$, однако для вычисления $\hat\sigma$ приходится пользоваться таблицами.

Проверим формулы для равномерного закона, честно посчитав (39) с матрицей $\boldsymbol V$:

Смысл по-пацански: для равномерного закона вся информация о сдвиге и масштабе сидит в двух крайних точках — минимуме и максимуме. Средние точки МНК «выбрасывает» сам: веса у них нулевые.

Проверь себя