Прогнозирование экстремальных значений по выборке

[Александрович]

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

[gris]

Александрович писал(а):

А когда говорят о "правиле 3-х сигм"... основанно на уверенности... чрезвычайно мала

Не уверенность, а холодный расчёт; не чрезвычайно мала, а равна ровно $1-2\Phi (3)$

[Александрович]

Совершенно верно! Но это и пугает, что "холодный расчет" построен на оценках параметров по выборке, а делается вид, что на знании самих параметров распределения. А равняется она ровно $1-2\Phi(3+ds(n))$, где $n-$ объем выборки.

[GAA]

Если параметры оцениваются по выборке, то, для указания интервала, в котором находится доля распределения не меньше, чем $\beta$, используется понятие толерантного интервала.

Def. Пусть случайная величина $Y$ не зависит от $X$; $X$ и $Y$ имеют функцию распределения $F_{\theta}$. Двусторонним толерантным интервалом уровня $\gamma$ называется такой интервал $[\underline{L}(X), \overline{L}(X)]$, что $\mathsf P_{\theta} \left\{ \mathsf P_{\theta} \{ \underline{L}(X) \le Y \le \overline{L}(X) | X \} \ge \beta \right\} \ge \gamma$.

Для нормального распределения толерантный интервал строят в виде $[M - \lambda S, M + \lambda S]$. Для некоторых значений объемов выборки, $\gamma$ и $\beta$ множители $\lambda$ приведены в [4, (с. 45, 237)]. О приближенном вычислении значений этих констант можно посмотреть в [7, §§20.37—20.40]. Множители для односторонних интервалов выражаются через квантили нецентрального распределения Стьюдента (см. [8, §10.4 Толерантные интервалы]). Квантили нецентрального распределения Стьюдента можно вычислить в Maple 12 (см. Quantile и NonCentralStudentT); в 11-ой не знаю, в 10-ой и более ранних версиях встроенной функции не было.

ref.
7. Кендалл М., Стьюарт А. Статистические выводы и связи. — М.: Наука, 1973.
8. Закс Ш. Теория статистических выводов. — М.: Мир, 1975.

[Александрович]

Уважаемый GAA вычислил что $t_{5,10,0.85}=6.779$, тогда значение стандартного нормального распределения $X = 3.032$, что соответствует $0.9988$ квантили распределения с.в. или $0.9879$ квантили распределения максимального значения (вместо $0.85$ для известных м.о. и ско). Теперь вопрос: Если исходная выборка принадлежит не нормальному, а иному виду распределения (также двухпараметрическому), получу ли я, подставив в формулу функции распределения вместо параметров их оценки, расчитанные по выборке, для $0.9988$ квантили значение - $X_{(n)}$?

[GAA]

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

0. Задача (2-я модификация). Пусть $X_i$, $i=1,\ldots, L$, $Y_i$, $i=1,\ldots, N$ — независимые, $(\mu, \sigma)$-нормально распределенные случайные величины. Требуется построить верхний $\gamma$-доверительный толерантный интервал для случайной величины $X_{(N)} = \max\limits_{i=1, \ldots, N} Y_i$ (на другом языке — доверительный интервал для квантили $\xi_p$ уровня $p$ распределения случайной величины $X_{(N)}$) в виде $(-\infty, M + \lambda S)$, т.е.
$\mathsf P_{\mu, \sigma} \left\{ \mathsf P_{\mu, \sigma} \{X_{(n)} < M + \lambda S| X_i, i=1,\ldots, L\} = p \right\} = \gamma$.
Здесь $\xi_p = F^{-1}(p)$, где $F$ — функция распределения случайной величины $X_{(N)}$.

1. Решение. Как уже выше в этой теме указывалось, $F(x) = \Phi^N (x)$, где $\Phi (x)$ — функция стандартного нормального распределения. Следовательно, $\xi_p = u(\sqrt [N] p)$, где $u(\beta)$ — квантиль уровня $\beta$ стандартного нормального распределения. Без ограничения общности можно считать, что $X_i$, $Y_i$ имеют $\mu=0$, $\sigma=1$. По сути повторяя рассуждения примера 10.9 [8], получим $\lambda = \frac {t_{\gamma} [L-1, \sqrt L u(\beta)]} {\sqrt L}$, где $ t_{\gamma} [L-1, \sqrt L u(\beta)]$ — квантиль уровня $\gamma$ нецентрального распределения Стьюдента c числом степеней свободы равным $\nu = L-1$ и параметром нецентральности $\delta = \sqrt L u(\beta)$.
Для $L=3$, $N = 2$, $p = 0.85$, $\gamma = 0.9$, имеем $\xi_{0.85} \approx 1.418341423$, $t_{0.9}[2, \sqrt 3 u(\sqrt [N] {0.85})] \approx 8.056133265$.
Расчеты выполнены в Maple 12:

> with(Statistics):
> L:= 3: N:= 2: p:= 0.85:
> epsilon:= 0.1: # epsilon = 1 - gamma
> u:= Quantile(Normal(0, 1), p^(1/N));
                   u := 1.418341423
> Quantile(NonCentralStudentT(L-1, sqrt(L)*u), 1-epsilon);
                   8.056133265

2. Демонстрационный пример. Пусть «составной» эксперимент включает:

1. Получение выборки объемом L; вычисление по подвыборке объема L частных значений оценок m, S.
2. i) Обнуления счетчика b оценки p.
ii) Выполнение серии из Length экспериментов, каждый из которых состоит из: генерации выборки объемом $N$; вычисления $x_{N}$; и увеличения на единицу b в случае выполнения неравенства $x_{N} < m + \frac {t}{\sqrt{L}} S$.
iii) По завершению серии экспериментов b/Length — оценка $p$ — сравнивается с заданным уровнем $p$.

Тогда, при достаточно большом числе повторений (Rep) составного эксперимента и достаточно большой длине Length серии экспериментов в каждом составном эксперименте, оценка gamma*=Count/Rep накрыть толерантным интервалом заданную долю p распределения будет приближенно равна заданному доверительному уровню $\gamma$. Для демонстрации этого был написан простой пример в среде Delphi 5.0. Для того, чтобы не загромождать текст, «генератор» нормально распределенных случайных величин вынесен в отдельный модуль NormRND [текст этого модуля (используемого и в предыдущем демонстрационном примере) помещен в приложении к этому сообщению].

uses NormRnd;
const
 p = 0.85;
 gamma = 0.9;
 t : Double = 8.056133265; { t_0.9 [2, sqrt(3)*u(0.85^(1/2))] = 8.056133265}
 L = 3;
 N = 2;
 Length = 100000;
 mu: Double = 1;
 sigma: Double  = 1;
 Rep = 100000;
var
 X: array[1..L+N] of Double;
 m, s, xN : Double;
 Count, b : LongInt;
 i, j, k  : LongInt;
 C        : Double;
begin
 Count:= 0; InitN1DRnd(mu, sigma); C:= t/sqrt(L);
 for j:= 1 to Rep
  do begin
      for i:= 1 to L do X[i]:= N1DRnd;
      m:= 0; for i:= 1 to L do m:= m + X[i]; m:= m/L;
      s:= 0; for i:= 1 to L do s:= s + sqr(X[i]-m); s:= sqrt(s/(L-1));
      b:= 0;
      for k:= 1 to Length
       do begin
           for i:= 1 to N do X[i]:= N1DRnd;
           xN := X[1]; for i:= 2 to N do if X[i] > xN then xN:= X[i];
           if xN < m + C*s then inc(b);
          end;
      if b/Length >= p then inc(Count);
     end;
 writeln('gamma* =',  Count / Rep, 'Pres kay "Enter"'); readln;
end.

При Rep = Length = 100000, получено $\gamma^* \approx 0.901$.

3. По поводу обобщения на произвольное двухпараметрическое распределение. Описанный способ построения верхнего толерантного интервала существенно опирается на то, что случайные величины $X_i$ являются независимыми и одинаково нормально распределенными. Для иных двухпараметрических распределений таким образом построить толерантный интервал нельзя.

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

Unit NormRND;

interface
procedure InitN1DRnd(Mean, StdDev : Double);
function N1DRnd: Double;

implementation
type
 TNRndParam  = record
               FNext: Boolean;
               Next : Double
              end;
 TN1DRndParam= record
                Mean: Double;
                StdDev: Double;
               end;
var
 NRndParam  : TNRndParam;
 N1DRndParam: TN1DRndParam;

procedure InitNRnd;
begin NRndParam.FNext := False; end;

function NRnd: Double;
var
  V1, V2, S, sqrtS : Double;
begin with NRndParam do begin
If Not FNext
 then
  begin
   repeat
    V1 := 2 * Random - 1; V2 := 2 * Random - 1;
    S := V1 * V1 + V2 * V2
   until s < 1;
    sqrtS := Sqrt(-2 * Ln(S) / S);
    NRnd := V1 * sqrtS; Next := V2 * sqrtS; FNext := True;
  end
 else
  begin
   NRnd := Next; FNext := False
  end
end {with} end; {function}

procedure InitN1DRnd;
begin N1DRndParam.Mean := Mean; N1DRndParam.StdDev := StdDev; InitNRnd end;

function N1DRnd: Double;
begin with N1DRndParam do begin N1DRnd := Mean + StdDev * NRnd end {with} end;

end. 

[Александрович]

Цитата:

Понять написанное Вами, Александрович, в трех предыдущих сообщениях — невозможно

Последнее сообщение никак не связано с двумя остальными. О чем я пытался сказать в последнем.
1. Вы показали как вычисляется $t_{5,10,0,85}$ и привели полученное значение $6.779$.
2. С учетом этого я расчитал вероятность нормального распределения, соответствующую $0,85$ квантили распределения наибольшей с.в. она получилась равной $0.9988$.
3. Теперь подставляя это значение в обратную функцию нормального распределения с выборочными оценка м.о. и с.к.о. я получу значения совпадающие с Вашими.
А теперь вопрос. Если функция распределения с.в. не нормальная, это будет работать?
Спасибо за модификацию, правда я ничего не понимаю в программировании.

[GAA]

Александрович писал(а):

2. С учетом этого я расчитал квантиль нормального распределения, соответствующая $0,85$ квантили распределения наибольшей с.в. она получилась равной $0.9988$.

Запишите подробно как вычисляли.

Добавлено спустя 5 минут 12 секунд:

Александрович писал(а):

Спасибо за модификацию, правда я ничего не понимаю в программировании.

Я знаю. Программы для тех, кто владеет "Паскалем". Читатели, которые им не владеют, многое не потеряют: суть изложена в пунктах 0 и 1.

[Александрович]

Я исправил ошибку, вместо $5$ подставлял $10$. Теперь стало понятно?

[GAA]

Александрович писал(а):

Теперь стало понятно?

Нет. Запишите подробно как вычисляли.

[Александрович]

$t_{5,10,0.85}=6.779$, -это для с.к.о. среднего, а для с.к.о. с.в. нужно разделить на корень из $5$, тогда $X = 3.032$, что после подстановки в стандартную нормальную ф.р. соответствует $0.9988$ вероятности с.в. или, если возвести в степень $10$ ($0.9879$) вероятности максимального значения (вместо $0.85$ для известных м.о. и ско).

[GAA]

Александрович писал(а):

$t_{5,10,0.85}=6.779$, -это для с.к.о. среднего, а для с.к.о. с.в. нужно разделить на корень из $5$, тогда $X = 3.032$, что после подстановки в стандартную нормальную ф.р. соответствует $0.9988$

На основании чего и с какой целью $\frac{t_{5,10,0.85}}{\sqrt 5}$ подставляется в функцию нормального распределения?

[Александрович]

Прошу меня великодушно извинить, думаю об одном, а пишу другое.
Ляпы, которые обнаружил, исправил. Теперь

Цитата:

На основании чего и с какой целью ... подставляется в функцию нормального распределения?

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

[GAA]

Александрович писал(а):

Определить значение вероятности для $0.85$ квантили наибольшего значения с.в.

Квантиль — детерминированная величина (не путать со случайной величиной принимающей постоянное значение). Поэтому: нельзя говорить о «значении вероятности для 0.85 квантили». Напишите четко и подробно: вероятность какого события Вы пытаетесь определить.

[Александрович]

Медиана распределения - детерминированная величина, а выборочная медиана - случайная. Поэтому вместо вероятности $0.5$ может получиться вероятность $0.7$ тогда я говорю, что выборочная квантиль $0.7$ соответствует квантили $0.5$.