C/C++/PLSQL. Pearson distribution (chi-square)

На днях возникла одна задача, касающаяся использования критерия Пирсона в экспертных оценках (для проверки гипотезы согласованности). Вроде бы ничего особенного: производишь расчеты, получаешь значения, а далее сравниваешь их с табличными значениями так называемых критических точек. Я немного был удивлен, что в Oracle отсутствует реализация функций получения значений как самого распределения, так и обратных значений. Так, например, в Excel данные функции реализованы в виде методов:
CHISQ.DIST.RT(x, deg_freedom) и CHISQ.INV.RT(probability, deg_freedom).

Рассмотрим подходы для реализации данных методов.
Распределение хи-квадрат характеризуется плотностью вероятности

и функциями распределения (x = χ²)

 Вычисление 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;