Минаючи нескінченність: t-тест своїми руками

У цьому пості мова піде про реалізацію процедури обчислення значення функції розподілу Стьюдента без використання яких-небудь спеціальних математичних бібліотек. Тільки 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_1D_2— дисперсії вибірок, а N_1N_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. Ймовірність, що значення випадково опиниться за межами цього діапазону становить: , де &\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, можна стверджувати, що t>T. Вірно і зворотне.
Формула самої функції F_n(t)буде розглянуто далі.
На стику нескінченностей
Для опису функції ймовірностей Стьюдента, нам знадобиться визначити дві додаткові функції:
  • гамма-функція Ейлера ;
  • гипергеометрическая функція важливо відзначити, що наведена формула коректна для значень .
Функція розподілу Стьюдента: .
При чисельній реалізації цієї формули для великих значень 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++.

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

Так як функція обчислюється через суму сходящегося ряду, нам необхідно в функцію передавати додатковий аргумент: необхідну кількість доданків.
// гіпергеометричний 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 доданків.
// гіпергеометричний 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++ коді просто міняємо сигнатуру першої функції:
// гіпергеометричний 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);
}

Після компіляції, можна порівняти обчислювані значення табличними. Для цього необхідно викликати функцію зі значенням статистики з таблиці, числом ступенів свободи з першого стовпця, отримане значення порівняти зі значенням у шапці таблиці для двостороннього тесту (попередньо відібравши це значення від одиниці).
І ще, необхідно пам'ятати, що наведена реалізація дає коректні результати для кількості ступенів свободи перевищує значення квадрата оцінюваної статистики (що і відбувається в реальних задачах з великими обсягами даних). Це пов'язано з обмеженням аргументу гипергеометрической функції .
В іншому випадку, можна просто відкидати нульову гіпотезу без розрахунку статистики, або шукати реалізацію, наприклад, з інтегруванням функції щільності розподілу Стьюдента.
Джерело: Хабрахабр

0 коментарів

Тільки зареєстровані та авторизовані користувачі можуть залишати коментарі.