// Test Chisqdist
System.out.println("Testing Chisqdist");
// Test special cases
df=10.0; chisq=1.0;
Chisqdist norm1=new Chisqdist(df);
localflag = abs(norm1.p(chisq)-pow(chisq,df/2.-1.)*exp(-chisq/2.)/pow(2.,df/2.)/factrl((int)(df/2.-1.))) > sbeps;
globalflag = globalflag || localflag;
if (localflag) {
fail("*** Chisqdist: Special case #1 failed");
}
df=10.0; chisq=2.0;
Chisqdist norm2=new Chisqdist(df);
localflag = abs(norm2.p(chisq)-pow(chisq,df/2.-1.)*exp(-chisq/2.)/pow(2.,df/2.)/factrl((int)(df/2.-1.))) > sbeps;
globalflag = globalflag || localflag;
if (localflag) {
fail("*** Chisqdist: Special case #2 failed");
}
df=10.0; chisq=3.0;
Chisqdist norm3=new Chisqdist(df);
localflag = abs(norm3.p(chisq)-pow(chisq,df/2.-1.)*exp(-chisq/2.)/pow(2.,df/2.)/factrl((int)(df/2.-1.))) > sbeps;
globalflag = globalflag || localflag;
if (localflag) {
fail("*** Chisqdist: Special case #3 failed");
}
df=6.0; chisq=1.0;
Chisqdist norm4=new Chisqdist(df);
localflag = abs(norm4.p(chisq)-pow(chisq,df/2.-1.)*exp(-chisq/2.)/pow(2.,df/2.)/factrl((int)(df/2.-1.))) > sbeps;
globalflag = globalflag || localflag;
if (localflag) {
fail("*** Chisqdist: Special case #4 failed");
}
// integral of distribution is one
sbeps=1.e-8;
df=100.;
func_Chisqdist dist = new func_Chisqdist(df);
Midpnt q2=new Midpnt(dist,0.0,4.0);
Midinf q3=new Midinf(dist,4.0,1.0e99);
integral=qromo(q2)+qromo(q3);
localflag = abs(1.0-integral) > sbeps;
// System.out.printf(setprecision(15) << 1.0-integral);
globalflag = globalflag || localflag;
if (localflag) {
fail("*** Chisqdist: Distribution is not normalized to 1.0");
}
// cdf agrees with incomplete integral
sbeps=1.e-7;
df=100.;
func_Chisqdist dist2 = new func_Chisqdist(df);
Chisqdist normcdf=new Chisqdist(df);
localflag=false;
for (i=0;i<N;i++) {
q2=new Midpnt(dist2,0.0,x[i]);
integral=qromo(q2);
c[i]=integral;
d[i]=normcdf.cdf(x[i]);
// System.out.printf(c[i]-d[i]);
localflag = localflag || abs(c[i]-d[i]) > sbeps;
}
globalflag = globalflag || localflag;
if (localflag) {
fail("*** Chisqdist: cdf does not agree with result of quadrature");
}
// inverse cdf agrees with cdf
df=100.;
Chisqdist normc=new Chisqdist(df);
Ran myran = new Ran(17);
sbeps=2.0e-12; // XXX 1.0e-12 not pass.
localflag=false;
for (i=0;i<1000;i++) {
chisq=df-3.0*sqrt(df)+6.0*sqrt(df)*myran.doub();
a=normc.cdf(chisq);
b=normc.invcdf(a);
if (abs(chisq-b) > sbeps) {
System.out.printf("%f %f %f\n",chisq, b ,abs(chisq-b));
}
localflag = localflag || abs(chisq-b) > sbeps;
}
globalflag = globalflag || localflag;
if (localflag) {
fail("*** Chisqdist: Inverse cdf does not accurately invert the cdf");
}
// Fingerprint test
df=100.;
Chisqdist normf=new Chisqdist(df);
for (i=0;i<N;i++) {
p[i]=normf.p(x[i]);
// System.out.printf(setprecision(17) << p[i] << " %f\n", pexp[i]);
}
// System.out.println("Chisqdist: Maximum discrepancy = %f\n", maxel(vecsub(p,pexp)));
localflag = maxel(vecsub(p,pexp)) > sbeps;
globalflag = globalflag || localflag;