Зависимость оценки параметров от числа факторов модели

8.2. Зависимость оценки параметров от числа факторов модели

Рассмотрим свойство адекватности линейных моделей. Для этого у модели (8.1) матрицу Х и вектор b разделим соответственно на две части, так что она принимает вид

у=Xb+e=[X₁, X₂]+ε

=X₁b₁+X₂b₂+ε.                                             (8.2.1)

Допустим, что при моделировании используется модель y=Х₁β₁+ε₁ с функцией Х₁β₁, но эта модель после проверки оказывается неадекватной. Статистическая проверка адекватности модели рассмотрена в разделе 6.4. В соответствии с критерием адекватности такой моделью является более разработанная модель y=Х₁β₁+Х₂β₂+ε с дополнительным членом Х₂β₂. Заметить, что адекватность модели зависит не только от числа используемых в ней факторов и параметров, но сейчас рассматриваются свойства оценки параметров только при изменении числа входящих в модель факторов.

Для неадекватной модели y=Х₁β₁+ε₁ оценка вектора β₁ параметров методом наименьших квадратов делается по формуле

=(Х₁^ТХ₁)^–1Х₁^Тy.

Но, если эта формула используется для адекватной модели с функцией Е(y)=Х₁β₁+Х₂β₂, то вектор ожидаемых значений найденного вектора ₁ оценки следующий

Е()= (Х₁^ТХ₁)^–1Х₁^ТЕ(y)

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

= (Х₁^ТХ₁)^–1Х₁^Т(Х₁β₁+Х₂β₂)

=β₁+Аβ₂,                                                        (8.2.2)

где А=(Х₁^ТХ₁)^–1Х₁^ТХ₂ - матрица совмещения векторов β₁ и β₂. Выражение (8.2.2) показывает, что до тех пор пока матрица А не будет нулевой (А=О), вектор  является результатом оценки не вектора β₁, а совмещённых векторов β₁ и β₂. Матрица А=О только, если Х₁^ТХ₂=О, то есть, когда все столбцы матрицы Х₁ ортогональны всем столбцам матрицы Х₂.

Пример 8.2.1. Таблица 8.2.1 содержит данные, полученные в эксперименте по изучению изменения процентного выхода (y) определенного красителя при изменении температуры реакции. Полагалось, что переменная (y) отклика достигает максимума в исследуемом диапазоне температур и правильная зависимость может быть представлена квадратичным (второго порядка) выражением относительно фактора температуры, в виде

y=β₀+β₁х+β₁₁х²+ε,

где х – нормированный фактор температуры.

Таблица 8.2.1. Выход красителя при разных температурах реакции, его оценочные значения и остатки

В матричном виде модель можно записать так

у=Xb+e,

где у=, Х=, b= и e=.

Оценка методом наименьших квадратов вектора параметров модели даёт вектор ^Т=[80,458  1,761  –11,618]. Следовательно, оценка ожидаемых значений переменных процентного выхода красителя делается по формуле

=80,458+1,761х–11,618х².

Экспериментальные данные и график функции модели второго порядка показаны на Рис.8.2.1 соответственно точками и кривой красного цвета. Результаты оценки ожидаемых значений переменных процентного выхода и остатков показаны в двух последних столбцах таблицы 8.2.1.

Рис. 8.2.1. Графики экспериментальных данных и функций первого (синим цветом) и второго (красным цветом) порядков.

Для моделирования представленных в таблице 8.2.1 данных эксперимента положим, что вместо модели y=β₀+β₁х+β₁₁х²+ε, используется модель y=β₀+β₁х+ε₁. Для этой модели с использованием метода наименьших квадратов получена функция =68,84+6,258х оценки ожидаемых значений переменных процентного выхода, показанная на Рис.8.2.1 прямой синей линией. Она моделирует данные очень плохо. Очевидно, что свободный член =68,84, являющийся значением  при х=0, и наклон =6,258 имеют математические ожидания соответственно параметры β₀ и β₁. Они совмещены с неучтённым параметром β₁₁ кривизны. Коэффициенты совмещения могут быть найдены вычислением матрицы А.

Так как для модели y=β₀+β₁х+ε₁ имеем Х₁=, β₁= и Х₂=, то матрица совмещения А=(Х₁^ТХ₁)^–1Х₁^ТХ₂=. Таким образом, имеем Е()=β₀+1β₁₁ и Е()=β₁–0,387β₁₁. Отсюда видно, что =68,84 является результатом оценки параметра β₀ модели с членом β₁₁х², то есть, результатом оценки совмещённых параметров β₀+1β₁₁. Аналогично,  - результат оценки не β₁, а совмещённых параметров β₁–0,387β₁₁.

В случаях, как в этом примере, когда возможно оценить по отдельности параметры правильной модели, то выше показанные выражения совмещения параметров будут справедливы и для оценки методом наименьших квадратов. Этот факт является прямым следствием теоремы Гаусса-Маркова [Box, Draper (2007) стр.77]. Например, в этом примере по формуле =80,458+1,761х–11,618х² делается оценка ожидаемых значений переменных процентного выхода с учётом х², а по формуле =68,84+6,258х оцениваются ожидаемые значения переменных процентного выхода без учёта х². Тогда легко убедиться, что в пределах ошибки округления, имеем

68,84=80,458–11,618 и 6,258≈1,761–0,387(–11,618).

□

Смещения результатов оценки для неадекватной модели

Теперь обсудим смещение и дисперсии результатов ,  и s² оценки для неадекватной модели. Рассмотрим сначала оценку вектора b₁ параметров неадекватной модели, когда матрица совмещения А≠О, то есть, когда Х₁^ТХ₂≠О. Запишем неадекватную модель в виде

у=X₁b₁^*+e^*,                                                   (8.2.3)

используя обозначение b₁^* для вектора параметров чтобы подчеркнуть, что эти параметры (и вектор  их оценки) отличаются от вектора b₁ (и его оценки ) адекватной модели (8.2.1), если столбцы матрицы Х не ортогональны (см. следствие 1 теоремы 8.2.1 и теорему 8.3). В следующей теореме обсуждается смещение вектора  оценки параметров модели (8.2.3) и даётся матрица ковариаций вектора .

Теорема 8.2.1. Если постулируется неадекватная модель у=X₁b₁^*+e^*, когда адекватной моделью является у=X₁b₁+X₂b₂+e, и соблюдается допущение С(у)=s²I, то вычисляемый методом наименьших квадратов вектор =(X₁^ТX₁)^–1X₁^Тy имеет следующие вектор математических ожиданий и матрицу ковариаций:

1. Е()=b₁+Аb₂, где А=(X₁^ТX₁)^–1X₁^ТX₂,                                                 (8.2.4)

2. С()=s²(X₁^ТX₁)^–1.                                                                                  (8.2.5)

Доказательство:

1. Е()=E[(X₁^ТX₁)^–1X₁^Тy]=(X₁^ТX₁)^–1X₁^ТЕ(y)

=(X₁^ТX₁)^–1X₁^Т(X₁b₁+X₂b₂)

=b₁+(X₁^ТX₁)^–1X₁^ТX₂b₂.

1. С()=С[(X₁^ТX₁)^–1X₁^Тy]

=(X₁^ТX₁)^–1X₁^Т(s²I)X₁(X₁^ТX₁)^–1           [в силу (3.6.10)]

=s²(X₁^ТX₁)^–1.

□

Таким образом, если постулирована неадекватная модель, то вектор  оценки смещён на (X₁^ТX₁)^–1X₁^ТX₂b₂. Матрица А=(X₁^ТX₁)^–1X₁^ТX₂ называется также матрицей смещения. Выражение (8.2.4) показывает, что, если матрица А≠О, то вектор  не результат оценки вектора b₁, а результат оценки совмещённых векторов b₁ и b₂.

Следствие 1. Если X₁^ТX₂=О, то есть, если столбцы матрицы X₁ ортогональны столбцам матрицы X₂, то результат  оценки остаётся несмещённым, то есть, Е()=b₁.

□

В следующих трёх теоремах обсуждаются свойства оценки ожидаемого значения переменной отклика, дисперсии и дисперсии результатов  оценки параметров при постулировании неадекватной модели. Это рассматривается также в [Rencher, Schaalje (2008) стр.170-173; Hocking (2003) стр.140-142].

Пусть вектор-строка x_0с= [1, x₀₁, x₀₂, ..., х_0(р–1)] содержит конкретные значения нормированных факторов модели для которых желательно оценить ожидаемое значение переменной отклика. Если, в соответствии с разделением матрицы X=[X₁, X₂] и вектора b^Т=[b₁^Т, b₂^Т], строку x_0с тоже разделить в виде [x_01с, x_02с], то для оценки ожидаемого значения переменной отклика можно использовать формулы =x_0с и =x_01с. В следующей теореме получено математическое ожидание результата  оценки ожидаемого значения переменной отклика.

Теорема 8.2.2. Пусть результат оценки =x_01с, где =(X₁^ТX₁)^–1X₁^Тy. Тогда, если вектор b₂≠0 и модель у=X₁b₁^*+e^* неадекватна, то математическое ожидание результата оценки  имеет вид

E()=E(x_01с)=x_01с(b₁+Аb₂),                              (8.2.6)

=x_0сb–(x_02с–x_01сА)b₂                                     (8.2.7)

≠x_0сb.

Доказательство: По теореме 8.2.1 имеем Е()=b₁+Аb₂, поэтому математическое ожидание E(x_01с)=x_01сE() и получается выражение (8.2.6). С учётом того, что x_0сb=[x_01с, x_02с]=x_01сb₁+x_02сb₂, выражение x_01с(b₁+Аb₂) можно записать в виде

x_01с(b₁+Аb₂)=x_01сb₁+x_01сАb₂

=x_0сb–x_02сb₂+x_01сАb₂

=x_0сb–(x_02с–x_01сА)b₂.

Оно показывает, что математическое ожидание E() смещено на (x_02с–А^Тx_01с)b₂ и не равно x_0сb.

□

Из теоремы 8.2.2 видно, что когда постулирована неадекватная модель, то для ожидаемого значения E(y₀)=x_0сb переменной отклика адекватной модели результат оценки =x_01с смещён. Когда же постулирована адекватная модель, то результат оценки x_0с не смещён, так как E(x_0с)=x_0сb=x_01сb₁+x_02сb₂, что даёт x_01сb₁, если b₂=0.

В следующей теореме сравниваются дисперсии результатов оценки  и , где  - элементы вектора  и  - элементы вектора . Сравниваются также дисперсии результатов оценки =x_01с и =x_0с ожидаемых значений переменной отклика.

Теорема 8.2.3. Пусть вектор =(X^ТX)^–1X^Тy оценки параметров адекватной модели разделён = в соответствии с разделением матрицы X=[X₁, X₂] и =(X₁^ТX₁)^–1X₁^Тy - вектор оценки параметров неадекватной модели. Тогда:

1. Дисперсии результатов оценки параметров адекватной модели больше дисперсий результатов оценки параметров неадекватной модели D()>D().
2. Дисперсия оценки ожидаемого значения переменной отклика адекватной модели не меньше дисперсии оценки ожидаемого значения переменной отклика неадекватной модели D(x_0с)≥D(x_01с).

Доказательство:

1. Используя матрицу X^ТX в разделённом виде, имеем ковариационную матрицу вектора

С()=С=s²(X^ТX)^–1=s²

=s²=s²,

где G_ij=X_i^ТX_j и G^ij - соответствующий блок разделённой обратной матрицы (X^ТX)^–1. Таким образом матрица ковариаций С()=s²G¹¹. В силу (П.5.6), G¹¹=G₁₁^–1+G₁₁^–1G₁₂B^–1G₂₁G₁₁^–1, где B=G₂₂–G₂₁G₁₁^–1G₁₂. Вектор оценки параметров неадекватной модели получается по формуле =(X₁^ТX₁)^–1X₁^Тy, а отдельный результат оценки параметра по формуле =s_jcy, где s_jc является j-й строкой матрицы (X₁^ТX₁)^–1X₁^Т. Вектор оценки параметров адекватной модели получается по формуле =(X^ТX)^–1X^Тy, а отдельный результат оценки параметра по формуле =r_jcy, где r_jc является j-й строкой матрицы

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

В этой матрице, в силу (П.5.6), G¹²=–G₁₁^–1G₁₂B^–1. Следовательно, матрица, содержащая сравниваемую со строкой s_jc строку r_jc, следующая

G₁₁^–1X₁^Т+G₁₁^–1G₁₂B^–1G₂₁G₁₁^–1X₁^Т–G₁₁^–1G₁₂B^–1X₂^Т=(X₁^ТX₁)^–1X₁^Т+АB^–1А^ТX₁^Т–АB^–1X₂^Т

=(X₁^ТX₁)^–1X₁^Т+АB^–1(А^ТX₁^Т–X₂^Т),

где А=G₁₁^–1G₁₂. В ней строка r_jc содержит строку s_jc матрицы (X₁^ТX₁)^–1X₁^Т и j-ю стороку ненулевой матрицы АB^–1(X₁А–X₂)^Т. При условии С(у)=s²I, дисперсии результатов оценки параметров адекватной модели D()=С(r_jcy)=s²r_jcr_jc^Т и неадекватной модели D() =С(s_jcy)=s²s_jcs_jc^Т. Строка r_jc содержит все элементы строки s_jc и дополнительные элементы j-й строки ненулевой матрицы АB^–1(X₁А–X₂)^Т. Поэтому r_jcr_jc^Т>s_jcs_jc^Т и D()>D().

1. Дисперсии оценки ожидаемых значений переменной отклика адекватной модели

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

=s²[x_01с, x_02с]

=s²(x_01сG¹¹x_01с^Т+x_01сG¹²x_02с^Т+x_02сG²¹x_01с^Т+x_02сG²²x_02с^Т)

и неадекватной модели

D(x_01с) =s²x_01с(X₁^ТX₁)^–1x_01с^Т

=s²x_01сG₁₁^–1x_01с^Т.

Используя (П.5.6), их разность получается в виде

D(x_0с)–D(x_01с)=s²(x_01сG¹¹x_01с^Т+x_01сG¹²x_02с^Т+x_02сG²¹x_01с^Т+x_02сG²²x_02с^Т–x_01сG₁₁^–1x_01с^Т)

=s²(x_01сG₁₁^–1G₁₂B^–1G₂₁G₁₁^–1x_01с^Т–x_01сG₁₁^–1G₁₂B^–1x_02с^Т

–x_02сB^–1G₂₁G₁₁^–1x_01с^Т+x_02сB^–1x_02с^Т)

=s²(x_01сАB^–1А^Тx_01с^Т–x_01сАB^–1x_02с^Т–x_02сB^–1А^Тx_01с^Т+x_02сB^–1x_02с^Т)

=s²(x_02сB^–1x_02с^Т–x_02сB^–1А^Тx_01с^Т–x_01сАB^–1x_02с^Т+x_01сАB^–1А^Тx_01с^Т)

=s²[x_02сB^–1(x_02с^Т–А^Тx_01с^Т)–x_01сАB^–1(x_02с^Т–А^Тx_01с^Т)]

=s²(x_02с–x_01сА)B^–1(x_02с–x_01сА)^Т

≥0,

так как матрица B^–1 положительно определённая. Это следует из того, что матрицу В можно представить в виде

B=G₂₂–G₂₁G₁₁^–1G₁₂

=X₂^ТX₂–X₂^ТX₁(X₁^ТX₁)^–1X₁^ТX₂

=X₂^Т[I–X₁(X₁^ТX₁)^–1X₁^Т]X₂

В ней матрица I–X₁(X₁^ТX₁)^–1X₁^Т идемпотентная. Следовательно, по теореме П.6.3 матрица В положительно определённая и по теореме П.6.5 матрица B^–1 положительно определённая. Таким образом, D(x_0с)≥D(x_01с).

□

По пункту 1 теоремы 8.2.3 дисперсии D() результатов оценки параметров адекватной модели больше дисперсий D() результатов оценки параметров неадекватной модели. Таким образом, для неадекватной модели дисперсии результатов оценки её параметров меньше, но сами результаты оценки параметров получаются смещёнными. С другой стороны, постулирование адекватной модели увеличивает дисперсии результатов оценки её параметров. По пункту 2 теоремы 8.2.3 дисперсия D() результата оценки ожидаемого значения переменной отклика для адекватной модели больше дисперсии D() результата этой оценки для неадекватной модели. И снова, для неадекватной модели дисперсия результата оценки ожидаемого значения переменной отклика меньше, но этот результат получается смещённым, а для адекватной модели дисперсия результата оценки ожидаемого значения переменной отклика получается больше.

Смещение результата оценки дисперсии

Рассмотрим теперь результаты оценки дисперсии для адекватной и неадекватной моделей. В случае адекватной модели у=Xb+e=X₁b₁+X₂b₂+e оценка дисперсии делается по формуле (7.3.7)

s²=(у–Х)^Т(у–Х)/(n–р).

По теореме 7.3.2 имеем E(s²)=s². Математическое ожидание результата оценки дисперсии в случае неадекватной модели даётся в следующей теореме.

Теорема 8.2.4. Если модель у=Xb+e=X₁b₁+X₂b₂+e адекватная, то для неадекватной модели у=X₁b₁^*+e^*, где X₁ – матрица размеров nxт и т<р, результат оценки дисперсии по формуле

s₁²=(у–Х₁)^Т(у–Х₁)/(n–т)                                 (8.2.8)

имеет математическое ожидание

E(s₁²)=s²+b₂^ТX₂^Т[I–X₁(X₁^ТX₁)^–1X₁^Т]X₂b₂/(n–т).    (8.2.9)

Доказательство: Запишем числитель правой части выражения (8.2.8) в виде квадратичной формы

S_E*=у^Ту–Х₁^Ту

=у^Ту–y^ТX₁(X₁^ТX₁)^–1Х₁^Ту

=y^Т[I–X₁(X₁^ТX₁)^–1Х₁^Т]y.

Для нахождения математического ожидания E(s₁²) оценочной дисперсии необходимо иметь математическое ожидание E(S_E*) квадратичной формы. По допущению Е(у)=Xb и теореме 5.2.1 имеем,

E(S_E*)=след{[I–X₁(X₁^ТX₁)^–1Х₁^Т]s²I}+b^ТX^Т[I–X₁(X₁^ТX₁)^–1Х₁^Т]Xb

=(n–т)s²+b^ТX^Т[I–X₁(X₁^ТX₁)^–1Х₁^Т]Xb,

где след[I–X₁(X₁^ТX₁)^–1Х₁^Т]=след(I)–след[Х₁^ТX₁(X₁^ТX₁)^–1]=п–т. В правой части этого выражения можно сделать также следующее преобразование

b^ТX^Т[I–X₁(X₁^ТX₁)^–1Х₁^Т]Xb=(b₁^ТX₁^Т+b₂^ТX₂^Т)[I–X₁(X₁^ТX₁)^–1Х₁^Т](X₁b₁+X₂b₂)

=b₂^ТX₂^Т[I–X₁(X₁^ТX₁)^–1Х₁^Т]X₂b₂.

Отсюда, в силу (8.2.8), получается искомое выражение (8.2.9)

E(s₁²)=E(S_E*)/(n–т)

=s²+b₂^ТX₂^Т[I–X₁(X₁^ТX₁)^–1Х₁^Т]X₂b₂/(n–т).

□

Поскольку квадратичная форма в выражении (8.2.9) неотрицательно определённая и, если постулирована неадекватная модель, а вектор b₂≠0, то результат оценки дисперсии для неадекватной модели смещён в большую сторону.

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

Пример 8.2.2. Положим, что используется модель у_i=b₀^*+b₁^*x_i+e_i^*, когда адекватной является модель у_i=b₀+b₁x_i+b₂x_i²+e_i. В этом случае ,  и s₁² будут смещены на величины, зависящие от выбора значений x_i фактора [см. (8.2.4) и (8.2.9)]. Для модели у_i=b₀^*+b₁^*x_i+e_i^* математическое ожидание ошибки e_i^* не является нулевым:

Е(e_i^*) =Е(у_i–b₀^*–b₁^*x_i)

=Е(у_i)–b₀^*–b₁^*x_i

=b₀+b₁x_i+b₂x_i²–b₀^*–b₁^*x_i

=b₀–b₀^*+(b₁–b₁^*)x_i+b₂x_i²

и поэтому допущение E(e_i^*)=0 не соблюдается.

□

В лекции "Торфяно-сланцевая промышленность" также много полезной информации.

Пример 8.2.3. Адекватной является модель у_i=b₀+b₁x_i+e_i, а используется модель у_i=b₁^*x_i+e_i^*. Для модели у_i=b₁^*x_i+e_i^* результат оценки параметра b₁^* методом наименьших квадратов

=(х^Tх)^–1х^Тy=.                                       (8.2.10)

Тогда, так как адекватной моделью является у_i=b₀+b₁x_i+e_i, то ожидаемое значение оценки  получаем в виде

E() ==

==.                        (8.2.11)

Таким образом, оценка  смещена на величину, зависящую от b₀ и значений x_i фактора.