package umontreal.iro.lecuyer.gof;
   import umontreal.iro.lecuyer.probdist.*;

import umontreal.iro.lecuyer.randvar.KernelDensityGen;


/**
 * This class provides methods to compute a kernel density estimator from a set
 * of n individual observations 
 * x0,…, xn-1, and returns its value
 * at m selected points. For details on how the kernel density is defined,
 * and how to select the kernel and the bandwidth h,
 * see the documentation of class 
 * {@link umontreal.iro.lecuyer.randvar.KernelDensityGen KernelDensityGen}
 * in package randvar.
 * 
 */
public class KernelDensity  {

   private static double estimate (EmpiricalDist dist,
                                   ContinuousDistribution kern,
                                   double h, double y) {
      // Computes and returns the kernel density estimate at $y$, where the 
      // kernel is the density kern.density(x), and the bandwidth is $h$.
      double z;
      double a = kern.getXinf();       // lower limit of density
      double b = kern.getXsup();       // upper limit of density
      double sum = 0;
      int n = dist.getN();
      for (int i = 0; i < n; i++) {
         z = (y - dist.getObs(i))/h;
         if ((z >= a) && (z <= b))
            sum += kern.density(z);
      }

      sum /= (h*n);
      return sum;
   }


   /**
    * Given the empirical distribution dist, this method computes the 
    * kernel density estimate at each of the m points Y[j],
    *  
    * j = 0, 1,…,(m - 1), where m is the length of Y, the kernel
    *  is kern.density(x),
    *  and the bandwidth is h. Returns the estimates as an array of m values.
    * 
    */
   public static double[] computeDensity (EmpiricalDist dist,
                                          ContinuousDistribution kern,
                                          double h, double[] Y) {
      int m = Y.length;
      double[] u = new double[m];
      for (int j = 0; j < m; j++)
         u[j] = estimate (dist, kern, h, Y[j]);
      return u;
   }


   /**
    * Similar to method  ,
    * but the bandwidth h is obtained from the method
    * {@link umontreal.iro.lecuyer.randvar.KernelDensityGen#getBaseBandwidth(EmpiricalDist) KernelDensityGen.getBaseBandwidth}(dist) in package randvar.
    * 
    */
   public static double[] computeDensity (EmpiricalDist dist,
                                          ContinuousDistribution kern,
                                          double[] Y) {
      double h = KernelDensityGen.getBaseBandwidth(dist);
      return computeDensity (dist, kern, h, Y);
   }

}