В этом посте речь пойдёт о реализации процедуры вычисления значения функции распределения Стьюдента без использования каких-либо специальных математических библиотек. Только Java (либо C/C++, код вполне универсален).
T-тест это...
Для старта, вспомним немного теории: t-тест, он же тест Стьдента, используется для проверки статистических гипотез о равенстве математических ожиданий двух выборок, либо о равенстве некоторому значению математического ожидания одной выборки.
T-тест бывает:
- одновыборочным или двухвыборочным (смотрите выше);
- односторонним или двухсторонним (оцениваемая величина строго положительна или может иметь отрицательные значения соотвественно, смотрите ниже );
- точным или ассимптотическим (оцениваемые выборки имеют нормальное распределение или лишь стремятся к нему, соответсвенно).
Для проведения теста нам необходимо вычислить значение t-статистики.
Для двух выборок:
Для одной выборки:
Рассчитанная величина t имеет распределение Стьюдента, если исходные данные распределены по Гауссу, или стремится к распределению Стьюдента в случае негауссовых выборок.
Распределение вероятностей описывается функцией 
Как и любое распределение вероятностей, 
Использование распределения рассмотрим на примере двусторонней оценки математического ожидания для одной выборки: предположим, что реальное значение математического ожидания равно m, тогда, рассчитанная статистика t должна иметь значение из окрестности 0. То есть, существует некоторое значение T, такое что: 
В случае одностороннего теста, рассматривается только один из «хвостов» распределения и 
Что касается значения T, то оно может быть взято из таблицы (например, здесь) или вычислено как значение обратной функции к 
Во многих случаях можно поступить проще: воспользоваться свойством монотонности функции распределения. Вычислив значение вероятности 
Формула самой функции 
На стыке бесконечноcтей
Для описания функции вероятностей Стьюдента, нам понадобится определить две дополнительные функции:
- гамма-функция Эйлера
;
- гипергеометрическая функция
, важно отметить, что приведённая формула корректна для значений
.
Функция распределения Стьюдента: 
При численной реализации этой формулы для больших значений n возникает проблема с множителем 
Для разрешения этой проблемы, воспользуемся следующим свойством гамма-функции:
С учётом этого свойства, можно переписать отношение гамма-функций следующим образом:
Предположим, что есть некоторое значение Q, такое, что отношение гамма-функций корректно вычисляется напрямую, если аргументы числителя и знаменателя не превышают Q (предполагаем, что мы работаем только с положительными действительными аргументами). Тогда рассмотрим функцию 
С учётом введённых определений, можно переписать функцию вероятностей:
и переходить к реализации.
Реализация
Все рассматриваемые далее процедуры у меня реализованы как статические метод специализированного класса в 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);
}После компиляции, можно сравнить вычисляемые значения с табличными. Для этого необходимо вызвать функцию со значением статистики из таблицы, числом степеней свободы из первого столбца, полученное значение сравнить со значением в шапке таблицы для двухстороннего теста (предварительно отняв это значение от единицы).
И ещё, необходимо помнить, что приведённая реализация даёт корректные результаты для количества степеней свободы превышающего значение квадрата оцениваемой статистики (что и происходит в реальных задачах с большими объемами данных). Это связано с ограничением аргумента гипергеометрической функции 
В противном случае, можно просто отвергать нулевую гипотезу без рассчета статистики, либо искать реализацию, например, с интегрированием функции плотности распределения Стьюдента.
Автор: JamaGava
