DEKernelDensity.java

package umontreal.ssj.stat.density;
 
import umontreal.ssj.probdist.ContinuousDistribution;

 
public class DEKernelDensity extends DensityEstimator {

   private ContinuousDistribution kernel;

   private double h;

   private double eps = 1.0E-10;

   public DEKernelDensity(double[] data) {
      this.data = data;
   }

   public DEKernelDensity(ContinuousDistribution kernel) {
      setKernel(kernel);
   }

   public DEKernelDensity(ContinuousDistribution kernel, double[] data) {
      setKernel(kernel);
      this.data = data;
   }

   public DEKernelDensity(double h) {
      setH(h);
   }

   public DEKernelDensity(double h, double[] data) {
      this(h);
      this.data = data;
   }

   public DEKernelDensity(ContinuousDistribution kernel, double h) {
      this(kernel);
      setH(h);
   }

   public DEKernelDensity(ContinuousDistribution kernel, double h, double[] data) {
      this(kernel, data);
      setH(h);
   }

   public void setData(double[] data) {
      this.data = data;
   }

   public void setH(double h) {
      this.h = h;
   }

   public double getH() {
      return h;
   }

   public ContinuousDistribution getKernel() {
      return kernel;
   }

   public void setKernel(ContinuousDistribution kernel) {
      this.kernel = kernel;
   }

   public double getEps() {
      return eps;
   }

   public void setEps(double eps) {
      this.eps = eps;
   }

   public double evalDensity(double x) {
      double density;
      int n = data.length;
      double invh = 1.0 / h;
      double invhn = invh / (double) n;
      double sum = 0.0;
      double term; // A term to be added to the sum that defines the density
                   // estimate.
      int imin = 0; // We know that the terms for i < imin do not contribute
                    // significantly.
 
      term = kernel.density((x - data[imin]) * invh);
      while ((term < eps) && (imin < n - 1) && (data[imin] < x))
         term = kernel.density((x - data[++imin]) * invh);
      sum = term; // The first significant term.
      for (int i = imin + 1; (i < n) && ((term > eps) || (data[i] < x)); i++) // Data
         // indexed
         // by i.
         sum += (term = kernel.density((x - data[i]) * invh));
      density = sum * invhn;
      return density;
   }

   @Override
   public double[] evalDensity(double[] evalPoints) {
      int k = evalPoints.length;
      double[] density = new double[k];
      int n = data.length;
      double invh = 1.0 / h;
      double invhn = invh / (double) n;
      double y;
      double sum = 0.0;
      double term; // A term to be added to the sum that defines the density
                   // estimate.
      int imin = 0; // We know that the terms for i < imin do not contribute
                    // significantly.
      for (int j = 0; j < k; j++) { // Evaluation points are indexed by j.
         y = evalPoints[j];
         term = kernel.density((y - data[imin]) * invh);
         while ((term < eps) && (imin < n - 1) && (data[imin] < y))
            term = kernel.density((y - data[++imin]) * invh);
         sum = term; // The first significant term.
         for (int i = imin + 1; (i < n) && ((term > eps) || (data[i] < y)); i++)
            // Data indexed by i.
            sum += (term = kernel.density((y - data[i]) * invh));
         density[j] = sum * invhn;
      }
 
      return density;
   }

   @Override
   public String toString() {
      String str = "KDE [h = " + h + ", Kernel: " + kernel.toString() + "]";
      return str;
   }
 
   // STATIC METHODS:

   public static double evalDensity(double x, ContinuousDistribution kernel, double h, double[] data, double eps) {
      double density;
      int n = data.length;
      double invh = 1.0 / h;
      double invhn = invh / (double) n;
      double sum = 0.0;
      double term; // A term to be added to the sum that defines the density
                   // estimate.
      int imin = 0; // We know that the terms for i < imin do not contribute
                    // significantly.
 
      term = kernel.density((x - data[imin]) * invh);
      while ((term < eps) && (imin < n - 1) && (data[imin] < x))
         term = kernel.density((x - data[++imin]) * invh);
      // System.out.println(imin);
      // pmin=imin;
      sum = term; // The first significant term.
      for (int i = imin + 1; (i < n) && ((term > eps) || (data[i] < x)); i++) // Data
         // indexed
         // by i.
         sum += (term = kernel.density((x - data[i]) * invh));
      density = sum * invhn;
      // System.out.println(density);
      return density;
   }

   public static double[] evalDensity(double[] evalPoints, ContinuousDistribution kernel, double h, double[] data,
         double eps) {
      int k = evalPoints.length;
      double[] density = new double[k];
      int n = data.length;
      double invh = 1.0 / h;
      double invhn = invh / (double) n;
      double y;
      double sum = 0.0;
      double term; // A term to be added to the sum that defines the density
                   // estimate.
      int imin = 0; // We know that the terms for i < imin do not contribute
                    // significantly.
      for (int j = 0; j < k; j++) { // Evaluation points are indexed by j.
         y = evalPoints[j];
         term = kernel.density((y - data[imin]) * invh);
         while ((term < eps) && (imin < n - 1) && (data[imin] < y))
            term = kernel.density((y - data[++imin]) * invh);
         sum = term; // The first significant term.
         for (int i = imin + 1; (i < n) && ((term > eps) || (data[i] < y)); i++)
            // Data indexed by i.
            sum += (term = kernel.density((y - data[i]) * invh));
         density[j] = sum * invhn;
      }
 
      return density;
   }

   public static double[][] evalDensity(double[] evalPoints, ContinuousDistribution kernel, double h, double[][] data,
         double eps) {
      int m = data.length;
      int k = evalPoints.length;
      double[][] density = new double[m][k];
      for (int rep = 0; rep < m; rep++)
         density[rep] = evalDensity(evalPoints, kernel, h, data[rep], eps);
      return density;
 
   }
 
}