Examples of com.nr.sf.KSdist

com.nr.sf.KSdist
Kolmogorov-Smirnov cumulative distribution functions and their inverses. Copyright (C) Numerical Recipes Software 1986-2007 Java translation Copyright (C) Huang Wen Hui 2012 @author hwh

   */
  public static void ksone(final double[] data, final UniVarRealValueFun func,
      final doubleW d, final doubleW prob) {
    int j,n=data.length;
    double dt,en,ff,fn,fo=0.0;
    KSdist ks = new KSdist();
    Arrays.sort(data);
    en=n;
    d.val=0.0;
    for (j=0;j<n;j++) {
      fn=(j+1)/en;
      ff=func.funk(data[j]);
      dt=max(abs(fo-ff),abs(fn-ff));
      if (dt > d.val) d.val=dt;
      fo=fn;
    }
    en=sqrt(en);
    prob.val=ks.qks((en+0.12+0.11/en)*d.val);
  }

View Full Code Here

   */
  public static void kstwo(final double[] data1, final double[] data2,
      final doubleW d, final doubleW prob) {
    int j1=0,j2=0,n1=data1.length,n2=data2.length;
    double d1,d2,dt,en1,en2,en,fn1=0.0,fn2=0.0;
    KSdist ks = new KSdist();
    Arrays.sort(data1);
    Arrays.sort(data2);
    en1=n1;
    en2=n2;
    d.val=0.0;
    while (j1 < n1 && j2 < n2) {
      if ((d1=data1[j1]) <= (d2=data2[j2]))
        do
          fn1=++j1/en1;
        while (j1 < n1 && d1 == data1[j1]);
      if (d2 <= d1)
        do
          fn2=++j2/en2;
        while (j2 < n2 && d2 == data2[j2]);
      if ((dt=abs(fn2-fn1)) > d.val) d.val=dt;
    }
    en=sqrt(en1*en2/(en1+en2));
    prob.val=ks.qks((en+0.12+0.11/en)*d.val);
  }

View Full Code Here

    doubleW ga = new doubleW(0);
    doubleW gb = new doubleW(0);
    doubleW gc = new doubleW(0);
    doubleW gd = new doubleW(0);


    KSdist ks = new KSdist();
    d1.val=0.0;
    for (j=0;j<n1;j++) {
      quadct(x1[j],y1[j],x1,y1,fa,fb,fc,fd);
      quadvl.quadvl(x1[j],y1[j],ga,gb,gc,gd);
      if (fa.val > ga.val) fa.val += 1.0/n1;
      if (fb.val > gb.val) fb.val += 1.0/n1;
      if (fc.val > gc.val) fc.val += 1.0/n1;
      if (fd.val > gd.val) fd.val += 1.0/n1;
      d1.val=max(d1.val,abs(fa.val-ga.val));
      d1.val=max(d1.val,abs(fb.val-gb.val));
      d1.val=max(d1.val,abs(fc.val-gc.val));
      d1.val=max(d1.val,abs(fd.val-gd.val));
    }
    pearsn(x1,y1,r1,dum,dumm);
    sqen=sqrt(n1);
    rr=sqrt(1.0-r1.val*r1.val);
    prob.val=ks.qks(d1.val*sqen/(1.0+rr*(0.25-0.75/sqen)));
  }

View Full Code Here

    doubleW gb = new doubleW(0);
    doubleW gc = new doubleW(0);
    doubleW gd = new doubleW(0);


    double d1,d2,rr,sqen;
    KSdist ks = new KSdist();
    d1=0.0;
    for (j=0;j<n1;j++) {
      quadct(x1[j],y1[j],x1,y1,fa,fb,fc,fd);
      quadct(x1[j],y1[j],x2,y2,ga,gb,gc,gd);
      if (fa.val > ga.val) fa.val += 1.0/n1;
      if (fb.val > gb.val) fb.val += 1.0/n1;
      if (fc.val > gc.val) fc.val += 1.0/n1;
      if (fd.val > gd.val) fd.val += 1.0/n1;
      d1=max(d1,abs(fa.val-ga.val));
      d1=max(d1,abs(fb.val-gb.val));
      d1=max(d1,abs(fc.val-gc.val));
      d1=max(d1,abs(fd.val-gd.val));
    }
    d2=0.0;
    for (j=0;j<n2;j++) {
      quadct(x2[j],y2[j],x1,y1,fa,fb,fc,fd);
      quadct(x2[j],y2[j],x2,y2,ga,gb,gc,gd);
      if (ga.val > fa.val) ga.val += 1.0/n1;
      if (gb.val > fb.val) gb.val += 1.0/n1;
      if (gc.val > fc.val) gc.val += 1.0/n1;
      if (gd.val > fd.val) gd.val += 1.0/n1;
      d2=max(d2,abs(fa.val-ga.val));
      d2=max(d2,abs(fb.val-gb.val));
      d2=max(d2,abs(fc.val-gc.val));
      d2=max(d2,abs(fd.val-gd.val));
    }
    d.val=0.5*(d1+d2);
    sqen=sqrt(n1*n2/(n1+n2));
    pearsn(x1,y1,r1,dum,dumm);
    pearsn(x2,y2,r2,dum,dumm);
    rr=sqrt(1.0-0.5*(r1.val*r1.val+r2.val*r2.val));
    prob.val=ks.qks(d.val*sqen/(1.0+rr*(0.25-0.75/sqen)));
  }

View Full Code Here

    // Test KSdist
    System.out.println("Testing KSdist");


    // Test special cases
    z=0.0;
    KSdist norm = new KSdist();
    localflag = abs(norm.pks(z)-0.0) > sbeps;
    globalflag = globalflag || localflag;
    if (localflag) {
      fail("*** KSdist: Special case #1 failed for Pks");
      
    }


    localflag = abs(norm.qks(z)-1.0) > sbeps;
    globalflag = globalflag || localflag;
    if (localflag) {
      fail("*** KSdist: Special case #1 failed for Qks");
      
    }
    
    z=1.0;
//    System.out.printf(setprecision(17) << norm.pks(z) << " %f\n", 1-2.0*exp(-2.0)+2.0*exp(-8.0)-2.0*exp(-18.0));
    localflag = abs(norm.pks(z)-1+2.0*exp(-2.0)-2.0*exp(-8.0)+2.0*exp(-18.0)) > sbeps;
    globalflag = globalflag || localflag;
    if (localflag) {
      fail("*** KSdist: Special case #2 failed for Pks");
      
    }


    localflag = abs(norm.qks(z)-2.0*exp(-2.0)+2.0*exp(-8.0)-2.0*exp(-18.0)) > sbeps;
    globalflag = globalflag || localflag;
    if (localflag) {
      fail("*** KSdist: Special case #2 failed for Qks");
      
    }


    z=100.0;
    localflag = abs(norm.pks(z)-1.0) > sbeps;
    globalflag = globalflag || localflag;
    if (localflag) {
      fail("*** KSdist: Special case #3 failed for Pks");
      
    }


    localflag = abs(norm.qks(z)-0.0) > sbeps;
    globalflag = globalflag || localflag;
    if (localflag) {
      fail("*** KSdist: Special case #3 failed for Qks");
      
    }


    // Pks correctly related to Qks
    localflag=false;
    for (i=0;i<N;i++)
      localflag = localflag || 1.0-norm.pks(x[i])-norm.qks(x[i]) > sbeps;
    globalflag = globalflag || localflag;
    if (localflag) {
      fail("*** KSdist: Pks and Qks do not add to 1.0 in all cases");
      
    }


    // inverse of Pks agrees with Pks
    Ran myran = new Ran(17);
    localflag=false;
    sbeps=6.e-10;  // XXX 5.e-10 not pass
    for (i=0;i<10;i++) {
      z=0.3+3.0*myran.doub();
      a=norm.pks(z);
      b=norm.invpks(a);
//      if (abs(z-b) > sbeps) {
//        System.out.printf(setprecision(15) << z << " %f\n", b << " %f\n", abs(z-b));
//      }
      localflag = localflag || abs(z-b) > sbeps;
    }
    globalflag = globalflag || localflag;
    if (localflag) {
      fail("*** KSdist: invPks does not accurately invert the Pks");
      
    }


    // inverse of Qks agrees with Qks
    localflag=false;
    sbeps=5.e-10;
    for (i=0;i<1000;i++) {
      z=0.3+3.0*myran.doub();
      a=norm.qks(z);
      b=norm.invqks(a);
//      if (abs(z-b) > sbeps) {
//        System.out.printf(setprecision(15) << z << " %f\n", b << " %f\n", abs(z-b));
//      }
      localflag = localflag || abs(z-b) > sbeps;
    }
    globalflag = globalflag || localflag;
    if (localflag) {
      fail("*** KSdist: invQks does not accurately invert the Qks");
      
    }


    // Fingerprint test for Pks
    sbeps=1.e-15;
    for (i=0;i<N;i++) {
      pks[i]=norm.pks(x[i]);
//      System.out.printf(setprecision(17) << pks[i]);
    }
//    System.out.println("KSdist: Maximum discrepancy = %f\n", maxel(vecsub(pks,pksexp)));
    localflag = maxel(vecsub(pks,pksexp)) > sbeps;
    globalflag = globalflag || localflag;
    if (localflag) {
      fail("*** KSdist: Fingerprint does not match expectations");
      
    }


    // Fingerprint test for Qks
    sbeps=1.e-15;
    for (i=0;i<N;i++) {
      qks[i]=norm.qks(x[i]);
//      System.out.printf(setprecision(17) << qks[i]);
    }
//    System.out.println("KSdist: Maximum discrepancy = %f\n", maxel(vecsub(qks,qksexp)));
    localflag = maxel(vecsub(qks,qksexp)) > sbeps;
    globalflag = globalflag || localflag;

View Full Code Here

TOP

Related Classes of com.nr.sf.KSdist

com.nr.stat.Stattests

com.nr.test.test_chapter6.Test_KSdist

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.