Code for chi-square distribution?

Help! Algorithm
1

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).
2

Why did you say that your values are slightly off? For example, by using (in R)

dchisq(8, df = 7)

i get the same result with your code by asking ?- chi2Density(8,7,V).
A small ad: some years ago I’ve written some statistical functions in prolog and published a package: plstat: Statistics using prolog (github.com). Maybe some of them can be of interest. Moreover, if you want to add also your code I’ll be happy to include it!

1 Like
3

@damiazz94 You are right! Checking against R actually gives the correct values.

With my predicate

?- chi2Density(8,7,D).
D = 0.08817913751079277.

And with R

dchisq(8, df = 7)
0.08817914

I was checking against Apples Numbers application with the function CHDIST (CHI2FÖRD in Swedish). With Numbers I get
CHI2FÖRD(8;7) = 0,332593902599308.

That’s odd… I’ll have to do some more investigating

That’s great. I had a look a while ago and checked it again now. The package has a really good list of algorithms. I’ll run some more tests on my chi square predicate and will get back to you.

Kind regards, JC

4

In R, 1 - pchisq(8, 7) yields 0.3325939, that’s your other result (the upper tail of the chi-square distribution). In statistical tests that is often the p value.

The density is not so often used (some Bayesian tests use it).

1 Like
5

@mgondan1 you are right, the p value is relevant instead. Thanks.