Оценивание

Для получения асимптотического распределения статистики $$Y_n$$ вида (2) можно применить метод линеаризации к асимптотически нормальному вектору выборочных моментов $$(M_{n1},M_{n2},...,M_{nq})$$ и функции $$g$$ из формулы (2).

В силу многомерной центральной предельной теоремы (см. лекцию 4) указанная асимптотическая нормальность имеет место, если, например,$$M|x_i|^{2q+1}<+\infty.$$

Это условие выполнено, в частности, для результатов измерений, распределения которых сосредоточены на ограниченных сверху и снизу интервалах.

При реализации намеченного плана для применения формулы (4) необходимо использовать асимптотические дисперсии и ковариации выборочных моментов, т.е. величины, обозначенные в формуле (4) как $$c_{rs}$$. Эти величины имеют вид [, с.388]:$$\begin{gathered} c_{rr}=\mu_{2r}-\mu_r^2-2r\mu_{r-1}\mu_{r+1}+r^2\mu_{r-1}^2\mu_2, \\ c_{rs}=\mu_{r+s}-\mu_r\mu_s+rs\mu_2\mu_{r-1}\mu_{s-1}-r\mu_{r-1}\mu_{s+1}-s\mu_{r+1}\mu_{s-1},r,s=1,2,...,\mu_0=0 \end{gathered}$$

Здесь $$\mu_r$$ - теоретический центральный момент порядка $$r$$, т.е.$$\mu_r=M(x_i-M(x_i))^r,r=1,2,...$$

Таким образом, для получения асимптотического распределения случайной величины $$Y_n$$ вида (2) достаточно знать теоретические центральные моменты результатов наблюдений и вид функции $$g$$. Отметим, что асимптотическим смещением оценок в рассматриваемом случае можно пренебречь, поскольку его вклад в средний квадрат ошибки статистической оценки - бесконечно малая величина более высокого порядка по сравнению с асимптотической дисперсией.

Однако моменты неизвестны. Их приходится оценивать. В соответствии с теоремами о наследовании сходимости для нахождения асимптотического распределения функции от выборочных моментов можно воспользоваться не теоретическими моментами, а их состоятельными оценками. Эти оценки можно получить разными способами. Можно непосредственно применить формулы (5), заменив теоретические моменты выборочными. Можно выразить моменты через параметры рассматриваемого распределения. Можно применять более сложные процедуры, например, на основе непараметрических устойчивых (робастных) оценок моментов типа урезанных средних Пуанкаре и др. (в первой в России книге по общей теории устойчивости [] проблематика робастных оценок рассмотрена в лекции 2).

Для оценивания параметров гамма-распределения воспользуемся известной формулой [, с.42], согласно которой для случайной величины $$X$$, имеющей гамма-распределение с параметрами формы $$a$$, масштаба $$b=1$$ и сдвига $$c=0$$,$$M(X^m)=\frac{\Gamma(a+m)}{\Gamma(a)}=a(a+1)...(a+m-1), m=1,2,...$$

Следовательно, $$M(X) = a, M(X^2) = a(a+1), D(X) = M(X^2) - (M(X))^2 = a(a+1) - a^2 = a$$. Найдем третий центральный момент $$M(X - M(X))^3$$. Справедливо равенство$$M(X - M(X))^3 = M(X^3) - 3 M(X^2) M(X) + 3 M(X) (M(X))^2 - (M(X))^3.$$

Из равенства (6) вытекает, что$$M(X - M(X))^3 = a(a+1)(a+2) - 3 a (a+1) a + 3 a a^2 - a^3 = 2a.$$

Если $$Y$$ - случайная величина, имеющая гамма-распределение с произвольными параметрами формы $$a$$, масштаба $$b$$ и сдвига $$c$$, то $$Y = bX + c$$. Следовательно,$$M(Y) = ab+c, D(Y) = ab^2, M(Y - M(Y))^3 = 2 a b^3.$$

Пример 2. Оценивание методом моментов параметров гамма-распределения в случае трех неизвестных параметров (строка 7 табл.6.1).

В соответствии с проведенными выше рассуждениями для оценивания трех параметров достаточно использовать три выборочных момента - выборочное среднее арифметическое$$\overline{x}=\frac{x_1+x_2+...+x_n}{n},$$ выборочную дисперсию$$s^2=\frac{1}{n-1}\sum_{i=1}^n(x_i-\overline{x})^2$$ и выборочный третий центральный момент$$m_3=\frac{1}{n}\sum_{i=1}^n(x_i-\overline{x})^3.$$

Приравнивая теоретические моменты, выраженные через параметры распределения, и выборочные моменты, получаем систему уравнений метода моментов:$$ab+c=\overline{x},ab^2=s^2,2ab^3=m_3.$$

Решая эту систему, находим оценки метода моментов. Подставляя второе уравнение в третье, получаем оценку метода моментов для параметра сдвига: 2s^2b=m_3,b*=\frac12\frac{m_3}{s^2}.

Подставляя эту оценку во второе уравнение, находим оценку метода моментов для параметра формы:$$a(b*)^2=a\left(\frac12\frac{m_3}{s^2}\right)^2=\frac{a}{4}\frac{m_3^2}{s^4}=s^2,\quad a^*=4\frac{s^6}{m_3^2}.$$

Наконец, из первого уравнения находим оценку для параметра сдвига:$$c*=\overline{x}-a*b*=\overline{x}-4\frac{s^6}{m_3^2}\frac12\frac{m_3}{s^2}=\overline{x}-2\frac{s^4}{m_3}.$$

Для реальных данных [], приведенных выше в табл.6.2, выборочное среднее арифметическое $$\overline{x}=57,88$$, выборочная дисперсия $$s^2=663,00$$, выборочный третий центральный момент $$m_3=14927,91$$. Согласно только что полученным формулам оценки метода моментов таковы: $$a*=5,23; b*=11,26, c*=-1,01$$.

Оценки параметров гамма-распределения, полученные методом моментов, являются функциями от выборочных моментов. В соответствии со сказанным выше они являются асимптотически нормальными случайными величинами. Их распределения аппроксимируются нормальными распределениями, математические ожидания которых равны соответствующим параметрам, а дисперсии находятся с помощью формулы (4) с учетом формул (5) и (6). В табл.6.3 приведены оценки метода моментов и их асимптотические дисперсии при различных вариантах сочетания известных и неизвестных параметров гамма-распределения.

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

Все оценки метода моментов, приведенные в ]. Они охватывают все постановки задач оценивания параметров гамма-распределения (см. ] разработаны специальные методы оценивания.

Поскольку асимптотическое распределение оценок метода моментов известно, то не представляет труда формулировка правил проверки статистических гипотез относительно значений параметров распределений, а также построение доверительных границ для параметров. Например, в вероятностной модели, когда все три параметра неизвестны, в соответствии с третьей строкой таблицы 3 нижняя доверительная граница для параметра а, соответствующая доверительной вероятности $$\gamma=0,95$$, в асимптотике имеет вид$$a_H=a*-1,96 \left\{ \frac{6a*}{n}([a*]^2+6a*+5) \right\}^{\frac12},$$ а верхняя доверительная граница для той же доверительной вероятности:$$a_B=a*+1,96 \left\{ \frac{6a*}{n}([a*]^2+6a*+5) \right\}^{\frac12},$$

где $$а*$$ - оценка метода моментов параметра формы (табл.6.3).

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

В работах, предназначенных для первоначального знакомства с математической статистикой, обычно рассматривают оценки максимального правдоподобия (сокращенно ОМП):$$\theta_0(n)=\theta_0(n;x_1,x_2,...,x_n)=Arg\min\limits_{\theta\in\Theta}\prod_{i=1}^nf(x_i,\theta).$$

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

Хорошо известно, что оценки максимального правдоподобия входят в класс наилучших асимптотически нормальных оценок (определение дано ниже). Однако при конечных объемах выборки в ряде задач ОМП недопустимы, так как они хуже (дисперсия и средний квадрат ошибки больше), чем другие оценки, в частности, несмещенные []. Именно поэтому в ГОСТ 11.010-81 для оценивания параметров отрицательного биномиального распределения используются несмещенные оценки, а не ОМП []. Из сказанного следует, что априорно предпочитать ОМП другим видам оценок можно - если можно - лишь на этапе изучения асимптотического поведения оценок.

В отдельных случаях ОМП находятся явно, в виде конкретных формул, пригодных для вычисления.

Пример 3. Найдем ОМП для выборки из нормального распределения, каждый элемент которой имеет плотность$$f(x,m,\sigma^2)=\frac{1}{\sigma\sqrt{2\pi}}\exp\left\{-\frac{(x-m)^2}{2\sigma^2}\right\}.$$

Таким образом, надо оценить двумерный параметр $$(m, \sigma^2)$$.

Произведение плотностей вероятностей для элементов выборки, т.е. функция правдоподобия, имеет вид$$H(m;\sigma^2)=\sigma^{-n}(2\pi)^{-n/2} \exp\left\{-\frac{1}{2\sigma^2}\sum_{i=1}^n(x_i-m)^2\right\}.$$

Требуется решить задачу оптимизации$$H(m;\sigma^2)\rightarrow\max.$$

Как и во многих иных случаях, задача оптимизации проще решается, если прологарифмировать функцию правдоподобия, т.е. перейти к функции$$h(m;\sigma^2)-\ln H(m;\sigma^2),$$ называемой логарифмической функцией правдоподобия. Для выборки из нормального распределения$$h(m;\sigma^2)=(-n)\ln\sigma+\left(-\frac{n}{2}\right)\ln(2\pi)-\frac{1}{2\sigma^2}\sum_{i=1}^n(x_i-m)^2.$$

Необходимым условием максимума является равенство 0 частных производных от логарифмической функции правдоподобия по параметрам, т.е.$$\frac{\partial h(m,\sigma^2)}{\partial m}=0,\frac{\partial h(m,\sigma^2)}{\partial(\sigma^2)}=0.$$

Система (10) называется системой уравнений максимального правдоподобия. В общем случае число уравнений равно числу неизвестных параметров, а каждое из уравнений выписывается путем приравнивания 0 частной производной логарифмической функции правдоподобия по тому или иному параметру.

При дифференцировании по $$m$$ первые два слагаемых в правой части формулы (9) обращаются в 0, а последнее слагаемое дает уравнение$$\frac{\partial}{\partial m}\sum_{i=1}^n(x_i-m)= \sum_{i=1}^n 2(x_i-m)(-1)=0,\sum_{i=1}^n x_i=nm.$$

Следовательно, оценкой $$m*$$ максимального правдоподобия параметра m является выборочное среднее арифметическое,$$m*=\overline{x}.$$

Для нахождения оценки дисперсии необходимо решить уравнение$$\frac{\partial}{\partial(\sigma^2)}h(m;\sigma^2)= \frac{\partial}{\partial(\sigma^2)}(-n)\ln\sqrt{(\sigma^2)}- \frac{\partial}{\partial(\sigma^2)}\frac{1}{2\sigma^2} \sum_{i=1}^n(x_i-m)^2=0.$$

Легко видеть, что$$\frac{\partial}{\partial(\sigma^2)}(-n)\ln\sqrt{(\sigma^2)}=\frac{(-n)}{2\sigma^2}, -\frac{\partial}{\partial(\sigma^2)}\frac{1}{2\sigma^2}\sum_{i=1}^n(x_i-m)^2= \frac{1}{2\sigma^4}\sum_{i=1}^n(x_i-m)^2.$$

Следовательно, оценкой $$(\sigma^2)*$$ максимального правдоподобия для дисперсии $$\sigma^2$$ с учетом найденной ранее оценки для параметра $$m$$ является выборочная дисперсия,$$(\sigma^2)*=\frac{1}{n}\sum_{i=1}^n(x_i-\overlina{x})^2.$$

Итак, система уравнений максимального правдоподобия решена аналитически, ОМП для математического ожидания и дисперсии нормального распределения - это выборочное среднее арифметическое и выборочная дисперсия. Отметим, что последняя оценка является смещенной.

Отметим, что в условиях примера 3 оценки метода максимального правдоподобия совпадают с оценками метода моментов. Причем вид оценок метода моментов очевиден и не требует проведения каких-либо рассуждений.

В большинстве случаев аналитических решений не существует, для нахождения ОМП необходимо применять численные методы. Так обстоит дело, например, с выборками из гамма-распределения или распределения Вейбулла-Гнеденко. Во многих работах каким-либо итерационным методом решают систему уравнений максимального правдоподобия ([] и др.) или впрямую максимизируют функцию правдоподобия типа (8) (см. [] и др.).

Однако применение численных методов порождает многочисленные проблемы. Сходимость итерационных методов требует обоснования. В ряде примеров функция правдоподобия имеет много локальных максимумов, а потому естественные итерационные процедуры не сходятся []. Для данных ВНИИ железнодорожного транспорта по усталостным испытаниям стали уравнение максимального правдоподобия имеет 11 корней []. Какой из одиннадцати использовать в качестве оценки параметра?

Как следствие осознания указанных трудностей, стали появляться работы по доказательству сходимости алгоритмов нахождения оценок максимального правдоподобия для конкретных вероятностных моделей и конкретных алгоритмов. Примером является статья [].

Однако теоретическое доказательство сходимости итерационного алгоритма - это еще не всё. Возникает вопрос об обоснованном выборе момента прекращения вычислений в связи с достижением требуемой точности. В большинстве случаев он не решен.

Но и это не все. Точность вычислений необходимо увязывать с объемом выборки - чем он больше, тем точнее надо находить оценки параметров, в противном случае нельзя говорить о состоятельности метода оценивания. Более того, при увеличении объема выборки необходимо увеличивать и количество используемых в компьютере разрядов, переходить от одинарной точности расчетов к двойной и далее - опять-таки ради достижения состоятельности оценок.

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

Нет необходимости абсолютизировать ОМП. Кроме них, существуют другие виды оценок, обладающих хорошими статистическими свойствами. Примером являются одношаговые оценки (ОШ-оценки).

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

6.2. Одношаговые оценки

Одношаговые оценки имеют столь же хорошие асимптотические свойства, что и оценки максимального правдоподобия, при тех же условиях регулярности, что и ОМП. Грубо говоря, они представляют собой результат первой итерации при решении системы уравнений максимального правдоподобия по методу Ньютона-Рафсона. Одношаговые оценки выписываются в виде явных формул, а потому требуют существенно меньше машинного времени, а также могут применяться при ручном счете (на калькуляторах). Снимаются вопросы о сходимости алгоритмов, о выборе момента прекращения вычислений, о влиянии округлений при вычислениях на окончательный результат. ОШ-оценки были использованы нами при разработке ГОСТ 11.011-83 вместо ОМП.

Как и раньше, рассмотрим выборку $$x_1, x_2,..., x_n$$ из распределения с плотностью $$f(x;\theta_0)$$, где $$f(x;\theta_0)$$ - элемент параметрического семейства плотностей распределения вероятностей $${f(x;\theta), \theta\in\Theta}$$. Здесь $$\Theta$$ - известное статистику k-мерное пространство параметров, являющееся подмножеством евклидова пространства Rk, а конкретное значение параметра ?0 неизвестно. Его и будем оценивать.

Обозначим $$\theta=(\theta^1,\theta^2,...,\theta^k)$$. Рассмотрим вектор-столбец частных производных логарифма плотности вероятности$$s(x,\theta)= \left|\left| \frac{\partial}{\partial\theta^{\alpha}}\ln f(x,\theta),\alpha=1,2,...,k \right|\right|$$ и матрицу частных производных второго порядка для той же функции$$b(x,\theta)= \left|\left| \frac{\partial}{\partial\theta^{\alpha}\partial\theta^{\beta}}\ln f(x,\theta),\alpha,\beta=1,2,...,k \right|\right|.$$

Положим$$s_n(\theta)=\frac{1}{n}\sum_{i=1}^n s(x_i,\theta),b_n(\theta)=\frac{1}{n}\sum_{i=1}^n b(x_i,\theta).$$

Пусть матрица информации Фишера $$I(\theta_0)=M[-b_n(\theta_0)]$$ положительно определена.

, с.269]. Оценку $$\theta(n)$$ параметра $$\theta_0$$ называют наилучшей асимптотически нормальной оценкой (сокращенно НАН-оценкой), если распределение случайного вектора $$\sqrt{n}(\theta(n)-\theta_0)$$ сходится при $$n\rightarrow\infty$$ к нормальному распределению с нулевым математическим ожиданием и ковариационной матрицей, равной $$I^{-1}(\theta_0)$$.

Определение 1 корректно: $$I^{-1}(\theta_0)$$ является нижней асимптотической границей для ковариационной матрицы случайного вектора $$\sqrt{n}(\theta*(n)-\theta_0)$$, где $$\theta*(n)$$ - произвольная оценка; ОМП - это НАН-оценки (см. [] и др.). Некоторые другие оценки также являются НАН-оценками, например, байесовские. Сказанное об ОМП и байесовских оценках справедливо при некоторых условиях регулярности (см., например, []). В ряде случаев несмещенные оценки являются НАН-оценками, более того, они лучше, чем ОМП (их дисперсия меньше), при конечных объемах выборки [].

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

Пусть $$\theta_1(n)$$ и $$I_n^{-1}$$ - некоторые оценки $$\theta_0$$ и $$I^{-1}(\theta_0)$$ соответственно.

Определение 2. Одношаговой оценкой (ОШ-оценкой, или ОШО) называется оценка$$\theta_2(n)=\theta_1(n)+I_{n}^{-1}S_n(\theta_1(n)).$$

Теорема 1 []. Пусть выполнены следующие условия.

(I) Распределение $$\sqrt{n}s_n(\theta_0)$$ сходится при $$n\rightarrow\infty$$ к нормальному распределению с математическим ожиданием 0 и ковариационной матрицей $$I(\theta_0)$$ и, кроме того, существует $$Mb_n(\theta_0)b'_n(\theta_0)$$.

(II) При некотором $$\varepsilon 0$$ и $$n\rightarrow\infty$$$$\sup_{\theta:0<|\theta-\theta_0|<\varepsilon} \frac{|s_n(\theta)-s_n(\theta_0)-b_n(\theta_0)(\theta-\theta_0)|}{|\theta-\theta_0|^2}=O_p(1).$$

(III) Для любого $$\varepsilon>0$$$$\lim_{n\rightarrow\infty}P\{n^{1/4}(|\theta_1(n)-\theta_0|+||I_n^{-1}-I^{-1}(\theta_0)||)>\varepsilon\}=0.$$

Тогда ОШ-оценка является НАН-оценкой.

Доказательство. Рассмотрим тождество$$\sqrt{n}(\theta_2(n)-\theta_0)=\sqrt{n}(\theta_1(n)-\theta_0)+\sqrt{n}I_n^{-1}S_n(\theta_1(n)).$$

В силу условия (II) теоремы$$\sqrt{n}I_n^{-1}S_n(\theta_1(n))=\sqrt{n}I_n^{-1}S_n(\theta_0)+ \sqrt{n}I_n^{-1}b_n(\theta_0)(\theta_1(n)-\theta_0)+ \sqrt{n}I_n^{-1}O_p(|\theta_1(n)-\theta_0|^2).$$

Из условия (I) теоремы следует, что первое слагаемое в правой части формулы (2) сходится при $$n\rightarrow\infty$$ по распределению к нормальному закону с математическим ожиданием 0 и ковариационной матрицей $$I^{-1}(\theta_0)$$. Согласно условию (III)$$\sqrt{n}|\theta_1(n)-\theta_0|^2\rightarrow 0$$ по вероятности. Кроме того, согласно тому же условию последовательность матриц $$I_n^{-1}$$ ограничена по вероятности. Поэтому третье слагаемое в правой части формулы (2) сходится к 0 по вероятности. Для завершения доказательства теоремы осталось показать, что$$\sqrt{n}(\theta_1(n)-\theta_0)+sqrt{n}I_n^{-1}b_n(\theta_0)(\theta_1(n)-\theta_0)\rightarrow 0$$ по вероятности. Левая часть формулы (3) преобразуется к виду$$(E+I_n^{-1}b_n(\theta_0))\sqrt{n}(\theta_1(n)-\theta_0),$$ где $$Е$$ - единичная матрица. Поскольку из условия (I) теоремы следует, что для $$b_n(\theta_0)$$ справедлива (многомерная) центральная предельная теорема, то$$b_n(\theta_0)=-I(\theta_0)+O_p(n^{-1/2}).$$

С учетом условия (III) теоремы заключаем, что$$E+I_n^{-1}b_n(\theta_0)=O_p(n^{-1/4}).$$

Из соотношений (4), (5) и условия (III) теоремы вытекает справедливость формулы (3). Теорема доказана.

Прокомментируем условия теоремы. Условия (I) и (II) обычно предполагаются справедливыми при рассмотрении оценок максимального правдоподобия []. Эти условия можно выразить в виде требований, наложенных непосредственно на плотность $$f(x;\theta)$$ из параметрического семейства, как это сделано, например, в []. Условие (III) теоремы, наложенное на исходные оценки, весьма слабое. Обычно используемые оценки $$\theta_1(n)$$ и $$I_n^{-1}$$ являются не $$n^{-1/4}$$ -состоятельными, а $$\sqrt{n}$$ -состоятельными, т.е. условие (III) заведомо выполняется.

Какие оценки годятся в качестве начальных? В качестве $$\theta_1(n)$$ можно использовать оценки метода моментов, как это сделано в ГОСТ 11.011-83 [], или, например, квантильные. В качестве $$I_n^{-1}$$ в теоретической работе [] предлагается использовать простейшую оценку$$I_n^{-1}=-b_n^{-1}(\theta_1(n)).$$

Для гамма-распределения с неизвестными параметрами формы, масштаба и сдвига ОШ-оценки применены в []. При этом оценка (6) оказалась непрактичной, поскольку с точностью до погрешностей измерений и вычислений $$\det(b_n) = 0$$ для реальных данных о наработке резцов до предельного состояния, приведенных выше в ] в качестве ОШ-оценки была применена непосредственно первая итерация метода Ньютона-Рафсона решения системы уравнений максимального правдоподобия, т.е. была использована оценка$$I_n^{-1}=I^{-1}(\theta_1(n)).$$

В формуле (7) непосредственно используется явный вид зависимости матрицы информации Фишера от неизвестных параметров распределения.

В других случаях выбор тех или иных начальных оценок, в частности, выбор между (6) и (7), может определяться, например, простотой вычислений. Можно использовать также устойчивые аналоги [] перечисленных выше оценок.

Необходимо отметить, что еще в 1925 г., т.е. непосредственно при разработке метода максимального правдоподобия, его создатель Р.Фишер считал, что первая итерация по методу Ньютона-Рафсона дает хорошую оценку вектору неизвестных параметров [, с.298]. Он однако рассматривал эту оценку как аппроксимацию ОМП. А.А. Боровков воспринимает ОШ-оценки как способ "приближенного вычисления оценок максимального правдоподобия" [, с.225] и показывает асимптотическую эквивалентность ОШ-оценок и ОМП (в более сильных предположениях, чем в теореме 1; другими словами, теорема 1 обобщает результаты А.А. Боровкова относительно ОШ-оценок). Мы же полагаем, что ОШ-оценки имеют самостоятельную ценность, причем не меньшую, а в ряде случаев большую, чем ОМП. По нашему мнению, ОМП целесообразно применять (на этапе асимптотики) только тогда, когда они находятся явно. Во всех остальных случаях следует использовать на этом этапе ОШ-оценки (или какие-либо иные, выбранные из дополнительных соображений).

С чем связана популярность оценок максимального правдоподобия? Из всех НАН-оценок они наиболее просто вводятся и ранее других предложены. Поэтому среди математиков сложилась устойчивая традиция рассматривать ОМП в курсах математической статистики. Однако при этом игнорируются вычислительные вопросы, а также отодвигаются в сторону многочисленные иные НАН оценки.

В прикладной статистике - иные приоритеты. На первом месте - ОШ-оценки, все остальные НАН-оценки, в том числе ОМП, рассматриваются в качестве дополнительных возможностей.

Пример 1. Найдем ОШ-оценки для гамма-распределения с плотностью$$f(x;a,b,c)= \left\{ \begin{aligned} \frac{1}{\Gamma(a)}(x-c)^{a-1}b^{-a}\exp\left[-\frac{x-c}{b}\right],x\ge c,\\ 0,\;x<c \end{aligned} \right.$$

Плотность вероятности в формуле (8) определяется тремя параметрами $$a, b, c$$, где $$a>0, b>0$$. При этом $$a$$ является параметром формы, $$b$$ - параметром масштаба и $$с$$ - параметром сдвига. Здесь $$\Gamma(а)$$ - одна из используемых в математике специальных функций, так называемая "гамма-функция", по которой названо и распределение, задаваемое формулой (8),$$\Gamma(a)=\int_0^{+\infty}x^{a-1}e^{-x}dx.$$

Как следует из явного вида плотности (8), логарифмическая функция правдоподобия имеет вид [ $$10$$, с.98]:$$L=\sum_{i=1}^n\ln f(x_i;a,b,c)=-n\ln\Gamma(a)-na\ln b+ (a-1)\sum_{i=1}^n\ln(x_i-c)-\frac{1}{b}\sum_{i=1}^n x_i+\frac{nc}{b},$$ а уравнения правдоподобия таковы:$$\begin{gathered} \frac{\partial L}{\partial a}=-n\Psi(a)+\sum_{i=1}^n \ln\left(\frac{x_i-c}{b}\right)=0, \\ \frac{\partial L}{\partial b}=-\frac{na}{b}+\frac{1}{b^2} \sum_{i=1}^n(x_i-c)=0, \\ \frac{\partial L}{\partial c}=-(a-1) \sum_{i=1}^n\frac{1}{x_i-c}+\frac{n}{b}=0, \end{gathered}$$ где$$\Psi(a)=\frac{d}{da}\ln\Gamma(a).$$ Ясно, что выписанная система нелинейных уравнений не имеет аналитического решения, в отличие от аналогичной системы для семейства нормальных распределений. Построим ОШ-оценки для задачи оценивания трех неизвестных параметров [ $$20$$ ].

В качестве начальных оценок $$\theta_1(n)$$ будем использовать оценки метода моментов (см. 6.1):$$a*=4\frac{s^6}{m_3^2}, b*=\frac12\frac{m_3}{s^2},c*=\overline{x}-a*b*$$ где \overline{x} - выборочное среднее арифметическое, $$s^2$$ - выборочная дисперсия, $$m^3$$ - выборочный третий центральный момент.

Матрица информации Фишера согласно [ $$10$$, с.98] при $$a > 2$$ имеет вид$$I(\theta)=I(a,b,c)= \begin{Vmatrix} \frac{d\Psi(a)}{da} \frac{1}{b} \frac{1}{b(a-1)} \\ \frac{1}{b} \frac{a}{b^2} \frac{1}{b^2} \\ \frac{1}{b(a-1)} \frac{1}{b^2} \frac{1}{b^2(a-2)} \end{Vmatrix}.$$

Вектор-столбец частных производных логарифма плотности вероятности$$s_(x;\theta)=s(x;a,b,c)=(s(1),s(2),s(3))'$$ имеет координаты$$\begin{gathered} s(1)=-\Psi(a)+\ln\left(\frac{x-c}{b}\right), \\ s(2)=-\frac{a}{b}+\frac{x-c}{b^2}, \\ s(3)=-\frac{a-1}{x-c}+\frac{1}{b}. \end{gathered}$$

Таким образом, для получения $$s_n(a*, b*, c*)$$ необходимо вычислить две суммы$$\sum_{i=1}^n\ln\left(\frac{x_i-x}{b}\right),\;\sum_{i=1}^n\frac{1}{x_i-c}$$ и произвести еще несколько арифметических действий, число которых не зависит от объема выборки.

Одношаговые оценки $$a_n, b_n, c_n$$ для параметров гамма-распределения вычисляют по формуле$$(a_n, b_n, c_n)=(a*,b*,c*)+I^{-1}(a*,b*,c*)s_n(a*,b*,c*),$$ где $$I^{-1}$$ - обратная матрица к матрице информации Фишера $$I$$, заданной формулой (9). Матрицу $$I^{-1}$$ нетрудно рассчитать аналитически. Формулы для нахождения одношаговых оценок расписаны в [ $$6$$ ]. Расчеты облегчает то обстоятельство, что для гамма-распределения вторая координата вектора $$s_n(a*, b*, c*)$$ тождественно равна 0, т.е. $$s_n^{(2)}(a*, b*, c*)\equiv0$$.

При $$n\rightarrow\infty$$ распределение вектора оценок $$a_n, b_n, c_n$$ приближается трехмерным нормальным распределением с математическим ожиданием, равным вектору истинных значений параметров $$(a, b, c)$$, и ковариационной матрицей $$I^{-1}(a_n, b_n, c_n)$$. На этом приближении основаны правила расчета доверительных границ для параметров гамма-распределения [6]. Дисперсии оценок неизвестны, но зато имеются известные статистику зависимости этих дисперсий от параметров гамма-распределения. Эти зависимости непрерывные. Они стоят на главной диагонали ковариационной матрицы $$I^{-1}(a_n, b_n, c_n)$$ ). Поэтому можно вместо неизвестных параметров подставить в них оценки этих параметров и на основе принципа наследования сходимости (см. лекцию 4 выше) получить состоятельные оценки дисперсий. Затем на основе оценок дисперсий обычным образом строятся доверительные интервалы для параметров гамма-распределения.

В табл.6.4 приведены результаты реализации описанной выше схемы расчетов - точечные и интервальные (при односторонней доверительной вероятности 0,95) оценки параметров гамма-распределения для данных, содержащихся в табл.6.2 предыдущего п. 6.1.

Приведенные в табл.6.4 данные получены на основе асимптотических формул. Из-за конечности объема выборки необходимо внести некоторые коррективы. Поскольку параметр формы всегда положителен, $$a > 0$$, то нижняя доверительная граница для этого параметра должна быть неотрицательна, т.е. следует положить $$a_H = 0$$. Поскольку плотность гамма-распределения положительна только правее параметра $$c$$, то, очевидно, $$c\le x_{\min} = 9,00$$, верхняя доверительная граница для параметра сдвига должна быть заменена на $$c_B=9,00$$.

Может ли параметр сдвига быть отрицательным в данной прикладной задаче? Отрицательность параметра сдвига означает, что с положительной вероятностью рассматриваемая случайная величина отрицательна, т.е. наработка резца до предельного состояния отрицательна. Ясно, что такого быть не может, хотя для специалиста по математической статистике отрицательность параметра сдвига вполне приемлема. Однако специалист по прикладной статистике должен признать неотрицательность параметра с при обработке данных, составляющих рассматриваемую выборку. Следовательно, нижнюю доверительную границу для параметра сдвига необходимо заменить на $$c_н = 0$$.

Как следует из проведенных выше рассуждений и выкладок (см. также [ $$10$$, с.98-100]), отношение дисперсий оценок метода моментов и ОШ-оценок имеет вид$$\frac{Da_n}{Da*}=\frac{\left\{(a-1)^3+\frac15(a-1)\right\}}{a(a+1)(a+5)}$$ при больших $$a$$. Это отношение, как и должно быть из общих соображений, всегда меньше 1. Отношение дисперсий возрастает при приближении к 0 коэффициента асимметрии распределения. Если $$a > 39,1$$ (коэффициент асимметрии меньше 0,102), то эффективность оценки метода моментов превышает 80%. При $$a = 20$$ (коэффициент асимметрии 0,20) она равна 65%. Напомним, что при безграничном росте параметра формы а гамма-распределение приближается к нормальному, для которого оценки метода моментов и ОМП совпадают, а потому имеют равные дисперсии. Поэтому вполне естественно, что отношение дисперсий в формуле (10) стремится к 1 при безграничном росте параметра формы $$a$$.

Хотя дисперсии оценок метода моментов, как правило, больше, чем дисперсии НАН-оценок, таких, как ОШО и ОМП, метод моментов играет большую роль в прикладной статистике. Во-первых, обычно их расчет проще (в частности, требует меньшего числа компьютерных операций), чем оценок других типов. К тому же оценки находятся с помощью выборочных моментов, которые, как правило, вычисляются на этапе описания статистических данных. Во-вторых, они служат основой для вычисления оценок других типов, например, ОШО. Для запуска итерационных методов нахождения ОМП также нужны начальные значения, и ими обычно являются оценки метода моментов. В-третьих, при учете погрешностей результатов наблюдений оценки метода моментов могут оказаться точнее ОМП и асимптотически эквивалентных им ОШО (см. лекцию 12 настоящего курса).

Методы оценивания параметров гамма-распределения и примеры расчетов для всех семи постановок, перечисленных в ]. Большинство из них основано на асимптотических (при $$n\rightarrow\infty$$ ) теоретических результатах прикладной статистики. Методом статистических испытаний (Монте-Карло) показано, что уже при $$n\ge 10$$ используемые приближения удовлетворительны. Другими словами, асимптотической нормальностью оценок и другими важными для проведенных выше рассуждений предельными результатами можно пользоваться уже при $$n\ge 10$$.

Алгоритмическое и программное обеспечение ОШ-оценок для распределения Вейбулла-Гнеденко и гамма-распределения рассмотрено в монографии []. История вопроса освещена в статье [].

6.3. Асимптотика решений экстремальных статистических задач

Если проанализировать приведенные выше (см. 5.5) постановки и результаты, касающиеся эмпирических и теоретических средних и законов больших чисел, то становится очевидной возможность их обобщения. Так, доказательства теорем практически не меняются, если считать, что функция $$f(x,y)$$ определена на декартовом произведении бикомпактных пространств $$X$$ и $$Y$$, а не на $$X^2$$. Тогда можно считать, что элементы выборки лежат в $$Х$$, а $$Y$$ - пространство параметров, подлежащих оценке.

Обобщения законов больших чисел. Пусть, например, выборка $$х_1=х_1(\omega), х_2=х_2(\omega), ..., х_n=х_n(\omega)$$ взята из распределения с плотностью $$p(x,y)$$, где $$у$$ - неизвестный параметр. Если положить$$f(x,y)=-\ln p(x,y),$$ то задача нахождения эмпирического среднего$$f_n(\omega,y)=\frac{1}{n}\sum_{k=1}^n f(x_k(\omega),y)\rightarrow\min$$ переходит в задачу оценивания неизвестного параметра y методом максимального правдоподобия$$\sum_{k=1}^n\ln p(x_k(\omega),y)\rightarrow\max.$$

Соответственно законы больших чисел переходят в утверждения о состоятельности этих оценок в случае пространств $$X$$ и $$Y$$ общего вида. При такой интерпретации функция $$f(x,y)$$ уже не является расстоянием или показателем различия. Однако для доказательства сходимости оценок к соответствующим значениям параметров это и не требуется. Достаточно непрерывности этой функции на декартовом произведении бикомпактных пространств $$X$$ и $$Y$$.

В случае функции $$f(x,y)$$ общего вида можно говорить об определении в пространствах произвольной природы оценок минимального контраста и их состоятельности. При этом при каждом конкретном значении параметра $$y$$ справедливо предельное соотношение$$f_n(\omega,y)=\frac{1}{n}\sum_{k=1}^n f(x_k(\omega),y)\rightarrow Mf(x_1(\omega),y)=g(y),$$ где $$f$$ - функция контраста. Тогда состоятельность оценок минимального контраста вытекает из справедливости предельного перехода$$Arg\min\left\{\frac{1}{n}\sum_{k=1}^n f(x_k(\omega),y)\right\}\rightarrow Arg\min\{Mf(x_1(\omega),y)\}.$$

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

Можно пойти и дальше в обобщении законов больших чисел. Пусть известно, что при каждом конкретном y при безграничном росте n имеет быть сходимость по вероятности$$f_n(\omega,y)\rightarrow f(y),$$ где $$f_n(\omega, y)$$ - последовательность случайных функций на пространстве $$Y$$, а $$f(y)$$ - некоторая функция на $$Y$$. В каких случаях и в каком смысле имеет место сходимость$$Arg\min\{f_n(\omega,y),y\in X\}\rightarrow Arg\min\{f(y),y\in X\}?$$

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

Причем под n здесь можно понимать натуральное число. А можно рассматривать сходимость по направленному множеству (см. 4.3), или же, что практически то же самое - "сходимость по фильтру" в смысле Картана и Бурбаки [, с.118]. В частности, можно описывать ситуацию вектором, координаты которого - объемы нескольких выборок, и все они безгранично растут. В классической математической статистике такие постановки рассматривать не любят.

Поскольку, как уже отмечалось, основные задачи прикладной статистики можно представить в виде оптимизационных, то ответ на поставленный вопрос о сходимости точек минимума дает возможность единообразного подхода к изучению асимптотики решений разнообразных экстремальных статистических задач. Одна из возможных формулировок, основанная на бикомпактности пространств $$X$$ и $$Y$$ и нацеленная на изучение оценок минимального контраста, дана и обоснована выше. Другой подход развит в работе []. Он основан на использовании понятий асимптотической равномерной разбиваемости и координатной асимптотической равномерной разбиваемости пространств. С помощью указанных подходов удается стандартным образом обосновывать состоятельность оценок характеристик и параметров в основных задачах прикладной статистики.

Рассматриваемую тематику можно развивать дальше, в частности, рассматривать аналоги законов больших чисел в случае пространств, не являющихся бикомпактными, а также изучать скорость сходимости $$Arg\min\{f_n(x(\omega), y),y\in X\}$$ к $$Arg\min\{f(y), y\in X\}$$.

Приведем примеры применения результатов о предельном поведении точек минимума.

Задача аппроксимации зависимости (параметрической регрессии). Пусть $$X$$ и $$Y$$ - некоторые пространства. Пусть имеются статистические данные - $$n$$ пар $$(x_k, y_k)$$, где $$x_k\in X, y_k\in Y, k=1,2,...,n$$. Задано параметрическое пространство $$\Theta$$ произвольной природы и семейство функций $$g(x,\theta):X\times\Theta\rightarrow Y$$. Требуется подобрать параметр $$\theta\in\Theta$$ так, чтобы $$g(x_k,\theta)$$ наилучшим образом приближали $$y_k, k=1,2,...,n$$. Пусть $$f_k$$ - последовательность показателей различия в $$Y$$. При сделанных предположениях параметр $$\theta$$ естественно оценивать путем решения экстремальной задачи:$$\theta_n=Arg\min_{\theta\in\Theta} f_k(g(x_k,\theta),y_k).$$

Часто, но не всегда, все $$f_k$$ совпадают. В классической постановке, когда $$X=R^k,Y=R^1$$, функции $$f_k$$ различны при неравноточных наблюдениях, например, когда число опытов меняется от одной точки $$x$$ проведения опытов к другой.

Если $$f_k(y_1,y_2) = f(y_1,y_2) = (y_1-y_2)^2$$, то получаем общую постановку метода наименьших квадратов (см. подробности в лекции 9):$$\theta_n=Arg\min_{\theta\in\Theta}\sum_{k=1}^n(g(x_k,\theta)-y_k)^2.$$

В рамках детерминированного анализа данных остается единственный теоретический вопрос - о существовании $$\theta_n$$. Если все участвующие в формулировке задачи (1) функции непрерывны, а минимум берется по бикомпакту, то $$\theta_n$$ существует. Есть и иные условия существования $$\theta_n$$ [, , ].

При появлении нового наблюдения $$x$$ в соответствии с методологией восстановления зависимости рекомендуется выбирать оценку соответствующего $$y$$ по правилу$$y*=g(x,\theta_n).$$

Обосновать такую рекомендацию в рамках детерминированного анализа данных невозможно. Это можно сделать только в вероятностной теории, равно как и изучить асимптотическое поведение $$\theta_n$$, доказать состоятельность этой оценки.

Как и в классическом случае, вероятностную теорию целесообразно строить для трех различных постановок.

Во всех трех случаях$$f_n(\omega,\theta)=\sum_{k=1}^n f_k(g(x_k,\theta),y_k),$$ однако случайность входит в правую часть по-разному в зависимости от постановки, от которой, в свою очередь, зависит и определение предельной функции $$f(\theta)$$.

Проще всего выглядит $$f(\theta)$$ в случае второй постановки при $$f_k\equiv f:f(\theta)=Mf(g(x_1(\omega),\theta),y_1(\omega))$$.

В случае первой постановки$$f(\theta)=\lim_{n\rightarrow\infty}\sum_{k=1}^n Mf_k(g(x_k,\theta),y_k(\omega))$$ в предположении существования указанного предела. Ситуация усложняется для третьей постановки:$$f(\theta)=\lim_{n\rightarrow\infty}\sum_{k=1}^n Mf_k(g(x_k,(\omega),\theta),y_k(\omega)).$$

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

] для случайной величины $$(\xi, \eta)$$ со значениями в $$X\times Y$$ регрессией $$\eta$$ на $$\xi$$ относительно меры близости $$f$$ естественно назвать решение задачи$$Mf(g(\xi),\eta)\rightarrow\min_{g},$$ где $$f:Y\times Y \rightarrow R^1, g:X\rightarrow Y$$, минимум берется по множеству всех измеримых функций.

Можно исходить и из другого определения. Для каждого $$x\in X$$ рассмотрим случайную величину $$\eta(x)$$, распределение которой является условным распределением $$\eta$$ при условии $$\xi=x$$. В соответствии с определением математического ожидания в пространстве общей природы назовем условным математическим ожиданием решение экстремальной задачи$$M(\eta|\xi=x)=Arg\min\{Mf(y,\eta(x)),y\in Y\}.$$

Оказывается, при обычных предположениях измеримости решение задачи (2) совпадает с $$M(\eta|\xi=x)$$. (Внутриматематические уточнения типа "равенство имеет место почти всюду" здесь опущены.)

Если заранее известно, что условное математическое ожидание $$M(\eta|\xi=x)$$ принадлежит некоторому параметрическому семейству $$g(x,\theta)$$, то задача нахождения регрессии сводится к оцениванию параметра $$\theta$$ в соответствии с рассмотренной выше второй постановкой вероятностной теории параметрической регрессии. Если же нет оснований считать, что регрессия принадлежит параметрическому семейству, то можно использовать непараметрические оценки регрессии. Они строятся с помощью непараметрических оценок плотности (см. лекцию 5).

Пусть $$\nu_1$$ - мера в $$X$$, $$\nu_2$$ - мера в $$Y$$, а их прямое произведение $$\nu=\nu_1\times\nu_2$$ - мера в $$X\times Y$$. Пусть $$g(x,y)$$ - плотность случайного элемента $$(\xi,\eta)$$ по мере $$\nu$$. Тогда условная плотность $$g(y|x)$$ распределения $$\eta$$ при условии $$\xi=х$$ имеет вид$$g(y|x)=\frac{g(x,y)}{\int\limits_Y g(x,y)\nu_2(dy)}$$ (в предположении, что интеграл в знаменателе отличен от 0). Следовательно,$$Mf(y,\eta(x))=\int\limits_Y f(y,a)g(a|x)\nu_2(da),$$ а потому$$M(\eta|\xi=x)=Arg\min_{y\in Y} Mf(y,\eta(x))= Arg\min_{y\in Y}\int\limits_Y f(y,a)g(a|x)\nu_2(da).$$

Заменяя $$g(x,y)$$ в (3) непараметрической оценкой плотности $$gn(x,y)$$, получаем оценку условной плотности$$g_n(y|x)=\frac{g_n(x,y)}{\int\limits_Y g_n(x,y)\nu_2(dy)}.$$

Если $$g_n(x,y)$$ - состоятельная оценка $$g(x,y)$$, то числитель (4) сходится к числителю (3). Сходимость знаменателя (4) к знаменателю (3) обосновывается с помощью предельной теории статистик интегрального типа (см. лекцию 7). В итоге получаем утверждение о состоятельности непараметрической оценки (4) условной плотности (3).

Непараметрическая оценка регрессии ищется как M_n(\eta|\xi=x)=Arg\min_{y\in Y}\int\limits_Y f(y,a)g_n(a|x)\nu_2(da).

Состоятельность этой оценки следует из приведенных выше общих результатов об асимптотическом поведении решений экстремальных статистических задач.

Применение к методу главных компонент. Исходные данные - набор векторов $$\xi_1,\xi_2,...,\xi_n$$, лежащих в евклидовом пространстве $$R^k$$ размерности $$k$$. Цель состоит в снижении размерности, т.е. в уменьшении числа рассматриваемых показателей. Для этого берут всевозможные линейные ортогональные нормированные центрированные комбинации исходных показателей, получают $$k$$ новых показателей, из них берут первые $$m$$, где $$m < k$$ (подробности см. в лекции 9). Матрицу преобразования $$C$$ выбирают так, чтобы максимизировать информационный функционал$$I_n(C)\frac{s^2(z(1))+s^2(z(2))+...+s^2(z(m))}{s^2(x(1))+s^2(x(2))+...+s^2(x(k))},$$ где $$x(i), i=1,2,...,k$$, - исходные показатели; исходные данные имеют вид $$\xi_j=(x_j(1),x_j(2), ..., x_j(k)), j=1,2,...,n$$ ; при этом $$z(\alpha), \alpha= 1,2,...,m$$, - комбинации исходных показателей, полученные с помощью матрицы $$C$$. Наконец, $$s^2(z(\alpha)), \alpha=1,2,...,m, s^2(x(i)), i=1,2,...,k$$, - выборочные дисперсии переменных, указанных в скобках.

Укажем подробнее, как новые показатели (главные компоненты) $$z(\alpha)$$ строятся по исходным показателям $$x(i)$$ с помощью матрицы $$C$$:$$z_j(\alpha)=\sum_{\beta=1}^k c_{\alpha\beta}(x_j(\beta)-\overline{x(\beta)}),\;\alpha=1,2,...,m,\;j=1,2,...,n,$$ где$$\overline{x(\beta)}=\frac{1}{n}\sum_{j=1}^n x_j(\beta).$$

Матрица $$C=||c_{\alpha\beta}||$$ порядка $$m\times k$$ такова, что$$\sum_{\beta=1}^k c_{\alpha\beta}^2=1, \alpha=1,2,...,m$$ (нормированность),$$\sum_{\beta=1}^k c_{\alpha\beta}c_\{\gamma\beta}=0, \alpha,\gamma=1,2,...,m, \alpha\ne\gamma$$ (ортогональность).

Решением основной задачи метода главных компонент является$$C_n=Arg\min(-I_n(C)),$$ где минимизируемая функция определена формулой (5), а минимизация проводится по всем матрицам $$C$$, удовлетворяющим условиям (6) и (7).

Вычисление матрицы $$C_n$$ - задача детерминированного анализа данных. Однако, как и в иных случаях, например, для медианы Кемени, возникает вопрос об асимптотическом поведении $$C_n$$. Является ли решение основной задачи метода главных компонент устойчивым, т.е. существует ли предел $$C_n$$ при $$n\rightarrow\infty$$? Чему равен этот предел?

Ответ, как обычно, может быть дан только в вероятностной теории. Пусть $$\xi_1,\xi_2,...,\xi_n$$ - независимые одинаково распределенные случайные векторы. Положим$$z_{\infty}=\sum_{\beta=1}^k c_{\alpha\beta}(x_1(\beta)-Mx_1(\beta)),\;\alpha=1,2,...,m\;,$$ где матрица $$C=||c_{\alpha\beta}||$$ удовлетворяет условиям (6) и (7). Введем функцию от матрицы$$I(C)=\frac{D(z_{\infty}(1))+D(z_{\infty}(2))+...+D(z_{\infty}(m))}{D(x(1))+D(x(2))+...+D(x(k))}.$$

Легко видеть, что при $$n\rightarrow\infty$$ и любом C$$I_n(C)\rightarrow I(C).$$

Рассмотрим решение предельной экстремальной задачи$$C_{\infty}=Arg\min(-I(C)).$$

Естественно ожидать, что$$\lim_{n\rightarrow\infty} C_n=C_{infty}.$$

Действительно, это соотношение вытекает из приведенных выше общих результатов об асимптотическом поведении решений экстремальных статистических задач.

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

6.4. Робастность статистических процедур

Термин "робастность" ( robustnes ) образован от англ. robust - крепкий, грубый. Сравните с названием одного из сортов кофе - robusta. Имеется в виду, что робастные статистические процедуры должны "выдерживать" ошибки, которые теми или иными способами могут попадать в исходные данные или искажать предпосылки используемых вероятностно-статистических моделей.

Термин "робастный" стал популярным в нашей стране в 1970-е годы. Сначала он использовался фактически как сужение термина "устойчивый" на алгоритмы статистического анализа данных классического типа (не включая теорию измерений, статистику нечисловых и интервальных данных). Затем реальная сфера его применения сузилась.

Пусть исходные данные - это выборка, т.е. совокупность независимых одинаково распределенных случайных величин с одной и той же функцией распределения $$F(x)$$. Наиболее простая модель изучения устойчивости - это модель засорения$$F(x)=(1-\varepsilon)F_0(x)+\varepsilon H(x).$$

Эта модель именуется также моделью Тьюки-Хубера. (Джон Тьюки - американский исследователь, П. Хубер или Хьюбер - швейцарский ученый.) Модель (1) показывает, что с близкой к 1 вероятностью, а именно, с вероятностью $$(1-\varepsilon)$$ наблюдения берутся из совокупности с функцией распределения $$F_0(x)$$ которая предполагается обладающей "хорошими" свойствами. Например, она имеет известный статистику вид (хотя бы с точностью до параметров), у нее существуют все моменты, и т.д. Но с малой вероятностью $$\varepsilon$$ появляются наблюдения из совокупности с "плохим" распределением, например, взятые из распределения Коши, не имеющего математического ожидания, резко выделяющиеся аномальные наблюдения, выбросы.

Актуальность модели (1) не вызывает сомнений. Наличие засорений (выбросов) может сильно исказить результаты эконометрического анализа данных. Ясно, что если функция распределения элементов выборки имеет вид (1), где первое слагаемое соответствует случайной величине с конечным математическим ожиданием, а второе - такой, для которого математического ожидания не существует (например, если $$H(x)$$ - функция распределения Коши), то для итоговой функции распределения (1) также не существует математического ожидания. Исследователя обычно интересуют характеристики первого слагаемого, но найти их, т.е. освободиться от влияния засорения, не так-то просто. Например, среднее арифметическое результатов наблюдений не будет иметь никакого предела (это - строгое математическое утверждение, вытекающее из того, что математическое ожидание не существует []).

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

Оценивать характеристики и параметры, проверять статистические гипотезы, вообще осуществлять статистический анализ данных все чаще рекомендуют на основе эмпирических квантилей (другими словами, порядковых статистик, членов вариационного ряда), отделенных от концов вариационного ряда. Речь идет об использовании статистик вида$$ax(0,1n)+bx(0,3n)+cx(0,5n)+dx(0,7n)+ex(0,9n),$$ где $$a, b, c, d, e$$ - заданные числа, $$x(0,1n), x(0,3n), x(0,5n), x(0,7n), x(0,9n)$$ - члены вариационного ряда с номерами, наиболее близкими к числам, указанным в скобках. Так ценой небольшой потери в эффективности избавляемся от засоренности, подобно описанной в модели (1).

Вариантом этого подхода является переход к сгруппированным данным. Отрезок прямой, содержащий основную часть наблюдений, разбивается на интервалы, и вместо количественных значений статистик подсчитывает лишь, сколько наблюдений попало в те или иные интервалы. Особое значение приобретают крайние интервалы - к ним относят все наблюдения, которые больше некоторого верхнего порога и меньше некоторого нижнего порога. Любым методам анализа сгруппированных данных резко выделяющиеся наблюдения не страшны.

Можно поставить под сомнение и саму опасность засорения. Дело в том, что практически все реальные величины ограничены. Все они лежат на каком-то интервале - от и до. Это совершенно ясно, если речь идет о физическом измерении - все результаты измерений укладываются в шкалу прибора. По-видимому, и для иных статистических измерений наибольшие сложности создают не сверхбольшие помехи, а те засорения, что находятся "на грани" между "интуитивно возможным" и "интуитивно невозможным".

Что же это означает для практики статистического анализа данных? Если элементы выборки по абсолютной величине не превосходят числа $$A$$, то все засорение может сдвинуть среднее арифметическое на величину $$\varepsilon A$$. Если засорение невелико, то и сдвиг мал.

Построена достаточно обширная и развитая теория, посвященная разработке и изучению методов анализа данных в модели (1). С ней можно познакомиться по монографиям [ $$25$$ - $$27$$ ]. К сожалению, в теории обычно предполагается известной степень засорения , а на практике эта величина неизвестна. Кроме того, теория обычно направлена на защиту от воздействий, якобы угрожающих из бесконечности (например, отсутствием математического ожидания), а на самом деле реальные данные финитны (сосредоточены на конечных отрезках). Все это объясняет, почему теория робастности, исходящая из модели (1), популярна среди теоретиков, но мало интересна тем, кто анализирует реальные технические, экономические, медицинские и иные статистические данные.

Рассмотрим несколько более сложную модель. Пусть наблюдаются реализации $$x_1,x_2,...,x_n$$ независимых случайных величин с функциями распределения $$F_1(x),F_2(x),...F_n(x)$$ соответственно. Эта модель соответствует гипотезе о том, что в процессе наблюдения (измерения) условия несколько менялись. Естественной представляется модель малых отклонений функций распределений наблюдаемых случайных величин от некоторой "базовой" функции распределения $$F_0(x)$$. Множество возможных значений функций распределений наблюдаемых случайных величин (т.е. совокупность допустимых отклонений согласно общей схеме устойчивости, рассмотренной в лекции 4) описывается следующим образом:$$E((F_1,F_2,...,F_n);\varepsilon)=\{(F_1,F_2,...,F_n):\sup_x|F_i(x)-F_0(x)|<\varepsilon, i=1,2,...,n\}.$$

Следующий тип моделей - это введение малой (т.е. слабой) зависимости между рассматриваемыми случайными величинами (см., например, монографию []). Ограничения на взаимную зависимость можно задать разными способами. Пусть $$F(x_1,x_2,...,x_n)$$ - совместная функция распределения $$n$$ -мерного случайного вектора, $$F_1(x_1), F_2(x_2),... , F_n(x_n)$$ - функции распределения его координат. Если все координаты независимы, то $$F(x_1,x_2,...,x_n) = F_1(x_1)F_2(x_2)...F_n(x_n)$$. Пусть $$\rho(i,j)$$ коэффициент корреляции между $$i$$ -ой и $$j$$ -ой случайными величинами - координатами вектора. Множество возможных совместных функций распределения (т.е. совокупность допустимых отклонений согласно общей схеме устойчивости, рассмотренной в лекции 4) описывается следующим образом:$$\begin{gathered} E(F(x_1,x_2,...,x_n);\varepsilon)=\{F(x_1,x_2,...,x_n):P(x_i(\omega)<x) =F_i(x),|\rho(i,j)|\le\varepsilon, \\ 1\le i<j\le n\}. \end{gathered}$$

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

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

Разработано много вариантов робастных методов анализа статистических данных (см. монографии [ $$14$$, 25-28]). Иногда говорят, что робастные методы позволяют использовать информацию о том, что реальные наблюдения лежат "около" тех или иных параметрических семейств, например, нормальных. В этом, дескать, их преимущество по сравнению с непараметрическими методами, которые предназначены для анализа данных, распределенных согласно произвольной непрерывной функции распределения. Однако количественных подтверждений этих уверений любителей робастных методов обычно не удается найти. В основном потому, что термин "около" трудно формализовать.

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

В каждой из этих моделей общая схема устойчивости (лекция 4) предлагает для решения целый спектр задач устойчивости. Кроме изучения свойств робастности известных статистических процедур можно в каждой из постановок находить оптимальные процедуры. Однако практическая ценность этих оптимальных процедур, как правило, невелика, поскольку в других постановках оптимальными будут уже другие процедуры.

Контрольные вопросы и задачи

Темы докладов, рефератов, исследовательских работ