package umontreal.iro.lecuyer.gof;


/**
 * This class provides methods to compute (or approximate)
 * the distribution functions of special types of
 * goodness-of-fit test statistics.
 * 
 */
public class FDist {
   private FDist() {}


   /**
    * Similar to
    * {@link umontreal.iro.lecuyer.probdist.KolmogorovSmirnovPlusDist KolmogorovSmirnovPlusDist}
    *  but for the case where the distribution function F has a jump of size
    *  a at a given point x0, is zero at the left of x0,
    *   and is continuous at the right of x0.
    * Restriction: 0 < a < 1.
    *   
    * @param N sample size
    * 
    *    @param a size of the jump
    * 
    *    @param x positive or negative Kolmogorov-Smirnov statistic
    * 
    *    @return the distribution function of the statistic evaluated at x
    * 
    */
   public static double kolmogorovSmirnovPlusJumpOne (int N, double a,
                                                      double x) {
      final double EPSILONLR = 1.E-15;
      final double EPSILON = 1.0E-290;
      final double NXAPARAM = 6.5;   // frontier: alternating series
      double LogCom;
      double q, p1, q1;
      double Sum = 0.0;
      double term;
      double Njreal;
      double jreal;
      int Sign;
      int j;
      int jmax;

      if (N < 1)
        throw new IllegalArgumentException (
                             "Calling kolmogorovSmirnovPlusJumpOne "+
                             "with N < 1");
      if (a >= 1.0 || a <= 0.0)
        throw new IllegalArgumentException (
                             "Calling kolmogorovSmirnovPlusJumpOne "+
                             "with a outside (0, 1)");
      if (x <= 0.0)
         return 0.0;
      if (x + a >= 1.0)
         return 1.0;
      LogCom = Math.log ((double)N);

      //--------------------------------------------------------------------
      // the alternating series is stable and fast for N*(x + a) very small
      //--------------------------------------------------------------------
      if (N*(x + a) < NXAPARAM && a + x < 0.5) {
         jmax = (int)(N*(x + a));
         for (j = 1; j <= jmax; j++) {
            jreal = j;
            Njreal = N - j;
            q = jreal/N - x;
            if ((q < 0.0 && ((j & 1) != 0)) ||
               ((q > 1.0) && (((N - j - 1) & 1) != 0)))
               Sign = -1;
            else
               Sign = 1;

            // we must avoid log (0.0)
            q1 = Math.abs (q);
            p1 = Math.abs (1.0 - q);
            if (q1 > EPSILON && p1 > EPSILON) {
               term = LogCom + jreal*Math.log (q1) +
                     (Njreal - 1.0)*Math.log (p1);
               Sum += Sign*Math.exp (term);
            }
            LogCom += Math.log (Njreal/(jreal + 1.0));
         }
         // add the term j = 0
         Sum += Math.exp ((N - 1)*Math.log (1.0 + x));
         return Sum*x;
      }

      //---------------------------------------------
      // For N (x + a) >= NxaParam or (a + x) > 0.5,
      // use the non-alternating series.
      //---------------------------------------------

      // EpsilonLR because the distribution has a jump
      jmax = (int)(N*(1.0 - a - x - EPSILONLR));
      for (j = 1; j <= jmax; j++) {
         jreal = j;
         Njreal = N - jreal;
         q = jreal/N + x;
         if (1.0 - q > EPSILON) {
            term = LogCom + (jreal - 1.0)*Math.log (q) + Njreal*Math.log (1.0 - q);
            Sum += Math.exp (term);
         }
         LogCom += Math.log (Njreal/(jreal + 1.0));
      }
      Sum *= x;

      // add the term j = 0
      if (1.0 - x > EPSILON)
         Sum += Math.exp (N*Math.log (1.0 - x));
      return 1.0 - Sum;
   }


   /**
    * Returns F(m), the distribution function of the scan statistic
    *   with parameters N and d, evaluated at m.
    *   For a description of this statistic and its distribution, see
    *   {@link FBar#scan(int,double,int) scan},
    *   which computes its complementary distribution
    *   
    * bar(F)(m) = 1 - F(m - 1).
    *  
    * @param N sample size ( >= 2)
    * 
    *    @param d length of the test interval (∈(0, 1))
    * 
    *    @param m scan statistic
    * 
    *    @return the distribution function of the statistic evaluated at m
    */
   public static double scan (int N, double d, int m) {
      return 1.0 - FBar.scan (N, d, m);
   }
}