Доверительные интервалы и интервалы предсказания

9.5. Доверительные интервалы и интервалы предсказания

В разделе 6.5 показано, как находить доверительные интервалы отдельных параметров простой модели. В этом разделе для модели с нормированными факторами находится доверительная область вектора b параметров, доверительные интервалы для отдельных параметров и линейной комбинации a^Tb параметров, для математического ожидания E(у) переменной отклика, для дисперсии s² и интервалы предсказания будущих наблюдений. Полагается, что во всём разделе вектор у имеет нормальное распределение N_n(Xb, s²I).

Доверительная область вектора параметров

Если в (9.3.14) матрица C равна единичной матрице I и вектор t равен b, то q=р и получаем величину , имеющую центральное распределение F(р, п–р). Зная, что s²=S_E/(п–р), можно сделать вероятностное утверждение

Рr[(–b)^ТX^ТX(–b)/(рs²)≤F_кр]=1–α,

где F_кр - критическое значение случайной переменной с распределением F(р, п–р) и интегральной вероятностью 1–α на интервале от 0 до F_кр. На основе этого утверждения 100(1–α)% совместная доверительная область для элементов b₀, b₁, ..., b_р–1 вектора b определяется состоящей из всех значений элементов этого вектора, удовлетворяющих условию

(–b)^ТX^ТX(–b)≤рs²F_кр.                                        (9.5.1)

Для р=2 эта область может быть построена в двух измерениях в виде эллипса. Для р>2 данная выражением (9.5.1) эллипсоидная область является громоздкой для интерпретирования и описания. Поэтому рассматриваются интервалы для отдельных параметров, как показано в разделе 6.5.

Доверительные интервалы для отдельных параметров

Рекомендуемые материалы

Если параметр b_j≠0, то в числителе выражения (9.3.6) можно вычесть b_j, так что получается статистика t_j=, имеющая центральное распределение t. Знаменатель этого выражения s часто называют среднеквадратической ошибкой результата  оценки параметра модели, то есть се()=s.

Считая, что вероятность Pr[–t_α/2(п–р)≤≤t_α/2(п–р)]=1–α, и преобразуя неравенство в скобках относительно b_j, получаем

Рr[–t_α/2(п–р)s≤b_j≤+t_α/2(п–р)s]=1–α

или

Рr[–t_α/2(п–р)се()≤b_j≤+t_α/2(п–р)се()]=1–α.

Следовательно, 100(1–α)% доверительный интервал для параметра b_j находится из выражения

±t_α/2(п–р)s.                                       (9.5.2)

Заметим, что доверительная вероятность 1–α имеет место только для одного доверительного интервала и одного из параметров. Доверительные интервалы для всех рассматриваемых одновременно р параметров с общей доверительной вероятностью 1–α даются, например, в книге [Rencher, Schaalje (2008) cтр. 215-217].

Пример 9.5.1. Вычислим 95% доверительные интервалы для всех параметров моделей из примера 8.3. Матрица (X^ТX)^–1 и вектор оценки параметров модели, матрица которой содержит неортогональные между собой вектор-столбцы нормированных значений факторов, имеют соответственно вид:

(X^ТX)^–1= и =.

Матрица (X_о^ТX_о)^–1 и вектор оценки параметров модели, матрица которой содержит все ортогональные между собой вектор-столбцы нормированных значений факторов, имеют соответственно вид:

(X_о^ТX_о)^–1= и =.

Результат s=144,711 оценки получается по формуле s=. Критическое значение t_α/2(п–р) находится по таблице [Box, Draper (2007) стр.745] или с использованием программы Mathcad. График функции плотности вероятности случайной переменной, имеющей распределение t c п–р=11 степенями свободы, показан на Рис.9.5.1.

Рис.9.5.1. График функции плотности вероятности случайной переменной с распределением t(11).

Как и функция плотности вероятности нормального распределения она симметрична относительно 0. Критические значения +t_α/2(п–р) и –t_α/2(п–р) вычисляются на основе того, что за их пределами площадь обеих хвостовых частей под кривой функции равна (1–А)/2. Тогда, если выбирается доверительная вероятность A=0,95, то 95% площади под кривой функции будет находиться между значениями –t_α/2(п–р) и t_α/2(п–р). Критическое значение –t_α/2(п–р)= –2,201 вычисляется с использованием встроенной в Mathcad обратной кумулятивной функции распределения t по формуле

–t_α/2(п–р)=qt[(1–А)/2, (п–р)].

И из-за симметрии кривой функции плотности вероятности значение t_α/2(п–р)=2,201.

Таким образом, по формуле (9.5.2) вычисляются 95% доверительные интервалы для всех параметров моделей. Результаты расчётов сведены в таблицу 9.5.1. Из этой таблицы видно, что для модели с матрицей, содержащей все ортогональные между собой столбцы, 95% доверительный интервал параметра b₁ меньше, чем доверительный интервал этого параметра для модели с матрицей, содержащей неортогональные между собой столбцы.

Таблица 9.5.1. Доверительные интервалы для параметров двух моделей.

□

Доверительный интервал для линейной комбинации параметров

Если линейная комбинация a^Tb≠0, то в числителе выражения (9.3.4) можно вычесть a^Tb из a^T и получить статистику

F_L=,

имеющую распределение F(1, п–р). Тогда, в силу (5.4.9), статистика

t_L=                                                    (9.5.3)

имеет распределение t(п–р) и 100(1–α)% доверительный интервал для одной линейной комбинации a^Tb получается в виде

a^T±t_α/2(п–р)s.                                   (9.5.4)

Доверительный интервал для математического ожидания Е(у)

Пусть вектор-строка x_с представляет конкретный набор значений элементов вектора-строки х_с= [1, x₁, x₂,..., х_(р–1)] факторов, являющейся одной из строк матрицы X. Пусть у - переменная отклика, соответствующая вектору-строке x_с значений влияющих на неё факторов. Тогда модель для этого опыта принимает вид

у=x_сb+e

и, если допустить, что функция модели правильна, так что Е(e)=0, то

E(у) =x_сb.                                                      (9.5.5)

Необходимо найти доверительный интервал для математического ожидания E(у) случайной переменной (у) отклика при заданных значениях элементов вектора x_с.

По следствию 1 теоремы 7.2.4 несмещенной оценкой математического ожидания E(у) с наименьшей дисперсией является

(у)=x_с.                                                     (9.5.6)

Так как (9.5.5) и (9.5.6) имеют соответственно вид x_сb и x_с, то в силу (9.5.4), 100(1–α)% доверительный интервал математического ожидания E(у) =x_сb получается в виде

x_с±t_α/2(п–р) s.                              (9.5.7)

Доверительная вероятность 1–α для интервала в (9.5.7) имеет место только для одного набора значений элементов вектора x_с.

Пример 9.5.2. Положим надо найти 95% доверительный интервал математического ожидания E(у) коэффициента усиления транзистора для данных примера 7.1 при x₁=225 и x₂=4,30 в опытах 9 и 12. В соответствии с выполненным в примере 8.3 нормированием факторов, нормированные их значения x₁=0,597 и x₂=–0,146, а вектор-строка x_с= [1  0,597  –0,146].

Для модели с неортогональными столбцами её матрицы оценка ожидаемого значения переменной отклика (у)=x_с=1361, а для модели с ортогональными столбцами её матрицы (у_о)=x_с=1358. Для расчёта доверительных интервалов используем матрицы (X^ТX)^–1 и (X_о^ТX_о)^–1 из примера 9.5.1, а также s=144,711 и t_α/2(п–р)=2,201. Для модели с неортогональными столбцами её матрицы математическое ожидание E(у)=x_сb и для модели с ортогональными столбцами её матрицы математическое ожидание E(у_о)=x_сb_о. В силу (9.5.7), 95% доверительные интервалы для них получаются соответственно в виде 1260≤x_сb≤1461 и 1258≤x_сb_о≤1458. Длина 95% доверительного интервала для модели с ортогональными столбцами её матрицы равна 199,845 и она меньше длины 95% доверительного интервала для модели с неортогональными столбцами её матрицы равной 201,148. Это ещё раз подтверждает отмеченную в разделе 6.5 предпочтительность ортогональных планов.

□

Интервал предсказания будущего наблюдения

В разделе 7.3 показано, что линейная модель может использоваться для предсказания значений переменной отклика при новых значениях факторов. Пусть строка x_0с= [1, x₀₁, x₀₂,..., х_0(р–1)] представляет новый набор значений элементов строки х_с=[1, x₁, x₂,..., х_(р–1)] факторов. Строка x_0с не является одной из строк матрицы X используемой для оценки параметров модели. Однако если значения элементов этой строки сильно отличаются от значений соответствующих элементов строк матрицы X, то предсказание с применением результата оценки параметров модели с матрицей Х может быть неверным.

Доверительный интервал предсказываемого значения переменной y₀ отклика, соответствующего строке x_0с значений факторов, называется интервалом предсказания. В этом случае говорят об интервале предсказания, а не о доверительном интервале, так как y₀ является случайной переменной отклика, а не параметром. Чтобы быть на 100(1–α)% уверенным, что интервал предсказания содержит предсказываемое значение переменной y₀, очевидно, что этот интервал должен быть больше доверительного интервала для ожидаемого значения E(y).

Для предсказания используется модель y₀=x_0сb+e и предсказываемое значение переменной y₀ находится по формуле =x_0с, что является также результатом оценки ожидаемого значения E(у₀)=x_0сb. Случайные величины y₀ и  независимы, так как значение переменной y₀ предсказываемого наблюдения считается независимым от п наблюдаемых значений переменных отклика, используемых для вычисления . Поэтому дисперсия разности y₀– получается в виде

D(y₀–)=D(y₀–x_0с)

=D(x_0сb+e –x_0с).

Так как x_0сb является постоянной величиной, то имеем

D(y₀–) =D(e)+D(x_0с)

=s²+s²x_0с(X^ТX)^–1x_0с^Т

=s²[1+x_0с(X^ТX)^–1x_0с^Т],                                   (9.5.8)

что можно оценить выражением s²[1+x_0с(X^ТX)^–1x_0с^Т]. Математическое ожидание Е(y₀–) =0, так как Е(e)=0 и Е()=b, а результат s² оценки дисперсии не зависит от y₀ и =x_0с. Таким образом, статистика t_у в виде

t_у=                                          (9.5.9)

имеет распределение t(п–р). Что с вероятностью 1–α её значения находится в интервале между –t_α/2(п–р) и t_α/2(п–р) можно записать в виде

Pr[–t_α/2(п–р)≤≤t_α/2(п–р)]=1–α.

Для получения 100(1–α)% интервала предсказания значений переменной y₀ неравенство в квадратных скобках преобразуем к виду

–t_α/2(п–р)s≤y₀≤+t_α/2(п–р)s

и, учитывая, что =x_0с, получаем выражение для этого интервала

x_0с±s t_α/2(п–р).                       (9.5.10)

Заметим, что доверительная вероятность 1–α для интервала предсказания (9.5.10) имеет место только для одного набора значений элементов строки x_0с.

В выражении 1+x_0с(X^ТX)^–1x_0с^Т второй член x_0с(X^ТX)^–1x_0с^Т, как правило, гораздо меньше единицы (при условии, что р гораздо меньше чем п), так как дисперсия оценки =x_0с гораздо меньше дисперсии переменной y₀. [Например, если бы X^ТX была диагональной и значения элементов x₀ были в области охватываемой значениями элементов строк матрицы X, то x_0с(X^ТX)^–1x_0с^Т была бы суммой из р слагаемых, каждое из которых вида x_0j²/, что порядка 1/п.] Таким образом, интервалы предсказания переменной y₀, как правило, значительно больше доверительных интервалов математического ожидания E(у)=x_сb.

Пример 9.5.3. Положим надо найти 95% интервал предсказания математического ожидания E(у₀) коэффициента усиления транзистора для данных примера 7.1 при x₀₁=225 и x₀₂=4,36. В соответствии с выполненным в примере 8.3 нормированием факторов, нормированные их значения соответственно x₀₁=0,597 и x₀₂=0,090, а строка x_0с= [1  0,597  0,090].

Для модели с неортогональными столбцами её матрицы результат оценки математического ожидания переменной отклика (у₀)=x_0с=1348, а для модели с ортогональными столбцами её матрицы (у_0о)=x_0с=1345. Для расчёта интервалов предсказания используем матрицы (X^ТX)^–1 и (X_о^ТX_о)^–1 примера 9.5.1, а также s=144,711 и t_α/2(п–р)=2,201. Для модели с неортогональными столбцами её матрицы математическое ожидание E(у₀) =x_0сb, а для модели с ортогональными столбцами её матрицы E(у₀) =x_0сb_о. В силу (9.5.10), они имеют 95% интервалы предсказания соответственно в виде: 1014≤x_0сb≤1681 и 1011≤x_0сb_о≤1679. Заметим, что эти интервалы предсказания для у₀ длинами соответственно 667,193 и 667,339 гораздо больше длин доверительных интервалов для Е(у) в примере 9.5.2.

□

Доверительный интервал для дисперсии

По пункту 2 теоремы 7.3.4 величина (п–р)s²/s² имеет распределение c²(п–р). Поэтому вероятностное высказывание, что её значения находятся в интервале между критическими значениями c²_1–α/2(п–р) и c²_α/2(п–р), можно представить в виде

Pr[c²_α/2(п–р)≤≤c²_1–α/2(п–р)]=1–α,                         (9.5.11)

где c²_α/2(п–р) - меньшее α/2-процентное критическое значение распределения хи-квадрат и c²_1–α/2(п–р) - большее 1–α/2-процентное критическое значение. Преобразование неравенства в квадратных скобках относительно s² приводит к неравенству

(п–р)s²/c²_1–α/2(п–р)≤s²≤(п–р)s²/c²_α/2(п–р),                           (9.5.12)

являющимся 100(1–α)% доверительным интервалом для дисперсии. Отсюда 100(1–α)% доверительный интервал стандартного отклонения s дается неравенством

[(п–р)s²/c²_1–α/2(п–р)]^1/2≤s≤[(п–р)s²/c²_α/2(п–р)]^1/2.                 (9.5.13)

Пример 9.5.4. В формулах расчёта доверительных интервалов (9.5.7) и интервалов предсказания (9.5.10) используется результат оценки стандартного отклонения s. Найдём доверительный интервал этого стандартного отклонения для данных из примеров 9.5.2 и 9.5.3.

На Рис.9.5.2 показано, что в отличие от распределения t распределение c²(11) несимметричное. Поэтому критические значения c²_1–α/2(п–р) и c²_α/2(п–р) необходимо вычислять отдельно. Для этого используется то, что площадь хвоста под кривой функции плотности вероятности слева от c²_α/2(п–р) равна (1–А)/2, а площадь под кривой функции плотности вероятности слева от c²_1–α/2(п–р) равно (1–А)/2+А=(1+А)/2. Например, если для рассматриваемых данных п–р=11 и доверительную вероятность А принять равной 0,95, то, используя встроенную в Mathcad обратную кумулятивную функцию распределения c², меньшее критическое значение c²_0,025(11)=3,816 рассчитывается по формуле

c²_0,025(11) =qchisq[(1–А)/2, 11],

а большее критическое значение c²_0,975(11)=21,920 рассчитывается по формуле

c²_0,975(11) =qchisq[(1+А)/2, 11].

Следовательно, выражение (9.5.13) вычисления доверительного интервала стандартного отклонения для рассматриваемых данных принимает вид

[11s²/c²_0,975(11)]^1/2≤s≤[11s²/c²_0,025(11)]^1/2

и в результате вычисления имеем 102,512≤s≤245,701.

Рис.9.5.2. График функции плотности вероятности случайной переменной с распределением c²(11).

□

Упражнения

9.1. Покажите, что H_c[I–Е/n]=H_c, как и в (9.1.4) пункта 1 теоремы 9.1.1, где матрица H_c=X_c(X_c^ТX_c)^–1X_c^Т.

9.2. Покажите, что среднее квадратичное S_Rс/(p–1) в таблице 9.1.1 имеет математическое ожидание E[S_Rс/(p–1)]=s²+b₁^ТX₁^ТX₁b₁/(p–1). Примените следующие два подхода:

(a) Используйте теорему 5.2.1.

(б) Используйте параметр не центральности в (5.19).

9.3. Разработайте процедуру проверки гипотезы H₀: b=0 для модели у=Xb+e, где вектор у имеет распределение N_n(Xb, s²I). (В начале раздела 9.1 было отмечено, что эта гипотеза не имеет большого практического интереса, потому что включает в себя b₀=0.) Используйте выражение y^Тy= (y^Тy–^ТX^Тy)+^ТX^Тy и действуйте следующим образом:

1. Покажите, что ^ТX^Тy=y^ТX(X^ТX)^–1X^Тy и y^Тy–^ТX^Тy=y^Т[I–X(X^ТX)^–1X^Т]у.

2. Пусть H=X(X^ТX)^–1X^Т. Покажите, что H и I–H идемпотентные и соответственно рангов р и п–р.

3. Покажите, что y^ТHy/s² имеет распределение c²(р, l₁), где l₁=b^ТX^ТXb/(2s²), и что y^Т(I–H)у/s² имеет распределение c²(п–р).

4. Покажите, что y^ТHy и y^Т(I–H)у независимы.

5. Покажите, что = имеет распределение F(р, п–р, l₁).

9.4. Покажите, что HH₁=H₁ и H₁H=H₁, как и в (9.2.9), где Н и H₁ являются такими, как определено в (9.2.5) и (9.2.6).

9.5. Покажите, что первое и второе условия следствия 1 теоремы 5.6.3 удовлетворяются для квадратичных форм выражения (9.2.6), как отмечено в доказательстве теоремы 9.2.2.

9.6. Найдите ожидаемое среднее квадратичное числителя данной в примере 9.2.2 выражением (9.2.14) статистики F_Н.

9.7. В доказательстве теоремы 9.2.4 покажите, что (^ТX₁^Т+^ТX₂^Т)(X₁^Т+X₂) –(^Т+^ТА^Т)X₁^ТX₁(+А)=^Т[X₂^ТX₂–X₂^ТX₁(X₁^ТX₁)^–1X₁^ТX₂] .

9.8. Покажите, что числитель ^ТX^Тy–^ТX₁^Тy выражения (9.3.3) равен /g_jj в (9.3.5) (для j=р–1), как указано после (9.3.5).

Информация в лекции "Испарение" поможет Вам.

9.9. Получите доверительный интервал для a^Tb в (9.5.4) по статистике t_L в (9.5.3).

9.10. Покажите, что t_у= в (9.5.9) имеет распределение t(п–р).

9.11.    1. Учитывая, что =y_0i/q представляет собой среднее q будущих наблюдений при х₀, покажите, что 100(1–α)% интервал предсказания для  задается в виде

х_0с±s t_α/2(п–р).

2. Покажите, что для простой линейной регрессии интервал предсказания  в пункте 1 сводится к виду +x₀±s t_α/2(п–2).

9.12. Получите доверительный интервал для s² в (9.5.12) из вероятностного высказывания (9.5.11).