Examples of com.nr.sf.Gamma

com.nr.sf.Gamma

    


    // Test Gamma
    System.out.println("Testing Gamma (gammp, gammq, gser, gcf, invgammp)");


    Gamma gam = new Gamma();


    // Test Gamma.gammp, Gamma.gammq (as well as gser, gcf)
    N=12;
    double a[]={0.5,0.5,0.5,1.0,1.0,1.0,3.0,3.0,3.0,10.0,10.0,10.0}; //,110.0,200.0
    double x[]={0.5,1.0,2.0,0.5,1.0,2.5,2.0,3.5,5.0,7.0,10.0,12.5}; //,100.0,220.0
    double y[]={6.826894921370857e-1,8.427007929497149e-1,9.544997361036406e-1,
      3.934693402873665e-1,6.321205588285578e-1,9.179150013761013e-1,
      3.233235838169363e-1,6.791528011378658e-1,8.753479805169189e-1,
      1.695040627613259e-1,5.420702855281473e-1,7.985688950544638e-1};
    double[] xx=buildVector(x),yy=buildVector(y),zz= new double[N],uu= new double[N],vv= new double[N],c=buildVector(N,1.);


    for (i=0;i<N;i++) {
      zz[i]=gam.gammp(a[i],x[i]);     // Test gammp
      uu[i]=gam.gammq(a[i],x[i]);     // Test gammq
      vv[i]=gam.invgammp(zz[i],a[i]);   // Test invgammp
    }


    System.out.printf("Gamma.gammp: Maximum discrepancy = %f\n", maxel(vecsub(zz,yy)));
    sbeps=5.e-15;
    localflag = maxel(vecsub(zz,yy)) > sbeps;
    globalflag = globalflag || localflag;
    if (localflag) {
      fail("*** Gamma.gammp: Incorrect function values");
      
    }


    System.out.printf("Gamma.gammq: Maximum discrepancy = %f\n", maxel(vecsub(vecadd(zz,uu),c)));
    sbeps=5.e-15;
    localflag = maxel(vecsub(vecadd(zz,uu),c)) > sbeps;
    globalflag = globalflag || localflag;
    if (localflag) {
      fail("*** Gamma.gammq: gammp and gammq do not sum to 1.0");
      
    }


    System.out.printf("Gamma.invgammp: Maximum discrepancy = %f\n", maxel(vecsub(xx,vv)));
    sbeps=5.e-14;
    localflag = maxel(vecsub(xx,vv)) > sbeps;
    globalflag = globalflag || localflag;
    if (localflag) {
      fail("*** Gamma.invgammp: Inverse does not return to original argument x[i]");
      
    }


    xx =new double[2];yy =new double[2];zz =new double[2];uu =new double[2];vv =new double[2];
    double[] aa = new double[2],cc=buildVector(2,1.0);
    aa[0]=110.0; aa[1]=200.0;
    xx[0]=100.0; xx[1]=210.0;
    yy[0]=1.705598979081085e-1; yy[1]=7.639696745011632e-1;
    zz[0]=gam.gammp(aa[0],xx[0]); zz[1]=gam.gammp(aa[1],xx[1]);
    uu[0]=gam.gammq(aa[0],xx[0]); uu[1]=gam.gammq(aa[1],xx[1]);
    vv[0]=gam.invgammp(zz[0],aa[0]); vv[1]=gam.invgammp(zz[1],aa[1]);


    System.out.printf("Gamma.gammp from gammpapprox: Maximum discrepancy = %f\n", maxel(vecsub(zz,yy)));
    sbeps=5.e-14;
    localflag = maxel(vecsub(zz,yy)) > sbeps;
    globalflag = globalflag || localflag;

View Full Code Here

        mychisq += 1.0*SQR(table[i][j]-mynij[i][j])/mynij[i][j];
      }
    }
    mycramerv=sqrt(mychisq/ntotal/(NTYPE-1));
    myconting=sqrt(mychisq/(mychisq+ntotal));
    Gamma gam = new Gamma();
    myprob=gam.gammq(mydf/2.0,mychisq/2.0);


    cntab(table,chisq,df,prob,cramrv,ccc);


  /*
    System.out.printf(fixed << setprecision(4));

View Full Code Here

    chsone(bins, ebins, df, chsq, prob, 1);
  }
  
  public static void chsone(final double[] bins, final double[] ebins,
      final doubleW df, final doubleW chsq, final doubleW prob, final int knstrn) {
    Gamma gam = new Gamma();
    int j,nbins=bins.length;
    double temp;
    df.val=nbins-knstrn;
    chsq.val=0.0;
    for (j=0;j<nbins;j++) {
      if (ebins[j]<0.0 || (ebins[j]==0. && bins[j]>0.))
        throw new IllegalArgumentException("Bad expected number in chsone");
      if (ebins[j]==0.0 && bins[j]==0.0) {
        --df.val;
      } else {
        temp=bins[j]-ebins[j];
        chsq.val += temp*temp/ebins[j];
      }
    }
    prob.val=gam.gammq(0.5*df.val,0.5*chsq.val);
  }

View Full Code Here

    chstwo(bins1, bins2, df, chsq,prob, 1);
  }
  
  public static void chstwo(final double[] bins1, final double[] bins2,
      final doubleW df, final doubleW chsq, final doubleW prob, final int knstrn) {
    Gamma gam = new Gamma();
    int j,nbins=bins1.length;
    double temp;
    df.val=nbins-knstrn;
    chsq.val=0.0;
    for (j=0;j<nbins;j++)
      if (bins1[j] == 0.0 && bins2[j] == 0.0)
        --df.val;
      else {
        temp=bins1[j]-bins2[j];
        chsq.val += temp*temp/(bins1[j]+bins2[j]);
      }
    prob.val=gam.gammq(0.5*df.val,0.5*chsq.val);
  }

View Full Code Here

   association, Cramer's V (cramrv) and the contingency coefficient C (ccc)
   */
  public static void cntab(final int[][] nn, final doubleW chisq,
      final doubleW df, final doubleW prob, final doubleW cramrv, final doubleW ccc) {
    final double TINY=1.0e-30;
    Gamma gam = new Gamma();
    int i,j,nnj,nni,minij,ni=nn.length,nj=nn[0].length;
    double sum=0.0,expctd,temp;
    double[] sumi= new double[ni],sumj = new double[nj];
    nni=ni;
    nnj=nj;
    for (i=0;i<ni;i++) {
      sumi[i]=0.0;
      for (j=0;j<nj;j++) {
        sumi[i] += nn[i][j];
        sum += nn[i][j];
      }
      if (sumi[i] == 0.0) --nni;
    }
    for (j=0;j<nj;j++) {
      sumj[j]=0.0;
      for (i=0;i<ni;i++) sumj[j] += nn[i][j];
      if (sumj[j] == 0.0) --nnj;
    }
    df.val=nni*nnj-nni-nnj+1;
    chisq.val=0.0;
    for (i=0;i<ni;i++) {
      for (j=0;j<nj;j++) {
        expctd=sumj[j]*sumi[i]/sum;
        temp=nn[i][j]-expctd;
        chisq.val += temp*temp/(expctd+TINY);
      }
    }
    prob.val=gam.gammq(0.5*df.val,0.5*chisq.val);
    minij = nni < nnj ? nni-1 : nnj-1;
    cramrv.val=sqrt(chisq.val/(sum*minij));
    ccc.val=sqrt(chisq.val/(chisq.val+sum));
  }

View Full Code Here

    sig = ssig;
    chi2 = 0.;
    q = 1.0;
    sigdat = 0.;
    
    Gamma gam = new Gamma();
    int i;
    double ss=0.,sx=0.,sy=0.,st2=0.,t,wt,sxoss;
    b=0.0;   // accumulate sums ...
    for (i=0;i<ndata;i++) {
      wt=1.0/SQR(sig[i]);   // .. with weights
      ss += wt;
      sx += x[i]*wt;
      sy += y[i]*wt;
    }
    sxoss=sx/ss;
    for (i=0;i<ndata;i++) {
      t=(x[i]-sxoss)/sig[i];
      st2 += t*t;
      b += t*y[i]/sig[i];
    }
    b /= st2;    // Solve for α, b, σ_α, and σ_b
    a=(sy-sx*b)/ss;
    siga=sqrt((1.0+sx*sx/(ss*st2))/ss);
    sigb=sqrt(1.0/st2);     // Calculate χ^2
    for (i=0;i<ndata;i++) chi2 += SQR((y[i]-a-b*x[i])/sig[i]);
    if (ndata>2) q=gam.gammq(0.5*(ndata-2),0.5*chi2);
  }

View Full Code Here

    sy = new double[ndat];
    ww = new double[ndat];
    
    final double POTN=1.571000,BIG=1.0e30,ACC=1.0e-6;
    final double PI=3.141592653589793238;
    Gamma gam = new Gamma();
    Brent brent = new Brent(ACC);
    Chixy chixy = new Chixy();
    int j;
    double amx,amn,scale,bmn,bmx,d1,d2,r2;
    double[] ang = new double[7];
    double[] ch = new double[7];
    doubleW varx = new doubleW(0);
    doubleW vary = new doubleW(0);
    doubleW dum1 = new doubleW(0);
    
    Moment.avevar(x,dum1,varx);
    Moment.avevar(y,dum1,vary);
    scale=sqrt(varx.val/vary.val);
    for (j=0;j<ndat;j++) {
      xx[j]=x[j];
      yy[j]=y[j]*scale;
      sx[j]=sigx[j];
      sy[j]=sigy[j]*scale;
      ww[j]=sqrt(SQR(sx[j])+SQR(sy[j]));
    }
    Fitab fit = new Fitab(xx,yy,ww);
    b = fit.b;
    offs=ang[0]=0.0;
    ang[1]=atan(b);
    ang[3]=0.0;
    ang[4]=ang[1];
    ang[5]=POTN;
    for (j=3;j<6;j++) ch[j]=chixy.get(ang[j]);
    brent.bracket(ang[0],ang[1],chixy);
    ang[0] = brent.ax; ang[1] = brent.bx; ang[2] = brent.cx;
    ch[0]  = brent.fa; ch[1]  = brent.fb; ch[2]  = brent.fc;
    b = brent.minimize(chixy);
    chi2=chixy.get(b);
    a=aa;
    q=gam.gammq(0.5*(ndat-2),chi2*0.5);
    r2=0.0;
    for (j=0;j<ndat;j++) r2 += ww[j];
    r2=1.0/r2;
    bmx=bmn=BIG;
    offs=chi2+1.0;

View Full Code Here

TOP

Related Classes of com.nr.sf.Gamma

com.nr.model.Fitab

com.nr.model.Fitexy

com.nr.stat.Stattests

com.nr.test.test_chapter14.Test_cntab

com.nr.test.test_chapter6.Test_Gamma

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.