Поправка Аспина-Уэлча

[AndreyL]

Если две выборки $x$ и $y$ взяты из одной генеральной совокупности, подчиняющейся нормальному распределению (т.е. истинные центры и дисперсии двух выборок равны), то статистика $t=\frac{\bar{x}-\bar{y}}{\sqrt{s^2 \frac{(n_x+n_y)}{n_x n_y}}}$ подчиняется распределению Стьюдента с $n_x+n_y-2$ степенями свободы. Здесь $n_x$ и $n_y$ - объемы соответствующих выборок, обобщенная дисперсия $s^2=\frac{s_x^2 (n_x-1) +s_y^2 (n_y-1)}{n_x+n_y-2}$. Это несложно проверить численно: моделируем две выборки из нормального распределения с одинаковыми центрами и стандартами, считаем статистику $t$ и так несколько тысяч раз - получаем выборку статистики с законом распределения Стьюдента с $n_x+n_y-2$ степенями свободы.
Если же выборки взяты из разных генеральных совокупностей с разными дисперсиями но одинаковыми центрами, то, в соответствии с поправкой Аспина-Уэлча, та-же самая статистика подчиняется распределению Стьюдента с $\frac{\left(\frac{s_x^2}{n_x}+\frac{s_y^2}{n_y} \right)^2}{\frac{s_x^4}{(n_x-1) n_x^2}+\frac{s_y^4}{(n_y-1) n_y^2}}$ степенями свободы. Как теперь это проверить численно? Сложность в том, что для каждой выборки (для каждой разницы выборочных дисперсий) свой закон распределения (свои степени свободы).

[GAA]

Как и всегда при проверке приближенного критерия, вначале следует исследовать вероятность ошибки первого рода. Например, можно выбрать несколько наиболее часто используемых уровней (0.20, 0.10, 0.05, 0.01) и найти методом статистических испытаний (Монте-Карло) оценки вероятностей первого рода, выбирая различные значения объемов выборок и дисперсий. Т.к. критерий Уэлча (Welch’s t test) является приближенным, то для сравнения имеет смысл (формально) применить критерий Стьюдента или другой уместный критерий.

[AndreyL]

Извиняюсь, по-видимому, не понятно задал вопрос. Проблема именно в генерировании выборок – точное значение выборочной дисперсии непредсказуемо, а отношение этих значений задает закон распределения. Получается, что для каждой пары выборок свой закон распределения.
На самом деле, огромное спасибо за ссылку! У меня была ошибка – статистика не «та-же самая», а совершенно другая. При расчете статистики я считаю дисперсии обеих выборок, а потом умножаю вторую выборку так, чтобы отношение дисперсий оставалось константой (истинные центры обеих выборок равны нулю). Правда, не знаю, несколько это правильно. Можно и не умножать – закон распределения статистики при количестве экспериментов 10 000 не сильно меняется.

[GAA]

Думаю, Вам не следует пытаться получить оценку функции распределения или плотности распределения. Посмотрите файл, на который ссылаются в статье из wiki. Содержимое этого файла типичный пример применения Монте-Карло для тестирования критерия (мне только не понравился выбор значений $\sigma_1$ и $\sigma_2$, я бы зафиксировал $\sigma_1$ и вычислял бы оценку вероятности ошибки первого рода как функцию отношения $\sigma_2/\sigma_1$).

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

Для фиксированных объемов выборки, одинаковых (для простоты равных нулю) значений математических ожиданий и некоторых значений $\sigma_1$ и $\sigma_2$. Выполняется серия экспериментов. В каждом эксперименте генерируется выборка $X$ объема $n_1$, с параметрами $\mu=0$, $\sigma_1$ и выборка $Y$ объема $n_2$, с параметрами, $\mu=0$, $\sigma_1$.
Для этих выборок находят оценки математического ожидания $\overline X$, $\overline Y$ (обозначим через $M_1$ и $M_2$) и несмещенные оценки дисперсий $S_1^2$ и $S_2^2$ (обозначим через $SS_1$ и $SS_2$). Далее вычисляется значение статистики $t$ и $\nu$, здесь $\nu$ округленное до целого значение приведенной Вами статистики. Проверяется условие $t > t_{\nu, \beta}$, где $t_{\nu, \beta}$ — квантиль уровня $\beta$ распределения Стьюдента с $\nu$ степенями свободы. Если условие выполнено, то увеличивается на 1 счетчик случаев, когда гипотеза была отвергнута. Обычно для каждого эксперимента выполняется несколько проверок с различными уровнями. После выполнения всех экспериментов значения счетчиков делят на число экспериментов в серии.
Пример программы на Borland Pascal (компилировалось и выполнялось в среде Delphi 5.0, но использованы средства языка, имеющиеся и в более ранних диалектах Borland Pascal).

 program welch;
uses NormRnd; {Модуль содержит генератор нормального распределения,
{в тексте это функции InitN1DRnd(0, 1) и N1DRnd}

const
 n1 = 5;
 n2 = 5;
 Rep = 1000000;
 sigma1 = 1;
 sigma2 = 16;
var
 quantile : array [1..100] of record l_80, l_90, l_95, l_99 : Extended; end;
 i, j : Integer;
 X: array[1..n1] of Extended;
 Y: array[1..n2] of Extended;
 m1, m2, ss1, ss2, t, nu:    Extended;
 scount_80, scount_90, scount_95, scount_99: Integer;
 wcount_80, wcount_90, wcount_95, wcount_99: Integer;
 txt: Text;

begin
Assign(txt, 'student.dat'); {Из файла 'student.dat' считываются квантили распределения Стюдента}
Reset(txt);                        {уровня 0.80, 0.90, 0.95, 0.99 и степенях свободы 1..100}
for i:= 1 to 100 do with quantile[i] do Readln(txt, l_80, l_90, l_95, l_99);
Close(txt);

scount_80:= 0; scount_90:= 0; scount_95:= 0; scount_99:= 0;
wcount_80:= 0; wcount_90:= 0; wcount_95:= 0; wcount_99:= 0;
InitN1DRnd(0, 1); {Инициализируем для получения стандартных нормальных случайных величин}
for j := 1 to Rep
 do begin
     m1:= 0; for i:= 1 to n1 do begin X[i]:= sigma1*N1DRnd; m1:= m1 + X[i]; end; m1:= m1/n1;
     m2:= 0; for i:= 1 to n2 do begin Y[i]:= sigma2*N1DRnd; m2:= m2 + Y[i]; end; m2:= m2/n2;
     ss1:=0; for i:= 1 to n1 do       ss1:= ss1 + sqr(X[i] - m1);   ss1:= ss1/(n1-1);
     ss2:=0; for i:= 1 to n2 do       ss2:= ss2 + sqr(Y[i] - m2);   ss2:= ss2/(n2-1);
     ss1:= ss1/n1; ss2:= ss2/n2;
     t := (m1-m2)/sqrt(ss1 + ss2);
     nu:= sqr(ss1 + ss2)/(sqr(ss1)/(n1-1) + sqr(ss2)/(n2-1));
     with Quantile[n1+n2-2]
      do begin
          if t > l_80 then inc(scount_80);
          if t > l_90 then inc(scount_90);
          if t > l_95 then inc(scount_95);
          if t > l_99 then inc(scount_99);
         end;
     with Quantile[Round(nu)]
      do begin
          if t > l_80 then inc(wcount_80);
          if t > l_90 then inc(wcount_90);
          if t > l_95 then inc(wcount_95);
          if t > l_99 then inc(wcount_99);
         end;
    end;
if ParamCount > 0 {если указан параметр командной строки}
 then begin       {то считаем, что задано имя файла и выводим в этот файл}
       Assign(txt, ParamStr(1));
       Rewrite(txt);
       Writeln(txt, 's: ', scount_80/Rep:8:3, scount_90/Rep:8:3, scount_99/Rep:8:3, scount_99/Rep:8:3);
       Writeln(txt, 'w: ', wcount_80/Rep:8:3, wcount_90/Rep:8:3, wcount_99/Rep:8:3, wcount_99/Rep:8:3);
       Close(txt);
      end
 else begin
       writeln('s: ', scount_80/Rep:8:3, scount_90/Rep:8:3, scount_95/Rep:8:3, scount_99/Rep:8:3);
       writeln('w: ', wcount_80/Rep:8:3, wcount_90/Rep:8:3, wcount_95/Rep:8:3, wcount_99/Rep:8:3);
      end;
end.

В качестве примера использования такого подхода

Для случая $\sigma_1 = 1$, $\sigma_2 =1$, $n_1=5$, $n_2=5$ получим значения
s: 0.200; 0.100; 0.050; 0.010
w: 0.197; 0.097; 0.047; 0.008

Для случая $\sigma_1 = 1$, $\sigma_2 =16$, $n_1=5$, $n_2=5$ получим значения
s: 0.212; 0.117; 0.068; 0.022
w: 0.200; 0.100; 0.050; 0.010

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

[AndreyL]

У меня немного другая задача - я просто численно проверяю критерии, которые использую в работе. Дело в том, что в современных справочниках иногда, к сожалению, содержатся ошибки. Поскольку ни по образованию, ни по профессии я не математик, то самый простой способ для меня - это численная проверка нескольких вариантов. Вот, благодаря Вам, понял, что пытался неправильно применить поправку Уэлча.
Еще один вопрос: существует ли аналог поправки Уэлча для нескольких выборок? Т.е. критерий равенства средних нескольких выборок при неравных дисперсиях. Что-нибудь похожее на дисперсионный анализ Фишера, но для случая, когда дисперсии выборок не равны. По идее та поправка Уэлча, которую мы обсудили, должна быть частным случаем, т.е. случаем, когда "несколько" равно двум.

[GAA]

В статье Sawilowsky (на которую ссылается статья в wiki) в связи с обобщением на случай многих выборок дается ссылка на работу
Welch, B. L. (1947). The generalization of Student’s problem when several different population variances are involved. Biometrika, 34, 28–35. (Cкан свободно доступен на нескольких ресурсах, например, тут: pdf. Если уберут, можно поискать по названию статьи.)

[AndreyL]

За статьи огромное спасибо! Вот только я не совсем понял формулу 27 в статье Уэлча. Утверждается, что в случае, если $k$ (количество выборок) не равно 2, то статистика $v=\frac{(y-\eta )}{\sqrt{\sum \lambda _i s_i^2}}$ аппроксимируется распределением Стьюдента с соответствующими степенями свободы. Однако не понятно, что в этом случае стоит в числителе. В начале статьи, сразу после уравнения 1, указано, что $ \eta =\alpha _1-\alpha _2$, а $ y=\bar{x}_1-\bar{x}_2$. Но тогда это опять же задача для двух выборок, а не обобщение на несколько выборок, или я чего-то не понимаю.

[GAA]

Выше я привел ссылку на слишком раннюю работу (не был времени посмотреть внимательно, только сегодня к вечеру освобожусь). В связи с Вашим вопросом, думаю, лучше посмотреть работу
Welch, B.L., On the comparison of several mean values: an alternative approach, Biometrika, 38, 330–336, 1951 (свободно доступный pdf).

[AndreyL]

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

[AndreyL]

Попробовал немного иначе. По идее, если взять выборку $Z$ случайной величины, подчиняющейся распределению с законом $F(z,p)$, где $p$ - параметры распределения, и достаточно большого объема $L$, то точки в координатах $\{ x,y \}= \{ rank(F(z_i, p_i))-\frac{1}{2 L}, F(z_i, p_i) \}$ должны лежать на прямой $y=x$, здесь $rank(F(z_i, p_i))$ - порядковый номер элемента $F(z_i, p_i)$ в отсортированной по возрастанию выборке $F(z, p)$. Т.е. рассматривается не сама статистика, а значение функции распределения от этой статистики. Тогда решается вопрос о разных степенях свободы. Правда не уверен, что это правильно.
Но, к сожалению, в области вероятностей больше 0.4 точки отклоняются от прямой вправо, т.е. большие значения статистики, заданной уравнением 27 в статье 1951 года, появляются чаще, чем предлагается распределением Фишера. Надо заметить, что при двух выборках эти точки лежат точно на прямой.