Зависимость оценки параметров от числа факторов модели
8.2. Зависимость оценки параметров от числа факторов модели
Рассмотрим свойство адекватности линейных моделей. Для этого у модели (8.1) матрицу Х и вектор b разделим соответственно на две части, так что она принимает вид
у=Xb+e=[X1, X2]+ε
=X1b1+X2b2+ε. (8.2.1)
Допустим, что при моделировании используется модель y=Х1β1+ε1 с функцией Х1β1, но эта модель после проверки оказывается неадекватной. Статистическая проверка адекватности модели рассмотрена в разделе 6.4. В соответствии с критерием адекватности такой моделью является более разработанная модель y=Х1β1+Х2β2+ε с дополнительным членом Х2β2. Заметить, что адекватность модели зависит не только от числа используемых в ней факторов и параметров, но сейчас рассматриваются свойства оценки параметров только при изменении числа входящих в модель факторов.
Для неадекватной модели y=Х1β1+ε1 оценка вектора β1 параметров методом наименьших квадратов делается по формуле
=(Х1ТХ1)–1Х1Тy.
Но, если эта формула используется для адекватной модели с функцией Е(y)=Х1β1+Х2β2, то вектор ожидаемых значений найденного вектора 1 оценки следующий
Е()= (Х1ТХ1)–1Х1ТЕ(y)
Рекомендуемые материалы
= (Х1ТХ1)–1Х1Т(Х1β1+Х2β2)
=β1+Аβ2, (8.2.2)
где А=(Х1ТХ1)–1Х1ТХ2 - матрица совмещения векторов β1 и β2. Выражение (8.2.2) показывает, что до тех пор пока матрица А не будет нулевой (А=О), вектор является результатом оценки не вектора β1, а совмещённых векторов β1 и β2. Матрица А=О только, если Х1ТХ2=О, то есть, когда все столбцы матрицы Х1 ортогональны всем столбцам матрицы Х2.
Пример 8.2.1. Таблица 8.2.1 содержит данные, полученные в эксперименте по изучению изменения процентного выхода (y) определенного красителя при изменении температуры реакции. Полагалось, что переменная (y) отклика достигает максимума в исследуемом диапазоне температур и правильная зависимость может быть представлена квадратичным (второго порядка) выражением относительно фактора температуры, в виде
y=β0+β1х+β11х2+ε,
где х – нормированный фактор температуры.
Таблица 8.2.1. Выход красителя при разных температурах реакции, его оценочные значения и остатки
Температура реакции Т (град. F) | х | х2 | Выход красителя (y%) | | Остатки |
56 | –1,648 | 2,717 | 45.9 | 45,984 | –0,084 |
60 | –0,33 | 0,109 | 79.8 | 78.615 | 1.185 |
61 | 0 | 0 | 78.9 | 80.458 | –1.558 |
63 | 0,659 | 0,435 | 77.1 | 76.568 | 0.532 |
65 | 1,319 | 1,739 | 62.5 | 62.575 | –0,075 |
В матричном виде модель можно записать так
у=Xb+e,
где у=, Х=
, b=
и e=
.
Оценка методом наименьших квадратов вектора параметров модели даёт вектор Т=[80,458 1,761 –11,618]. Следовательно, оценка ожидаемых значений переменных процентного выхода красителя делается по формуле
=80,458+1,761х–11,618х2.
Экспериментальные данные и график функции модели второго порядка показаны на Рис.8.2.1 соответственно точками и кривой красного цвета. Результаты оценки ожидаемых значений переменных процентного выхода и остатков показаны в двух последних столбцах таблицы 8.2.1.
Рис. 8.2.1. Графики экспериментальных данных и функций первого (синим цветом) и второго (красным цветом) порядков.
Для моделирования представленных в таблице 8.2.1 данных эксперимента положим, что вместо модели y=β0+β1х+β11х2+ε, используется модель y=β0+β1х+ε1. Для этой модели с использованием метода наименьших квадратов получена функция =68,84+6,258х оценки ожидаемых значений переменных процентного выхода, показанная на Рис.8.2.1 прямой синей линией. Она моделирует данные очень плохо. Очевидно, что свободный член
=68,84, являющийся значением
при х=0, и наклон
=6,258 имеют математические ожидания соответственно параметры β0 и β1. Они совмещены с неучтённым параметром β11 кривизны. Коэффициенты совмещения могут быть найдены вычислением матрицы А.
Так как для модели y=β0+β1х+ε1 имеем Х1=, β1=
и Х2=
, то матрица совмещения А=(Х1ТХ1)–1Х1ТХ2=
. Таким образом, имеем Е(
)=β0+1β11 и Е(
)=β1–0,387β11. Отсюда видно, что
=68,84 является результатом оценки параметра β0 модели с членом β11х2, то есть, результатом оценки совмещённых параметров β0+1β11. Аналогично,
- результат оценки не β1, а совмещённых параметров β1–0,387β11.
В случаях, как в этом примере, когда возможно оценить по отдельности параметры правильной модели, то выше показанные выражения совмещения параметров будут справедливы и для оценки методом наименьших квадратов. Этот факт является прямым следствием теоремы Гаусса-Маркова [Box, Draper (2007) стр.77]. Например, в этом примере по формуле =80,458+1,761х–11,618х2 делается оценка ожидаемых значений переменных процентного выхода с учётом х2, а по формуле
=68,84+6,258х оцениваются ожидаемые значения переменных процентного выхода без учёта х2. Тогда легко убедиться, что в пределах ошибки округления, имеем
68,84=80,458–11,618 и 6,258≈1,761–0,387(–11,618).
□
Смещения результатов оценки для неадекватной модели
Теперь обсудим смещение и дисперсии результатов ,
и s2 оценки для неадекватной модели. Рассмотрим сначала оценку вектора b1 параметров неадекватной модели, когда матрица совмещения А≠О, то есть, когда Х1ТХ2≠О. Запишем неадекватную модель в виде
у=X1b1*+e*, (8.2.3)
используя обозначение b1* для вектора параметров чтобы подчеркнуть, что эти параметры (и вектор их оценки) отличаются от вектора b1 (и его оценки
) адекватной модели (8.2.1), если столбцы матрицы Х не ортогональны (см. следствие 1 теоремы 8.2.1 и теорему 8.3). В следующей теореме обсуждается смещение вектора
оценки параметров модели (8.2.3) и даётся матрица ковариаций вектора
.
Теорема 8.2.1. Если постулируется неадекватная модель у=X1b1*+e*, когда адекватной моделью является у=X1b1+X2b2+e, и соблюдается допущение С(у)=s2I, то вычисляемый методом наименьших квадратов вектор =(X1ТX1)–1X1Тy имеет следующие вектор математических ожиданий и матрицу ковариаций:
1. Е()=b1+Аb2, где А=(X1ТX1)–1X1ТX2, (8.2.4)
2. С()=s2(X1ТX1)–1. (8.2.5)
Доказательство:
- Е(
)=E[(X1ТX1)–1X1Тy]=(X1ТX1)–1X1ТЕ(y)
=(X1ТX1)–1X1Т(X1b1+X2b2)
=b1+(X1ТX1)–1X1ТX2b2.
- С(
)=С[(X1ТX1)–1X1Тy]
=(X1ТX1)–1X1Т(s2I)X1(X1ТX1)–1 [в силу (3.6.10)]
=s2(X1ТX1)–1.
□
Таким образом, если постулирована неадекватная модель, то вектор оценки смещён на (X1ТX1)–1X1ТX2b2. Матрица А=(X1ТX1)–1X1ТX2 называется также матрицей смещения. Выражение (8.2.4) показывает, что, если матрица А≠О, то вектор
не результат оценки вектора b1, а результат оценки совмещённых векторов b1 и b2.
Следствие 1. Если X1ТX2=О, то есть, если столбцы матрицы X1 ортогональны столбцам матрицы X2, то результат оценки остаётся несмещённым, то есть, Е(
)=b1.
□
В следующих трёх теоремах обсуждаются свойства оценки ожидаемого значения переменной отклика, дисперсии и дисперсии результатов оценки параметров при постулировании неадекватной модели. Это рассматривается также в [Rencher, Schaalje (2008) стр.170-173; Hocking (2003) стр.140-142].
Пусть вектор-строка x0с= [1, x01, x02, ..., х0(р–1)] содержит конкретные значения нормированных факторов модели для которых желательно оценить ожидаемое значение переменной отклика. Если, в соответствии с разделением матрицы X=[X1, X2] и вектора bТ=[b1Т, b2Т], строку x0с тоже разделить в виде [x01с, x02с], то для оценки ожидаемого значения переменной отклика можно использовать формулы =x0с
и
=x01с
. В следующей теореме получено математическое ожидание результата
оценки ожидаемого значения переменной отклика.
Теорема 8.2.2. Пусть результат оценки =x01с
, где
=(X1ТX1)–1X1Тy. Тогда, если вектор b2≠0 и модель у=X1b1*+e* неадекватна, то математическое ожидание результата оценки
имеет вид
E()=E(x01с
)=x01с(b1+Аb2), (8.2.6)
=x0сb–(x02с–x01сА)b2 (8.2.7)
≠x0сb.
Доказательство: По теореме 8.2.1 имеем Е()=b1+Аb2, поэтому математическое ожидание E(x01с
)=x01сE(
) и получается выражение (8.2.6). С учётом того, что x0сb=[x01с, x02с]
=x01сb1+x02сb2, выражение x01с(b1+Аb2) можно записать в виде
x01с(b1+Аb2)=x01сb1+x01сАb2
=x0сb–x02сb2+x01сАb2
=x0сb–(x02с–x01сА)b2.
Оно показывает, что математическое ожидание E() смещено на (x02с–АТx01с)b2 и не равно x0сb.
□
Из теоремы 8.2.2 видно, что когда постулирована неадекватная модель, то для ожидаемого значения E(y0)=x0сb переменной отклика адекватной модели результат оценки =x01с
смещён. Когда же постулирована адекватная модель, то результат оценки x0с
не смещён, так как E(x0с
)=x0сb=x01сb1+x02сb2, что даёт x01сb1, если b2=0.
В следующей теореме сравниваются дисперсии результатов оценки и
, где
- элементы вектора
и
- элементы вектора
. Сравниваются также дисперсии результатов оценки
=x01с
и
=x0с
ожидаемых значений переменной отклика.
Теорема 8.2.3. Пусть вектор =(XТX)–1XТy оценки параметров адекватной модели разделён
=
в соответствии с разделением матрицы X=[X1, X2] и
=(X1ТX1)–1X1Тy - вектор оценки параметров неадекватной модели. Тогда:
- Дисперсии результатов оценки параметров адекватной модели больше дисперсий результатов оценки параметров неадекватной модели D(
)>D(
).
- Дисперсия оценки ожидаемого значения переменной отклика адекватной модели не меньше дисперсии оценки ожидаемого значения переменной отклика неадекватной модели D(x0с
)≥D(x01с
).
Доказательство:
- Используя матрицу XТX в разделённом виде, имеем ковариационную матрицу вектора
С()=С
=s2(XТX)–1=s2
=s2=s2
,
где Gij=XiТXj и Gij - соответствующий блок разделённой обратной матрицы (XТX)–1. Таким образом матрица ковариаций С()=s2G11. В силу (П.5.6), G11=G11–1+G11–1G12B–1G21G11–1, где B=G22–G21G11–1G12. Вектор оценки параметров неадекватной модели получается по формуле
=(X1ТX1)–1X1Тy, а отдельный результат оценки параметра по формуле
=sjcy, где sjc является j-й строкой матрицы (X1ТX1)–1X1Т. Вектор оценки параметров адекватной модели получается по формуле
=(XТX)–1XТy, а отдельный результат оценки параметра по формуле
=rjcy, где rjc является j-й строкой матрицы
(XТX)–1XТ==
.
В этой матрице, в силу (П.5.6), G12=–G11–1G12B–1. Следовательно, матрица, содержащая сравниваемую со строкой sjc строку rjc, следующая
G11–1X1Т+G11–1G12B–1G21G11–1X1Т–G11–1G12B–1X2Т=(X1ТX1)–1X1Т+АB–1АТX1Т–АB–1X2Т
=(X1ТX1)–1X1Т+АB–1(АТX1Т–X2Т),
где А=G11–1G12. В ней строка rjc содержит строку sjc матрицы (X1ТX1)–1X1Т и j-ю стороку ненулевой матрицы АB–1(X1А–X2)Т. При условии С(у)=s2I, дисперсии результатов оценки параметров адекватной модели D()=С(rjcy)=s2rjcrjcТ и неадекватной модели D(
) =С(sjcy)=s2sjcsjcТ. Строка rjc содержит все элементы строки sjc и дополнительные элементы j-й строки ненулевой матрицы АB–1(X1А–X2)Т. Поэтому rjcrjcТ>sjcsjcТ и D(
)>D(
).
- Дисперсии оценки ожидаемых значений переменной отклика адекватной модели
D(x0с) =s2x0с(XТX)–1x0сТ
=s2[x01с, x02с]
=s2(x01сG11x01сТ+x01сG12x02сТ+x02сG21x01сТ+x02сG22x02сТ)
и неадекватной модели
D(x01с) =s2x01с(X1ТX1)–1x01сТ
=s2x01сG11–1x01сТ.
Используя (П.5.6), их разность получается в виде
D(x0с)–D(x01с
)=s2(x01сG11x01сТ+x01сG12x02сТ+x02сG21x01сТ+x02сG22x02сТ–x01сG11–1x01сТ)
=s2(x01сG11–1G12B–1G21G11–1x01сТ–x01сG11–1G12B–1x02сТ
–x02сB–1G21G11–1x01сТ+x02сB–1x02сТ)
=s2(x01сАB–1АТx01сТ–x01сАB–1x02сТ–x02сB–1АТx01сТ+x02сB–1x02сТ)
=s2(x02сB–1x02сТ–x02сB–1АТx01сТ–x01сАB–1x02сТ+x01сАB–1АТx01сТ)
=s2[x02сB–1(x02сТ–АТx01сТ)–x01сАB–1(x02сТ–АТx01сТ)]
=s2(x02с–x01сА)B–1(x02с–x01сА)Т
≥0,
так как матрица B–1 положительно определённая. Это следует из того, что матрицу В можно представить в виде
B=G22–G21G11–1G12
=X2ТX2–X2ТX1(X1ТX1)–1X1ТX2
=X2Т[I–X1(X1ТX1)–1X1Т]X2
В ней матрица I–X1(X1ТX1)–1X1Т идемпотентная. Следовательно, по теореме П.6.3 матрица В положительно определённая и по теореме П.6.5 матрица B–1 положительно определённая. Таким образом, D(x0с)≥D(x01с
).
□
По пункту 1 теоремы 8.2.3 дисперсии D() результатов оценки параметров адекватной модели больше дисперсий D(
) результатов оценки параметров неадекватной модели. Таким образом, для неадекватной модели дисперсии результатов оценки её параметров меньше, но сами результаты оценки параметров получаются смещёнными. С другой стороны, постулирование адекватной модели увеличивает дисперсии результатов оценки её параметров. По пункту 2 теоремы 8.2.3 дисперсия D(
) результата оценки ожидаемого значения переменной отклика для адекватной модели больше дисперсии D(
) результата этой оценки для неадекватной модели. И снова, для неадекватной модели дисперсия результата оценки ожидаемого значения переменной отклика меньше, но этот результат получается смещённым, а для адекватной модели дисперсия результата оценки ожидаемого значения переменной отклика получается больше.
Смещение результата оценки дисперсии
Рассмотрим теперь результаты оценки дисперсии для адекватной и неадекватной моделей. В случае адекватной модели у=Xb+e=X1b1+X2b2+e оценка дисперсии делается по формуле (7.3.7)
s2=(у–Х)Т(у–Х
)/(n–р).
По теореме 7.3.2 имеем E(s2)=s2. Математическое ожидание результата оценки дисперсии в случае неадекватной модели даётся в следующей теореме.
Теорема 8.2.4. Если модель у=Xb+e=X1b1+X2b2+e адекватная, то для неадекватной модели у=X1b1*+e*, где X1 – матрица размеров nxт и т<р, результат оценки дисперсии по формуле
s12=(у–Х1)Т(у–Х1
)/(n–т) (8.2.8)
имеет математическое ожидание
E(s12)=s2+b2ТX2Т[I–X1(X1ТX1)–1X1Т]X2b2/(n–т). (8.2.9)
Доказательство: Запишем числитель правой части выражения (8.2.8) в виде квадратичной формы
SE*=уТу–Х1Ту
=уТу–yТX1(X1ТX1)–1Х1Ту
=yТ[I–X1(X1ТX1)–1Х1Т]y.
Для нахождения математического ожидания E(s12) оценочной дисперсии необходимо иметь математическое ожидание E(SE*) квадратичной формы. По допущению Е(у)=Xb и теореме 5.2.1 имеем,
E(SE*)=след{[I–X1(X1ТX1)–1Х1Т]s2I}+bТXТ[I–X1(X1ТX1)–1Х1Т]Xb
=(n–т)s2+bТXТ[I–X1(X1ТX1)–1Х1Т]Xb,
где след[I–X1(X1ТX1)–1Х1Т]=след(I)–след[Х1ТX1(X1ТX1)–1]=п–т. В правой части этого выражения можно сделать также следующее преобразование
bТXТ[I–X1(X1ТX1)–1Х1Т]Xb=(b1ТX1Т+b2ТX2Т)[I–X1(X1ТX1)–1Х1Т](X1b1+X2b2)
=b2ТX2Т[I–X1(X1ТX1)–1Х1Т]X2b2.
Отсюда, в силу (8.2.8), получается искомое выражение (8.2.9)
E(s12)=E(SE*)/(n–т)
=s2+b2ТX2Т[I–X1(X1ТX1)–1Х1Т]X2b2/(n–т).
□
Поскольку квадратичная форма в выражении (8.2.9) неотрицательно определённая и, если постулирована неадекватная модель, а вектор b2≠0, то результат оценки дисперсии для неадекватной модели смещён в большую сторону.
Подводя итоги данного раздела можно сказать, что постулирование неадекватной модели, то есть, модели с недостаточным числом факторов, ведёт к смещению результатов оценки её параметров, результата
оценки ожидаемого значения переменной отклика и смещению результата оценки дисперсии. Для достижения адекватности постулирование модели с дополнительными факторами ведёт к увеличению дисперсий результатов
оценки параметров и результата
оценки ожидаемого значения переменной отклика. Таким образом, получение адекватной модели вынуждает искать такое число её факторов, которое обеспечивает адекватность и минимальные дисперсии результатов оценки. Следовательно, разработка адекватной модели состоит в поиске минимально необходимого числа используемых в её функции факторов.
Пример 8.2.2. Положим, что используется модель уi=b0*+b1*xi+ei*, когда адекватной является модель уi=b0+b1xi+b2xi2+ei. В этом случае ,
и s12 будут смещены на величины, зависящие от выбора значений xi фактора [см. (8.2.4) и (8.2.9)]. Для модели уi=b0*+b1*xi+ei* математическое ожидание ошибки ei* не является нулевым:
Е(ei*) =Е(уi–b0*–b1*xi)
=Е(уi)–b0*–b1*xi
=b0+b1xi+b2xi2–b0*–b1*xi
=b0–b0*+(b1–b1*)xi+b2xi2
и поэтому допущение E(ei*)=0 не соблюдается.
□
В лекции "Торфяно-сланцевая промышленность" также много полезной информации.
Пример 8.2.3. Адекватной является модель уi=b0+b1xi+ei, а используется модель уi=b1*xi+ei*. Для модели уi=b1*xi+ei* результат оценки параметра b1* методом наименьших квадратов
=(хTх)–1хТy=
. (8.2.10)
Тогда, так как адекватной моделью является уi=b0+b1xi+ei, то ожидаемое значение оценки получаем в виде
E() =
=
==
. (8.2.11)
Таким образом, оценка смещена на величину, зависящую от b0 и значений xi фактора.