Ортогональное разложение функции модели

8.3. Ортогональное разложение функции модели

Пусть в модели (8.1) с р факторами, значения которых представлены столбцами матрицы Х, вектор b параметров разделён на два подвектора b₁ и b₂ параметров для р₁ и р₂ факторов со значениями в столбцах соответственно матриц Х₁ и Х₂, так что функция модели Xb представляется в виде

Xb=Х₁b₁+Х₂b₂.                                            (8.3.1)

В ней более простая часть Х₁b₁ может быть использована для моделирования изменений переменных отклика в эксперименте, а часть Х₂b₂ представляет новые члены модели, которые должны быть добавлены, если модель с функцией Х₁b₁ окажется неадекватной. Далее положим, что матрицы Х₁ и Х₂ между собой ортогональны, так что Х₁^ТХ₂=О. Тогда формулу =(Х^ТХ)^–1Х^Тy, дающую вектор  оценки вектора b методом наименьших квадратов, тоже можно разделить на две части

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

Более того, получаемая в результате регрессии сумма квадратов (–b)^ТХ^ТХ(–b) для данной выражением (8.3.1) функции модели может быть записана тоже в разделённом виде. Это получается, если вектор (b–) и матрицу Х представить так

(b–)^Т=[(b₁–₁), (b₂–₂)] и Х=[Х₁, Х₂].

Тогда, так как Х^ТХ=[Х₁, Х₂]==, то

(b–)^ТХ^ТХ(b–)=[(b₁–₁), (b₂–₂)]

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

=(b₁–₁)^ТХ₁^ТХ₁(b₁–₁)+(b₂–₂)^ТХ₂^ТХ₂(b₂–₂)

и соответствующий дисперсионный анализ модели с разделённой функцией показан в таблице 8.3.1.

Таблица 8.3.1. Дисперсионный анализ модели с разделенной на взаимно ортогональные части функцией, то есть, при Х₁^ТХ₂=О

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

Теорема 8.3. Если X₁^ТX₂=O, то результат оценки вектора b₁ параметров адекватной модели у=X₁b₁+X₂b₂+e является таким же, как и результат оценки вектора b₁^* параметров неадекватной модели у=X₁b₁^*+e^*.

Доказательство: Вектор =(X₁^ТX₁)^–1X₁^Тy является вектором оценки вектора b₁^* параметров методом наименьших квадратов. Для оценки вектора b₁ адекватной модели разделим вектор =(X^ТX)^–1X^Тy его оценки чтобы получить

=.

Используя обозначения в доказательстве теоремы 8.2.3, это выражение преобразуется к виду

=

=.

В силу (П.5.6), получаем

=G¹¹X₁^Ту+G¹²X₂^Ту

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

где B=G₂₂–G₂₁G₁₁^–1G₁₂. Если G₁₂=X₁^ТX₂=O, то  сводится к

=G₁₁^–1X₁^Ту

=(X₁^ТX₁)^–1X₁^Тy,

что является таким же, как для .

□

Ортогонализация двух частей матрицы модели

В случае, когда подматрицы Х₁ и Х₂ разделённой на части матрицы Х=[Х₁, Х₂] модели взаимно не ортогональны, что получается, если значения факторов при планировании эксперимента выбираются произвольно, то для того чтобы матрицы Х₁ и Х₂ были между собой ортогональны необходимо осуществить их взаимную ортогонализацию.

Рассмотрим модель у=X₁b₁+e. Оценка параметров этой модели выполняется по формуле ₁=(Х₁^ТХ₁)^–1Х₁^Тy и вектор оценки ожидаемых значений переменных отклика находится в виде =Х₁₁. Следующим из нормальных уравнений метода наименьших квадратов важным следствием оценки является то, что вектор остатков

y–=y–Х₁₁

ортогонален каждому столбцу матрицы Х₁. Это даёт общий метод нахождения такого составляющего вектора у_о для вектора y, который ортогонален имеющему р₁ измерений векторному пространству Ст(Х₁) чьими базисными векторами являются столбцы матрицы Х₁. Ортогональный пространству Ст(Х₁) вектор y_о находится так

y_о=y–=y–Х₁₁=(I–Р₁)y,

где матрица Р₁=Х₁(Х₁^ТХ₁)^–1Х₁^Т.

Теперь, вместо вектора y рассмотрим матрицу Х₂ в той же роли. В этом случае модель принимает матричный вид Х₂=Х₁A+E, где A - матрица параметров и E - матрица ошибок. Считая Е(E)=О, методом наименьших квадратов результат оценки матрицы A параметров модели получается по формуле А=(Х₁^ТХ₁)^–1Х₁^ТХ₂, а результат оценки ожидаемых значений матрицы Х₂ находится в виде =Х₁А. Тогда, искомая Х_О составляющая матрицы Х₂ ортогональная матрице Х₁ находится в виде

Х_О=Х₂–.

Матрица Х_О содержит р₂ столбцов, каждый из которых ортогонален каждому из р₁ столбцов матрицы Х₁ и может быть получена также по формуле

Х_О=Х₂–Х₁А,

где А - матрица оценки параметров размеров р₁хр₂, полученная посредством регрессии каждого столбца матрицы Х₂ во все столбцы матрицы Х₁. Как было отмечено ранее, матрицу А называют также матрицей совмещения или смещения.

В случае, когда матрицы Х₁ и Х₂ между собой не ортогональны, модель с функцией (8.3.1) можно преобразовать следующим образом

y=Х₁b₁+Х₁Аb₂–Х₁Аb₂+Х₂b₂+ε

=Х₁(b₁+Аb₂)+(Х₂–Х₁А)b₂+ε

=Х₁b+Х_Оb₂+ε,

где b=b₁+Аb₂. Здесь матрицы Х₁ и Х_О между собой ортогональны. Соответствующий дисперсионный анализ модели с ортогональными между собой частями её матрицы показан в таблице 8.3.2.

Таблица 8.3.2. Дисперсионный анализ модели с ортогональными между собой частями её матрицы

В результате ортогонализации между собой частей матрицы модели должны быть отмечены следующие пункты:

1. Полученный для модели y=Х₁b+Х_Оb₂+ε вектор  оценки такой же, как и для модели y=Х₁b+ε₁, которая содержит только матрицу Х₁.

2. Если модель y=Х₁b+ε₁ неадекватна и необходимо включить в неё дополнительные факторы, то вектор ₂ оценки дополнительных параметров будет тем же, что получается методом наименьших квадратов для модели с функцией Xb=Х₁b₁+Х₂b₂ и ортогональными между собой матрицами Х₁ и Х₂, или для модели y=Х_Оb₂+ε₂.

(Выполнение пунктов 1 и 2 обеспечивается тем, что матрицы Х₁ и Х_О ортогональны между собой, то есть Х₁^ТХ_О=О, но сами по себе они могут быть не ортогональными.)

3. Если необходимы дополнительные факторы для превращения неадекватной модели y=Х₁b₁+ε в адекватную, где b₂≠0, то вектор  будет результатом смещённой оценки вектора b₁, а вернее результатом оценки совмещения b₁+Аb₂.

4. Необходимость присутствия части Х₂b₂ в функции модели будет показана величиной суммы квадратов “Дополнительной для Х₂” с р₂ степенями свободы, когда в этой сумме квадратов вектор b₂=0.

5. “Дополнительная для Х₂” сумма квадратов может быть получена посредством использования сначала более простой модели y=Х₁b₁+ε, а затем более разработанной модели y=Х₁b₁+Х₂b₂+ε и вычислением разности сумм квадратов остатков для этих моделей.

Ортогонализация столбцов матрицы модели

Кроме взаимной ортогонализации двух частей матрицы модели, состоящих из нескольких столбцов, можно сделать ортогональными между собой все столбцы этой матрицы. Есть процедура последовательного преобразования столбцов матрицы модели, так что каждый новый преобразованный столбец ортогонален ранее преобразованным столбцам [Draper, Smith (1998) стр.383]. Преобразование столбцов осуществляется по формуле

x_jо=(I–X(X^TX)^–1X^T)x_j                                                           (8.3.2)

где X - матрица уже преобразованных столбцов, x_j - следующий подлежащий преобразованию вектор-столбец матрицы и х_jо - преобразованный вектор-столбец ортогональный всем столбцам матрицы Х. В этой процедуре вектор х_jо в действительности является вектором остатков для вектора x_j после того как он был спроецирован в пространство вектор-столбцов матрицы Х.

Исходной матрицей этого преобразования может быть матрица, состоящая из столбца 1 значений фактора х₀ и столбца нормированных значений фактора х₁, так как само нормирование делает его ортогональным столбцу 1. Если существует линейная зависимость между столбцами формируемой матрицы Х и вновь присоединяемым столбцом, то этот столбец получается состоящим из очень малых чисел или возможно нулей. Такой столбец присоединять нет смысла, так как фактор с нулевыми или почти нулевыми значениями оказывает приблизительно нулевое воздействие на переменную отклика. Добавление такого столбца приводит к линейной или почти линейной зависимости столбцов матрицы Х. Это, в свою очередь, приводит к вырожденной или плохо обусловленной матрице Х^ТХ и её обратную (Х^ТХ)^–1 нельзя получить.

Пример 8.3. Выполним нормирование рассматриваемых в примере 7.1  факторов по формуле x=(ξ–)/S, где = и S=. Тогда данные таблицы 7.1 с нормированными значениями факторов получаются, как показано в таблице 8.3.3.

В результате такого нормирования факторов x₁ и x₂ получены нормированные их значения. Векторы нормированных значений этих факторов можно считать ортогональными вектору 1 матрицы модели, так как 1^Тx₁=1,477х10^–14 и 1^Тx₂=3,889х10^–14. Однако эти векторы не ортогональны между собой, так как x₁^Тx₂=1,22. Чтобы сделать их ортогональными воспользуемся формулой (8.3.2). В этой формуле матрица Х=[1, x₁], а преобразуемым вектором является x₂. Полученные в результате преобразования значения элементов преобразованного вектора х_2о показаны в таблице 8.3.3. Этот вектор ортогонален векторам 1 и x₁, так как 1^Тx_2о=0 и x₁^Тx_2о=0. Заметим, что если необходима строгая ортогональность столбцов матрицы модели, то начинать надо с ортогонализации векторов 1 и x₁. При этом исходной матрицей будет Х=1.

Таблица 8.3.3. Значения нормированных факторов, коэффициента (у) усиления и оценок их ожидаемых значений, остатков и других данных анализа

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

=1250+172,5х₁–54,3х₂,

а для модели с ортогональными между собой векторами нормированных факторов - по формуле

=1250+167,8х₁–54,3х₂,

что подтверждает результаты, полученные в разделе 6.2 в процессе ортогонализации столбцов матрицы модели и в разделе 8.3 в процессе ортогонализации двух частей матрицы модели.

□

Оценка переменных отклика при ортогонализации матрицы модели

В разделе 7.5 показано, что при нормировании факторов модели вектор оценки ожидаемых значений переменных отклика остаётся тем же, что и для модели с ненормированными факторами. Покажем, что это справедливо и при ортогонализации столбцов матрицы модели у=Хb+e.

Начиная со столбца x₁ ортогональное преобразование столбцов матрицы Х=[1, x₁, x₂, …, x_p–1] выполняется по формулам

u₁= [I–1(1^T1)^–11^T]x₁

u₂= [I–X₁(X₁^TX₁)^–1X₁^T]x₂, где X₁= [1, u₁]

u₃= [I–X₂(X₂^TX₂)^–1X₂^T]x₃, где X₂= [1, u₁, u₂]

…

u_p–1= [I–X_p–2(X_p–2^TX_p–2)^–1X_p–2^T]x_p–1, где X_p–2=[1, u₁, u₂,..., u_p–2]

Следовательно, матрица U модели с взаимно ортогональными столбцами имеет вид

U=[1, u₁, u₂, …, u_p–1]

={1, [I–1(1^T1)^–11^T]x₁, [I–X₁(X₁^TX₁)^–1X₁^T]x₂, …, [I–X_p–2(X_p–2^TX_p–2)^–1X_p–2^T]x_p–1}

Она получается в результате произведения разделённой по столбцам матрицы Х=[1, x₁, x₂, …, x_p–1] на диагональную матрицу

D=диаг{I, [I–1(1^T1)^–11^T], [I–X₁(X₁^TX₁)^–1X₁^T],..., [I–X_p–2(X_p–2^TX_p–2)^–1X_p–2^T]}.

Вектор оценки параметров модели, имеющей матрицу U, находится по формуле =(U^ТU)^–1U^Тy. Подставляя в эту формулу ХD вместо U получаем

=[(ХD)^Т(ХD)]^–1(ХD)^Тy

=(ХD)^–1[(ХD)^Т]^–1(ХD)^Тy

=D^–1Х^–1(D^ТХ^Т)^–1(ХD)^Тy

=D^–1Х^–1(Х^Т)^–1(D^Т)^–1D^ТХ^Тy [в силу (П.5.5)]

=D^–1(Х^ТХ)^–1Х^Тy

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

=D^–1,

где  - вектор оценки параметров, вычисляемый по формуле (8.1.4).

Преобразование i-й строки матрицы X в i-ю строку матрицы U выполняется по формуле u_iс=x_iсD, а преобразование u_iс в x_iс в виде u_iсD^–1=x_iс. Следовательно, результат  оценки ожидаемого значения i-й переменной отклика находится так

=u_iс=u_iсD^–1

=x_iс.

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