Package com.nr.sf

Examples of com.nr.sf.Chisqdist


    // 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;
View Full Code Here


  class func_Chisqdist implements UniVarRealValueFun{
    double df;
    Chisqdist normi;
    func_Chisqdist(double ddf) {
      df = ddf;
      normi = new Chisqdist(df);
    }
View Full Code Here

TOP

Related Classes of com.nr.sf.Chisqdist

Copyright © 2018 www.massapicom. All rights reserved.
All source code are property of their respective owners. Java is a trademark of Sun Microsystems, Inc and owned by ORACLE Inc. Contact coftware#gmail.com.