/* *
 * Title:        

 * Description:  

 * Copyright:    Copyright (c) 

 * Company:      

 * @author
 * @version 1.0
 */

package umontreal.iro.lecuyer.probdist;

import umontreal.iro.lecuyer.util.*;
import umontreal.iro.lecuyer.functions.MathFunction;
 

/**
 * Extends the class {@link ContinuousDistribution} for the
 * Cramér-von Mises distribution (see).
 * Given a sample of n independent uniforms Ui over [0, 1],
 * the Cramér-von Mises statistic Wn2 is defined by
 * 

 * 
 * 

 * Wn2 = 1/12n + ∑j=1n(U(j) - (j-0.5)/n)2,
 * 

 * where the U(j) are the Ui sorted in increasing order.
 *   The  distribution function (the cumulative probabilities)
 *   is defined as 
 * Fn(x) = P[Wn2 <= x].
 * 
 */
public class CramerVonMisesDist extends ContinuousDistribution {
   protected int n;

   private static class Function implements MathFunction {
      protected int n;
      protected double u;

      public Function (int n, double u) {
         this.n = n;
         this.u = u;
      }

      public double evaluate (double x) {
         return u - cdf(n,x);
      }
   }



   /**
    * Constructs a Cramér-von Mises distribution for a sample of size n.
    * 
    */
   public CramerVonMisesDist (int n) {
      setN (n);
   }


   public double density (double x) {
      return density (n, x);
   }

   public double cdf (double x) {
      return cdf (n, x);
   }

   public double barF (double x) {
      return barF (n, x);
   }

   public double inverseF (double u) {
      return inverseF (n, u);
   }

   public double getMean() {
      return CramerVonMisesDist.getMean (n);
   }

   public double getVariance() {
      return CramerVonMisesDist.getVariance (n);
   }

   public double getStandardDeviation() {
      return CramerVonMisesDist.getStandardDeviation (n);
   }

   /**
    * Computes the density function
    *   for a Cramér-von Mises distribution with parameter n.
    * 
    */
   public static double density (int n, double x) {
      if (n <= 0)
        throw new IllegalArgumentException ("n <= 0");

      if (x <= 1.0/(12.0*n) || x >= n/3.0)
         return 0.0;

      if (n == 1)
         return 1.0 / Math.sqrt (x - 1.0/12.0);

      if (x <= 0.002 || x > 3.95)
         return 0.0;

      throw new UnsupportedOperationException("density not implemented.");
   }


   /**
    * Computes the Cramér-von Mises distribution function with parameter n.
    *   Returns an approximation of 
    * P[Wn2 <= x], where Wn2 is the
    *   Cramér von Mises  statistic (see).
    *   The approximation is based on the distribution function of 
    * W2 = limn -> ∞Wn2, which has the following series expansion derived
    *   by Anderson and Darling:
    *    

    * 
    * 

    *         P(W2 <= x)  =  ∑j=0∞(- 1)j(4j+1)1/2    exp{ - }K1/4([tex2html_wrap_indisplay246]),
    * 

    * where 
    * Kν is the  modified Bessel function of the 
    *   second kind.
    *   To correct for the deviation between 
    * P(Wn2 <= x) and 
    * P(W2 <= x),
    *   we add a correction in 1/n, obtained empirically by simulation.
    *   For n = 10, 20, 40, the error is less than
    *   0.002, 0.001, and 0.0005, respectively, while for
    *   n >= 100 it is less than 0.0005.
    *   For 
    * n -> ∞, we estimate that the method returns
    *   at least 6 decimal digits of precision.
    *   For n = 1, the method uses the exact distribution:
    *   
    * P(W12 <= x) = 2(x - 1/12)1/2 for 
    * 1/12 <= x <= 1/3.
    * 
    */
   public static double cdf (int n, double x) {
      if (n <= 0)
         throw new IllegalArgumentException ("n <= 0");

      if (n == 1) {
         if (x <= 1.0/12.0)
            return 0.0;
         if (x >= 1.0/3.0)
            return 1.0;
         return 2.0*Math.sqrt (x - 1.0/12.0);
      }

      if (x <= 1.0/(12.0*n))
         return 0.0;

      if (x <= (n + 3.0)/(12.0*n*n)) {
         double t = Num.lnFactorial(n) - Num.lnGamma (1.0 + 0.5*n) +
            0.5*n*Math.log (Math.PI*(x - 1.0/(12.0*n)));
         return Math.exp(t);
      }

      if (x <= 0.002)
         return 0.0;
      if (x > 3.95 || x >= n/3.0)
         return 1.0;

      final double EPSILON = Num.DBL_EPSILON;
      final int JMAX = 20;
      int j = 0;
      double Cor, Res, arg;
      double termX, termS, termJ;

      termX = 0.0625/x;            // 1 / (16x)
      Res = 0.0;

      final double A[] = {
         1.0,
         1.11803398875,
         1.125,
         1.12673477358,
         1.1274116945,
         1.12774323743,
         1.1279296875,
         1.12804477649,
         1.12812074678,
         1.12817350091
      };

      do {
         termJ = 4*j + 1;
         arg = termJ*termJ*termX;
         termS = A[j]*Math.exp (-arg)*Num.besselK025 (arg);
         Res += termS;
         ++j;
      } while (!(termS < EPSILON || j > JMAX));

      if (j > JMAX)
         System.err.println ("cramerVonMises: iterations have not converged");
      Res /= Math.PI*Math.sqrt (x);

      // Empirical correction in 1/n
      if (x < 0.0092)
         Cor = 0.0;
      else if (x < 0.03)
         Cor = -0.0121763 + x*(2.56672 - 132.571*x);
      else if (x < 0.06)
         Cor = 0.108688 + x*(-7.14677 + 58.0662*x);
      else if (x < 0.19)
         Cor = -0.0539444 + x*(-2.22024 + x*(25.0407 - 64.9233*x));
      else if (x < 0.5)
         Cor = -0.251455 + x*(2.46087 + x*(-8.92836 + x*(14.0988 -
                  x*(5.5204 + 4.61784*x))));
      else if (x <= 1.1)
         Cor = 0.0782122 + x*(-0.519924 + x*(1.75148 +
               x*(-2.72035 + x*(1.94487 - 0.524911*x))));
      else
         Cor = Math.exp (-0.244889 - 4.26506*x);

      Res += Cor/n;
      // This empirical correction is not very precise, so ...
      if (Res <= 1.0)
         return Res;
      else
         return 1.0;
   }


   /**
    * Computes the complementary distribution function 
    * bar(F)n(x)
    *   with parameter n.
    * 
    */
   public static double barF (int n, double x) {
      return 1.0 - cdf(n,x);
   }


   /**
    * Computes 
    * x = Fn-1(u), where Fn is the 
    *   Cramér-von Mises distribution with parameter n.
    * 
    */
   public static double inverseF (int n, double u) {
      if (n <= 0)
         throw new IllegalArgumentException ("n <= 0");
      if (u < 0.0 || u > 1.0)
         throw new IllegalArgumentException ("u must be in [0,1]");

      if (u >= 1.0)
         return n/3.0;
      if (u <= 0.0)
         return 1.0/(12.0*n);

      if (n == 1)
         return 1.0/12.0 + 0.25*u*u;

      Function f = new Function (n,u);
      return RootFinder.brentDekker (0.0, 10.0, f, 1e-6);
   }


   /**
    * Returns the mean of the distribution with parameter n.
    * 
    * @return the mean
    * 
    */
   public static double getMean (int n) {
      return 1.0 / 6.0;
   }


   /**
    * Returns the variance of the distribution with parameter n.
    * 
    * @return variance
    * 
    */
   public static double getVariance (int n) {
      return (4.0 * n - 3.0) / (180.0 * n ); 
   }


   /**
    * Returns the standard deviation of the distribution with
    *    parameter n.
    * 
    * @return the standard deviation
    * 
    */
   public static double getStandardDeviation (int n) {
     return Math.sqrt(getVariance(n));
   }


   /**
    * Returns the parameter n of this object.
    * 
    */
   public int getN() {
      return n;
   }


   /**
    * Sets the parameter n of this object.
    * 
    */
   public void setN (int n) {
      if (n <= 0)
         throw new IllegalArgumentException ("n <= 0");
      this.n = n;
      supportA = 1.0/(12.0*n);
      supportB = n/3.0;
   }


   /**
    * Return an array containing the parameter n of this object.
    * 
    * 
    */
   public double[] getParams () {
      double[] retour = {n};
      return retour;
   }


   public String toString () {
      return getClass().getSimpleName() + " : n = " + n;
   }

}