На днях возникла одна задача, касающаяся использования критерия Пирсона в экспертных оценках (
для проверки гипотезы согласованности). Вроде бы ничего особенного: производишь расчеты, получаешь значения, а далее сравниваешь их с табличными значениями
так называемых критических точек. Я немного был удивлен, что в Oracle отсутствует реализация функций получения значений как самого распределения, так и обратных значений. Так, например, в Excel данные функции реализованы в виде методов:
CHISQ.DIST.RT(x, deg_freedom) и
CHISQ.INV.RT(probability, deg_freedom).
Рассмотрим подходы для реализации данных методов.
Распределение хи-квадрат характеризуется плотностью вероятности
и функциями распределения (x = χ
2)
Вычисление P(x, v) выполняется разложением в ряд:
Также для реализации необходимо вычисление гамма-функции Г(v/2). В общем виде гамма-функция имеет следующий вид
и обладает следующими свойствами
В случае, если вещественная часть аргумента гамма-функции относится к натуральным числам, тогда гамма-функция может быть записана следующим образом
Также стоит рассмотреть особый случай, когда 0 < z < 1, причем в силу того, что число степеней свободы является целым числом, то можно ограничиться точкой z = 0.5
Ниже приведены примеры расчетов гамма-функции
Реализация функций на C++:
#include <iostream>
#include <math.h>
// specific factorial function
float factorial(float n) {
if (n == 0) { return 1; }
if (n < 0) { return sqrt(M_PI); }
return n * factorial(n - 1);
}
// gamma function
float gamma(float n) {
return factorial(n - 1);
}
// chi square distribution function (v.1)
float chi_square_distribution(float x, float v) {
float a = v / 2, b = x / 2;
// count the Gamma function
float g = gamma(a);
// count of the probability density
float f = pow(x, a - 1) / g / exp(b) / pow(2, a);
// count the sum of series
float s = 0, z = 0;
int k = 0, w = 0;
g = 1;
do {
k = k + 1;
w = w + 2;
z = s;
g = g * (v + w);
s = s + (pow(x, k) / g);
} while (s != z);
// count of probabilities
float p = 2 * x * f * (1 + s) / v;
float q = 1 - p;
return q;
}
// chi square distribution function (v.2)
float _chi_square_distribution(float x, float v) {
float a = v / 2, b = x / 2, g;
// count the Gamma function...
// for even degree
if (int(v) % 2 == 0) {
g = 1;
for (int i = 1; i < a; i++) {
g = g * i;
}
}
// for odd degree
if (int(v) % 2 != 0) {
g = sqrt(M_PI);
for (float i = 0.5; i < a; i++) {
g = g * i;
}
}
// count of the probability density
float f = pow(x, a - 1) / g / exp(b) / pow(2, a);
// count the sum of series
float s = 0, z = 0;
int k = 0, w = 0;
g = 1;
do {
k = k + 1;
w = w + 2;
z = s;
g = g * (v + w);
s = s + (pow(x, k) / g);
} while (s != z);
// count of probabilities
float p = 2 * x * f * (1 + s) / v;
float q = 1 - p;
return q;
}
// recursio search chi square inverse
float search_chi_square(float x, float h, float p, float v) {
float z = x;
x = x + h;
float y = chi_square_distribution(x, v);
if (z == x || y == p) {
return x;
}
if (y < p) {
return search_chi_square(x - h, h / 2, p, v);
}
if (y > p) {
return search_chi_square(x, h, p, v);
}
return 0; // warning control reaches end of non-void function -wreturn-type
}
// chi shuare inverse
float chi_square_inverse(float p, float v) {
return search_chi_square(0, 0.1, p, v);
}
int main(int argc, char **argv) {
std::cout << "chi_square(1, 3) = " << chi_square_distribution(1, 3) << std::endl; // 0.801252
std::cout << "chi_square_inverse(0, 2) = " << chi_square_inverse(0, 2) << std::endl; // 0
std::cout << "chi_square_inverse(0.5, 2) = " << chi_square_inverse(0.5, 2) << std::endl; // 1.38629
std::cout << "chi_square_inverse(0.05, 2) = " << chi_square_inverse(0.05, 2) << std::endl; // 5.99147
std::cout << "chi_square_inverse(0.7, 10) = " << chi_square_inverse(0.7, 10) << std::endl; // 7.26722
std::cout << "chi_square_inverse(0.1, 20) = " << chi_square_inverse(0.1, 20) << std::endl; // 28.4122
return 0;
}
Реализация функций на PL/SQL:
create or replace package stat_pkg is
-- chi square distribution
function chi_square_distribution(x in number, v in number) return number;
-- chi square inverse
function chi_square_inverse(p in number, v in number) return number;
end stat_pkg;
create or replace package body stat_pkg is
/*
* chi square distribution
*/
function chi_square_distribution(x in number, v in number) return number
is
a number;
b number;
g number;
f number;
s number;
z number;
k number;
w number;
p number;
q number;
-- factorial function
function factorial(n in number) return number
is
begin
if n = 0 then return 1; end if;
if n < 0 then return sqrt(asin(1) * 2); end if;
return n * factorial(n - 1);
end;
-- gamma function
function gamma(n in number) return number
is
begin
return factorial(n - 1);
end;
begin
a := v / 2;
b := x / 2;
g := gamma(a);
f := power(x, a - 1) / g / exp(b) / power(2, a);
s := 0; z := 0; k := 0; w := 0;
g := 1;
loop
k := k + 1;
w := w + 2;
z := s;
g := g * (v + w);
s := s + (power(x, k) / g);
exit when s = z;
end loop;
p := 2 * x * f * (1 + s) / v;
q := 1 - p;
return q;
end;
/*
* chi square inverse
*/
function chi_square_inverse(p in number, v in number) return number
is
-- recursio search chi square inverse
function search_chi_square(x in number, h in number, p in number, v in number) return number
is
y number;
begin
y := chi_square_distribution(x + h, v);
if x = (x + h) or y = p then return x; end if;
if y < p then return search_chi_square(x, h / 2, p , v); end if;
if y > p then return search_chi_square(x + h, h, p, v); end if;
end;
begin
return search_chi_square(0, 0.1, p, v);
end;
begin
-- init
null;
end stat_pkg;
