## Wednesday, June 5, 2013

### How to calculate PValue from ChiSquare in Java

Several days ago I faced a problem: I needed to calculate pValue for a given ChiSquare with a given Degrees of Freedom. I didn't have time to do the research and to implement the algorithm from the scratch. In fact I had just a couple of hours to solve the problem. So I started searching...

I was unable to find PValue calculator written in Java but was able to find one written in JavaScript there. So I converted the code into Java and want to share it to you (just call pochisq(chiSquareValue, degreesOfFreedom)):

```public class ChiSquareUtils {
private static final double LOG_SQRT_PI = log(sqrt(PI));
private static final double I_SQRT_PI = 1 / sqrt(PI);
public static final int MAX_X = 20; // max value to represent exp(x)

/* POCHISQ -- probability of chi-square value
Hill, I. D. and Pike, M. C. Algorithm 299
Collected Algorithms for the CACM 1967 p. 243
Updated for rounding errors based on remark in
ACM TOMS June 1985, page 185
*/
public static double pochisq(double x, int df) {
double a, s;
double e, c, z;

if (x <= 0.0 || df < 1) {
return 1.0;
}
a = 0.5 * x;
boolean even = (df & 1) == 0;
double y = 0;
if (df > 1) {
y = ex(-a);
}
s = (even ? y : (2.0 * poz(-sqrt(x))));
if (df > 2) {
x = 0.5 * (df - 1.0);
z = (even ? 1.0 : 0.5);
if (a > MAX_X) {
e = (even ? 0.0 : LOG_SQRT_PI);
c = log(a);
while (z <= x) {
e = log(z) + e;
s += ex(c * z - a - e);
z += 1.0;
}
return s;
} else {
e = (even ? 1.0 : (I_SQRT_PI / sqrt(a)));
c = 0.0;
while (z <= x) {
e = e * (a / z);
c = c + e;
z += 1.0;
}
return c * y + s;
}
} else {
return s;
}
}

public static double poz(double z) {
double y, x, w;
double Z_MAX = 6.0; // Maximum meaningful z value
if (z == 0.0) {
x = 0.0;
} else {
y = 0.5 * Math.abs(z);
if (y >= (Z_MAX * 0.5)) {
x = 1.0;
} else if (y < 1.0) {
w = y * y;
x = ((((((((0.000124818987 * w
- 0.001075204047) * w + 0.005198775019) * w
- 0.019198292004) * w + 0.059054035642) * w
- 0.151968751364) * w + 0.319152932694) * w
- 0.531923007300) * w + 0.797884560593) * y * 2.0;
} else {
y -= 2.0;
x = (((((((((((((-0.000045255659 * y
+ 0.000152529290) * y - 0.000019538132) * y
- 0.000676904986) * y + 0.001390604284) * y
- 0.000794620820) * y - 0.002034254874) * y
+ 0.006549791214) * y - 0.010557625006) * y
+ 0.011630447319) * y - 0.009279453341) * y
+ 0.005353579108) * y - 0.002141268741) * y
+ 0.000535310849) * y + 0.999936657524;
}
}
return z > 0.0 ? ((x + 1.0) * 0.5) : ((1.0 - x) * 0.5);
}

public static double ex(double x) {
return (x < -MAX_X) ? 0.0 : Math.exp(x);
}
}
```
And here is a test for this PValue from ChiSquare implementation:
```
/*
Test values are obtained from
- http://www.fourmilab.ch/rpkp/experiments/analysis/chiCalc.html
- http://www.danielsoper.com/statcalc3/calc.aspx?id=11
*/
@Test
public void testPiValueOfChi2() {
assertThat(getFirstDigits(pochisq(2.314, 1)), equalTo(1282));
assertThat(getFirstDigits(pochisq(2.314, 2)), equalTo(3144));
assertThat(getFirstDigits(pochisq(17, 1)), equalTo(0));
assertThat(getFirstDigits(pochisq(1, 1)), equalTo(3173));
assertThat(getFirstDigits(pochisq(3, 2)), equalTo(2231));
assertThat(getFirstDigits(pochisq(12, 1)), equalTo(5));
}

```