Hi!

I want to write a predicate that performs a McNemar test, a non-parametric test for comparing the correctness of two classification models. For the significance level I need a predicate that computes the density D of the chi-square distribution for different values of the test statistic X and degrees of freedom V.

I wrote the predicate below based on this formula for the density function and this code for the gamma function. As far as I can tell the gamma/2 predicate works correctly, but I get density values D that are slightly off.

*Does anybody know where I might find Prolog code for computing the density of the chi-square distribution?*

Aleph has a chi square predicate but if computes the chi square value from density and degrees of freedom. I also know I could use R to do this but I would like to use Prolog only.

Kind regards, JC

```
chi2Density(X, V, D):-
V1 is V / 2,
gamma(V1, Gamma),
D is 2^(-V / 2) * exp(-X / 2) * X^(V / 2 - 1) / Gamma.
% Gamma
% From https://rosettacode.org/wiki/Gamma_function#Prolog
% Seems correct
gamma_coefficients(
[ 1.00000000000000000000000, 0.57721566490153286060651, -0.65587807152025388107701,
-0.04200263503409523552900, 0.16653861138229148950170, -0.04219773455554433674820,
-0.00962197152787697356211, 0.00721894324666309954239, -0.00116516759185906511211,
-0.00021524167411495097281, 0.00012805028238811618615, -0.00002013485478078823865,
-0.00000125049348214267065, 0.00000113302723198169588, -0.00000020563384169776071,
0.00000000611609510448141, 0.00000000500200764446922, -0.00000000118127457048702,
0.00000000010434267116911, 0.00000000000778226343990, -0.00000000000369680561864,
0.00000000000051003702874, -0.00000000000002058326053, -0.00000000000000534812253,
0.00000000000000122677862, -0.00000000000000011812593, 0.00000000000000000118669,
0.00000000000000000141238, -0.00000000000000000022987, 0.00000000000000000001714
]).
tolerance(1e-17).
gamma(X, _):-
X =< 0.0,
!,
fail.
gamma(X, Y):-
X < 1.0,
small_gamma(X, Y),
!.
gamma(1, 1):- !.
gamma(1.0, 1):- !.
gamma(X, Y):-
X1 is X - 1,
gamma(X1, Y1),
Y is X1 * Y1.
small_gamma(X, Y):-
gamma_coefficients(Cs),
recip_gamma(X, 1.0, Cs, 1.0, 0.0, Y0),
Y is 1 / Y0.
recip_gamma(_, _, [], _, Y, Y) :- !.
recip_gamma(_, _, [], X0, X1, Y) :- tolerance(Tol), abs(X1 - X0) < Tol, Y = X1, !. % early exit
recip_gamma(X, PrevPow, [C|Cs], _, X1, Y) :-
Power is PrevPow * X,
X2 is X1 + C*Power,
recip_gamma(X, Power, Cs, X1, X2, Y).
```