Статистика интервальных данных

Необходимость учета погрешностей измерений. Положим$$v=f(x_1,x_2,...,x_n)=\frac{1}{n}\sum_{1\le i\le n}\ln\left(\frac{\overline{x}}{x_i}\right).$$ Из свойств функции $$H(ullet)$$ следует [, с.14], что при малых $$v$$$$a^*\approx\frac{1}{2v}.$$

В силу состоятельности оценки максимального правдоподобия $$a^*$$ из формулы (13) следует, что $$v\rightarrow 0$$ по вероятности при $$a\rightarrow\infty$$.

Согласно модели статистики интервальных данных результатами наблюдений являются не $$x_i$$, а $$y_i$$, вместо $$v$$ по реальным данным рассчитывают$$w=f(y_1,y_2,...,y_n)=\frac{1}{n}\sum_{1\le i\le n}\ln\left(\frac{\overline{y}}{y_i}\right).$$

Имеем$$w-v=\ln\left(\frac{\overline{y}}{\overline{x}}\right) -\frac{1}{n}\sum_{1\le i\le n}\ln\left(1+\frac{\varepsilon_i}{x_i}\right).$$

В силу закона больших чисел при достаточно малой погрешности $$\varepsilon$$, обеспечивающей возможность приближения $$\ln(1+\alpha)\approx\alpha$$ для слагаемых в формуле (14), или, что эквивалентно, при достаточно малых предельной абсолютной погрешности $$\Delta$$ в формуле (1) или достаточно малой предельной относительной погрешности $$\delta$$ имеем при $$n\rightarrow\infty$$$$w-v\rightarrow\frac{M(\varepsilon_i)}{M(x_i)}-M\left(\frac{\varepsilon_i}{x_i}\right)=c$$ по вероятности (в предположении, что все погрешности одинаково распределены). Таким образом, наличие погрешностей вносит сдвиг, вообще говоря, не исчезающий при росте объема выборки. Следовательно, если $$c\ne 0$$, то оценка максимального правдоподобия не является состоятельной. Имеем$$a^*(y)-a^*\approx-\frac{c}{2v^2},$$ где величина $$a^*(y)$$ определена по формуле (12) с заменой $$x_i$$ на $$y_i, i=1,2,...,n$$. Из формулы (13) следует [], что$$a^*(y)-a\approx-2\left(a^*\right)^2 c,$$ т.е. влияние погрешностей измерений увеличивается по мере роста $$a$$.

Из формул для $$v$$ и $$w$$ следует, что с точностью до бесконечно малых более высокого порядка$$w-v\approx\sum_{1\le i\le n}\frac{\partial f}{\partial x_i}\varepsilon_i= \frac{1}{n}\sum_{1\le i\le n}\ln\left(\frac{1}{\overline{x}}-\frac{1}{x_i}\right)\varepsilon_i.$$

С целью нахождения асимптотического распределения $$w$$ выделим, используя формулу (16) и формулу для $$v$$, главные члены в соответствующих слагаемых$$w=\ln M(x_1)+\frac{1}{n}\sum_{1\le i\le n} \left\{ \frac{x_i=M(x_1)}{M(x_i)}-\ln x_i+\left(\frac{1}{M(x_1)}-\frac{1}{x_i}\right)\varepsilon_i \right\} +O_p\left(\frac{1}{n}\right).$$

Таким образом, величина w представлена в виде суммы независимых одинаково распределенных случайных величин (с точностью до зависящего от случая остаточного члена порядка $$1/n$$ ). В каждом слагаемом выделяются две части - одна, соответствующая $$v$$, и вторая, в которую входят $$\varepsilon_i$$. На основе представления (17) можно показать, что при $$n\rightarrow\infty, \varepsilon\rightarrow 0$$ распределения случайных величин $$v$$ и $$w$$ асимптотически нормальны, причем M(w)\approx M(v)+c, D(w)\approx D(v).

Из асимптотического совпадения дисперсий $$v$$ и $$w$$, вида параметров асимптотического распределения (при $$a\rightarrow\infty$$ ) оценки максимального правдоподобия $$a^*$$ и формулы (15) вытекает одно из основных соотношений статистики интервальных данных о квадрате средней ошибки$$M\left(a^*(y)-a\right)^2\approx 4a^4c^2+\frac{a(2a-1)}{n}.$$

Соотношение (18) уточняет утверждение о несостоятельности $$a^*$$. Из него следует также, что не имеет смысла безгранично увеличивать объем выборки $$n$$ с целью повышения точности оценивания параметра $$a$$, поскольку при этом уменьшается только второе слагаемое в (18), а первое остается постоянным.

В соответствии с общим подходом статистики интервальных данных в стандарте [] предлагается определять рациональный объем выборки nrat из условия "уравнивания погрешностей" (это условие было впервые предложено в монографии []) различных видов в формуле (18), т.е. из условия$$4a^4c^2=\frac{a(2a-1)}{n_{rat}}.$$

Упрощая это уравнение в предположении $$a\rightarrow\infty$$, получаем, что$$n_{rat}=frac{1}{2a^2c^2}.$$

Согласно сказанному выше, целесообразно использовать лишь выборки с объемами $$n\le n_{rat}$$. Превышение рационального объема выборки $$n_{rat}$$ не дает существенного повышения точности оценивания.

Применение методов теории устойчивости. Найдем асимптотическую нотну. Как следует из вида главного линейного члена в формуле (17), решение оптимизационной задачи$$w-v\rightarrow\max,|\varepsilon_i|\le\Delta,$$ соответствующей ограничениям на абсолютные погрешности, имеет вид$$\varepsilon_i=\left\{ \begin{aligned} \Delta,\frac{1}{\overline{x}}-\frac{1}{x_i}\ge 0, \\ -\Delta,\frac{1}{\overline{x}}-\frac{1}{x_i}< 0 \end{aligned} \right. .$$

Однако при этом пары $$(x_i,\varepsilon_i)$$ не образуют простую случайную выборку, так как в выражения для $$\varepsilon_i$$ входит $$\overline{x}$$. Однако при $$n\rightarrow\infty$$ можно заменить $$\overline{x}$$ на $$М(х_1)$$. Тогда получаем, что$$w-v\approx A\Delta$$ при $$a>1$$, где$$A=M\left|\frac{1}{M(x_1)}-frac{1}{x_1}\right|=\int\limits_0^{\infty} \left|\frac{1}{ab}-\frac{1}{x}\right|f(x;a,b)dx.$$

Таким образом, с точностью до бесконечно малых более высокого порядка нотна имеет вид$$N_{a^*}(y)=2(a^*)^2c,\quad c=A\Delta.$$

Применим полученные результаты к построению доверительных интервалов. В постановке классической математической статистики (т.е. при $$\varepsilon=0$$ ) асимптотический (при $$a\rightarrow\infty$$ ) доверительный интервал для параметра формы $$a$$, соответствующий доверительной вероятности , имеет вид []:$$\left[ a^*-u\left(\frac{1+\gamma}{2}\right)\sigma^*(a^*); a^*+u\left(\frac{1+\gamma}{2}\right)\sigma^*(a^*) \right],$$ где $$u\left(\frac{1+\gamma}{2}$$ - квантиль порядка $$\frac{1+\gamma}{2}$$ стандартного нормального распределения с математическим ожиданием 0 и дисперсией 1,$$[\sigma^*(x^*)]^2= \frac{a}{n(a^*\psi'(x^*)-1)}, \psi(a)=\left.\frac{d\Gamma(a)}{da}right/\Gamma(a).$$

В постановке статистики интервальных данных (т.е. при $$\varepsilon\ne 0$$ ) следует рассматривать доверительный интервал$$\left[ a^*-2(a^*)^2|c|-u\left(\frac{1+\gamma}{2}\right)\sigma^*(a^*); a^*+2(a^*)^2|c|-u\left(\frac{1+\gamma}{2}\right)\sigma^*(a^*) \right],$$ где$$c=\frac{M(\varepsilon_i)}{M(x_i)}-M\left(\frac{\varepsilon_i}{x_i}\right)$$ в вероятностной постановке (пары (x_i,\varepsilon_i) образуют простую случайную выборку) и $$c=A\Delta$$ в оптимизационной постановке. Как в вероятностной, так и в оптимизационной постановках длина доверительного интервала не стремится к 0 при $$n\rightarrow\infty$$.

Если ограничения наложены на предельную относительную погрешность, задана величина $$\delta$$, то значение с можно найти с помощью следующих правил приближенных вычислений [, с.142].

(I)Относительная погрешность суммы заключена между наибольшей и наименьшей из относительных погрешностей слагаемых.

(II) Относительная погрешность произведения и частного равна сумме относительных погрешностей сомножителей или, соответственно, делимого и делителя.

Можно показать, что в рамках статистики интервальных данных с ограничениями на относительную погрешность правила (I) и (II) являются строгими утверждениями при $$\delta\rightarrow 0$$.

Обозначим относительную погрешность некоторой величины $$t$$ через $$\text{ОП}(t)$$, абсолютную погрешность - через $$\text{АП}(t)$$.

Из правила (I) следует, что $$\text{ОП}(\overline{x})=\delta$$, а из правила (II) - что$$\text{ОП}\left(\frac{\overline{x}}{x_i}\right)=2\delta.$$

Поскольку рассмотрения ведутся при $$a\rightarrow\infty$$ то в силу неравенства Чебышева$$\frac{\overline{x}}{x_i}\rightarrow 1$$ по вероятности при $$a\rightarrow\infty$$ поскольку и числитель, и знаменатель в (19) с близкой к 1 вероятностью лежат в промежутке $$[ab-db\sqrt{a};ab+db\sqrt{a}]$$, где константа $$d$$ может быть определена с помощью упомянутого неравенства Чебышева.

Поскольку при справедливости (19) с точностью до бесконечно малых более высокого порядка$$\ln\left(\frac{\overline{x}}{x_i}\right)\approx\frac{\overline{x}}{x_i}-1,$$ то с помощью трех последних соотношений имеем$$\textit{ОП}\left(\frac{\overline{x}}{x_i}\right)= \textit{АП}\left(\frac{\overline{x}}{x_i}\right)= \textit{АП}\left(\ln\left(\frac{\overline{x}}{x_i}\right)\right) =2\delta.$$

Применим еще одно правило приближенных вычислений [, с.142].

(III) Предельная абсолютная погрешность суммы равна сумме предельных абсолютных погрешностей слагаемых.

Из (20) и правила (III) следует, что $$\textit{АП}(v)=2\delta$$.

Из (15) и (21) вытекает [, с.44, форм. (18)], что $$\textit{АП}(a^*)=4a^2\delta$$, откуда в соответствии с ранее полученной формулой для рационального объема выборки с заменой $$c=2\delta$$ получаем, что$$n_{rat}=\frac{1}{8a^2\delta^2}.$$

В частности, при $$a=5,00, \delta=0,01$$ получаем $$n_{rat}=50$$, т.е. в ситуации, в которой были получены данные о наработке резцов до предельного состояния [, с.29], проводить более 50 наблюдений нерационально.

В соответствии с ранее проведенными рассмотрениями асимптотический доверительный интервал для a, соответствующий доверительной вероятности $$\gamma=0,95$$, имеет вид$$\left[ a^*-4(a^*)^2\delta-1,96\sqrt{\frac{a^*(2a^*-1)}{n}}; a^*+4(a^*)^2\delta-1,96\sqrt{\frac{a^*(2a^*-1)}{n}} \right].$$

В частности, при $$a^*=5,00, \delta=0,01, n=50$$ имеем асимптотический доверительный интервал $$[2,12; 7,86]$$ вместо $$[3,14; 6,86]$$ при $$\delta=0$$.

При больших $$a$$ в силу соображений, приведенных при выводе формулы (19), можно связать между собой относительную и абсолютную погрешности результатов наблюдений $$x_i$$:$$\delta=\frac{\Delta}{M(x_1)}=\frac{\Delta}{ab}.$$

Следовательно, при больших $$a$$ имеем$$c=2\delta=A\Delta, A=\frac{2\delta}{\Delta}=\frac{2}{ab}.$$

Таким образом, проведенные рассуждения дали возможность вычислить асимптотику интеграла, задающего величину $$A$$.

Сравнение методов оценивания. Изучим влияние погрешностей измерений (с ограничениями на абсолютную погрешность) на оценку $$\widehat{a}$$ метода моментов. Имеем $$\textit{АП}(\overline{x})=\Delta, \; \textit{АП}((\overline{x})^2)\approx 2\overline{x}\Delta\approx 2ab\Delta$$.

Погрешность $$s^2$$ зависит от способа вычисления $$s^2$$. Если используется формула$$s^2=\frac{1}{n-1}\sum_{1\le i\le n}(x_i-\overline{x})^2,$$ то необходимо использовать соотношения$$\textit{АП}(x_i-\overline{x})^2=2\Delta, \textit{АП}\left[(x_i-\overline{x})^2\right]\approx 2|x_i-\overline{x}|\Delta.$$

По сравнению с анализом влияния погрешностей на оценку $$а^*$$ здесь возникает новый момент - необходимость учета погрешностей в случайной составляющей отклонения оценки $$\widehat{a}$$ от оцениваемого параметра, в то время как при рассмотрении оценки максимального правдоподобия погрешности давали лишь смещение. Примем в соответствии с неравенством Чебышева$$|x_i-\overline{x}|\approx\sqrt{D(x_1)},$$ тогда$$\textit{АП}[(x_i-\overline{x})^2]\approx 2b\sqrt{a}\Delta,\; \textit{АП}(s^2)\approx 2b\sqrt{a}\Delta.$$

Если вычислять $$s^2$$ по формуле$$s^2=\frac{1}{n-1}\sum_{1\le i\le n} x_i^2-\frac{n}{n-1}(\overline{x})^2,$$ то аналогичные вычисления дают, что$$\textit{АП}(s^2)\approx 4ab\Delta,$$ т.е. погрешность при больших а существенно больше. Хотя правые части формул (22) и (24) тождественно равны, погрешности вычислений по этим формулам весьма отличаются. Связано это с тем, что в формуле (24) последняя операция - нахождение разности двух больших чисел, примерно равных по величине (для выборки из гамма-распределения при большом значении параметра формы).

Из полученных результатов следует, что$$\textit{АП}(\widehat{a})=\textit{АП}\left(\frac{(\overline{x})^2}{s^2}\right)\approx\frac{2\Delta}{b}(1+sqrt{a}).$$ При выводе этой формулы использована линеаризация влияния погрешностей (выделение главного линейного члена). Используя связь (21) между абсолютной и относительной погрешностями, можно записать$$\textit{АП}(\widehat{a})\approx 2a(1+\sqrt{a})\delta.$$

Эта формула отличается от приведенной в [, с.44, форм. (19)]$$\textit{АП}(\widehat{a})\approx 2a(1+3\sqrt{a})\delta,$$ поскольку в [] вместо (23) использовалась оценка$$|x_i-\overline{x}|<3\sqrt{D(x_1)}.$$

Используя соотношение (23), мы характеризуем влияние погрешностей "в среднем".

Доверительный интервал, соответствующий доверительной вероятности 0,95, имеет вид \left[ \widehat{a}-2\widehat{a}(1+\sqrt{\widehat{a}})\delta-1,96\sqrt{\frac{2\widehat{a}(\widehat{a}+1)}{n}}; \widehat{a}+2\widehat{a}(1+\sqrt{\widehat{a}})\delta+1,96\sqrt{\frac{2\widehat{a}(\widehat{a}+1)}{n}} \right].

Если $$\widehat{a} = 5,00, \delta = 0,01, n = 50$$, то получаем доверительный интервал $$[2,54; 7,46]$$ вместо $$[2,86; 7,14]$$ при $$\delta=0$$. Хотя при $$\delta=0$$ доверительный интервал для $$a$$ при использовании оценки метода моментов $$\widehat{a}$$ шире, чем при использовании оценки максимального правдоподобия $$а^*$$, при $$\delta = 0,01$$ результат сравнения длин интервалов противоположен.

Необходимо выбрать способ сравнения двух методов оценивания параметра $$a$$, поскольку в длины доверительных интервалов входят две составляющие - зависящая от доверительной вероятности и не зависящая от нее. Выберем $$\delta = 0,68$$, т.е. $$u\left(\frac{1+\gamma}{2}\right)$$. Тогда оценке максимального правдоподобия $$a^*$$ соответствует полудлина доверительного интервала$$v(a^*)=4a^2\delta+\sqrt{\frac{a(2a-1)}{n}},$$ а оценке $$\widehat{a}$$ метода моментов соответствует полудлина доверительного интервала$$v(\widehat{a})=2a(1+\sqrt{a})\delta+\sqrt{\frac{2a(a+1)}{n}}.$$

Ясно, что больших $$a$$ или больших $$n$$ справедливо неравенство $$v(a^*)>v(\widehat{a})$$, т.е. метод моментов лучше метода максимального правдоподобия, вопреки классическим результатам Р.Фишера при $$\delta=0$$ [,с.99].

Из (25) и (26) элементарными преобразованиями получаем следующее правило принятия решений. Если$$\delta\sqrt{n}\ge\frac{\sqrt{2a(a+1)}-\sqrt{2a(a-1)}}{4a^2-2a(1+\sqrt{a})}=B(a),$$ то $$v(a^*)\ge v(\widehat{a})$$ и следует использовать $$\widehat{a}$$ ; а если $$\delta\sqrt{n}<B(a)$$, то $$v(a^*)<v(\widehat{a})$$ и надо применять $$а^*$$. Для выбора метода оценивания при обработке реальных данных целесообразно использовать $$B(\widehat{a})$$ (см. раздел 5 в ГОСТ 11.011-83 [, с.10-11]).

Пример анализа реальных данных опубликован в [].

На основе рассмотрения проблем оценивания параметров гамма-распределения можно сделать некоторые общие выводы. Если в классической теории математической статистики:

а) существуют состоятельные оценки $$a_n$$ параметра $$a$$,$$\lim_{n\rightarrow\infty}M(a_n-a)^2=0;$$

б) для повышения точности оценивания объем выборки целесообразно безгранично увеличивать;

в) оценки максимального правдоподобия лучше оценок метода моментов,

то в статистике интервальных данных, учитывающей погрешности измерений, соответственно:

а) не существует состоятельных оценок: для любой оценки $$a_n$$ существует константа $$c$$ такая, что$$\lim_{n\rightarrow\infty}M(a_n-a)^2\ge c>0;$$

б) не имеет смысла рассматривать объемы выборок, большие "рационального объема выборки" $$n_{rat}$$ ;

в) оценки метода моментов в обширной области параметров $$(a,n,\delta)$$ лучше оценок максимального правдоподобия, в частности, при $$a\rightarrow\infty$$ и при $$n\rightarrow\infty$$.

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

Метрологические, методические, статистические и вычислительные погрешности. Целесообразно выделить ряд видов погрешностей статистических данных. Погрешности, вызванные неточностью измерения исходных данных, называем метрологическими. Их максимальное значение можно оценить с помощью нотны. Впрочем, выше на примере оценивания параметров гамма-распределения показано, что переход от максимального отклонения к реально имеющемуся в вероятностно-статистической модели не меняет выводы (с точностью до умножения предельных значений погрешностей $$\Delta$$ или $$\delta$$ на константы). Как правило, метрологические погрешности не убывают с ростом объема выборки.

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

Статистическая погрешность - это та погрешность, которая традиционно рассматривается в математической статистике. Ее характеристики - дисперсия оценки, дополнение до 1 мощности критерия при фиксированной альтернативе и т.д. Как правило, статистическая погрешность стремится к 0 при росте объема выборки.

, с.51-52].

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

Вычислительную погрешность здесь подробно не рассматриваем. Ряд интересных результатов о ее роли в статистике получили Н.Н. Ляшенко и М.С. Никулин [].

Проведем сравнение методов оценивания параметров в более общей постановке.

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

Оценки минимального контраста. Пусть $$X$$ - пространство, в котором лежат независимые одинаково распределенные случайные элементы $$x_1, x_2, ..., x_n, ..$$.. Будем оценивать элемент пространства параметров $$\Theta$$ с помощью функции контраста $$f:X\times\Theta\rightarrow R^1$$. Оценкой минимального контраста называется$$\theta_n=Arg\min\left\{\sum_{1\le i\le n}f(x_i,\theta),\theta\in\Theta\right\}.$$

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

Оценками минимального контраста являются, в частности, многие робастные статистики [,]. Эти оценки широко используются в статистике объектов нечисловой природы [,], поскольку при $$X=\Theta$$ переходят в эмпирические средние, а если $$X=\Theta$$ - пространство бинарных отношений - в медиану Кемени.

Пусть в $$X$$ имеется мера $$\mu$$ (заданная на той же $$\sigma$$ -алгебре, что участвует в определении случайных элементов $$x_i$$ ), и $$p(x;\theta)$$ - плотность распределения $$x_i$$ по мере $$\mu$$. Если$$f(x;\theta)=-\ln p(x;\theta),$$ то оценка минимального контраста переходит в оценку максимального правдоподобия.

Асимптотическое поведение оценок минимального контраста в случае пространств $$X$$ и $$\Theta$$ общего вида хорошо изучено [], в частности, известны условия состоятельности оценок. Здесь ограничимся случаем $$X = R^1$$, но при этом введя погрешности измерений $$\varepsilon_i$$. Примем также, что $$\Theta=(\theta_{min},\theta_{max})\subseteq R^1$$.

В рассматриваемой математической модели предполагается, что статистику известны лишь искаженные значения $$y_i=x_i+\varepsilon_i,i=1,2,...,n$$. Поэтому вместо $$\theta_n$$ он вычисляет$$\theta_n^*=Arg\min\left\{\sum_{1\le i\le n}f(y_i,\theta),\theta\in\Theta\right\}.$$

Будем изучать величину $$\theta_n^*-\theta_n$$ в предположении, что погрешности измерений $$\varepsilon_i$$ малы. Цель этого изучения - продемонстрировать идеи статистики интервальных данных при достаточно простых предположениях. Поэтому естественно следовать условиям и ходу рассуждений, которые обычно принимаются при изучении оценок максимального правдоподобия [, п.33.3].

Пусть $$\theta_0$$ - истинное значение параметра, функция f(x;\theta) трижды дифференцируема по $$\theta$$, причем$$\left|\frac{\partial^3 f(x;\theta)}{\partial \theta^3}\right|<H(x)$$ при всех $$x, \theta$$. Тогда$$\frac{\partial f(x;\theta)}{\partial\theta}=\frac{\partial f(x;\theta_0)}{\partial\theta}= \frac{\partial^2 f(x;\theta_0)}{\partial\theta^2} (\theta-\theta_0)+\frac12\alpha(x)H(x)(\theta-\theta_0)^2,$$ где $$|\alpha(x)|<1$$.

Используя обозначения векторов $$x=(x_1,x_2,...,x_n),y=(y_1,y_2,...,y_n)$$, введем суммы$$B_0(x)=\frac{1}{n}\sum_{1\le i\le n}\frac{\partial f(x;\theta_0)}{\partial\theta}, B_1(x)=\frac{1}{n}\sum_{1\le i\le n}\frac{\partial^2 f(x_i;\theta_0)}{\partial\theta^2}, R(x)=\frac{1}{n}\sum_{1\le i\le n}H(x_i).$$

Аналогичным образом введем функции $$B_0(y), B_1(y), R(y)$$, в которых вместо $$x_i$$ стоят $$y_i, i=1,2,...,n$$.

Поскольку в соответствии с теоремой Ферма оценка минимального контраста $$\theta_n$$ удовлетворяет уравнению$$\sum_{1\le i\le n}\frac{\partial f(x_i;\theta_n)}{\partial\theta}=0,$$ то, подставляя в (27) $$x_i$$ вместо $$x$$ и суммируя по $$i = 1,2,...,n$$, получаем, что$$0=B_0(x)+B_1(x)(\theta_n-\theta_0)+\frac{\beta R(x)}{2}(\theta_n-\theta_0)^2,|\beta|<1,$$ откуда$$\theta_n-\theta_0=\frac{-B_0(x)}{B_1(x)+\frac{\beta R(x)}{2}}(\theta_n-\theta_0).$$

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

Теорема 1. Пусть для любого $$x$$ выполнено соотношение (27). Пусть для случайной величины $$x_1$$ с распределением, соответствующим значению параметра $$\theta=\theta_0$$, существуют математические ожидания$$M\frac{\partial f(x_1;\theta_0)}{\partial\theta_0}=0, M\frac{\partial^2 f(x_1;\theta_0)}{\partial\theta_0^2}=A\ne 0, MH(x_1)=M<+\infty.$$

Тогда существуют оценки минимального контраста $$\theta_n$$ такие, что $$\theta_n\rightarrow\theta_0$$ при $$n\rightarrow\infty$$ (в смысле сходимости по вероятности).

Доказательство. Возьмем $$\varepsilon>0$$ и $$\delta>0$$. В силу закона больших чисел (теорема Хинчина) существует $$n(\varepsilon,\delta)$$ такое, что для любого $$n>n(\varepsilon,\delta)$$ справедливы неравенства$$P\{|B_0|\ge\delta^2\}<\varepsilon/3, P\{|B_1|<|A|/2\}<\varepsilon/3, P\{R(x)>2M\}<\varepsilon/3.$$

Тогда с вероятностью не менее $$1-\varepsilon$$ одновременно выполняются соотношения$$|B_0|\le\delta^2,|B_1|\ge|A|/2,R(x)\le 2M.$$

При $$\theta\in[\theta_0-\delta;\theta_0+\delta]$$ рассмотрим многочлен второй степени$$y(\theta)=B_0(x)+B_1(x)(\theta-\theta_0)+\frac{\betaR(x)}{2}(\theta-\theta_0)^2$$ (см. формулу (29)). С вероятностью не менее $$1-\varepsilon$$ выполнены соотношения$$\left|B_0+\frac{\beta R(x)}{2}(\theta-\theta_0)^2\right|\le|B_0|+ \frac{R(x)\delta^2}{2}\le\delta^2(M+1),|B_1\delta|\ge\frac{|A|\delta}{2}.$$

Если $$0<2(M+1)\delta<|A|$$, то знак $$y(\theta)$$ в точках $$\theta_1=\theta_0-\delta$$ и $$\theta_2=\theta_0+\delta$$ определяется знаком линейного члена $$B_1(\theta_i-\theta_0),i=1,2$$, следовательно, знаки $$y(\theta_1)$$ и $$y(\theta_2)$$ различны, а потому существует $$\theta_n\in[\theta_0-\delta;\theta_0+\delta]$$ такое, что $$y(\theta_n)=0$$, что и требовалось доказать.

Теорема 2. Пусть выполнены условия теоремы 1 и, кроме того, для случайной величины $$x_1$$, распределение которой соответствует значению параметра $$\theta=\theta_0$$, существует математическое ожидание$$M\left(\frac{\partial f(x_1,\theta_0)}{\partial\theta}\right)^2=\sigma^2.$$

Тогда оценка минимального контраста имеет асимптотически нормальное распределение:$$\lim_{n\rightarrow\infty}P\left\{\sqrt{n}\frac{|A|}{\sigma}(\theta_n-\theta_0)<x\right}=\Phi(x)$$

для любого $$x$$, где $$\Phi(x)$$ - функция стандартного нормального распределения с математическим ожиданием 0 и дисперсией 1.

, $$\S$$ 2.4]).

Нотна оценки минимального контраста. Аналогично (30) нетрудно получить, что$$\theta_n^*-\theta_0=\frac{-B_0(y)}{B_1(y)+\frac{\beta(y)R(y)}{2}}(\theta_n^*-\theta_0), |\beta(y)|<.$$

Следовательно, $$\theta_n^*-\theta_n$$ есть разность правых частей формул (30) и (34). Найдем максимально возможное значение (т.е. нотну) величины при ограничениях (1) на абсолютные погрешности результатов измерений.

Покажем, что при $$\Delta\rightarrow 0$$ для некоторого $$C>0$$ нотна имеет вид$$N_{\theta_n}(x)=\sup_{\{\varepsilon\}}|\theta_n^*-\theta_n|=C\Delta(1+o(1)).$$

Поскольку $$\theta_n^*-\theta_n=(\theta_n^*-\theta_0)+(\theta_0-\theta_n)$$, то из (33) и (35) следует, что$$\sup_{\{\varepsilon\}}M(\theta_n^*-\theta_n)^2= \left(C^2\Delta^2+\frac{\sigma^2}{A^2 n}\right)(1+o(2)).$$

Можно сказать, что наличие погрешностей $$\varepsilon_i$$ приводит к появлению систематической ошибки (смещения) у оценки метода максимального правдоподобия, и нотна является максимально возможным значением этой систематической ошибки.

В правой части (36) первое слагаемое - квадрат асимптотической нотны, второе соответствует статистической ошибке. Приравнивая их, получаем рациональный объем выборки$$n_{rat}=\left(\frac{\sigma}{CA\Delta}\right)^2.$$

Остается доказать соотношение (35) и вычислить $$C$$. Укажем сначала условия, при которых $$\theta_n^*\rightarrow\theta_0$$ (по вероятности) при $$n\rightarrow\infty$$ одновременно с $$\Delta\rightarrow 0$$.

Теорема 3. Пусть существуют константа $$\Delta_0$$ и функции $$g_1(x), g_2(x), g_3(x)$$ такие, что при $$0\le\Delta\le\Delta_0$$ и $$-1\le\gamma\le 1$$ выполнены неравенства (ср. формулу (27))$$\begin{aligned} \left|\frac{\partial f(x;\theta_0)}{\partial \theta}- \frac{\partial f(x+\gamma\Delta;\theta_0)}{\partial \theta}\right|\le g_1(x)\Delta, \\ \left|\frac{\partial^2 f(x;\theta_0)}{\partial \theta^2}- \frac{\partial^2 f(x+\gamma\Delta;\theta_0)}{\partial \theta^2}\right|\le g_2(x)\Delta, \\ |H(x)-H(x+\gamma\Delta)|\le g_3(x)\Delta \end{aligned}$$ при всех $$x$$. Пусть для случайной величины $$x_1$$, распределение которой соответствует $$\theta=\theta_0$$, существуют $$m_1 = Mg_1(x_1), m_2 = Mg_2(x_1)$$ и $$m_3 = Mg_3(x_1)$$. Пусть выполнены условия теоремы 1. Тогда $$\theta_n^*\rightarrow\theta_0$$ (по вероятности) при $$\Delta\rightarrow 0,n\rightarrow\infty$$.

Доказательство проведем по схеме доказательства теоремы 1. Из неравенств (37) вытекает, что$$\begin{aligned} |B_0(y)-B_0(x)|\le\Delta\left(\frac{1}{n}\sum_{1\le i\le n}g_1(x_i)\right), \\ |B_1(y)-B_1(x)|\le\Delta\left(\frac{1}{n}\sum_{1\le i\le n}g_2(x_i)\right), \\ |R(y)-R(x)|\le\Delta\left(\frac{1}{n}\sum_{1\le i\le n}g_3(x_i)\right). \end{aligned}$$

Возьмем $$\varepsilon>0$$ и $$\delta>0$$. В силу закона больших чисел (теорема Хинчина) существует $$n(\varepsilon,\delta)$$ такое, что для любого $$n>n(\varepsilon,\delta)$$ справедливы неравенства$$\begin{aligned} P\left\{|B_0|\ge\frac{\delta^2}{2}\right\}<\frac{\varepsilon}{6}, \left\{|B_1|<\frac{3|A|}{4}\right\}<\frac{\varepsilon}{6}, P\left\{R(x)>\frac{3M}{2}\right\}<\frac{\varepsilon}{6}, \\ P\left\{\frac{1}{n}\sum_{1\le i\le n}g_j(x_i)>2m_j\right\}<\frac{\varepsilon}{6},j=1,2,3. \end{aligned}$$

Тогда с вероятностью не менее $$1-\varepsilon$$ одновременно выполняются соотношения$$|B_0|<\frac12\delta^2,|B_1\ge\frac{3|A|}{4}|,R(x)\le\frac{3M}{2}, \frac{1}{n}\sum_{1\le i\le n}g_j(x_i)\le 2m_j, j=1,2,3.$$

В силу (38) при этом$$|B_0(y)|<\frac12\delta^2+2\Delta m_1, |B_1(y)|\ge\frac{3|A|}{4}-2\Delta m_2, R(y)\le\frac{3M}{2}+2\Delta m_3.$$

Пусть$$0\le\Delta\le\min\left\{\frac14\frac{\delta^2}{m_1};\frac18\frac{|A|}{m_2};\frac14\frac{M}{m_3}\right\}.$$

Тогда с вероятностью не менее $$1-\varepsilon$$ одновременно выполняются соотношения (ср. (32))$$|B_0(y)|\le\delta^2,|B_1(y)|\ge|A|/2,R(y)\le 2M.$$

Завершается доказательство дословным повторением такового в теореме 1, с единственным отличием - заменой в обозначениях $$x$$ на $$y$$.

Теорема 4. Пусть выполнены условия теоремы 3 и, кроме того, существуют математические ожидания (при $$\theta=\theta_0$$ )$$M\left|\frac{\partial^2 f(x_1,\theta_0)}{\partial x\partial\theta}\right|,\; M\left|\frac{\partial^3 f(x_1,\theta_0)}{\partial x\partial\theta^2}\right|.$$

Тогда выполнено соотношение (35) с$$C=\frac{1}{|A|}M\left|\frac{\partial^2 f(x_1,\theta_0)}{\partial x\partial\theta}\right|.$$

Доказательство. Воспользуемся следующим элементарным соотношением. Пусть $$a$$ и $$b$$ - бесконечно малые по сравнению с $$Z$$ и $$B$$ соответственно. Тогда с точностью до бесконечно малых более высокого порядка$$\frac{Z+a}{B+b}-\frac{Z}{B}=\frac{aB-bZ}{B^2}.$$

Чтобы применить это соотношение к анализу $$\theta_n^*-\theta_n$$ в соответствии с (30), (34) и теоремой 2, положим$$Z=B_0(x),\;a=B_0(y)-B_0(x),\;B=B_1(x),\;b=(B_1(y)-B_1(x))+\frac{\beta(y)R(y)}{2}(\theta_n^*-\theta_0).$$

В силу условий теоремы 4 при малых $$\varepsilon_i$$ с точностью до членов более высокого порядка$$B_0(y)-B_0(x)=\frac{1}{n}\sum_{1\le i\le n}\frac{\partial^2 f(x_i,\theta_0)}{\partial x_i\partial\theta_0}\varepsilon_i, B_1(y)-B_1(x)=\frac{1}{n}\sum_{1\le i\le n}\frac{\partial^3 f(x_i,\theta_0)}{\partial x_i\partial\theta_0^2}\varepsilon_i.$$

При $$\Delta\rightarrow 0$$ эти величины бесконечно малы, а потому с учетом сходимости $$B_1(x)$$ к $$A$$ и теоремы 3$$\theta_n^*-\theta_n=\frac{1}{A_2}{(B_0(y)-B_0(x))A-(B_1(y)-B_1(x))B_0(x)}=\frac{1}{A^2}\sum_{1\le i\le n}\gamma_i\varepsilon_i$$ с точностью до бесконечно малых более высокого порядка, где$$\gamma_i=\frac{\partial^2 f(x_i,\theta_0)}{\partial x_i\partial\theta_0}A- \frac{\partial^3 f(x_i,\theta_0)}{\partial x_i\partial\theta_0^2}B_0(x).$$

Ясно, что задача оптимизации$$\left\{ \begin{aligned} \sum_{1\le i\le n}\gamma_i\varepsilon_i\rightarrow\max \\ |\varepsilon_i|\le\Delta, i=1,2,...,, \end{aligned} \right.$$ имеет решение$$\varepsilon_i= \left\{ \begin{aligned} \Delta,\gamma_i\ge 0, \\ -\Delta,\gamma_i<0, \end{aligned} \right.$$ при этом максимальное значение линейной формы есть $$\Delta\sum_{1\le i\le n}|\gamma_i|$$. Поэтому$$\sup_{\{\varepsilon\}}|\theta_n^*-\theta_n|=\frac{\Delta}{A^2n}\sum_{1\le i\le n}|\gamma_i|.$$

С целью упрощения правой части (42) воспользуемся тем, что$$\frac{1}{n}\sum_{1\le i\le n}|\gamma_i|=\frac{|A|}{n}\sum_{1\le i\le n}\left|\frac{\partial^2 f(x_i;\theta_0)}{\partial x\partial\theta_0}\right|+ \alpha\frac{|B_0(x)|}{n}\sum_{1\le i\le n}\left|\frac{\partial^3 f(x_i;\theta_0)}{\partial x\partial\theta_0^2}\right|,$$ где $$|\alpha|\le 1$$.

Поскольку при $$m\rightarrow\infty$$$$\frac{1}{n}\sum_{1\le i\le n}\left|\frac{\partial^3 f(x_i;\theta_0)}{\partial x\partial\theta_0^2}\right|\rightarrow M\left|\frac{\partial^3 f(x_1;\theta_0)}{\partial x\partial\theta_0^2}\right|<+\infty,\;B_0(x)\rightarrow 0$$ по вероятности, то второе слагаемое в (43) сходится к 0, а первое в силу закона больших чисел с учетом (39) сходится к $$CA^2$$, где $$C$$ определено в (40). Теорема 4 доказана.

Оценки метода моментов. Пусть $$g:R^k\rightarrowR^1,h_j:R^1\rightarrow R^1, j=1,2,...,k$$, - некоторые функции. Рассмотрим аналоги выборочных моментов$$m_j=\frac{1}{n}\sum_{1\le i\le n}h_j(x_i),\;j=1,2,...,k.$$

Оценки метода моментов имеют вид$$\widehat{\theta}_n=g(m_1,m_2,...,m_k)$$ (функции $$g$$ и $$h_j$$ должны удовлетворять некоторым дополнительным условиям [, с.80], которые здесь не приводим). Очевидно, что$$\theta_n(y)-\theta_n(x)=\sum_{1\le j\le k}\frac{\partial g}{\partial m_j}(m_j(y)-m_j(x)),$$ $$m_j(y)-m_j(x)=\frac{1}{n}\sum_{1\le j\le n}\frac{dh_j(x_i)}{dx_i}\varepsilon_i,j=1,2,...,k,$$ с точностью до бесконечно малых более высокого порядка, а потому с той же точностью$$\theta_n(y)-\theta_n(x)=\frac{1}{n}\sum_{1\le j\le n}\left(\sum_{1\le j\le k}\frac{\partial g}{\partial m_j}\frac{dh_j(x_i)}{dx_i}\right)\varepsilon_i.$$

Теорема 5. Пусть при $$\theta=\theta_0$$ существуют математические ожидания$$M_j=Mm_j=Mh_j(x_1), M\left(\frac{dh_j(x_1)}{dx_1}\right), j=1,2,...,n,$$ функция $$g$$ дважды непрерывно дифференцируема в некоторой окрестности точки $$(M_1,M_2,...,M_k)$$. Пусть существует функция $$t:R^1\rightarrow R^1$$ такая, что$$\sup_{|x-y|\le\Delta}\left|h_j(y)-h_j(x)-\frac{dh_j(x)}{dx}(y-x)\right|\let(x)\Delta^2,\;j=1,2,...,k,$$ причем $$Mt(x_1)$$ существует. Тогда$$\sup_{\{\varepsilon\}}|\widehat{\theta}_n(y)-\widehat{\theta}_n(x)|=C_1\Delta$$ с точностью до бесконечно малых более высокого порядка, причем$$C_1=M\left|\sum_{1\le j\le k}\frac{\partial g(M_1,M_2,...,M_k)}{\partial m_j}\frac{dh_j(x_1)}{dx_1}\right|.$$

Доказательство теоремы 5 сводится к обоснованию проведенных ранее рассуждений, позволивших получить формулу (45). В условиях теоремы 5 собраны предположения, достаточные для такого обоснования. Так, условие (46) дает возможность обосновать соотношения (44); существование $$M\left(\frac{dh_j(x_1)}{dx_1}\right)$$ обеспечивает существование $$C_1$$, и т.д. Завершает доказательство ссылка на решение задачи оптимизации (41) и применение закона больших чисел.

Полученные в теоремах 4 и 5 нотны оценок минимального контраста и метода моментов, асимптотические дисперсии этих оценок (см. теорему 2 и [] соответственно) позволяют находить рациональные объемы выборок, строить доверительные интервалы с учетом погрешностей измерений, а также сравнивать оценки по среднему квадрату ошибки (36). Подобное сравнение было проведено для оценок максимального правдоподобия и метода моментов параметров гамма-распределения. Установлено, что классический вывод о преимуществе оценок максимального правдоподобия [, с.99-100] неверен в случае $$\Delta>0$$.

12.3. Интервальные данные в задачах проверки гипотез

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

Рассмотрим подходы статистики интервальных данных в задачах проверки статистических гипотез. Пусть принятие решения основано на сравнении рассчитанного по выборке значения статистики критерия $$f=f(y_1,y_2,...,y_n)$$ с граничным значением $$C$$: если $$f>C$$, то гипотеза отвергается, если же $$f\le C$$, то принимается. С учетом погрешностей измерений выборочное значение статистики критерия может принимать любое значение в интервале $$[f(y)-N_f(y);\;f(y)+N_f(y)]$$. Это означает, что "истинное" значение порога, соответствующее реально используемому критерию, находится между $$C-N_f(y)$$ и $$C+N_f(y)$$, а потому уровень значимости описанного правила (критерия) лежит между $$1-P(C+N_f(y))$$ и $$1-P(C-N_f(y))$$, где $$P(Z)=P(f^lt;Z)$$.

Пример 1. Пусть $$x_1,x_2,...,x_n$$ - выборка из нормального распределения с математическим ожиданием $$a$$ и единичной дисперсией. Необходимо проверить гипотезу $$H_0:a=0$$ при альтернативе $$H_0:a\ne 0$$

Как известно из любого учебного курса математической статистики, следует использовать статистику $$f=\sqrt{n}|\overline{y}|$$ и порог C=\Phi(1-\alpha/2), где $$\alpha$$ - уровень значимости, $$\Phi(ullet)$$ - функция стандартного нормального распределения с математическим ожиданием 0 и дисперсией 1. В частности, $$C=1,96$$ при $$\alpha=0,05$$.

При ограничениях (1) на абсолютную погрешность $$N_f(y)=\sqrt{n}\Delta$$. Например, если $$\Delta = 0,1$$, а $$n = 100$$, то $$N_f(y) = 1,0$$. Это означает, что истинное значение порога лежит между 0,96 и 2,96, а истинный уровень значимости - между 0,003 и 0,34. Можно сделать и другой вывод: нулевую гипотезу $$H_0$$ допустимо отклонить на уровне значимости 0,05 лишь тогда, когда $$f > 2,96$$.

Если же $$n = 400$$ при $$\Delta = 0,1$$, то $$N_f(y)=2,0$$ и $$C-N_f(y)=-0,04$$, в то время как $$C+N_f(y)=3,96$$. Таким образом, даже в случае $$x=0$$ гипотеза $$H_0$$ может быть отвергнута только из-за погрешностей измерений результатов наблюдений.

Вернемся к общему случаю проверки гипотез. С учетом погрешностей измерений граничное значение $$C_{\alpha}$$ в статистике интервальных данных целесообразно заменить на $$C_{\alpha}+N_f(y)$$. Такая замена дает гарантию, что вероятность отклонения нулевой гипотезы $$H_0$$, когда она верна, не более $$\alpha$$. При проверке гипотез аналогом статистической погрешности, рассмотренной выше в задачах оценивания, является $$C_{\alpha}$$. Суммарная погрешность имеет вид $$C_{\alpha}+N_f(y)$$. Исходя из принципа уравнивания погрешностей [], целесообразно определять рациональный объем выборки из условия$$C_{\alpha}=N_f(y)$$

Если $$f=|f_1|$$, где $$f_1$$ при справедливости $$H_0$$ имеет асимптотически нормальное распределение с математическим ожиданием 0 и дисперсией $$\sigma_2/n$$ то$$C_{\alpha}=u\left(1-\frac{\alpha}{2}\right)\frac{\sigma}{\sqrt{n}}$$ при больших $$n$$, где $$u(1-\alpha/2)$$ - квантиль порядка $$1-\alpha/2$$ стандартного нормального распределения с математическим ожиданием 0 и дисперсией 1. Из (47) вытекает, что в рассматриваемом случае$$n_{rat}=\frac{3,84}{\Delta^2}=384$$

Пример 2. Рассмотрим статистику одновыборочного критерия Стьюдента$$t=\sqrt{n}\frac{\overline(y)}{s(y)}=\frac{\sqrt{n}}{v},$$ где $$v$$ - выборочный коэффициент вариации. Тогда с точностью до бесконечно малых более высокого порядка нотна для $$t$$ имеет вид$$N_t(y)=\frac{\sqrt{n}}{v^2}N_v(y),$$ где $$N_v(y)$$ - рассмотренная ранее нотна для выборочного коэффициента вариации. Поскольку распределение статистики Стьюдента $$t$$ сходится к стандартному нормальному, то небольшое изменение предыдущих рассуждений дает$$n_{rat}=\frac{v^4u^2(1-\alpha/2)}{N_v^2(y)}.$$

]. Статистика этого критерия имеет вид$$D_{mn}=\sup_x|F_m(x)-G_n(x)|,$$ где $$F_m(x)$$ - эмпирическая функция распределения, построенная по первой выборке объема $$m$$, извлеченной из генеральной совокупности с функцией распределения $$F(x)$$, а $$G_n(x)$$ - эмпирическая функция распределения, построенная по второй выборке объема $$n$$, извлеченной из генеральной совокупности с функцией распределения $$G(x)$$. Нулевая гипотеза имеет вид $$H_0:F(x)\equiv G(x)$$, альтернативная состоит в ее отрицании: $$H_1:F(x)\ne G(x)$$ при некотором $$x$$. Значение статистики сравнивают с порогом $$D(\alpha,m,n)$$ зависящим от уровня значимости $$\alpha$$ и объемов выборок $$m$$ и $$n$$. Если значение статистики не превосходит порога, то принимают нулевую гипотезу, если больше порога - альтернативную. Пороговые значения $$D(\alpha,m,n)$$ берут из таблиц []. Описанный критерий иногда неправильно называют критерием Колмогорова-Смирнова. История вопроса описана в [].

При ограничениях (1) на абсолютные погрешности и справедливости нулевой гипотезы $$H_0:F(x)\equiv G(x)$$ нотна имеет вид (при больших объемах выборок)$$N_D=\sup_x|F(x+\Delta)-F(x-\Delta)|.$$

Если $$F(x)=G(x)=x$$ при $$0\le x\le 1$$, то $$N_D=2\Delta$$. С помощью условия $$C_\alpha=N_f(y)$$ при уровне значимости $$\alpha=0,05$$ и достаточно больших объемах выборок (т.е. используя асимптотическое выражение для порога согласно []) получаем, что выборки имеет смысл увеличивать, если$$\frac{mn}{m+n}\le\frac{0,46}{\Delta^2}.$$

Правая часть этой формулы при $$\Delta=0,1$$ равна 46. Если $$m = n$$, то последнее неравенство переходит в $$n\le 92$$.

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

12.4. Линейный регрессионный анализ интервальных данных

Перейдем к многомерному статистическому анализу. Сначала с позиций асимптотической математической статистики интервальных данных рассмотрим оценки метода наименьших квадратов (МНК).

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

Если какая-либо группа объектов характеризуется переменными $$X_1, X_2, ..., X_m$$ и проведен эксперимент, состоящий из $$n$$ опытов, где в каждом опыте эти переменные измеряются один раз, то экспериментатор получает набор чисел: $$X_{1j}, X_{2j},...,X_{mj} (j = 1,...,n)$$.

Однако процесс измерения, какой бы физической природы он ни был, обычно не дает однозначный результат. Реально результатом измерения какой-либо величины $$X$$ являются два числа: $$X_{\textit{Н}}$$ - нижняя граница и $$X_{\textit{В}}$$ - верхняя граница. Причем $$X_{\textit{ИСТ}}\in [X_{\textit{Н}},X_{\textit{В}}]$$, где $$X_{\textit{ИСТ}$$ - истинное значение измеряемой величины. Результат измерения можно записать как $$X:[X_{\textit{Н}},X_{\textit{В}}]$$. Интервальное число $$X$$ может быть представлено другим способом, а именно, $$X:[Х_m,\Delta_x]$$, где $$X_{\textit{Н}} = X_m - \Delta_x, X_{\textit{В}} = X_m + \Delta_x$$. Здесь $$X_m$$ - центр интервала (как правило, не совпадающий с $$X_{\textit{ИСТ}}$$ ), а $$\Delta_x$$ - максимально возможная погрешность измерения.

Метод наименьших квадратов для интервальных данных. Пусть математическая модель задана следующим образом:$$y=Q(x,b)+\varepsilon,$$ где $$x=(x_1,x_2,...,x_m)$$ - вектор влияющих переменных (факторов), поддающихся измерению; $$b=(b_1,b_2,...,b_r)$$ - вектор оцениваемых параметров модели; $$y$$ - отклик модели (скаляр); $$Q(x,b)$$ - скалярная функция векторов $$x$$ и $$b$$ ; наконец, $$\varepsilon$$ - случайная ошибка (невязка, погрешность).

Пусть проведено $$n$$ опытов, причем в каждом опыте измерены (один раз) значения отклика $$(y)$$ и вектора факторов $$(x)$$. Результаты измерений могут быть представлены в следующем виде:$$X=\{х_{ij};i=1,...,n;j=1,...,m\},Y=(y_1,y_2,...,y_n),Е=(\varepsilon_1,\varepsilon_2,...,\varepsilon_n),$$ где $$X$$ - матрица значений измеренного вектора $$(x)$$ в $$n$$ опытах; $$Y$$ - вектор значений измеренного отклика в $$n$$ опытах; $$E$$ - вектор случайных ошибок. Тогда выполняется матричное соотношение:$$Y=Q(X,b)+E,$$ где $$Q(X,b)=(Q(x_1,b),Q(x_2,b),...,Q(x_n,b))^T$$, причем $$x_1,x_2,...,x_n$$ - $$m$$ -мерные вектора, которые составляют матрицу $$X=(x_1,x_2,...,x_n)T$$.

Введем меру близости $$d(Y,Q)$$ между векторами $$Y$$ и $$Q$$. В МНК в качестве $$d(Y,Q)$$ берется квадратичная форма взвешенных квадратов $$\varepsilon_i^2$$ невязок $$\varepsilon_i=y_i-Q(x_i,b)$$, т.е.$$d(Y,Q) = [Y-Q(X,b)]^T W[Y-Q(X,b)],$$ где $$W=\{w_ij, i,j=1,...,n\}$$ - матрица весов, не зависящая от $$b$$. Тогда в качестве оценки $$b$$ можно выбрать такое $$b^*$$, при котором мера близости $$d(Y,Q)$$ принимает минимальное значение, т.е.$$b^*=\{b:d(Y,Q)\rightarrow\min_{\{b\}}\}.$$

В общем случае решение этой экстремальной задачи может быть не единственным. Поэтому в дальнейшем будем иметь в виду одно из этих решений. Оно может быть выражено в виде $$b^* = f(X,Y)$$, где $$f(X,Y) = (f_1(X,Y), f_2(X,Y),..., f_m(X,Y))^T$$, причем $$f_i(X,Y)$$ непрерывны и дифференцируемы по $$(X,Y)\in Z$$, где $$Z$$ - область определения функции $$f(X,Y)$$. Эти свойства функции $$f(X,Y)$$ дают возможность использовать подходы статистики интервальных данных.

Преимущество метода наименьших квадратов заключается в сравнительной простоте и универсальности вычислительных процедур. Однако не всегда оценка МНК является состоятельной (при функции $$Q(X,b)$$, не являющейся линейной по векторному параметру b), что ограничивает его применение на практике.

Важным частным случаем является линейный МНК, когда $$Q(x,b)$$ есть линейная функция от $$b$$:$$y=b_0x_0+b_1x_1+...+b_mx_m+\varepsilon=bx^T+\varepsilon,$$ где, возможно, $$x_0 = 1$$, а $$b_0$$ - свободный член линейной комбинации. Как известно, в этом случае МНК-оценка имеет вид:$$b^*(X^TWX)^{-1}X^TWY.$$

Если матрица $$X^TWX$$ не вырождена, то эта оценка является единственной. Если матрица весов $$W$$ единичная, то$$b^*=(X^TX)^{-1}X^TY.$$

Пусть выполняются следующие предположения относительно распределения ошибок $$\varepsilon_i$$:

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

Как и в других постановках асимптотической математической статистики интервальных данных, при использовании МНК измеренные величины отличаются от истинных значений из-за наличия погрешностей измерения. Запишем истинные данные в следующей форме:$$X_R=\{x_{ij}^R;i=\overline{1,n};j=\overline{1,m}\},Y_R=(y_1^R,y_2^R,...,y_n^R),$$ где $$R$$ - индекс, указывающий на то, что значение истинное. Истинные и измеренные данные связаны следующим образом:$$X=X_R+\Delta X, Y=Y_R+\Delta Y,$$ где $$\Delta X=\{\Delta x_{ij};i=\overline{1,n};j=\overline{1,m}\},\Delta Y=(\Delta y_1,\Delta y_2,...,\Delta y_n)$$. Предположим, что погрешности измерения отвечают граничным условиям$$|\Delta x_{ij}|\le\Delta_j^x(i=1,2,...,n,\;j=1,2,...,m), |\Delta y_i|\le\Delta^y(i=1,2,...,n),$$ аналогичным ограничениям (1).

Пусть множество $$W$$ возможных значений $$(X_R,Y_R)$$ входит в $$Z$$ - область определения функции $$f(X,Y)$$. Рассмотрим $$b^{*R}$$ - оценку МНК, рассчитанную по истинным значениям факторов и отклика, и $$b^*$$ - оценку МНК, найденную по искаженным погрешностями данным. Тогда$$\Delta b^*=b^{*R}-b^*=f(X_R,Y_R)-f(X,Y).$$

Ввести понятие нотны придется несколько иначе, чем это было сделано выше, поскольку оценивается не одномерный параметр, а вектор. Положим:$$n(1)=(\sup\Delta b_1^*,\sup\Delta b_2^*,...,\sup\Delta b_r^*)^T, n(2)=-(\inf\Delta b_1^*,\inf\Delta b_2^*,...,\inf\Delta b_r^*)^T.$$

Будем называть $$n(1)$$ - нижней нотной, а $$n(2)$$ - верхней нотной. Предположим, что при безграничном возрастании числа измерений $$n$$, т.е. при $$n\rightarrow\infty$$, векторы $$n(1), n(2)$$ стремятся к постоянным значениям $$N(1), N(2)$$ соответственно. Тогда $$N(1)$$ будем называть нижней асимптотической нотной, а $$N(2)$$ - верхней асимптотической нотной.

Рассмотрим доверительное множество $$B_{\alpha}=B_{\alpha}(n,b^{*R})$$ для вектора параметров $$b$$, т.е. замкнутое связное множество точек в $$r$$ -мерном евклидовом пространстве такое, что $$P(b\in B_{\alpha})=\alpha$$ где $$\alpha$$ - доверительная вероятность, соответствующая $$B_{\alpha}(\alpha\approx 1)$$. Другими словами, $$B_{\alpha}(n,b^{*R})$$ есть область рассеивания (аналог эллипсоида рассеивания) случайного вектора $$b^{*R}$$ с доверительной вероятностью $$\alpha$$ и числом опытов $$n$$.

Из определения верхней и нижней нотн следует, что всегда $$b^{*R}\in [b^*-n(1);b^*+n(2)]$$.. В соответствии с определением нижней асимптотической нотны и верхней асимптотической нотны можно считать, что $$b^{*R}\in [b^*-N(1);b^*+N(2)]$$. при достаточно большом числе наблюдений $$n$$. Этот многомерный интервал описывает $$r$$ -мерный гиперпараллелепипед $$P$$.

Каким-либо образом разобьем $$P$$ на $$L$$ гиперпараллелепипедов. Пусть $$b_k$$ - внутренняя точка $$k$$ -го гиперпараллелепипеда. Учитывая свойства доверительного множества и устремляя $$L$$ к бесконечности, можно утверждать, что $$P(b\in C)\ge\alpha$$ где$$C=\lim_{L\rightarrow\infty}\bigcup_{1\le k\le L} B_{\alpha}(n,b_k).$$

Таким образом, множество $$C$$ характеризует неопределенность при оценивании вектора параметров $$b$$. Его можно назвать доверительным множеством в статистике интервальных данных.

Введем некоторую меру $$M(X)$$, характеризующую "величину" множества $$X\subseteq R^r$$. По определению меры она удовлетворяет условию: если $$X=Z\cup Y$$ и $$Z\cap Y=\varnothing$$, то $$M(X)=M(Z)+M(Y)$$. Примерами такой меры являются площадь для $$r = 2$$ и объем для $$r = 3$$. Тогда:$$M(C)=M(P)+M(F),$$ где $$F=C\P$$. Здесь $$M(F)$$ характеризует меру статистической неопределенности, в большинстве случаев она убывает при увеличении числа опытов $$n$$. В то же время $$M(P)$$ характеризует меру интервальной (метрологической) неопределенности, и, как правило, $$M(P)$$ стремится к некоторой постоянной величине при увеличении числа опытов $$n$$. Пусть теперь требуется найти то число опытов, при котором статистическая неопределенность составляет $$\delta$$ -ю часть общей неопределенности, т.е.$$M(F)=\delta M(C),$$ где $$\delta<1$$. Тогда, подставив соотношение (50) в равенство (49) и решив уравнение относительно $$n$$, получим искомое число опытов. В асимптотической математической статистике интервальных данных оно называется "рациональным объемом выборки". При этом $$\delta$$ есть "степень малости" статистической неопределенности $$М(P)$$ относительно всей неопределенности. Она выбирается из практических соображений. При использовании "принципа уравнивания погрешностей" согласно [] имеем $$\delta = 1/2$$.

Метод наименьших квадратов для линейной модели. Рассмотрим наиболее важный для практики частный случай МНК, когда модель описывается линейным уравнением (см. выше).

Для простоты описания преобразований пронормируем переменные $$х_{ij},у_i$$. следующим образом:$$x_{ij}^0=(x_{ij}-\overline{x}_j)/s(x_j),\; y_i^0=(y_i-\overline{y})/s(y),$$ где$$\overline{x}_j=\frac{1}{n}\sum_{1\le i\le n}x_{ij}, \; s^2(x_j)=\frac{1}{n}\sum_{1\le i\le n}(x_{ij}-\overline{x}_j), \; \overline{y}=\frac{1}{n}\sum_{1\le i\le n}y_i, \; s^2(y)=\frac{1}{n}\sum_{1\le i\le n}(y_i-\overline{y})^2.$$

Тогда$$\overline{x}_j^0=0s^2(x_j^0)=\frac{1}{n}\sum_{1\le i\le n}(x_{ij^0}-\overline{x}_j^0)=1,\; \overline{y}^0=0s^2(y^0)=\frac{1}{n}\sum_{1\le i\le n}(y_i^0-\overline{y}^0)^2=1,\;j=1,2,...,m.$$

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

1. Для рассматриваемых переменных существуют следующие пределы:$$\lim_{n\rightarrow\infty}\frac{1}{n}\sum_{1\le i\le n}\x_{ij}x_{ik}=0,\; j,k=1,2,...,m.$$

2. Количество опытов n таково, что можно пользоваться асимптотическими результатами, полученными при $$n\rightarrow\infty$$.

3. Погрешности измерения удовлетворяют одному из следующих типов ограничений:

тип 1. Абсолютные погрешности измерения ограничены согласно (48);

$$тип 2$$. Относительные погрешности измерения ограничены:$$|\Delta x_{ij}|\le\delta_j^x|x_{ij}|(i=1,2,....,n,\;j=1,2,...,m),\;|\Delta y_i|\le\delta^y|y_i|(i=1,2,...,n);$$

тип 3. Ограничения наложены на сумму погрешностей:$$\sum_{j=1}^m|\Delta_{ij}|\le\alpha_x(i=1,2,...,n,\;j=1,2,...,m),\;|\Delta y_i|\le\alpha_y(i=1,2,...,n)$$ (поскольку все переменные отнормированы, т.е. представляют собой относительные величины, то различие в размерности исходных переменных не влияет на возможность сложения погрешностей).

Перейдем к вычислению нотны оценки МНК. Справедливо равенство:$$\begin{aligned} \Delta b^*=b^{*R}-b^*=(X_R^T X_R)^{-1}X_R^T Y_R-(X^TX)^{-1}X^TY=\\ =(X_R^TX_R)^{-1}X_R^TY_R-((X_R+\Delta X)^T(X_R+\Delta X))^{-1}(X_R+\Delta X)(Y_R+\Delta Y). \end{aligned}$$

Воспользуемся следующей теоремой из теории матриц [].

Теорема. Если функция $$f(\lambda)$$ разлагается в степенной ряд в круге сходимости $$|\lambda - \lambda_0| < r$$, т.е.$$f(\lambda)=\sum_{k=0}^{\infty}\alpha_k(\lambda - \lambda_0)^k,$$ то это разложение сохраняет силу, если скалярный аргумент заменить любой матрицей $$A$$, характеристические числа которой $$\lambda_k, k = 1,...,n$$, лежат внутри круга сходимости.

Из этой теоремы вытекает, что:$$(E-A)^{-1}=\sum_{P=0}^{\infty}A^P,\quad\textit{если}\quad|\lambda_k|<1;\;k=1,...,n.$$

Легко убедиться, что:$$\begin{aligned} ((X_R+\Delta X)^T(X_R+\Delta X))^{-1}=-Z(E-\Delta\cdot Z)^{-1}, \\ \textit{где }Z=-(X_R^TX_R)^{-1},\Delta=X_R^T\Delta X+\Delta X^TX_R+\Delta X^T\Delta X. \end{aligned}$$

Это вытекает из последовательности равенств:$$\begin{aligned} ((X_R+\Delta X)^T(X_R+\Delta X))^{-1}=(X_R^TX_R+X_R^T\Delta X+\Delta X^T\Delta X)^{-1}=(X_R^TX_R+\Delta)^{-1}= \\ =(E+\Delta(X_R^TX_R)^{-1}X_R^TX_R)^{-1}=(X_R^TX_R)^{-1}(E+\Delta(X_R^TX_R)^{-1})^{-1}=-Z(E-\Delta\cdot Z)^{-1}. \end{aligned}$$

Применим приведенную выше теорему из теории матриц, полагая $$A = \Delta Z$$ и принимая, что собственные числа этой матрицы удовлетворяют неравенству $$|\lambda_k|<1$$. Тогда получим:$$((X_R+\Delta X)^T(X_R+\Delta X))^{-1}=-Z\sum_{P=0}^{\infty}(\Delta\cdot Z)^P=(X_R^TX_R)^{-1}\sum_{P=0}^{\infty}(-\Delta\cdot(X_R^TX_R)^{-1})^P.$$

Подставив последнее соотношение в заключение упомянутой теоремы, получим:$$\begin{aligned} \Delta b^*=(X_R^TX_R)^{-1}X_R^TY_R-((X_R^TX_R)^{-1} \sum_P^{\infty}(-\Delta\cdot(X_R^TX_R)^{-1})^P)(X_R+\Delta X)^T(Y_R+\Delta Y)= \\ =(X_R^TX_R)^{-1}X_R^TY_R-((X_R^TX_R)^{-1}\sum_P^{\infty}(-\Delta\cdot(X_R^TX_R)^{-1})^P) (X_R^TY_R+\Delta X^TY_R+ X_R^T\Delta Y+\Delta X^T\Delta Y). \end{aligned}$$

Для дальнейшего анализа понадобится вспомогательное утверждение. Исходя из предположений 1-3, докажем, что:$$(X_R^TX_R)^{-1}\approx\frac{1}{n}E.$$ Доказательство. Справедливо равенство$$X_R^TX_R=n \begin{pmatrix} D^*(x_1)\cdots\text{cov}^*(x_1,x_m) \\ \cdots\cdots\cdots \\ \text{cov}^*(x_1,x_m)\cdotsD^*(x_m) \end{pmatrix} =n\cov^*(x),$$ где $$D^*(x_i)\text{cov}^*(x_i,x_j)$$ - состоятельные и несмещенные оценки дисперсий и коэффициентов ковариации. Следовательно,$$D^*(x_i)=D(x_i)+o(1/n),\;\text{cov}^*(x_i,x_j)-\text{cov}^*(x_i,x_j)+o(1/n),$$ тогда$$XX_R^TX_R=n\text{cov}^*(x)=n(||\text{cov(x_i,x_j)}||+o(1/n)),$$ где$$o(1/n)=\{a_{ij}=o(1/n)\}(i=\overline{1,n},j=\overline{1,m}).$$

Другими словами, каждый элемент матрицы, обозначенной как $$o(1/n)$$, есть бесконечно малая величина порядка $$1/n$$. Для рассматриваемого случая $$\text{cov}(x)=E$$, поэтому$$X_R^TX_R=n\text{cov}^*(x)=n(E+o(1/n)).$$

Предположим, что $$n$$ достаточно велико и можно считать, что собственные числа матрицы $$o(1/n)$$ меньше единицы по модулю, тогда$$(X_R^TX_R)^{-1}=\frac{1}{n}\cdot(E+o(1/n))^{-1}\approx\frac{1}{n}(E+o(1/n))=\frac{1}{n}E+o(1/n^2)\approx\frac{1}{n}E,$$ что и требовалось доказать.

Подставим доказанное асимптотическое соотношение в формулу для приращения $$b^*$$, получим$$\begin{aligned} \Delta b^*=b^{*R}-\frac{1}{n}\sum_{P=0}^{\infty}\left(-\Delta\cdot\frac{1}{n}\right)^P (nb^*R+\Delta X^TY_R+X_R^T\Delta Y+\Delta X^T\Delta Y)= \\ =b^{*R}-\frac{1}{n}\sum_{P=0}^{\infty}(-(X_R^T\Delta X+\Delta X^TX_R+\Delta X^T\Delta X)\cdot\left(\frac{1}{n}\right)^P(nb^{*R}+\Delta X^TY_R+X_R^T\Delta Y+\Delta X^T\Delta Y)= \\ b^{*R}-\frac{1}{n}\left(E-(X_R^T\Delta X+\Delta X^TX_R+\Delta X^T\Delta X)\frac{1}{n}+(X_R^T\Delta X+\Delta X^TX_R+\Delta X^T\Delta X)^2\left(\frac{1}{n}\right)^2\right)\cdot \\ \cdot(nb^{*R}+\Delta X^TY_R+X_R^T\Delta Y +\Delta X^T\Delta Y). \end{aligned}$$

Перейдем от матричной к скалярной форме, опуская индекс (R):$$\begin{aligned} \Delta b_k^*=\frac{1}{n}\{\sum_j^m\sum_i^n(x_{ik}\Delta x_{ij}+\Delta x_{ik}x_{ij})b_j^*-\sum_i^n(\Delta x_{ik}y_i+x_{ik}\Delta y_i)\}; \\ \Delta b_k^*=\frac{1}{n}\{2\sum_i^nx_{ik}\Delta x_{ik}b_k^*+\sum_{j\ne k}^m \sum_i^n [(x_{ik}\Delta x_{ij}+\Delta x_{ik}x_{ij})b_j^*-\sum_{i}^n(\Delta x_{ik}y_i+x_{ik}\Delta y_i)\}= \\ =\frac{1}{n}\{2\sum_i^nx_{ik}\Delta x_{ik}b_k^*+\sum_{j\ne k}^m \sum_i^n [(x_{ik}\Delta x_{ij}+\Delta x_{ik}x_{ij})b_j^*-\frac{1}{m-1}\Delta x_{ik}y_i]-\sum_i^n x_{ik}\Delta y_i\}= \\ =\frac{1}{n}\{\sum_{j\ne k}^m \sum_i^n[\frac{2}{m-1}x_{ik}\Delta x_{ik}b_k^*+(x_{ik}\Delta x_{ij}+\Delta x_{ik}x_{ij})b_j^*-\frac{1}{m-1}\Delta x_{ik}y_i]-\sum_i^n x_{ik}\Delta y_i\}= \\ =\frac{1}{n}\{\sum_{j\ne k}^m \sum_i^n[(\frac{2}{m-1}x_{ik}b_k^*+x_{ij}b_j^*-\frac{1}{m-1}y_i)\Delta x_{ik}-x_{ik}b_j^*\Delta x_{ij}]-\sum_i^n x_{ik}\Delta y_i\} \end{aligned}$$

Будем искать $$\max(|\Delta b_k^*|)$$ по $$\Delta x_{ij}$$ и $$\Delta y_i (i=1,...,n; j=1,...,m)$$. Для этого рассмотрим все три ранее введенных типа ограничений на ошибки измерения.

Тип 1 (абсолютные погрешности измерения ограничены). Тогда:$$\max_{\Delta x,\Delta y}(|\Delta b_k^*|)=\frac{1}{n} \left\{ \sum_{j\ne k}^m \sum_i^n \left[\left|\left( \frac{2}{m-1}x_{ik}b_k^*+x_{ij}b_j^*-\frac{1}{m-1}y_i \right) \right| \Delta_k^x+|x_{ik}b_j^*|\Delta_j^x \right] -\sum_i^n |x_{ik}|\Delta y \right\}.$$

Тип 2 (относительные погрешности измерения ограничены). Аналогично получим:$$\sum_{j=1}^m|\Delta x_{ij}|<\alpha_x\;(i=1,2,...,n,\;j=1,2,...,m),\;|\Delta y_i|<\alpha_y\;(i=1,2,...,n).$$

Тип З (ограничения наложены на сумму погрешностей). Предположим, что $$|\Delta b_k^*|$$ достигает максимального значения при таких значениях погрешностей $$\Delta x_{ij}$$ и $$\Delta y_i$$, которые мы обозначим как:$$\{\Delta x_{ij}^*,\;i=\overline{1,2,...,n};j=1,2,...,m\},\;\{\Delta y_i^*,\;i=1,2,...,n\}.$$ тогда:$$\max_{\Delta x,\Delta y}(|\Delta b_k^*|)=\frac{1}{n} \left\{ \sum_{j\ne k}^m \sum_i^n \left[\left( \frac{2}{m-1}x_{ik}b_k^*+x_{ij}b_j^*-\frac{1}{m-1}y_i \right) x_{ik}^*+x_{ik}b_j^*x_{ij}^* \right] -\sum_i^n x_{ik}y_i^* \right\}.$$

Ввиду линейности последнего выражения и выполнения ограничения типа 3:$$\begin{aligned} \max_{\Delta x,\Delta y}(|\Delta b_k^*|)=\frac{1}{n} \left\{ \sum_{j\ne k}^m \sum_i^n \left[\left| \frac{2}{m-1}x_{ik}b_k^*+x_{ij}b_j^*-\frac{1}{m-1}y_i \right|\cdot |\Delta x_{ik}^*|+|x_{ik}b_j^*|\cdot|\Delta x_{ij}^*| \right] -\sum_i^n |x_{ik}|\cdot|\Delta y_i^*| \right\}, \\ \sum_j^m|\Delta x_{ij}^*|-\alpha_x\quad(j=1,2,...,m),\quad |\Delta y_i^*|=\alpha_y. \end{aligned}$$

Для простоты записей выкладок сделаем следующие замены:$$\begin{aligned} |\Delta x_{ij}|=\alpha_{ij}\ge 0,\;C_k=n\sum_i^n|x_{ik}|\cdot|\Delta y_i^*|\ge 0, \\ K_i^k=\sum_{j\ne k}^m\left|\frac{2}{m-1}x_{ik}b_k^*+x_{ij}b_j^*-\frac{1}{m-1}y_i\right|\ge 0, \\ |x_{ik}b_j^*|=R_{ij}^k\ge 0. \end{aligned}$$

Теперь для достижения поставленной цели можно сформулировать следующую задачу, которая разделяется на $$m$$ типовых задач оптимизации:$$f_k(\{\alpha_{ij}\})\rightarrow\max_{\alpha_{ij}}(i=1,2,...,n;j=1,2,...,m; k=1,2,...,m),$$ где$$f_k(\{\alpha_{ij}\})=\frac{1}{n}\left\{\sum_i^n K_i^k\alpha_{ik}+\sum_{j\ne m}^m\sum_i^n R_{ij}^k\alpha_{ij}\right\}+C_k,$$ при ограничениях$$\sum_j^m\alpha_{ij}=\alpha_x\quad (j=1,2,...,m).$$

Перепишем минимизируемые функции в следующем виде:$$f_k=\frac{1}{n}\sum_i^n\left(K_i^k\alpha_{ik}+\sum_{j\ne m}^m R_{ij}^k\alpha_{ij}\right)+C_k=\frac{1}{n}\sum_i^n f_i^k+C_k.$$

Очевидно, что $$f_{ik} > 0$$.

Легко видеть, что$$\begin{aligned} n\cdot\max_{\alpha_{ij}}(f_k)=\max_{\alpha_{i1}}(f_1^k)+\max_{\alpha_{i2}}(f_2^k)+... + \max_{\alpha_{in}}(f_n^k)+C_k=\sum_i^n\max_{\alpha_{ii}}(f_i^k)+C_k, \\ \textin{где } i=1,2,...,n;j=1,2,...,m. \end{aligned}$$

Следовательно, необходимо решить $$nm$$ задач$$\{f_i^k\}\rightarrow\max_{a_{ij}}(i=1,2,...,n;\;j=1,2,...,m;\;k=1,2,...,m)$$ при ограничениях "типа равенства":$$\begin{aligned} \sum_j^m\alpha_{ij}=\alpha_x\;(i=1,2,...,n),\\ \textit{где }f_i^k=K_i^k\alpha_{ik}+\sum_{j\ne m}^m R_{ij}^k\alpha_{ij}=\sum_j^m S_{ij}^k\alpha_{ij}, \\ \textit{причем }S_{ij}^k= \left\{ \begin{aligned} k_i^k,\textit{ если } j=k,\\ R_{ij}^k,\textit{ если } j\ne k. \end{aligned} \right. \end{aligned}$$

Сформулирована типовая задача поиска экстремума функции. Она легко решается. Поскольку$$\max_{a_{ij}}(f_l^k)=\max_j(S_{ij}^k)\cdot\alpha_x,$$ то максимальное отклонение МНК-оценки k-ого параметра равно$$\max_{\Delta X,\Delta Y}(|\Delta \widehat{b}_k|)=\max_{a_{ij}}(f_k)=\frac{1}{n} \alpha_x\sum_i^n\max_j(S_{ij}^k)+\frac{1}{n}C_k,\;(i=1,2,...,n;\;j=1,2,...,m).$$

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

Оценивание линейной корреляционной связи. В качестве примера рассмотрим оценивание линейной корреляционной связи случайных величин $$y$$ и $$x_1, x_2,..., x_m$$ с нулевыми математическими ожиданиями. Пусть эта связь описывается соотношением:$$y=\sum_{j=1}^m b_jx_j+e,$$ где $$b_1, b_2,...,b_m$$ - постоянные, а случайная величина $$e$$ некоррелирована с $$x_1, x_2,..., х_m$$. Допустим, необходимо оценить неизвестные параметры $$b_1, b_2,..., b_m$$ по серии независимых испытаний:$$y_i=\sum_{j=1}^m b_jx_{ij}+e_i,\;(i=1,2,...,n).$$

Здесь при каждом $$i = 1,2,...,n$$ имеем новую независимую реализацию рассматриваемых случайных величин. В этой частной схеме оценки наименьших квадратов $$b_1^{*R}, b_2^{*R}, ..., b_m^{*R}$$ параметров $$b_1, b_2, ..., b_m$$ являются, как известно, состоятельными [].

Пусть величины $$x_1, x_2, ..., x_m$$ в дополнение к попарной независимости имеют единичные дисперсии. Тогда из закона больших чисел [] следует существование следующих пределов (ср. предположение 1 выше):$$\begin{aligned} \lim_{n\rightarrow\infty}\left\{\frac{1}{n}\sum_i^nx_{ij}^R\right\}=M\{x_j\}=0\quad(j=\overline{1,m}), \\ \lim_{n\rightarrow\infty}\left\{\frac{1}{n}\sum_i^n(x_{ij}^R-M\{x_j\})^2\right\}=D\{x_j\}=1 \quad(j=\overline{1,m}), \\ \lim_{n\rightarrow\infty}\left\{\frac{1}{n}\sum_i^n(x_{ij}^R-M\{x_j\}) (x_{ik}^R-M\{x_k\})\right\}=0 \quad(j,k=\overline{1,m}), \\ \lim_{n\rightarrow\infty}\left\{\frac{1}{n}\sum_i^n y_i^R\right\}=M\{y\}=b_1 M\{x_1\}+...+b_m M\{x_m\} + M\{e\}=0, \\ \lim_{n\rightarrow\infty}\left\{\frac{1}{n}\sum_i^n(y_i^R-M\{y\})^2\right\}=D\{y\}=b_1^2+...+ b_m^2+\sigma^2, \end{aligned}$$ где $$\sigma$$ - среднее квадратическое отклонение случайной величины $$е$$.

Пусть измерения производятся с погрешностями, удовлетворяющими ограничениям типа 1, тогда максимальное приращение величины $$|\Delta b^*_k|$$, как показано выше, равно:$$\max_{\Delta x,\Delta y}(|\Delta b_k^*|)=\frac{1}{n} \left\{ \sum_{j\ne k}^m\sum_i^n \left[\left| \frac{2}{m-1}x_{ik}^R b_k^*+x_{ij}^R b_j^*-\frac{1}{m-1}y_i^r \right| \cdot\Delta_k^x+|x_{ik}^R b_j^*|\cdot\Delta_j^x \right] +\sum_i^n|x_{ik}^R|\cdot\Delta y \right\}.$$

Перейдем к предельному случаю и выпишем выражение для нотны:$$\begin{aligned} N_k=\lim_{n\rightarrow\infty}\{\max_{\Delta x, \Delta y}(|\Delta b_k^*|)\}= \\ =\sum_{j\ne k}^m[\{|\frac{2}{m-1}x_k b_k+x_j b_j-\frac{1}{m-1}y|\}\cdot\Delta_k^x+M\{|x_k b_j|\}\cdot\Delta_j^x+M\{|x_k|\}\cdot\Delta y. \end{aligned}$$

В качестве примера рассмотрим случай $$m = 2$$. Тогда$$\begin{aligned} N_1=M\{|2x_1b_1+x_2b_2-y|\}\Delta_1^x+M\{b_2x_1\}\Delta_2^x+M\{|x_1|\}\Delta y, N_2=M\{|2x_2b_2+x_1b_1-y|\}\Delta_2^x+M\{b_1x_2}\Delta_1^x+M\{|x_2|\}\Delta y. \end{aligned}$$

Приведенное выше выражение для максимального приращения метрологической погрешности не может быть использовано в случае $$m=1$$. Для $$m=1$$ выведем выражение для нотны, исходя из соотношения:$$\Delta b_k^*=\frac{1}{n} \left\{ \sum_j^m\sum_i^n(x_{ik}\Delta x_{ij}+\Delta x_{ik}x_{ij}),\;b_j^*-\sum_i^n(\Delta x_{ik}y_i+x_{ik}\Delta y_i) \right\}.$$

Подставив $$m = 1$$, получим:$$\Delta b^*=\frac{1}{n} \left\{ \sum_i^n(2x_i\Delta x_i)b^*-\sum_i^n(\Delta x_i y_i+x_i\Delta y_i) \right\} =\frac{1}{n} \left\{ \sum_i^n((2x_ib^*-y_i)\Delta x_i+x_i\Delta y_i) \right\}.$$

Следовательно, нотна выглядит так:$$N_f=M\{|2xb^* - y|\}\Delta x+M\{|x|\}\Delta y.$$

Для нахождения рационального объема выборки необходимо сделать следующее.

Этап 1. Выразить зависимость размеров и меры области рассеивания $$B_\alpha(n,b)$$ от числа опытов $$n$$ (см. выше).

Этап 2. Ввести меру неопределенности и записать соотношение между статистической и интервальной неопределенностями.

Этап 3. По результатам этапов 1 и 2 получить выражение для рационального объема выборки.

Для выполнения этапа 1 определим область рассеивания следующим образом. Пусть доверительным множеством $$B_{\alpha}(n,b)$$ является $$m$$ -мерный куб со сторонами длиною $$2K$$, для которого$$P(b\in B_{\alpha}(n,b^{*R}))=\alpha.$$

Исследуем случайный вектор $$b^*$$ и$$\begin{gathered} b^*R=(X_R^TX_R)^{-1}X_R^TY_R=(X_R^TX_R)^{-1}X_R^T(X_Rb+e)= \\ =(X_R^TX_R)^{-1}X_R^TX_Rb+(X_R^TX_R)^{-1}X_R^Te=b+(X_R^TX_R)^{-1}X_R^Te. \end{gathered}$$

Как известно, если элементы матрицы $$A = \{a_{ij}\}$$ - случайные, т.е. $$A$$ - случайная матрица, то ее математическим ожиданием является матрица, составленная из математических ожиданий ее элементов, т.е. $$M\{A\}=\{M\{a_{ij}\}\}$$.

Утверждение 1. Пусть $$A=\{a_{ij}\}$$ и $$B=\{b_{ij}\}$$ - случайные матрицы порядка $$(m \times n)$$ и $$(n \times r)$$ соответственно, причем любая пара их элементов $$(а_{ij},b_{kl})$$ состоит из независимых случайных величин. Тогда математическое ожидание произведения матриц равно произведению математических ожиданий сомножителей, т.е. $$M\{AB\} = M\{A\} M\{B\}$$.

Доказательство. На основании определения математического ожидания матрицы заключаем, что$$A\cdot B= \left\{\sum_k^n a_i{ik}\cdot b_{kj}\right\}\rightarrow M\{A\cdot B\}= \left\{M\left\{\sum_k^n a_i{ik}\cdot b_{kj}\right\}\right\}= \left\{\sum_k^n M\{a_{ik}\cdot b_{kj}\}\right\},$$ но так как случайные величины $$a_{ik}, b_{kj}$$ независимы, то$$M\{A\cdot B\}=\left\{\sum_k^n M\{a_{ik}\}\cdot M\{b_{kj}\}\right\}=M\{A\}\cdot M\{B\},$$ что и требовалось доказать.

Утверждение 2. Пусть $$A=\{a_{ij}\}$$ и $$B=\{b_{ij}\}$$ - случайные матрицы порядка $$(m\times n)$$ и $$(n\times r)$$ соответственно. Тогда математическое ожидание суммы матриц равно сумме математических ожиданий слагаемых, т.е. $$M\{А+В\} = M\{А\} + M\{В\}$$.

Доказательство. На основании определения математического ожидания матрицы заключаем, что$$M\{А+В\} = \{M\{а_{ij}+b_{ij}\}\} = \{M\{а_{ij}\} + M\{b_{ij}\}\} =M\{А\} + M\{В\},$$ что и требовалось доказать.

Найдем математическое ожидание и ковариационную матрицу вектора $$b^*$$ с помощью утверждений 1, 2 и выражения для $$b^{*R},$$ приведенного выше. Имеем$$M\{b^*R\}=b+M\{(X_R^TX_R)^{-1}X_R^Te\}=b+M\{(X_R^TX_R)^{-1}X_R^T\}\cdot M\{e\}.$$

Но так как $$M\{e\}=0$$, то $$M\{b^{*R}\}=b$$. Это означает, что оценка МНК является несмещенной.

Найдем ковариационную матрицу:$$D\left\{b^{*R}\right\}=M\left\{(b^{*R}-b)(b^{*R}-b)^r\right\}= M\left\{(X_R^TX_R)^{-1}X_R^T\cdot e\cdot e^T\cdot X_R(X_R^TX_R)^{-1}\right\}.$$

Можно доказать, что$$D\left\{b^{*R}\right\}=M\left\{(X_R^TX_R)^{-1}X_R^T\cdot M\left\{e\cdot e^T\right\}\cdot X_R(X_R^TX_R)^{-1}\right\},$$ но$$M\left\{e\cdot e^T\right\}=D\{e\}=\sigma^2 E,$$ поэтому$$D\left\{b^{*R}\right\}=M\left\{(X_T^RX_R)^{-1}X_R^T(\sigma^2 E)X_R(X_R^TX_R)^{-1}\right\}=\sigma^2 M\{(X_R^TX_R)^{-1}\}.$$

Как выяснено ранее, для достаточно большого количества опытов $$n$$ выполняется приближенное равенство$$(X_R^TX_R)^{-1}\approx\frac{1}{n}E,$$ $$D\{b^{*R}\}=\frac{\sigma^2}{n}E.$$ Осталось определить вид распределения вектора $$b^{*R}$$. Из выражения для $$b^{*R}$$, приведенного выше, и асимптотического соотношения (51) следует, что$$b^{*R}=b+\frac{1}{n}X_R^Te.$$

Можно утверждать, что вектор $$b^{*R}$$ имеет асимптотически нормальное распределение, т.е.$$b^{*R}\in N(b,\frac{\sigma^2}{n}E).$$

Тогда совместная функция плотности распределения вероятностей случайных величин $$b^{*R}_1, b^{*R}_2,..., b^{*R}_m$$ будет иметь вид:$$f(b^{*R})=\frac{1}{(2\pi)^{m/2}\cdot(\det C)^{1/2}}\cdot\exp[-\frac12(b^{*R}-b)^T\cdot C^{-1}\cdot(b^{*R}-b)],$$ где$$C=D(b^{*R})=\frac{\sigma^2}{n}E.$$

Тогда справедливы соотношения$$C^{-1}=\frac{n}{\sigma^2}E,\;\det C=\det\left(\frac{n}{\sigma^2}E\right)=\left(\frac{\sigma^2}{n}\right)^m.$$

Подставим в формулу (52), получим$$\begin{aligned} f(b^{*R})=\frac{1}{(2\pi)^{m/2}\cdot(\sigma^2/n)^{m/2}}\cdot\exp[-\frac{n}{2\sigma^2}(b^{*R}-b)^T\cdot C^{-1}\cdot(b^{*R}-b)]= \\ =\frac{1}{(\sigma\sqrt{2\pi/n})^m}\exp[-\frac{n}{2\sigma^2}(b^{*R}-b)^T\cdot C^{-1}\cdot(b^{*R}-b)]= \\ =\frac{1}{(\sigma\sqrt{2\pi/n})^m}\exp[-\frac{n}{2\sigma^2}(\beta_1^2+\beta_2^2+...+\beta_m^2)], \end{aligned}$$ где$$\beta_i=b_i^{*R}-b_i,\;i=1,2,...,m.$$

Вычислим асимптотическую вероятность попадания описывающего реальность вектора параметров $$b$$ в $$m$$ -мерный куб с длиной стороны, равной $$2k$$, и с центром $$b^{*R}$$.$$\begin{gathered} P(-k<\beta_1<k,-k<\beta_2<k,...,-k<\beta_m<k)= \\ =\frac{1}{(\sigma\sqrt{2\pi/n})^m}\{\int\limits_{-k}^k..\int\limits_{-k}^k\exp[-\frac{n}{2\sigma^2}(\beta_1^2+\beta_2^2+...+\beta_m^2)]\cdot d\beta_1...d\beta_m\}= \\ =\frac{1}{(\sigma\sqrt{2\pi/n})^m}\{\int\limits_{-k}^k\exp[-\frac{n}{2\sigma^2}\beta_1^2]d\beta_1...\int\limits_{-k}^k\exp[-\frac{n}{2\sigma^2}\beta_i^2]d\beta_i\}. \end{gathered}$$

Сделаем замену$$t_i=\sqrt{n/2}\cdot\frac{1}{\sigma}\beta_i,\;i=1,2,...,m.$$

Тогда$$\begin{aligned} P=P(-k<\beta_1<k,-k<\bela_2<k,...,-k<\beta_m<k)= \\ =\frac{(\sigma\sqrt{2/n})^m}{\sigma\sqrt{2\pi/n})^m}[\int\limits_{-T}^T e^{-t^2}dt]^m = [(1/\sqrt{\pi})\int\limits_{-T}^T e^{-t^2}dt]^m]=[\Phi_0(T)]^m, \end{aligned}$$ где $$T = (n/2)^{1/2}(k/\sigma)$$, а $$\Phi_0(T)$$ - интеграл Лапласа,$$\Phi_0(T)=\Phi(\sqrt{2}T)-\Phi(-\sqrt{2}T),$$ где $$\Phi(t)$$ - функция стандартного нормального распределения с математическим ожиданием 0 и дисперсией 1. Из последнего соотношения получаем$$T=\Phi^{-1}(P^{1/m}),$$ где $$\Phi^{-1}(P)$$ - обратная функция Лапласа. Отсюда следует, что$$k=\sigma(2/n)^{1/2}\Phi^{-1}(P^{1/m}).$$

Напомним, что доверительная область $$B_{\alpha}(n,b)$$ - это $$m$$ -мерный куб, длина стороны которого равна $$K$$, т.е.$$P(b\in B_{\alpha}(n,b))=P(-K<\beta_1<K,-K<\beta_2<K,...,-K<\beta_m<K)=\alpha.$$

Подставляя $$P=\alpha$$ в формулу (53), получим$$K=k=\sigma(2/n)^{1/2}\Phi^{-1}(\alpha^{1/m}).$$

Соотношение (54) выражает зависимость размеров доверительной области (т.е. длины ребра куба $$K$$ ) от числа опытов $$n$$, среднего квадратического отклонения $$\sigma$$ ошибки $$e$$ и доверительной вероятности $$\alpha$$. Это соотношение понадобится для определения рационального объема выборки.

Переходим к этапу 2. Необходимо ввести меру разброса (неопределенности) и установить соотношение между статистической и интервальной (метрологической) неопределенностями в соответствии с ранее сформулированным общим подходом.

Пусть $$A$$ - некоторое измеримое множество точек в $$m$$ -мерном евклидовом пространстве, характеризующее неопределенность задания вектора $$a\in A$$. Тогда необходимо ввести некую меру $$M(A)$$, измеряющую степень неопределенности. Такой мерой может служить $$m$$ -мерный объем $$V(A)$$ множества $$A$$ (т.е. его мера Лебега или Жордана), $$M(A)=V(A)$$.

Пусть $$P$$ - $$m$$ -мерный параллелепипед, характеризующий интервальную неопределенность. Длины его сторон равны значениям нотн $$2N_1, 2N_2,..., 2N_m$$, а центр $$a$$ (точка пересечений диагоналей параллелепипеда) находится в точке $$b^{*R}$$. Пусть $$C$$ - измеримое множество точек, характеризующее общую неопределенность. В рассматриваемом случае это $$m$$ -мерный параллелепипед, длины сторон которого равны $$2(N_1+K),2(N_2+K),...,2(N_m+K)$$, а центр находится в точке $$b^{*R}$$. Тогда$$M(P)=V(P)=2^m N_1 N_2... N_m,$$ $$M(C)=V(C)=2^m(N_1+K)(N_2+K)...(N_m+K).$$

Справедливо соотношение (49), согласно которому $$M(C)=M(P)+M(F)$$, где множество $$F=C\P$$ характеризует статистическую неопределенность.

На этапе 3 получаем по результатам этапов 1 и 2 выражение для рационального объема выборки. Найдем то число опытов, при котором статистическая неопределенность составит $$\delta$$ 100% от общей неопределенности, т.е. согласно правилу (50)$$M(F)=M(C)-M(P)=\delta M(C),$$ где $$0<\delta<1$$. Подставив (55) и (56) в (57), получим$$2^m\Pi_{i=1}^m(N_i+K)-2^m\Pi_{i=1}^m(N_i)=2^m\delta\Pi_{i=1}^m(N_i+K).$$

Следовательно, $$(1-\delta)\Pi_{i=1}^m(N_i+K)/\Pi_{i=1}^m(N_i)=1$$.

Преобразуем эту формулу:$$(1-\delta)\Pi_{i=1}^m(1+K/N_i)=1,$$ откуда$$\Pi_i^m(1+K/N_i)=1/(1-\delta).$$

Если статистическая погрешность мала относительно метрологической, т.е. величины $$K/N_i$$ малы, то$$\Pi_i^m(1+K/N_i)\approx 1+\sum_i^m(K/N_i).$$

При $$m = 1$$ эта формула является точной. Из нее следует, что для дальнейших расчетов можно использовать соотношение$$1+\sum_i^m(K/N_i)=1/(1-\delta).$$

Отсюда нетрудно найти $$K$$:$$K=\frac{\delta}{1-\delta}(1/\sun_{i=1}^m(1/N_i)).$$

Подставив в формулу (58) зависимость $$K=K(n)$$, полученную в формуле (54), находим приближенное (асимптотическое) выражение для рационального объема выборки:$$n_{\textit{рац}}=2\left(\frac{1-\delta}{\delta}\sigma\sum_{i=1}^m(1/N_i)\cdot\Phi^{-1}(\alpha^{1/m})\right)^2.$$

При $$m = 1$$ эта формула также справедлива, более того, является точной.

Переход от произведения к сумме является обоснованным при достаточно малом $$\delta$$, т.е. при достаточно малой статистической неопределенности по сравнению с метрологической. В общем случае можно находить $$K$$ и затем рациональный объем выборки тем или иным численным методом.

Пример 1. Представляет интерес определение nрац для случая, когда $$m = 2$$, поскольку простейшая линейная регрессия с $$m =2$$ широко применяется. В этом случае базовое соотношение имеет вид$$(1 + K/N_1)(1 + K/N_2) = 1/(1 - \delta).$$

Решая это уравнение относительно $$K$$, получаем$$K=0.5\{-(N_1 + N_2) + [(N_1 + N_2)^2 + 4 N_{1N2} (\delta/(1 - \delta)]^{1/2}\}.$$

Далее, подставив в формулу (54), получим уравнение для рационального объема выборки в случае $$m = 2$$:$$\sigma(2/n)^{1/2}\Phi^{-1}(\alpha^{1/2})=0.5\{-(N_1+N_2)+[(N_1+N_2)^2+4N_1N_2(\delta/(1-\delta)]^{1/2}\}.$$

Следовательно,$$n_{rat}=\frac{8\left\{\Phi^{-1}(\sqrt{\alpha})\right\}^2}{\left\{-\frac{N_1}{\sigma}-\frac{N_2}{\sigma}+\sqrt{\left(\frac{N_1}{\sigma}+\frac{N_2}{\sigma}\right)^2+4\frac{N_1N_2\delta}{\sigma^2(1-\delta)}}\right\}^2}.$$

При использовании "принципа уравнивания погрешностей" согласно [] $$\delta=1/2$$. При доверительной вероятности $$\alpha=0,95$$ имеем $$\sqrt{\alpha}=0,9747$$ и согласно [] $$\Phi^{-1}(\sqrt{\alpha})=1,954$$. Для этих численных значений$$n_{rat}=\frac{30,545}{\left\{-\frac{N_1}{\sigma}-\frac{N_2}{\sigma}+\sqrt{\left(\frac{N_1}{\sigma}+\frac{N_2}{\sigma}\right)^2+4\frac{N_1N_2}{\sigma^2}\right\}^2}.$$

Если $$\frac{N_1}{\sigma}=\frac{N_2}{\sigma}=0,1$$ то $$n_{rat}=4455$$. Если же $$\frac{N_1}{\sigma}=\frac{N_2}{\sigma}=0,5$$ то $$n_{rat}=178$$. Если первое из этих чисел превышает обычно используемые объемы выборок, то второе находится в "рабочей зоне" регрессионного анализа.

Парная регрессия. Наиболее простой и одновременно наиболее широко применяемый частный случай парной регрессии рассмотрим подробнее. Модель имеет вид$$y_i=ax_i+b+\epsilon_i,\;i=1,2,...,n.$$

Здесь $$x_i$$ - значения фактора (независимой переменной), $$y_i$$ - значения отклика (зависимой переменной), $$\varepsilon_i$$ - статистические погрешности, $$a, b$$ - неизвестные параметры, оцениваемые методом наименьших квадратов. Она переходит в модель (используем альтернативную запись линейной модели)$$y=X\beta+\varepsilon,$$ если положить$$y= \begin{pmatrix} y_1 \\ \ldots \\ y_n \end{pmatrix},\quad X= \begin{pmatrix} x_1 1 \\ \ldots \ldots \\ x_n 1 \end{pmatrix},\quad \varepsilon= \begin{pmatrix} \varepsilon_1 \\ \ldots \\ \varepsilon_n \end{pmatrix},\quad \beta= \begin{pmatrix} a \\ b \end{pmatrix}.$$

Естественно принять, что погрешности факторов описываются матрицей$$\Delta_x=\alpha= \begin{pmatrix} \Delta x_1 0 \\ \ldots \ldots \\ \Delta x_n 0 \end{pmatrix}= \begin{pmatrix} \alpha_1 0 \\ \ldots \ldots \\ \alpha_n 0 \end{pmatrix}.$$

В рассматриваемой модели интервального метода наименьших квадратов$$X=X_R+\alpha,\;y=y_R+\gamma=y_R+ \begin{pmatrix} \gamma_1 \\ \ldots \\ \gamma_n \end{pmatrix}$$ где $$X, y$$ - наблюдаемые (т.е. известные статистику) значения фактора и отклика, $$X_R, y_R$$ - истинные значения переменных, $$\alpha,\gamma$$ - погрешности измерений переменных. Пусть $$\beta^*$$ - оценка метода наименьших квадратов, вычисленная по наблюдаемым значениям переменных, $$\beta_R^*$$ - аналогичная оценка, найденная по истинным значениям. В соответствии с ранее проведенными рассуждениями$$\beta^*-\beta=[-(X_0^TX_0)^{-1}\Delta(X_0^TX_0)^{-1}X_0^T+(X_0^TX_0)^{-1}\alpha^T]y_0+(X_0^TX_0)^{-1}X_0^T\gamma$$ с точностью до бесконечно малых более высокого порядка по $$|\alpha|$$ и $$|\gamma|$$. В формуле (59) использовано обозначение $$\Delta=X_0^T\alpha+\alpha^TX_0$$. Вычислим правую часть в (59), выделим главный линейный член и найдем нотну.

Легко видеть, что$$X^TX= \begin{pmatrix} \sum x_i^2 \sum x_i \\ \sum x_i n \end{pmatrix},$$

где суммирование проводится от 1 до $$n$$. Для упрощения обозначений в дальнейшем до конца настоящего параграфа не будем указывать эти пределы суммирования. Из (60) вытекает, что$$(X^TX)^{-1}= \left. \begin{pmatrix} n -\sum x_i \\ -\sum x_i \sum x_i^2 \end{pmatrix} \right/ \left[ n\sum x_i^2-\left(\sum x_i\right)^2 \right].$$

Легко подсчитать, что$$X^T\alpha+\alpha^TX= \begin{pmatrix} 2\sum x_i\alpha_i \sum\alpha_i \\ \sum\alpha_i n \end{pmatrix}.$$

Положим$$S_0^2=\frac{1}{n}\sum(x_i-\overline{x})^2.$$

Тогда знаменатель в (61) равен $$n^2S_0^2$$. Из (61) и (62) следует, что$$(X^TX)^{-1}(X^T\alpha+\alpha^TX)=\dfrac{1}{n^2S_0^2} \begin{pmatrix} 2n\sum x\alpha-\sum x\sum\alpha n\sum\alpha \\ -2\sum x\sum x\alpha+\sum x^2\sum\alpha -\sum x\sum\alpha \end{pmatrix}.$$

Здесь и далее опустим индекс $$i$$, по которому проводится суммирование. Это не может привести к недоразумению, поскольку всюду суммирование проводится по индексу $$i$$ в интервале от 1 до $$n$$. Из (61) и (63) следует, что$$(X^TX)^{-1}(X^T\alpha+\alpha^T X)(X^TX)^{-1}=\dfrac{1}{n^4S_0^4}\begin{pmatrix}A B \\ C D \end{pmatrix},$$ где$$\begin{gathered} A=2n^2\sum x\alpha-2n\sum x\sum\alpha, \\ B=C=-2n\sum x\sum x\alpha + \left(\sum x\right)^2\sum\alpha+n\sum\alpha\sum x^2, \\ D=2\left(\sum x\right)^2\sum x\alpha-2\sum\alpha\sum x\sum x^2. \end{gathered}$$

Наконец, вычисляем основной множитель в (59)$$(X^TX)^{-1}(X^T\alpha+\alpha^TX)(X^TX)^{-1}X^T=\dfrac{1}{n^4S_0^4} \begin{pmatrix} z_{11}z_{12}\ldotsz_{1i}\ldotsz_{1n} \\ z_{21}z_{22}\ldotsz_{2i}\ldotsz_{2n} \end{pmatrix},$$ где$$z_{1i}=Ax_i+B, z_{2i}=Cx_i+D,\;i=1,2,...,n.$$

Перейдем к вычислению второго члена с $$\alpha$$ в (59). Имеем$$(X^TX)^{-1}\alpha^T=\dfrac{1}{n^2S_0^2} \begin{pmatrix} w_{11} \ldots w_{1i} \ldots w_{1n} \\ w_{21} \ldots w_{2i} \ldots w_{2n} \end{pmatrix},$$ где$$w_{1i}=n\alpha_i,\;w_{2i}=-\alpha\sum x,\;i=1,2,...,n.$$

Складывая правые части (65) и (67) и умножая на $$y$$, получим окончательный вид члена с $$\alpha$$ в (59):$$\{(X^TX)^{-1}(X^T\alpha+\alpha^TX)(X^TX)^{-1}X^T+(X^TX)^{-1}\alpha^T\}y={F\choose G},$$ где$$\begin{aligned} G=\left(\sum xy\right)\left(-2n\sum x\sum x\alpha+n\sum\alpha\sum x^2+\sum\alpha\left(\sum x\right)^2\right)/n^4S_0^4- \\ -\left(\sum y\alpha\right)\left(\sum x\right)/n^2S_0^2+\left(\sum y\right)\left(2\sum x\alpha\left(\sum x\right)^2-2\sum\alpha\sum x\sum x^2\right)/n^4S_0^4. \end{aligned}$$

Для вычисления нотны выделим главный линейный член. Сначала найдем частные производные. Имеем$$\begin{aligned} \frac{\partial F}{\partial\alpha_j}=\left(\sum xy\right)(2n^2x_j-2n\sum x)/n^4S_0^4+y_j/nS_0^2+ \\ +\left(\sum y\right)\left(n\sum x^2+\left(\sum x\right)^2-2n\left(\sum x\right)x_j\right)/n^4S_0^4; \end{aligned}$$ $$\begin{aligned} \frac{\partial G}{\partial\alpha_j}=\left(\sum xy\right)\left(-2n\left(\sum x\right)x_j+n\sum x^2+\left(\sum x\right)^2\right)/n^4S_0^4- \\ -y_i\left(\sum x\right)/n^2S_0^2+\left(\sum y\right)\left(2x_j\left(\sum x\right)^2-2\sum x\sum x^2\right)/n^4S_0^4. \end{aligned}$$

Если ограничения имеют вид$$|\alpha_j|\le\Delta,\;j=1,2,...,n,$$ то максимально возможное отклонение оценки $$a^*$$ параметра $$a$$ из-за погрешностей $$\alpha_j$$ таково:$$N_a(x)=\sum_{1\le j\le n}\left|\frac{\partial F}{\partial\alpha_j}\right|\Delta+O(\Delta^2),$$ где производные заданы формулой (70).

Пример 2. Пусть вектор $$(x,y)$$ имеет двумерное нормальное распределение с нулевыми математическими ожиданиями, единичными дисперсиями и коэффициентом корреляции $$\rho$$. Тогда$$\lim_{\Delta\rightarrow 0}\lim_{n\leftarrow\infty}\frac{N_a(x)}{\Delta}=\lim_{n\rightarrow\infty}\sum_{1\le j\le n}\left|\frac{\partial F}{\partial\alpha_j}\right|=M|2\rho x+y|=\sqrt{\frac{2(1+8\rho^2)}{\pi}}.$$

При этом$$\lim_{n\rightarrow\infty}\frac{\partial G}{\partial\alpha_j}=\rho,$$ следовательно, максимально возможному изменению параметра $$b^*$$ соответствует сдвиг всех $$x_i$$ в одну сторону, т.е. наличие систематической ошибки при определении $$x$$ -ов. В то же время согласно (71) значения $$\alpha_j$$ в асимптотике выбираются по правилу$$\alpha_j= \left\{ \begin{aligned} \Delta,\;2\rho x_j+y_j>0, \\ -\Delta,\;2\rho x_j+y_j\le 0. \end{aligned} \right\}$$

Таким образом, максимальному изменению $$a^*$$ соответствуют не те $$\alpha_j$$, что максимальному изменению $$b^*$$. В этом - новое по сравнению с одномерным случаем. В зависимости от вида ограничений на возможные отклонения, в частности, от вида метрики в пространстве параметров, будут "согласовываться" отклонения по отдельным параметрам. Ситуация аналогична той, что возникает в классической математической статистике в связи с оптимальным оцениванием параметров. Если параметр одномерен, то ситуация с оцениванием достаточно прозрачна - есть понятие эффективных оценок, показателем качества оценки является средний квадрат ошибки, а при ее несмещенности - дисперсия. В случае нескольких параметров возникает необходимость соизмерить точность оценивания по разным параметрам. Есть много критериев оптимальности (см., например, []), но нет признанных правил выбора среди них.

Вернемся к формуле (59). Интересно, что отклонения вектора параметров, вызванные отклонениями значений факторов $$\alpha$$ и отклика $$\gamma$$, входят в (59) аддитивно. Хотя$$\begin{aligned} \sup_{\alpha,\gamma}||\beta^*-\beta||\ne\sup_{\alpha}|\{-(X_0^TX_0)^{-1} \Delta(X_0^TX_0)^{-1}X_0^T+(X_0^TX_0)^{-1}\alpha^T\}y_0|+ \\ +\sup_\gamma|(X_0^TX_0)^{-1}X_0^T\gamma|, \end{aligned}$$ но для отдельных компонент (не векторов!) имеет место равенство.

В случае парной регрессии$$(X_0^TX_0)^{-1}X_0^T\gamma=\frac{1}{n^2S_0^2} \left(\sum\gamma_i\left(nx_i-\sum x\right);\;\sum\gamma_i\left(-x_i\sum x+\sum x^2\right)\right)^T.$$

Из формул (68), (69) и (72) следует, что$$\beta^*-\beta={{a^*(X,y)-a^*(X_0,y_0)}\choose{b^*(X,y)-b^*(X_0,y_0)}}={{F+F_1}\choose{G+G_1}},$$ где $$F$$ и $$G$$ определены в (69), а$$F_1=\frac{1}{n^2S_0^2}\left(n\sum\gamma x-\sum x\sum\gamma\right),\; G_1=\frac{1}{n^2S_0^2}\left(\sum\gamma\sum x^2-\sum\gamma x\sum x\right).$$

Итак, продемонстрирована возможность применения основных подходов статистики интервальных данных в регрессионном анализе. Пример практического применения этих подходов при оценивании зависимости затрат от объема выпуска продукции дан в статье: Гуськова Е.А., Орлов А.И. Интервальная линейная парная регрессия (обобщающая статья). - Журнал "Заводская лаборатория". 2005. Т.71. No.3. С.57-63.

12.5. Интервальный дискриминантный анализ

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

В линейном дискриминантном анализе правило принятия решений основано на линейной функции $$f(x)$$ от распознаваемого вектора $$x\in R^k$$. Рассмотрим для простоты случай двух классов. Правило принятия решений определяется константой $$C$$ - при $$f(x)>C$$ распознаваемый объект относится к первому классу, при $$f(x)<C$$ - ко второму.

В первоначальной вероятностной модели Р.Фишера предполагается, что классы заданы обучающими выборками объемов $$N_1$$ и $$N_2$$ соответственно из многомерных нормальных распределений с разными математическими ожиданиями, но одинаковыми ковариационными матрицами. В соответствии с леммой Неймана-Пирсона, дающей правило принятия решений при поверке статистических гипотез, дискриминантная функция является линейной. Для ее практического использования теоретические характеристики распределения необходимо заменить на выборочные. Тогда дискриминантная функция приобретает следующий вид$$f(x)=\left(x-\frac12(\overline{x}_1+\overline{x}_2)\right)S^{-1}(\overline{x}_1-\overline{x}_2).$$

Здесь $$\overline{x}_1$$ - выборочное среднее арифметическое по первой выборке $$x_{\alpha}^{(1)}, \; \alpha=1,2,...,N_1$$ а $$\overline{x}_2$$ - выборочное среднее арифметическое по второй выборке $$x_{\beta}^{(2)}, \; \beta=1,2,...,N_2$$. В роли $$S$$ может выступать любая состоятельная оценка общей для выборок ковариационной матрицы. Обычно используют следующую оценку, естественным образом сконструированную на основе выборочных ковариационных матриц:$$S=\frac{\sum\limits_{\alpha=1}^{N_1}(x_{\alpha}^{(1)}-\overline{x}_1)(x_{\alpha}^{(1)}-\overline{x}_1)^T + \sum\limits_{\beta=1}^{N_2}(x_{\beta}^{(2)}-\overline{x}_2)(x_{\beta}^{(2)}-\overline{x}_2)^T}{N_1+N_2-2}.$$

В соответствии с подходом статистики интервальных данных считаем, что специалисту по анализу данных известны лишь значения с погрешностями$$y_{\alpha}^{(1)}=x_{\alpha}^{(1)}+\varepsilon_{\alpha}^{(1)},\;\alpha=1,2,...,N_1,\; y_{\beta}^{(2)}=x_{\beta}^{(2)}+\varepsilon_{\beta}^{(2)},\;\beta=1,2,...,N_2.$$

Таким образом, вместо $$f(x)$$ статистик делает выводы на основе искаженной линейной дискриминантной функции $$f_1(x)$$, в которой коэффициенты рассчитаны не по исходным данным $$x_{\alpha}^{(1)},x_{\beta}^{(2)}$$, а по искаженным погрешностями значениям $$y_{\alpha}^{(1)},y_{\beta}^{(2)}$$.

Это - модель с искаженными параметрами дискриминантной функции. Следующая модель - такая, в которой распознаваемый вектор $$x$$ также известен с ошибкой. Далее, константа $$C$$ может появляться в модели различными способами: задаваться априори абсолютно точно; задаваться с какой-то ошибкой, не связанной с ошибками, вызванными конечностью обучающих выборок; рассчитываться по обучающим выборкам, например, с целью уравнять ошибки классификации, т.е. провести плоскость дискриминации через середину отрезка, соединяющего центры классов. Итак - целый спектр моделей ошибок.

На какие статистические процедуры влияют ошибки в исходных данных? Здесь тоже много постановок. Можно изучать влияние погрешностей измерений на значения дискриминантной функции $$f$$, например, в той точке, куда попадает вновь поступающий объект $$x$$. Очевидно, случайная величина $$f(x)$$ имеет некоторое распределение, определяемое распределениями обучающих выборок. Выше описана модель Р.Фишера с нормально распределенными совокупностями. Однако реальные данные, как правило, не подчиняются нормальному распределению []. Тем не менее линейный статистический анализ имеет смысл и для распределений, не являющихся нормальными (при этом вместо свойств многомерного нормального распределения приходится опираться на многомерную Центральную предельную теорему и теорему о наследовании сходимости []). В частности, приравняв метрологическую ошибку, вызванную погрешностями исходных данных, и статистическую ошибку, получим условие, определяющее рациональность объемов выборок. Здесь два объема выборок, а не один, как в большинстве рассмотренных постановок статистики интервальных данных. С подобным мы сталкивались ранее при рассмотрении двухвыборочного критерия Смирнова.

Естественно изучать влияние погрешностей исходных данных не при конкретном $$x$$, а для правила принятия решений в целом. Может представлять интерес изучение характеристик этого правила по всем $$x$$ или по какой-либо области возможных значений $$x$$. Более интересно рассмотреть показатель качества классификации, связанный с пересчетом на модель линейного дискриминантного анализа [].

Математический аппарат изучения перечисленных моделей развит выше в предыдущих пунктах настоящей главы. Некоторые результаты приведены в []. Из-за большого объема выкладок ограничимся приведенными здесь замечаниями.

12.6. Интервальный кластер-анализ

Кластер-анализ, как известно [], имеет целью разбиение совокупности объектов на группы сходных между собой. Многие методы кластер-анализа основаны на использовании расстояний между объектами. (Степень близости между объектами может измеряться также с помощью мер близости и показателей различия, для которых неравенство треугольника выполнено не всегда.) Рассмотрим влияние погрешностей измерения на расстояния между объектами и на результаты работы алгоритмов кластер-анализа.

С ростом размерности $$p$$ евклидова пространства диагональ единичного куба растет как $$\sqrt{p}$$. А какова погрешность определения евклидова расстояния? Пусть двум рассматриваемым векторам соответствуют $$X_0=(x_1,x_2,...,x_p)$$ и $$Y_0=(y_1,y_2,...,y_p)$$ - векторы размерности $$p$$. Они известны с погрешностями $$\varepsilon=(\varepsilon_1,\varepsilon_2,...,\varepsilon_p)$$ и $$\delta=(\delta_1,\delta_2,...,\delta_p)$$, т.е. статистику доступны лишь векторы $$X=X_0+\varepsilon,\; Y=Y_0+\delta$$ Легко видеть, что$$\rho^2(X,Y)=\rho^2(X_0,Y_0)+2\sum_{1\le i\le p}(x_i-y_i)(\varepsilon_i-\delta_i)+\sum_{1\le i\le p}(\varepsilon_i-\delta_i)^2.$$

Пусть ограничения на абсолютные погрешности имеют вид$$|\varepsilon_i|\le\Delta,\;|\delta_i|\le\Delta,\;i=1,2,...,n.$$

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

Пусть $$p\Delta^2\rightarrow 0$$. Тогда последнее слагаемое в (73) не превосходит $$4p\Delta^2$$ поэтому им можно пренебречь. Тогда из (73) следует, что нотна евклидова расстояния имеет вид$$N_{\rho^2}(X_0,Y_0)=4\sum_{1\le i\le p}|x_i-y_i|\Delta$$ с точностью до бесконечно малых более высокого порядка. Если случайные величины $$|x_i-y_i|$$ имеют одинаковые математические ожидания и для них справедлив закон больших чисел (эти предположения естественны, если переменные перед применением кластер-анализа стандартизованы), то существует константа $$C$$ такая, что$$N_{\rho^2}(X_0,Y_0)=Cp\Delta$$ с точностью до бесконечно малых более высокого порядка при малых \Delta больших $$p$$ и $$p\Delta^2\rightarrow 0$$.

Из рассмотрений настоящего пункта вытекает, что$$\rho(X,Y)=\rho(X_0,Y_0)+\theta\frac{Cp\Delta}{2\rho(X_0,Y_0)}$$ при некотором $$\theta$$ таком, что $$|\theta|<1$$.

Какое минимальное расстояние является различимым? По аналогии с определением рационального объема выборки при проверке гипотез предлагается уравнять слагаемые в (74), т.е. определять минимально различимое расстояние $$\rho_{min}$$ из условия$$\rho_{min}=\frac{Cp\Delta}{2\rho_{min}},\;\rho_{min}=\sqrt{\frac{Cp\Delta}{2}}.$$

Естественно принять, что расстояния, меньшие \rho_{min}, не отличаются от 0, т.е. точки, лежащие на расстоянии $$\rho\le\rho_{min}$$, не различаются между собой.

Каков порядок величины $$C$$? Если $$x_i$$ и $$y_i$$ независимы и имеют стандартное нормальное распределение с математическим ожиданием 0 и дисперсией 1, то, как легко подсчитать, $$M|x_i-y_i|=2/\sqrt{\pi}=1,13$$ и соответственно $$C = 4,51$$. Следовательно, в этой модели$$\rho_{min}=1,5\sqrt{p\Delta}.$$

Формула (75) показывает, что хотя с ростом размерности пространства $$p$$ растет диаметр (длина диагонали) единичного куба - естественной области расположения значений переменных, с той же скоростью растет и естественное квантование расстояния с помощью порога неразличимости , т.е. увеличение размерности (вовлечение новых переменных), вообще говоря, не улучшает возможности кластер-анализа.

Можно сделать выводы и для конкретных алгоритмов. В дендрограммах (например, результатах работы иерархических агломеративных алгоритмах ближнего соседа, дальнего соседа, средней связи) можно порекомендовать склеивать (т.е. объединять) уровни, отличающиеся менее чем на $$\rho_{min}$$. Если все уровни склеятся, то можно сделать вывод, что у данных нет кластерной структуры, они однородны. В алгоритмах типа "Форель" центр тяжести текущего кластера определяется с точностью $$\pm\Delta$$ по каждой координате, а порог для включения точки в кластер (радиус шара $$R$$ ) из-за погрешностей исходных данных может измениться согласно (74) на$$\pm\frac{2,25}{R}p\Delta.$$

Поэтому кроме расчетов с $$R$$ рекомендуется провести также расчеты с радиусами $$R_1$$ и $$R_2$$, где$$R_1=R\left(1-\frac{2,25}{R^2}p\Delta\right),\;R_2=R\left(1+\frac{2,25}{R^2}p\Delta\right),$$ и сравнить полученные разбиения. Быть адекватными реальности могут только выводы, общие для всех трех расчетов. Эти рекомендации развивают общую идею [] о целесообразности проведения расчетов при различных значениях параметров алгоритмов с целью выделения выводов, инвариантных по отношению к выбору конкретного алгоритма.

12.7. Статистика интервальных данных и оценки погрешностей характеристик финансовых потоков инвестиционных проектов

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

Основная идея формулируется так. Все знают, что любое инженерное измерение проводится с некоторой погрешностью. Эту погрешность обычно приводят в документации и учитывают при принятии решений. Ясно, что и любое экономическое измерение также проводится с погрешностью. А вот какова она? Необходимо уметь ее оценивать, поскольку ошибки при принятии экономических решений обходятся дорого.

Например, как принимать решение о выгодности или невыгодности инвестиционного проекта? Как сравнивать инвестиционные проекты между собой? Как известно, для решения этих задач используют такие экономические характеристики, как NPV (Net Present Value) - чистая текущая стоимость (этот термин переводится с английского также как чистый дисконтированный доход, чистое приведенное значение и др.), внутренняя норма доходности, срок окупаемости, показатели рентабельности и др.

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

В настоящее время широко используются различные теоретические подходы к сравнению инвестиционных проектов и облегчающие расчеты компьютерные системы, в частности, Project Expert, COMFAR, PROPSIN, Альт-Инвест, ТЭО-ИНВЕСТ. Однако ряд важных моментов в них не учтен.

Введем основные понятия. Дисконт-функция как функция от времени показывает, сколько стоит для фирмы 1 руб. в заданный момент времени, если его привести к начальному моменту. Если дисконт-функция - константа для разных отраслей, товаров и проектов, то эта константа называется дисконт-фактором, или просто дисконтом. Дисконт-функция определяется совместным действием различных факторов, в частности, реальной процентной ставки и индекса инфляции. Реальная процентная ставка описывает "нормальный" рост экономики (т.е. без инфляции). В стабильной ситуации доходность от вложения средств в различные отрасли, в частности, в банковские депозиты, примерно одинакова. Сейчас она, по оценке ряда экспертов, составляет около 12%. Итак, нынешний 1 руб. превращается в 1,12 руб. через год, а потому 1 руб. через год соответствует 1/1,12 = 0,89 руб. Сейчас - это и есть максимум дисконта.

Обозначим дисконт буквой $$C$$. Если $$q$$ - банковский процент (плата за депозит), т.е. вложив в начале года в банк 1 руб., в конце года получим $$(1+ q)$$ руб., то дисконт определяется по формуле $$C=1/(1+q)$$. При таком подходе полагают, что банковские проценты одинаковы во всех банках. Более правильно было бы считать $$q$$, а потому и $$C$$, нечисловыми величинами, а именно, интервалами $$[q1; q2]$$ и $$[C1; C2]$$. Следовательно, экономические выводы должны быть исследованы на устойчивость (применяют и термин "чувствительность" ) по отношению к возможным отклонениям.

Как функцию времени $$t$$ дисконт-функцию обозначим $$C(t)$$. При постоянстве дисконт-фактора имеем $$C(t) = Сt$$. Если $$q = 0,12, C = 0,89$$, то 1 руб. за 2 года превращается в 1,122 = 1,2544, через 3 - в 1,4049. Итак, 1 руб., получаемый через 2 года, соответствует 1/1,2544=0,7972 руб., т.е. 79,72 коп. сейчас, а 1 руб., обещанный через 3 года, соответствует 0,71 руб. сейчас. Другими словами, $$C(2) = 0,80$$, а $$C(3) = 0,71$$. Если дисконт-фактор зависит от времени, в первый год равен $$C1$$, во второй - $$C2$$, в третий - $$C3$$,..., в $$t$$ -ый год - $$Ct$$, то $$C(t)=C_1C_2C_3...C_t$$.

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

Примитивный способ расчета срока окупаемости состоит в делении объема вложений $$A$$ на ожидаемый ежегодный доход $$B$$. Тогда срок окупаемости равен $$A/B$$. Этот способ некорректен. Если дисконт-фактор равен $$C$$, то максимально возможный суммарный доход равен$$BC+BC^2+BC^3+BC^4+BC^5+...=BC(1+C+C^2+C^3+C^4+...) = BC / (1-C).$$

Если $$A/B$$ меньше $$C/(1-C)$$, то можно рассчитать срок окупаемости проекта, но он будет больше, чем $$A/B$$. Если же $$A/B$$ больше или равно $$C/(1-C)$$, то проект не окупится никогда. Поскольку максимум $$C$$ равен 0,89, то проект не окупится никогда, если $$A/B$$ не меньше 8,09.

Пусть вложения равны 1 млн руб., ежегодная прибыль составляет 500 тыс., т.е. $$A/B=2$$, дисконт-фактор $$C=0,8$$. При примитивном подходе (при $$C=1$$ ) срок окупаемости равен 2 годам. А на самом деле? За $$k$$ лет будет возвращено$$BC(1+C+C^2+C^3+C^4+...+C^k)=BC(1-C^{k+1})/(1-C).$$

Срок окупаемости $$k$$ получаем из уравнения $$1=0,5\times 0,8(1-0,8^{k+1})/(1 - 0,8)$$, откуда $$k= 2,11$$. Он оказался равным 2,11 лет, т.е. увеличился примерно на 6 недель. Это немного. Однако если $$B = 0,2$$, то имеем уравнение $$1=0,2\times 0,8(1-0,8^{k+1})/ (1-0,8)$$. У этого уравнения нет корней, поскольку $$A/B=5>C/(1-C)=0,8/(1-0,8)=4$$. Проект не окупится никогда. Прибыль можно ожидать лишь при $$A/B<4$$. Рассмотрим промежуточный случай, $$B=0,33$$, с "примитивным" сроком окупаемости 3 года. Тогда имеем уравнение $$1=0,33\times 0,8(1-0,8^{k+1})/(1-0,8)$$, откуда $$k = 5,40$$.

Рассмотрим финансовый поток $$a(0), a(1), a(2), a(3), ... , a(t), ..$$.. (для простоты примем, что платежи или поступления происходят раз в год). Выше рассмотрен поток с одним платежом $$a(0)=(-A)$$ и дальнейшими поступлениями $$a(1) = a(2) = a(3) = ... = a(t) = .... = B$$. Чистая текущая стоимость (Net Present Value, сокращенно NPV ), рассчитывается для финансового потока путем приведения затрат и поступлений к начальному моменту времени:$$NPV = a(0) +a(1)C(1)+a(2)C(2)+a(3)C(3)+...+ a(t)C(t) + ...,$$ где $$C(t)$$ - дисконт-функция. В простейшем случае, когда дисконт-фактор не меняется год от года и имеет вид $$C=1/(1+q)$$, формула для NPV конкретизируется:$$NPV=NPV(q)=a(0)+a(1)/(1+q) + a(2)/(1+q)2 +a(3)/(1+q)3 +...+a(t)/(1+q)t + ...$$

Пусть, например, $$a(0)= -10, a(1)=3, a(2)=4, a(3)=5$$. Пусть $$q=0,12$$, тогда$$NPV(0,12)=-10+3\times 0,89+4\tomes 0,80+5\times 0,71=-10+2,67 + 3,20+3,55= -0,58.$$

Итак, проект невыгоден для вложения капитала, поскольку NPV (0,12) отрицательно. При отсутствии дисконтирования (при $$C=1, q=0$$ ) вывод иной:$$NPV(0) = - 10 + 3 + 4 + 5 = 2,$$ проект выгоден.

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

Приведем пример исследования NPV на устойчивость (чувствительность) к малым отклонениям значений дисконт-функции. Для этого надо найти максимально возможное отклонение NPV при допустимых отклонениях значений дисконт-функции (или, если угодно, значений банковских процентов). В качестве примера рассмотрим$$\begin{gathered} NPV = NPV (a(0), a(1), C(1), a(2), C(2), a(3), C(3))= \\ = a(0) + a(1)C(1) + a(2)C(2) + a(3)C(3). \end{gathered}$$

Предположим, что изучается устойчивость (чувствительность) для ранее рассмотренных значений$$a(0)=-10, a(1)=3, a(2)=4, a(3)=5, C(1)=0,89, C(2)=0,80, C(3)=0,71.$$

Пусть максимально возможные отклонения $$C(1), C(2), C(3)$$ равны $$\pm 0,05$$. Тогда, максимум значений NPV равен$$NPVmax = -10+3\times 0,94+4\times 0,85+5\times 0,76 = -10+ 2,82 + 3,40 + 3,80 = 0,02,$$ в то время как минимум значений NPV есть$$NPVmin = -10+3\times 0,84+4\times 0,75+5\times 0,66 = -10 +2,52 +3,00+3,30 = -1,18.$$

Для NPV получаем интервал от (-1,18) до (+0,02). В нем есть и положительные, и отрицательные значения. Следовательно, нет однозначного заключения - проект убыточен или выгоден. Для принятия решения не обойтись без экспертов.

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

12.8. Место статистики интервальных данных (СИД) в прикладной статистике

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

Нечеткость и СИД. С формальной точки зрения описание нечеткости интервалом - это частный случай описания ее нечетким множеством. В СИД функция принадлежности нечеткого множества имеет специфический вид - она равна 1 в некотором интервале и 0 вне него. Такая функция принадлежности описывается всего двумя параметрами (границами интервала). Эта простота описания делает математический аппарат СИД гораздо более прозрачным, чем аппарат теории нечеткости в общем случае. Это, в свою очередь, позволяет продвинуться дальше, чем при использовании функций принадлежности произвольного вида.

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

Математическая статистика и СИД. Как уже отмечалось, математическая статистика и СИД отличаются тем, в каком порядке делаются предельные переходы $$n\rightarrow\infty$$ и $$\Delta\rightarrow 0$$. При этом СИД переходит в математическую статистику при $$\Delta=0$$. Правда, тогда исчезают основные особенности СИД: нотна становится равной 0, а рациональный объем выборки - бесконечности. Рассмотренные выше методы СИД разработаны в предположении, что погрешности малы (но не исчезают) и объем выборки велик. СИД расширяет классическую математическую статистику тем, что в исходных статистических данных каждое число заменяет интервалом. С другой стороны, можно считать СИД новым этапом развития математической статистики.

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

]. К этой модели СИД не имеет отношения.

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

], показали, что это подозрение не подтверждается. Влияние погрешностей измерений по порядку такое же, только вместо максимально возможного отклонения (нотны) приходится рассматривать математическое ожидание соответствующего отклонения (см. выше). Подчеркнем, что применение в СИД вероятностно-статистических моделей погрешностей не менее перспективно, чем минимаксных.

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

Анализ чувствительности и СИД. При анализе чувствительности, как и в СИД, рассчитывают производные по используемым переменным, или непосредственно находят изменения при отклонении переменной на +10% от базового значения. Однако этот анализ делают по каждой переменной отдельно. В СИД все переменные рассматриваются совместно, и находится максимально возможное отклонение (нотна). При малых погрешностях удается на основе главного члена разложения функции в многомерный ряд Тейлора получить удобную формулу для нотны. Можно сказать, что СИД - это многомерный анализ чувствительности.

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

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