Минуя бесконечность: t-тест своими руками16.08.2016 00:33

В этом посте речь пойдёт о реализации процедуры вычисления значения функции распределения Стьюдента без использования каких-либо специальных математических библиотек. Только Java (либо C/C++, код вполне универсален).

T-тест это…

Для старта, вспомним немного теории: t-тест, он же тест Стьдента, используется для проверки статистических гипотез о равенстве математических ожиданий двух выборок, либо о равенстве некоторому значению математического ожидания одной выборки.

T-тест бывает:

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

Для проведения теста нам необходимо вычислить значение t-статистики.

Для двух выборок:
$t = \frac{\overline X_1 - \overline X_2}{\sqrt{\frac{D_1}{n_1}+\frac{D_2}{n_2}}}$ , где $\overline X_1$ и $\overline X_2$ — расчётные математические ожидания выборок, $D_1$ и $D_2$ — дисперсии выборок, а $N_1$ и $N_2$ — количество элементов в первой и второй выборке соответственно.

Для одной выборки:
$t = \frac{\overline X - m}{\sqrt{D/N}}$ , где $\overline X$ и D- рассчитанные математическое ожидание и дисперсия выборки, N — количество элементов, m — величина, равенство которой математического ожидания проверяется в тесте.

Рассчитанная величина t имеет распределение Стьюдента, если исходные данные распределены по Гауссу, или стремится к распределению Стьюдента в случае негауссовых выборок.

Распределение вероятностей описывается функцией $F_n(t)={P(x<t)}$ , где n — количество степеней свободы: $n =N-1$ в случае одной выборки и $n = N_1 + N_2 - 2$ в случае двух выборок.

Как и любое распределение вероятностей, $F_n$ монотонно стремится к 0 или к 1 при аргументе стремящемся к отрицательной или положительной бесконечности соответственно.

Использование распределения рассмотрим на примере двусторонней оценки математического ожидания для одной выборки: предположим, что реальное значение математического ожидания равно m, тогда, рассчитанная статистика t должна иметь значение из окрестности 0. То есть, существует некоторое значение T, такое что: $P_0(-T<t<T)->1» />. Вероятность, что значение случайно окажется за пределами этого диапазона составляет: <img src=$ , где $&\alpha&$ — величина, называемая уровнем значимости. Если рассчитанное значение статистики t находится за пределами интервала (-T; T), то с вероятностью $P_1$ можно утверждать, что $\overline X$ отлично от m.

В случае одностороннего теста, рассматривается только один из «хвостов» распределения и $P_1=1-&\alpha&$ .

Что касается значения T, то оно может быть взято из таблицы (например, здесь) или вычислено как значение обратной функции к $F_n$ для требуемого $P_1$ . Табличные значения есть не для всех возможных уровней значимости и количеств степеней свободы, а вычисление обратной функции весьма трудоёмко.

Во многих случаях можно поступить проще: воспользоваться свойством монотонности функции распределения. Вычислив значение вероятности $P_t=F_n(t)$ , можно сравнить его с $P_1$ . Если окажется, что $P_t>P_1» />, можно утверждать, что <img src=$ будет рассмотрена далее.

На стыке бесконечноcтей

Для описания функции вероятностей Стьюдента, нам понадобится определить две дополнительные функции:

гамма-функция Эйлера $%5CGamma(z)%3D%5Cint%5Climits_0%5E1%20(-$ ;
гипергеометрическая функция $H(a%2Cb%3Bc%3Bz)%20%3D%201%2B%20%5Csum%5$ , важно отметить, что приведённая формула корректна для значений $%5Cleft%7Cz%5Cright%7C%3C1$ .

Функция распределения Стьюдента: $F_n(t)%3D%5Cfrac%7B1%7D%7B2%7D%2B%5Cfrac$ .

При численной реализации этой формулы для больших значений n возникает проблема с множителем $%5Cfrac%7B%5CGamma%5Cleft(%5Cfrac%7Bn%2B$ . Не смотря на то, что само отношение имеет достаточно небольшие значения, числитель и знаменатель по отдельности могут принимать значения, для которых не хватит разрядности типа с плавающей точкой. Численно, мы получаем неопределённость вида $%5Cfrac%7B%5Cinfty%7D%7B%5Cinfty%7D$ .

Для разрешения этой проблемы, воспользуемся следующим свойством гамма-функции:

$%5CGamma(x)%3D%5Cfrac%7B%5CGamma%5Cleft($

С учётом этого свойства, можно переписать отношение гамма-функций следующим образом:

$%5Cfrac%7B%5CGamma(x)%7D%7B%5CGamma(y)%7$

Предположим, что есть некоторое значение Q, такое, что отношение гамма-функций корректно вычисляется напрямую, если аргументы числителя и знаменателя не превышают Q (предполагаем, что мы работаем только с положительными действительными аргументами). Тогда рассмотрим функцию $R%5Cleft(x%2C%20y%5Cright)%3D%5Cfrac%7B%$ , которая вычисляется рекурсивно:

$R%5Cleft(%20x%2C%20y%20%5Cright)%3D%5Cle$

С учётом введённых определений, можно переписать функцию вероятностей:

$F_n(t)%3D%5Cfrac%7B1%7D%7B2%7D%2B%5Cfrac$

и переходить к реализации.

Реализация

Все рассматриваемые далее процедуры у меня реализованы как статические метод специализированного класса в Java, чтобы использовать их без создания экземпляров. Для реализации в C/C++ спецификатор «static» не актуален (хотя и не криминален).

Гамма-функция Эйлера

Начнём с гамма-функции Эйлера: на её основе рассчитываются многие статистические величины. Предполагаем, что гамма-функция будет вычисляться для действительных положительных аргументов.

Формулу для численного расчёта можно подсмотреть, например, здесь. Программный код на Java выглядит следующим образом:

    // Euler's gamma-function
static double gamma(double x) {       
    double tmp = (x - 0.5) * Math.log(x + 4.5) - (x + 4.5);
    double ser = 1.0 + 
                    76.18009173   / (x + 0.0) -  86.50532033   / (x + 1.0) +
                    24.01409822   / (x + 2.0) -  1.231739516   / (x + 3.0) +
                    0.00120858003 / (x + 4.0) -  0.00000536382 / (x + 5.0);    
     return Math.exp(tmp + Math.log(ser * Math.sqrt(2 * Math.PI)));
}

Пару лёгких движений по удалению всех «Math.» и мы получим код на C/C++.

Отношение гамма-функций

Функция принимает на вход два аргумента: первый — аргумент числителя, второй — знаменателя. Если максимальный из аргументов меньше либо равен 100, то вычисляем отношение напрямую. Иначе, вычисляем рекурсивно по частям, деля аргументы надвое. Число 100 выбрано эмпирически.

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

       // gamma(x) / gamma(y)
static double gammaRatio(double x, double y) {          
    double m = Math.abs(Math.max(x, y));          
    if (m <= 100.0)
        return gamma(x) / gamma(y);          
    return Math.pow(2, x - y)  *  
                   gammaRatio(x * 0.5, y * 0.5)  * 
                   gammaRatio(x * 0.5 + 0.5, y * 0.5 + 0.5);
}

Поиск максимального из аргументов выполняется для большей универсальности функции. При вычислении функции распределения Стьюдента, аргумент числителя всегда больше аргумента знаменателя. Без «Math.» мы получим код на C/C++.

Гипергеометрическая функция

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

    // hypergeometric function
static double hyperGeom(double a, double b, double c, double z, int deep) {
    double S = 1.0, M, d;    
    for (int i = 1; i <= deep; i++) {     
        M = Math.pow(z, (double)i);
            for (int j = 0; j < i; j++) {
                    d = (double)j;
                    M *= (a + d) * (b + d) / ( (1.0 + d) * (c + d) );
            }
            S += M;           
    }      
    return S;       
}

По умолчанию будем использовать 20 слагаемых.

    // hypergeometric function with default deep value= 20
static double hyperGeom(double a, double b, double c, double z) {
    return hyperGeom(a, b, c, z, 20);
 }

При использовании в C/C++ коде просто меняем сигнатуру первой функции:

    // hypergeometric function
double hyperGeom(double a, double b, double c, double z, int deep = 20) {       
    ...

Вторая функция уже не нужна.

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

Функция распределения Стьюдента

Наконец, собираем всё вместе:

    // Student's distribution
static double tDist(double x, int n) {
    double dN = (double)n;
    return 0.5 + x * gammaRatio(0.5 * (dN + 1.0), 0.5 * dN) * 
               hyperGeom(0.5, 0.5 * (dN + 1.0), 1.5, -x*x / dN) / Math.sqrt(Math.PI * dN);
}

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

И ещё, необходимо помнить, что приведённая реализация даёт корректные результаты для количества степеней свободы превышающего значение квадрата оцениваемой статистики (что и происходит в реальных задачах с большими объемами данных). Это связано с ограничением аргумента гипергеометрической функции $%5Cleft%7Cz%5Cright%7C%3C1$ .

В противном случае, можно просто отвергать нулевую гипотезу без рассчета статистики, либо искать реализацию, например, с интегрированием функции плотности распределения Стьюдента.

Комментарии (4)

Psychopompe

15 августа 2016 в 18:44

+1

↑

↓

Отличная штука. Планируется продолжение по матфункциям?
- JamaGava
  
  15 августа 2016 в 18:59
  
  0
  
  ↑
  
  ↓
  
  Спасибо. Да планирую продолжать по статистике.
Zulfat

15 августа 2016 в 18:58

+1

↑

↓

Мне кажется у вас небольшая опечатка в формуле, в которой вы расписываете функцию распределения Стъюдента через гамма-функцию и гипергеометрическую функцию, там вместо $F{n}(t)$ должно быть $F{n}(x)$. Ну или я что-то недопонял.
- JamaGava
  
  15 августа 2016 в 18:59
  
  0
  
  ↑
  
  ↓
  
  Спасибо большое. Уже исправил.