package stats;

/**
*
*This class contains Java translations of FORTRAN
*routines written by Don Amos and S.L. Daniel to calculate the
*noncentral T cumulative distribution function.
*
*@version .5 --- February 17, 2002
*
*/

public class CDF_nct_Amos extends Object implements Gaus8_fcn, Fzero_methods {

      double h2o2,p_inv,delta_inv,dof_inv;

/**
*
*This method calculates the noncentral T cumulative distribution function.
*It is a Java translation of the public domain FORTRAN
*routine fnct written by Don Amos and S.L. Daniel.
*Steve Verrill translated the FORTRAN code (a write date of September, 1972,
*and a modification date of July, 1991)
*into Java.  This translation was performed on February 17, 2001.
*
*@param   x     The method returns the value of the noncentral T
*               cumulative distribution function at x.
*@param   delta The noncentrality parameter
*@param   dof   The degrees of freedom
*
*@version .5 --- February 17, 2002
*
*/


   public static double fnct(double x, double delta, double dof) {

/*

Here is a copy of the documentation in the FORTRAN code:

      FUNCTION FNCT(X,DELTA,F)
C 
C     WRITTEN BY D.E. AMOS AND S.L. DANIEL, SEPTEMBER, 1972.
C     MODIFIED, JULY, 1991
C 
C     REFERENCES
C         A SURVEY OF PROPERTIES AND APPLICATIONS OF THE NON-CENTRAL T
C         DISTRIBUTION BY D.B. OWEN, TECHNOMETRICS, VOL. 10, NO. 3, 
C         AUGUST, 1968, PP. 464-465.
C 
C         SLA-73-0333 
C 
C     ABSTRACT
C         FNCT COMPUTES THE CUMULATIVE NON-CENTRAL T-DISTRIBUTION FOR 
C         T=U/SQRT(V/F) WHERE U IS N(DELTA,1) AND V IS CHI-SQUARE 
C         WITH F DEGREES OF FREEDOM. FINITE SUMS WHOSE FIRST
C         TERMS CONTAIN THE CUMULATIVE NORMAL AND/OR THE T(H,A) 
C         INTEGRAL OF OWEN ARE EVALUATED FOR EVEN AND ODD F.
C 
C     DESCRIPTION OF ARGUMENTS
C 
C         INPUT 
C           X      - T-DISTRIBUTION VARIABLE,  UNRESTRICTED
C           DELTA  - NON-CENTRALITY PARAMETER, UNRESTRICTED 
C           F      - DEGREES OF FREEDOM,       F.GE.1.0E0
C 
C         OUTPUT
C           FNCT   - CUMULATIVE DISTRIBUTION AT X, FNCT=PR(T.LE.X)
C 
C     ERROR CONDITIONS
C         IMPROPER INPUT PARAMETER IS A FATAL ERROR 
C 
C     FNCT   USES COMMON BLOCKS MLTHA
c
c                               x    x   x    x      x      x      x
c   need to compile with xerror.f and machcon.f
c
C     FNCT   USES SUBROUTINES FNORM,THA,FTHA,RERF,I1MACH,R1MACH,XERROR
C-----------------------------------------------------------------------

*/

      final double srtp = 2.50662827463100E0;
      final double osrtp = 3.98942280401432E-1;

      double rel[] = new double[1];

      double t2[] = new double[1];

      double p,tol,a,b,srb,x1,xn,t1,asrb,x2,x3,g2,
             gprime,cm0,cm1,cm2,fk,ak,ss,cmk,y1,arg;

      int nz[] = new int[1];

      int ndof,kode,n,j,nt;

      if (dof < 1.0) {

         p = -1.0;

         System.out.print("\nError in call to fnct:" +
         " The degrees of freedom" +
         " was less than 1.\n\n");

         return p;

      }

/*

From the FORTRAN code:

C-----------------------------------------------------------------------
C     SET PARAMETERS RELATED TO MACHINE CONSTANTS.
C     TOL IS THE MAXIMUM OF UNIT ROUNDOFF, 1.0E-15
C     ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT.
C     K1 = I1MACH(12)
C     K2 = I1MACH(13)
C     K = MIN0(IABS(K1),IABS(K2))
C     ELIM = 2.303E0*(FLOAT(K)*R1MACH(5)-3.0E0)
C-----------------------------------------------------------------------

In the FORTRAN code "unit roundoff" is taken to be
SLATEC's r1mach(4).  It appears that for IEEE 754
double precision aritmetic this should be replaced
by 2^{-53} x 2 = 2.22e-16.  The maximum of this and
1.0e-15 is 1.0e-15, so:

*/

      tol = 1.0e-15;
      rel[0] = tol*10.0;

/*

SLATEC's i1mach(12) = the smallest exponent of single precision reals
         i1mach(13) = the largest exponent of single precision reals
         r1mach(5) = log, base 10, of 2 = .3010299957

For IEEE 754 doubles, the smallest exponent is -1022.
The largest exponent is 1023.  So:

      elim = 2.303*(1022*.3010299957 - 3.0);

*/

      a = x/Math.sqrt(dof);
      b = dof/(dof + x*x);
      srb = Math.sqrt(b);
      ndof = (int)dof;

      if (ndof % 2 == 1) {

//   dof is odd

         x1 = delta*srb;
         xn = -x1;
         kode = 1;

         t1 = CDF_Normal_Amos.fnorm(xn,kode,nz);

         CDF_nct_Amos.tha(x1,a,rel,t2);        

         if (ndof == 1) {

            p = t1 + t2[0] + t2[0];

            return p;

         }

// 102

         asrb = a*srb;
         x2 = x1*a;
         x3 = delta*a;

         kode = 1;

         g2 = CDF_Normal_Amos.fnorm(x2,kode,nz);

         gprime = osrtp*Math.exp(-x1*x1/2.0);

         cm0 = asrb*gprime*g2;

         gprime = osrtp*Math.exp(-delta*delta/2.0);

         cm1 = b*(x3*cm0 + (a/srtp)*gprime);

         cm2 = 0.5*b*(x3*cm1 + cm0);

         fk = 3.0;

         ak = 1.0;

         n = (ndof - 3)/2;

         if (n - 1 < 0) {

// 103

            p = t1 + t2[0] + t2[0] + cm1 + cm1;

            return p;

         }

         if (n - 1 == 0) {

// 104

            x1 = b/3.0;

            x2 = cm1 + (x1 + x1)*(x3*cm2 + cm1);

            p = t1 + t2[0] + t2[0] + x2 + x2;

            return p;

         }

//   n - 1 > 0
//   105

         ss = cm1;

         for (j = 1; j <= n; j++) {

            cmk = ((fk - 1.0)/fk)*b*(ak*x3*cm2 + cm1);
            ss += cmk;
            fk++;

            ak = 1.0/((fk - 2.0)*ak);

//   update and skip even values

            cm1 = cm2;
            cm2 = cmk;
            cmk = ((fk - 1.0)/fk)*b*(ak*x3*cm2 + cm1);
            fk++;

            ak = 1.0/((fk - 2.0)*ak);

            cm1 = cm2;
            cm2 = cmk;

         }

         p = t1 + t2[0] + t2[0] + ss + ss;    

         return p;

      } else {

//   dof is even

         y1 = -delta;

         kode = 1;

         t1 = CDF_Normal_Amos.fnorm(y1,kode,nz);

         x1 = a*srb;
         x2 = delta*x1;
         x3 = delta*a;

         kode = 1;

         g2 = CDF_Normal_Amos.fnorm(x2,kode,nz);

         arg = -delta*delta*srb*srb/2.0;

         gprime = osrtp*Math.exp(arg);

         cm0 = x1*gprime*g2;

         ndof = (int)dof;

         n = ndof/2;

         if (n == 1) {

            p = t1 + srtp*cm0;

            return p;

         }

         gprime = osrtp*Math.exp(-delta*delta/2.0);

         cm1 = b*(x3*cm0 + (a/srtp)*gprime);

         cm2 = (b/2.0)*(x3*cm1 + cm0);

         if (n == 2) {

            p = t1 + srtp*(cm0 + cm2);

            return p;

         }

         fk = 3.0;
         ak = 1.0;
         ss = cm0 + cm2;
         nt = n - 2;

         for (j = 1; j <= nt; j++) {

            cmk = ((fk - 1.0)/fk)*b*(ak*x3*cm2 + cm1);

            fk++;

            ak = 1.0/((fk - 2.0)*ak);

            cm1 = cm2;

            cm2 = cmk;

            cmk = ((fk - 1.0)/fk)*b*(ak*x3*cm2 + cm1);

            ss += cmk;

            fk++;

            ak = 1.0/((fk - 2.0)*ak);

            cm1 = cm2;

            cm2 = cmk;

         }

         p = t1 + srtp*ss;

         return p;        

      }

   }








/**
*
*This method calculates the inverse of the
*noncentral T cumulative distribution function.
*It uses the fzero zero-finding routine and the
*fnct noncentral T cumulative distribution function.
*It was written by a statistician rather than
*a numerical analyst.  When a noncentral T inverse
*function written in Java by a numerical analyst becomes
*available, it should be used rather than this routine.
*
*
*@param   p     The method returns the x value
*               for which CDF_nct(delta,dof)(x) = p.
*@param   delta The noncentrality parameter
*@param   dof   The degrees of freedom
*
*@version .5 --- February 26, 2002
*
*/


   public static double fnct_inv(double p, double delta, double dof) {

      double norminv,add,add0,start,pdelta,pp,plow,phigh,xx;

      double b[] = new double[2];
      double c[] = new double[2];
      double r,re,ae;

      int iflag[] = new int[2];

      CDF_nct_Amos cdf_nct_amos = new CDF_nct_Amos();

      cdf_nct_amos.p_inv = p;
      cdf_nct_amos.delta_inv = delta;
      cdf_nct_amos.dof_inv = dof;

      norminv = CDF_Normal.xnormi(p);
      add = Math.sqrt(delta*delta/(2.0*dof) + 1.0);
      add0 = norminv*add;

// an initial estimate of the inverse:

      xx = delta + add0;

      pp = CDF_nct_Amos.fnct(xx,delta,dof);

// obtaining a region in which the inverse lies:

      if (pp < p) {

         plow = pp;
         b[1] = xx;

         while (true) {

            xx += add;

            phigh = CDF_nct_Amos.fnct(xx,delta,dof);

            if (phigh >= p) break;

         }

         c[1] = xx;

      } else {

         phigh = pp;
         c[1] = xx;

         while (true) {

            xx -= add;

            plow = CDF_nct_Amos.fnct(xx,delta,dof);

            if (plow <= p) break;

         }

         b[1] = xx;

      }

      r = b[1] + ((p - plow)/(phigh - plow))*(c[1] - b[1]);

      re = .000001;
      ae = .00001;

      Fzero.fzero(cdf_nct_amos,b,c,r,re,ae,iflag);

      return b[1];

   }



/**
*
*This method calculates the t(h,a) integral of Owen.
*It is a Java translation of the public domain FORTRAN
*routine tha (dated January, 1975) written by Don Amos and S.L. Daniel.
*Steve Verrill translated the FORTRAN code 
*into Java in February of 2002.
*
*@param   h     
*@param   a
*@param   rel[1]   Relative error parameter for gaus8
*                  rel = 5.0e-(s+1) for s significant digits
*@param   ans[1]
*
*@version .5 --- February 17, 2002
*
*/


   public static void tha(double h, double a, double rel[], double ans[]) {


/*

Here is a copy of the documentation in the FORTRAN code:

      SUBROUTINE THA(H,A,REL,ANS) 
C 
C     WRITTEN BY D.E. AMOS AND S.L. DANIEL, JANUARY, 1975.
C 
C     REFERENCE SLA-73-0333 
C 
C     ABSTRACT
C         SUBROUTINE THA COMPUTES THE T(H,A) INTEGRAL OF OWEN,
C         T(H,A)=INTEGRAL ON (0,A) OF EXP(-H*H*(1+X*X)/2)/(1+X*X).
C         THE IDENTITIES
C                     T(-H,A)=T(H,A) ,  T(H,-A)=-T(H,A) 
C         T(H,A)=G(U)+G(V)-4.0*G(U)*G(V)-T(V,1/A)   ,   A.GT.1
C                   G(X)=0.25*(1+ERF(X/SQRT(2)) 
C                           U=H , V=A*H 
C         AND THEIR VARIANTS ARE USED TO KEEP THE VARIABLES POSITIVE
C         AND LIMIT THE INTERVAL OF INTEGRATION TO A SUBSET OF (0,1). 
C         THE INTEGRATION IS PERFORMED BY GAUS8 ON FUNCTION FTHA.  THE
C         CASES H=0 AND A=0 ARE EVALUATED EXPLICITLY. 
C 
C     DESCRIPTION OF ARGUMENTS
C 
C         INPUT 
C           H      - ARGUMENT, UNRESTRICTED 
C           A      - ARGUMENT, UNRESTRICTED 
C           REL    - RELATIVE ERROR PARAMETER FOR GAUS8 
C                    REL=5.E-(S+1) FOR S SIGNIFICANT DIGITS 
C         OUTPUT
C           ANS    - A VALUE FOR T(H,A) 
C 
C     ERROR CONDITIONS
C         ERROR MESSAGE FROM GAUS8 - A FATAL ERROR
C 
C     THA    USES COMMON BLOCKS MLTHA 
C     THA    USES SUBROUTINES RERF,FTHA,GAUS8,I1MACH,R1MACH,XERROR
C-----------------------------------------------------------------------

*/

      final double otpi = 1.59154943091895E-1;
      final double srt2 = 1.41421356237310E0;

      double aa,hh,ab,elim,h2,c1,c2,t1,t2,t3,t4,t5;

      int ierr[] = new int[1];
      int nr[] = new int[1];

      CDF_nct_Amos cdf_nct_amos = new CDF_nct_Amos();


      aa = a;
      hh = h;

      if (aa == 0.0) {

         ans[0] = 0.0;

         return;

      }

      ab = Math.abs(aa);

      if (hh == 0.0) {

         ans[0] = otpi*Math.atan(ab);

         if (aa < 0.0) ans[0] = -ans[0];

         return;

      }

      if (hh < 0.0) hh = Math.abs(hh);

//      kk = Math.min(Math.abs(i1mach[15]),Math.abs(i1mach[16]);
//  For IEEE 754 arithmetic, the double precision smallest exponent
//  is -1022, the double precision largest exponent is 1023.
//  This yields kk = 1022.  zzz = log_10(2) = .3010299957.
//  Presumably we don't need many digits here, but ...

//      elim = 2.303*(kk*zzz - 3.0);

      elim = 2.303*(1022.0*.3010299957 - 3.0);

      if (ab <= 1.0) {

         h2 = hh*hh;

         cdf_nct_amos.h2o2 = 0.5*h2;

         if (cdf_nct_amos.h2o2 >= elim) {

            ans[0] = 0.0;

            return;

         }

         ans[0] = Gaus8.gaus8(cdf_nct_amos,0.0,ab,rel,ierr);

         if (ierr[0] == 2) {

            System.out.print("\n\nERROR: In tha, gaus8 returned ierr[0]" +
            " equal to 2.\n\n");

            return;

         }

         ans[0] *= otpi;

         if (aa < 0.0) ans[0] = -ans[0];

         return;

      } else {

// 14

         c1 = hh/srt2;

         hh *= ab;

         c2 = hh/srt2;

         h2 = hh*hh;

         cdf_nct_amos.h2o2 = 0.5*h2;

         ab = 1.0/ab;

         if (cdf_nct_amos.h2o2 < elim) {

            t3 = Gaus8.gaus8(cdf_nct_amos,0.0,ab,rel,ierr);

            if (ierr[0] == 2) {

               System.out.print("\n\nERROR: In tha, gaus8 returned ierr[0]" +
               " equal to 2.\n\n");

               return;

            }

         } else {

            t3 = 0.0;

         }

// 17

         if (c1 < 2.0 || c2 < 2.0) {

// Note that we are using an error function different from
// the one included in Amos and Daniel's fnct file.

            t1 = Derf.derf(c1);
            t2 = Derf.derf(c2);

            ans[0] = 0.25*(2.0 + t1 + t2 
                   - (1.0 + t1)*(1.0 + t2)) - t3*otpi;

            if (aa < 0.0) ans[0] = -ans[0];

            return;

         }

// 30

         t4 = Derf.derfc(c1);
         t5 = Derf.derfc(c2);

         ans[0] = 0.25*(t4 + t5 - t4*t5) - t3*otpi;

         if (aa < 0.0) ans[0] = -ans[0];

         return;

      }

   }



   public double f_to_integrate(double x) {

//   for tha

      double opx2,ftha;

      opx2 = 1.0 + x*x;

      ftha = Math.exp(-h2o2*opx2)/opx2;

      return ftha;

   }



   public double f_to_zero(double x) {

      double diff;

      diff = p_inv - CDF_nct_Amos.fnct(x,delta_inv,dof_inv);

      return diff;         

   }
      





}
