Оценка параметров модели и дисперсии

6.2. Оценка параметров модели и дисперсии

Результаты оценки параметров

Для оценки параметров q₀, q₁ и дисперсии s² используется случайная выборка у₁, у₂, ..., у_п наблюдений переменных отклика в п опытах эксперимента и соответствующие им числовые значения x₁, x₂, ..., x_п влияющей на отклик переменной x. Для получения результатов  и  оценки параметров используется метод наименьших квадратов, который не требует никаких допущений о распределении случайных переменных модели.

Методом наименьших квадратов ищутся такие результаты  и  оценки параметров q₀ и q₁, которые делают минимальной сумму квадратов разностей

=                                   (6.2.1)

наблюдаемых в п опытах значений у_i переменных отклика и их ожидаемых значений Е(у_i)=q₀+q₁x_i в соответствии с постулируемой моделью (6.1.1). Эта сумма квадратов является функцией параметров q₀ и q₁ и поэтому обозначается S(q₀, q₁). Для нахождения  и , при которых функция S(q₀, q₁)= принимает минимальное значение, берутся её первые частные производные по q₀ и q₁ и приравниваются нулю:

=–2=0,                                    (6.2.2)

=–2=0.                     (6.2.3)

Эти уравнения можно записать соответственно в виде

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

=0 и =0.

Решение первого уравнения относительно q₀ даёт

q₀=

=–q₁,                                                      (6.2.4)

где = и =. Подставляя его во второе уравнение, получаем

=0

и отсюда находим формулу для вычисления

==

=           .           [в силу (1.4.3)]          (6.2.5)

Подставляя в (6.2.4) вместо  его найденное значение , получаем формулу расчёта

=–.                                                 (6.2.6)

Чтобы убедиться, что при результатах  и  оценки параметров функция S(q₀, q₁) принимает минимальное значение, нужно рассмотреть вторые производные этой функции по q₀ и q₁. Вторые производные функции S(q₀, q₁) получаются соответственно =2п и =2. Значения этих производных положительные, поэтому при q₀= и q₁= функция S(q₀, q₁) принимает минимальное значение [Выгодский (2006) стр. 429; Khuri (2003) стр. 114].

Заметим, что при выполнении оценки по формулам (6.2.5) и (6.2.6) не используются допущения раздела 6.1. Получаемые из уравнения =+x_i результаты оценки ожидаемых значений переменных отклика оценивают их значения, полученные в результате эксперимента, и которые на самом деле возможно должны моделироваться нелинейной функцией. Тем не менее, если три допущения раздела 6.1 соблюдаются, то полученные методом наименьших квадратов  и  являются несмещёнными результатами оценки и имеют наименьшую дисперсию среди всех возможных линейных несмещенных результатов оценки.

Уравнения (6.1.1) моделей для п опытов эксперимента можно представить в матричном или векторном видах

у=Xq+e

=q₀1+q₁x+e,                                                  (6.2.7)

где у= - вектор полученных в опытах значений переменных отклика, X= - матрица модели, q= - вектор параметров модели и e=- вектор ошибок. В этом случае оценка параметров модели может быть сделана тоже с использованием матриц. Так, произведение матрицы модели на себя принимает вид X^ТX=, а обратная этого произведения (X^ТX)^–1=. Произведение матрицы модели на вектор значений переменных отклика имеет вид X^Тy=. В теореме 7.2.1 следующей главы доказано, что оценки  и  получаются в результате решения нормальных уравнений X^ТXq=X^Тy по формуле =(X^ТX)^–1X^Тy, то есть

==.    (6.2.8)

Результаты  и  оценки здесь такие же, как по формулам (6.2.5) и (6.2.6). Это можно показать так

===,

что, как и в (6.2.5). Для  начнём с выражения (6.2.6)

=–=–

=–

=.

Второй и четвёртый члены числителя сокращаются и, как в (6.2.7), получаем

=.

Пример 6.2.1. В интегральной микросхеме коэффициент (у) усиления транзистора между эмиттером и коллектором зависит от двух контролируемых в процессе напыления переменных: эмиттерной дозы (x в единицах по 10¹⁴ ионов) и времени (x₁ в мин.) разгонки примеси эмиттера. Здесь рассмотрим часть данных для 10 образцов после напыления при x₁=225, сведённых в таблицу 6.2.1 [Myers и соавт. (2016) стр.17-18].

Таблица 6.2.1. Значения коэффициента (у) усиления транзистора и переменной x

По формулам (6.2.5) и (6.2.6) находятся =2201,7 и =–197,6. Таким образом, уравнение оценки ожидаемых значений коэффициента усиления в зависимости от переменной x получается в виде

=2201,7–197,6x.

На Рис. 6.1 показаны график зависимости  от x в виде прямой линии синим цветом вместе с 10 точками с координатами (x, у). Из рисунка, очевидно, что наклон  является скоростью изменения  при изменении x, а значение  равно значению  при x=0.

Кажущаяся линейной зависимость на Рис. 6.1 не устанавливает причинно следственной зависимости коэффициента усиления от эмиттерной дозы (выводы, которые здесь можно сделать, см. в разделе 6.3). Допущение D(e_i)=s² (постоянной дисперсии) для всех i=1, 2, ..., 10 представляется разумным.

Рис. 6.1. Линия регрессии и данные эксперимента для коэффициента усиления и эмиттерной дозы.

□

Объяснение результатов и условие их раздельной оценки

Обсудим теперь смысл результатов оценки параметров, использованных в уравнении =2201,7–197,6x оценки ожидаемых значений переменных отклика из примера 6.2.1. Заметим, что на Рис.6.1 результат 2201,7 оценки параметра q₀ равный  при x=0, не показана. За пределами интервала от 4,00 до 4,72 единиц переменной x, при которых коэффициент усиления транзистора действительно измерялся, прямая линия на Рис.6.1 не показана, так как нет данных, чтобы проверить её обоснованность за пределами этого интервала. В частности, результат 2201,7 должен расцениваться просто как точка, через которую проходит прямая линия в диапазоне значений переменной x опытов эксперимента.

Второй результат –197,6 оценки параметра q₁ в уравнении оценки ожидаемых значений переменных отклика определяет наклон линии в используемых единицах измерений. Таким образом, по уравнению =2201,7–197,6x величина  уменьшается на 197,6 единиц при изменении на единицу переменной x.

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

Рассмотрим это условие на примере [Box, Draper (2007) стр. 35]. В нём имеется уравнение линейной модели y=θ₁x₁+θ₂x₂+e, где переменной (y) отклика является скорость протекания химической реакции, а x₁ и x₂ - процентные содержания двух катализаторов А и В, влияющих на скорость реакции. Полагается, что выбран такой план эксперимента, в котором значения x₁ и x₂ оказались пропорциональны один другому, так что для каждого опыта x₂=δx₁. Тогда, например, при δ=2 в каждом опыте процентное содержание катализатора В будет в два раза больше, чем катализатора А. В этом случае матрица модели имеет, например, вид X=. Тогда уравнение y=θ₁x₁+θ₂x₂+e модели может быть записано в виде

y=θ₁x₁+θ₂δx₁+e

=(θ₁+δθ₂)x₁+e

=δ^–1(θ₁+δθ₂)x₂+e.

Методом наименьших квадратов могут быть найдены нормальные уравнения для оценки параметров θ₁ и θ₂, но они не обеспечивают единственности их оценки. Эти параметры не могут быть оценены раздельно. В этом случае можно оценить только их линейную комбинацию θ₁+δθ₂. Причина этого в том, что когда x₂=δx₁, то влияние на переменную отклика переменной x₁ (катализатор А) полностью неразличимо от влияния переменной x₂ (катализатор В). Равенство x₂=δx₁ означает, что x₂–δx₁=0. В общем, это происходит всегда, когда линейная зависимость вида α₁x₁+α₂x₂=0 (для данного примера α₁=–δ, α₂=1) связывает линейно зависимые столбцы матрицы X.

В начале раздела П.4 приложения даётся определение линейно независимых векторов и столбцов матрицы. Следовательно, для раздельной оценки параметров модели вектор-столбцы матрицы модели должны быть линейно независимы. Это условие соблюдается для столбцов матрицы модели в примере 6.2.1.

Математические ожидания и дисперсии результатов оценки

Результаты  и  оценки параметров модели являются линейными функциями значений у₁, у₂, ..., у_п переменных отклика. Используя три допущения раздела 6.1, можно получить следующие математические ожидания и дисперсии для  и .

В числителе правой части выражения (6.2.5) имеем

=

=

=.

Тогда формула (6.2.5) принимает вид

=.

Теперь, используя первое допущение Е(у_i)=q₀+q₁x_i раздела 6.1, получаем

Е()==

==

=q₁.                                                                                         (6.2.9)

Математическое ожидание для  находится следующим образом

E()=E(–)=E()–E()

=–q₁=–q₁

=q₀+q₁–q₁=q₀.                                          (6.2.10)

Таким образом, математические ожидания для  и  равны самим оцениваемым параметрам и поэтому их результаты оценки являются несмещёнными. В векторном виде для модели (6.2.7) это можно записать так

E()=q.                                                         (6.2.11)

Дисперсия  определяется с использованием его выражения по формуле (6.2.5), а также второго D(у_i)=s² и третьего C(у_i, у_j)=0 допущений раздела 6.1. В силу (3.2.8), имеем

D()==

=.                                            (6.2.12)

По формуле (6.2.6) для  можно записать

=–=–.

Тогда дисперсия  находится следующим образом

D()=D=

=s²

=s²

=s².                                          (6.2.13)

Обратим внимание, что при нахождении математического ожидания Е() и дисперсии D() рассматриваются случайные изменения  от выборки к выборке значений случайных переменных у_i. Полагается, что n значений x₁, x₂, ..., x_п влияющей на отклик переменной x остаются теми же в опытах эксперимента при получении выборочных значений случайных переменных у_i, так что дисперсии D() и D() постоянны.

Матрица дисперсий и ковариаций вектора  оценки параметров модели находится в виде

D()=E{[–E()][–E()]^T}=(X^ТX)^–1X^ТE{[y–E(y)][y^T–E(y^T)]}X(X^ТX)^–l

= (X^ТX)^–lX^TE(ee^T)X(X^ТX)^–l

= (X^ТX)^–ls².                                                                                      (6.2.14)

Условия оценки параметров с минимальной дисперсией

Из выражения (6.2.12) видно, что дисперсия D() становится минимальной, когда сумма  максимальна. Если значения x_i влияющей на отклик переменной находятся в интервале а≤x_i≤b, где а и b - крайние числа интервала, то при четном п сумма  становится максимальной, если в опытах эксперимента одна половина значений переменной x выбирается равной а, а другая половина равной b. Это можно показать следующим образом.

Пусть р значений переменной x равны а, а оставшиеся п–р значений равны b. Тогда усреднённое значение этой переменной =[рa+(n–р)b]/n. В этом случае  можно представить в виде

=р{a–[рa+(n–р)b]/n}²+(n–р){b–[рa+(n–р)b]/n}²

=р{[na–рa–nb+рb]/n}²+(n–р){[nb–рa–nb+рb]/n}²

=р{[n(a–b)–р(a–b)]/n}²+(n–р)[–р(a–b)/n]²

=р(n–р)²(a–b)²/n²+(n–р)р²(a–b)²/n²

=[р(n–р)²+р²(n–р)](a–b)²/n²

=р(n–р)(n–р+р)(a–b)²/n²

=р(n–р)(a–b)²/n.

Далее возмём первую производную от полученного выражения по р и приравняем её нулю

=(n–2р)(a–b)²/n=0.

Отсюда получаем р=п/2. А если взять вторую производную, то получаем

=–2(a–b)²/n.

Вторая производная получается отрицательной, следовательно при р=п/2 достигается максимум суммы .

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

Кроме этого, в силу (6.2.13), очевидно, что дисперсия D() становится минимальной, когда =0. Для этого при обработке результатов двухуровневых факторных экспериментов каждый влияющий на отклик фактор нормируется по формуле (2.6.4) чтобы усреднённое нормированного фактора было равно нулю.

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

При соблюдении первого допущения раздела 6.1, математическое ожидание вектора случайных переменных отклика модели (6.2.7) имеет вид Е(у)=Xq. Если в это выражение вместо вектора q подставить вектор =(X^ТX)^–1X^Тy его оценки, то получается вектор оценки ожидаемых значений случайных переменных ≡=X. Разность векторов у и  даёт вектор остатков или остаточных ошибок

е=у–=у–X(X^ТX)^–1X^Тy

=[I–X(X^ТX)^–1X^Т]y                                         (6.2.15)

Произведение этого вектора и матрицы X даёт нулевой вектор

X^Те=X^Т[I–X(X^ТX)^–1X^Т]y=[X^Т–X^Т]y=0.

По определению произведения матрицы на вектор это значит, что произведение вектора е с любым вектор-столбцом матрицы X даёт нулевой результат.

Если при планировании эксперимента векторы столбцы матрицы X не сделаны ортогональными, то на практике они обычно получаются не ортогональными. Это видно из примера 6.2.1, где первый и второй столбцы не ортогональны, то есть 1^Тx≠0. Однако можно найти составляющий вектор x_о вектора x, который ортогонален вектору 1, и переписать функцию модели с использованием ортогональных векторов. Для нахождения вектора x_о, являющегося составляющим вектора x и ортогонального вектору 1, воспользуемся тем, что вектор е остатков ортогонален векторам 1 и x. Временно считая x вектором переменных отклика и 1 вектором значений влияющей на отклик переменной, методом наименьших квадратов получим оценку вектора ожидаемых значений отклика x в виде =4,352x1. В силу (6.2.15), вектор остатков для данных примера 6.2.1находится так x_о= x–= x–4.352x1, что в численном выражении имеет вид

x_о^Т=[–0,352  0,248  –0,152  –0,252  0,248  –0,052  –0,352  0,348  –0,052  0,368].

Теперь перепишем функцию модели (6.2.7) в виде

Е(у)=θ₀1+θ₁x+4,352θ₁1–4,352θ₁1

=(θ₀+4,352θ₁)1+θ₁(x–4,352x1),

откуда получаем

Е(у)=θ1+θ₁x_о,

где θ=θ₀+4,352θ₁. Для такой функции модели имеем матрицу X_о=[1, x_о] с использованием которой вычисляем

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

И последним найдём также функцию Е(у)=θ₀1 модели, где параметр θ₀ тоже оценивается методом наименьших квадратов по формуле =(1^Т1)^–11^Тy=1341,5.

Теперь можно сравнить три выражения, полученные для оценки вектора ожидаемых значений переменных отклика:

Ø Для модели с одним параметром                                            =1341,5x1

Ø Для модели с двумя параметрами                                          =2201,7x1–197,6x

Ø Для модели с ортогональными столбцами её матрицы         =1341,5x1–197,6x_о

Из сравнения делаем следующие заключения:

1. Так как векторы 1 и x_о ортогональны, то коэффициент перед вектором 1 в модели с ортогональными столбцами её матрицы и двумя параметрами является тем же, что и коэффициент перед 1 в модели с одним параметром.

2. Коэффициент перед вектором x_о в модели с ортогональными столбцами её матрицы является тем же, что и коэффициент перед вектором x в модели с двумя параметрами и не ортогональными столбцами её матрицы.

Полученные выше выражения оценки представлены графически на Рис.6.2.1.

Рис. 6.2.1. Плоскость оценки векторов ожидаемых значений переменных отклика тремя моделями.

Вектор 1 и ортогональный ему вектор x_о, как и векторы 1 и x, могут использоваться для задания плоскости, где расположен вектор  оценки ожидаемых значений переменных отклика. При не ортогональности базисных векторов 1 и x конец вектора  имеет координаты 2201,7x1 и –197,6x. При ортогональности базисных векторов 1 и x_о, тот же конец вектора  имеет координаты 1341,5x1 и –197,6x_о.

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

Оценка дисперсии

Методом наименьших квадратов невозможно оценить дисперсию D(у_i)=s². Нахождение минимума функции S(q₀, q₁) дает только результаты  и  оценки параметров модели. Для оценки дисперсии используется выражение (3.2.2), то есть,

D(у_i)=E[у_i–E(у_i)]².

В опытах эксперимента по второму допущению раздела 6.1 дисперсия s² считается одинаковой для всех переменных у_i (i =1, 2, ..., п) отклика. Используя обозначение  для результата оценки ожидаемого значения E(у_i) случайной переменной отклика, дисперсия s² оценивается выражением

s²===,                     (6.2.16)

где  и  даются формулами (6.2.5) и (6.2.6), а S_E=. Результат разности у_i– = называют остатком от значения у_i случайной переменной отклика или остаточной ошибкой и S_E называют суммой квадратов остатков или суммой квадратов остаточных ошибок.

Так как знаменатель в выражении (6.2.16) принят равным n–2, то статистика s² является несмещенным результатом оценки дисперсии s². Это можно показать следующим образом. Сумму S_E квадратов остатков можно записать в виде

S_E ====

=–2+

=–2+

=–

=–.

Усреднённое значение выборки переменных отклика принимает вид

===q₀+q₁+=q₀+q₁+

Математическое ожидание суммы  получается

Е[]=Е[]=Е{}

=Е{}

=

=,

так как по аналогии с (5.1.2), имеем

=–n

=e^ТIe–e^Т(Е/n)e

=e^Т(I–Е/n)e

и по теореме 5.2.1 математическое ожидание этой суммы

=след[(I–Е/n)( s²I)]+е1^Т(I–Е/n)е1

=s²след(I–Е/n)+е²(1^Т1–1^Т11^Т1/n)  [так как Е=11^Т]

=s²(n–n/n)+е²(n–n²/n)                     [так как 1^Т1=n]

=s²(n–1)+0.

На основе выражения (3.2.4) можно получить Е()=D()+[Е()]²=+q₁². Следовательно, математическое ожидание суммы S_E

E(S_E) =––q₁²

=(n–2)s²

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

E(s²)===s²,                                  (6.2.17)

что и указывает на несмещённость статистики s².

Интуитивно, в выражении (6.2.16) делителем взята разность n–2 вместо n–1, как в формуле (1.4.3) выборочной дисперсии, так как функция модели Е(у_i) =q₀+q₁x_i имеет два оцениваемых параметра и поэтому оценка ожидаемого значения E(у_i) должна быть лучше, чем усреднённого . Таким образом, можно ожидать, что S_E= меньше суммы . В самом деле, используя формулы (6.2.5) и (6.2.6), можно записать числитель выражения (6.2.16) в виде

S_E==

==

=–2+

=–2+

Лекция "20.1 Борьба за власть в руководстве страны после смерти И.В. Сталина" также может быть Вам полезна.

=–,                                                                   (6.2.18)

чем и подтверждается, что  меньше суммы .

Пример 6.2.2. Для данных примера 6.2.1 оценим дисперсию по формуле (6.2.16)

s²==2,6149х10⁴.

□

Оценка дисперсии по формуле (6.2.16) зависит от числа параметров модели. Но если в эксперименте поставлены повторные опыты, то есть опыты при одних и тех же значениях переменной x, то можно оценить дисперсию независимо от числа параметров модели. Такая оценка дисперсии получена в конце следующего раздела, но вначале рассмотрим проверку гипотез о параметрах модели.