Параметрические методы оценивания параметра 0 предполагают соответствие вида предполагаемого распределения g(x, 0) неизвестному истинному. Получаемая при этом по выборке независимых значений

максимально правдоподобная оценка в виде векторазначений параметров (аргумента), обеспечивающего максимальное значение функции правдоподобия

обладает минимально возможной дисперсией, т.е. является эффективной оценкой параметра 0 при условии равенства предполагаемого распределения истинному. Отличиеотобусловливает снижение эффективности оценок. Это отличие, в частности, может быть обусловлено присутствием в выборке "посторонних включений" – наблюдений из совокупности, описываемой другими законами распределения. Оценки максимального правдоподобия могут быть менее эффективными по сравнению с оценками, не лучшими в идеальных условиях, но выигрывающими в эффективности в реальных ситуациях статистического оценивания. Такие оценки благодаря Хыоберу получили название робастных.

Мерой относительной эффективности оценок выступает отношение их погрешностей. В качестве погрешности как меры точности измерения широко используются средняя абсолютная ошибка

и средняя квадратическая ошибка

где – г-е значение случайной величины; – среднее значение случайной величины.

Выбор этих и других мер погрешности относится к категории предпочтений того или иного критерия точности оценивания. При выборе достаточно общего байесовского критерия минимума среднего риска ошибок мера погрешности (7.6) соответствует линейной функции стоимости ошибки, а мера (7.7) – квадратичной.

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

Задача оценивания параметра сдвига симметричного распределения является одной из важнейших статистических задач, имеющих прикладное значение. Примерами таких распределений могут служить распределение Лапласа

и нормальное распределение

где– параметр сдвига распределения относительно нуля, определяющий положение центра симметрии.

Зависимость симметричного распределения от параметра сдвига можно представить в виде

Оценку максимального правдоподобия (7.5) для параметра сдвига для случая нормального распределения признака (7.8) можно получить путем дифференцирования плотности вероятности функции правдоподобия или монотонно связанного с ней ее логарифма (что намного удобнее)

по параметру сдвига р и приравнивания результата к нулю. В результате для распределения (7.8) с точностью до постоянного множителя, не зависящего от р, получаем уравнение

левая часть которого представляет собой сумму так называемых оценочных функций (score functions)

Оценочная функция может иметь вид, отличный от выражения (7.10).

Оценочную функцию можно использовать для определения весовой функции , если она существует:

Выразив оценочную функцию через весовую из формулы (7.11) и подставив ее в уравнение (7.9), убедимся в том, что весовая функция соответствует своему названию в смысле определения веса каждого наблюдения в формировании оценки параметра сдвига:

Для оценочной функции (7.10) все наблюдения х, равноправны в формировании оценки Д. Для случая отсутствия в выборке "посторонних" объектов это логично. Однако наличие аномальных наблюдений может существенно исказить оценку параметра сдвига нормальной совокупности. Избежать этого можно путем выявления аномалий и их исключения из выборки подобно извлечению одного или нескольких лезвий из складного ножа (jackknife ). Этот принцип лежит в основе джекнайф- процедур оценивания. Их недостатком является отсечение в явном или неявном виде не истинно аномальных наблюдений, а наблюдений, признаваемых аномальными или "подозрительными" на основе выбранного решающего правила, что может привести к искажениям и информационным потерям.

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

Для нормально распределенной генеральной совокупности с плотностью вероятности (7.8) средняя арифметическая величина является несмещенной, состоятельной и эффективной оценкой параметра сдвига в виде математического ожидания р. Однако эффективность ее падает с утяжелением "хвостов" распределения, т.е. наличием достаточно большого числа наблюдений, значительно удаленных от среднего значения. Дж. Тыоки исследовал влияние выбросов на эффективность оценки генерального среднего (параметра сдвига). В качестве модели распределения, полагаемого при оценивании нормальным, он использовал смесь двух нормальных распределений, в которой к основному распределению добавлено с весомраспределение с тем же параметром сдвига, но втрое большей дисперсией :

Величина е определяет вероятность попадания аномальных наблюдений в нормальную выборку с единичной дисперсией, и она, как правило, невелика. "Гьюки показал, что при таком засорении оценки методом максимального правдоподобия неустойчивы: их эффективность резко снижается и оказывается худшей, чем оценка усеченного среднего

где – наблюдения , для которых модуль отклонения от р меньше некоторого порога k. Функция веса всех наблюдений при определении среднего значения приведена на рис. 7.6.

Рис. 7.6.

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

Хьюбер в качестве функции, описывающей "засорения", рассматривал произвольную симметричную функциюс нулевым математическим ожиданием. Оценочную функцию необходимо выбрать таким образом, чтобы при наихудшем засорении оценка обладала минимальным средним квадратом отклонения от истинного значения параметра сдвига:

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

где – производная оценочной функции по параметру сдвига ц.

Правая часть этого равенства представляет собой отношение средних значений оценочной функции и ее производной.

Асимптотическая дисперсия оценкисоставит

Согласно теореме Хыобера (6.2.5)

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

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

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

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

у =Xq +e

Q 0 1 +q 1 x +e , (6.2.7)

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

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

===,

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

=–

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

=.

Пример 6.2.1 . В интегральной микросхеме коэффициент (у) усиления транзистора между эмиттером и коллектором зависит от двух контролируемых в процессе напыления переменных: эмиттерной дозы (x в единицах по 10 14 ионов) и времени (x 1 в мин.) разгонки примеси эмиттера. Здесь рассмотрим часть данных для 10 образцов после напыления при x 1 =225, сведённых в таблицу 6.2.1 .

Таблица 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 2 (постоянной дисперсии) для всех i =1, 2, ..., 10 представляется разумным.

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

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

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

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

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

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

y=θ 1 x 1 +θ 2 δ x 1 +e

=(θ 1 +δ θ 2)x 1 +e

=δ –1 (θ 1 +δ θ 2)x 2 +e.

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

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

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

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

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

=

=

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

=.

Теперь, используя первое допущение Е (у i )=q 0 +q 1 x i раздела 6.1, получаем

Е ()==

==

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

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

=–q 1 =–q 1

Q 0 +q 1 –q 1 =q 0 . (6.2.10)

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

E ()=q . (6.2.11)

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

D ()==

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

=–=–.

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

D ()=D =

=s 2

S 2

S 2 . (6.2.13)

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

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

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

= (X Т X ) –l X T E (ee T)X (X Т X ) –l

= (X Т X ) –l s 2 . (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 } 2 +(n –р ){b –[рa +(n –р )b ]/n } 2

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

X Т е =X Т [I –X (X Т X ) –1 X Т ]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) в виде

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

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

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

Е (у )=θ1 +θ 1 x о,

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

(X о Т X о) –1 =, X о Т y = и ==.

И последним найдём также функцию Е (у )=θ 0 1 модели, где параметр θ 0 тоже оценивается методом наименьших квадратов по формуле =(1 Т 1 ) –1 1 Т 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 2 . Нахождение минимума функции S (q 0 , q 1) дает только результаты и оценки параметров модели. Для оценки дисперсии используется выражение (3.2.2), то есть,

D (у i )=E [у i –E (у i )] 2 .

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