DEDerivativeGaussian.java

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

 
public class DEDerivativeGaussian extends DensityDerivativeEstimator {

   public DEDerivativeGaussian(int order) {
      setOrder(order);
   }

   public DEDerivativeGaussian(int order, double[] data) {
      this(order);
      this.data = data;
   }

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

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

 
   public DEDerivativeGaussian(int order, double h) {
      this(order);
      setH(h);
   }

   public DEDerivativeGaussian(int order, double h, double[] data) {
      this(order, h);
      this.data = data;
 
   }

   @Override
   public double evalDensity(double x) {
      double z;
      double sign = getOrder() % 2 == 0 ? 1.0 : -1.0;
      double density = 0.0;
      int n = data.length;
      double nInv = 1.0 / (double) n;
      double hInv = 1.0 / getH();
 
      double norma = sign * nInv * Math.pow(hInv, getOrder() + 1.0);
 
      for (int i = 0; i < n; i++) {
         z = (x - data[i]) * hInv;
 
         density += NormalDist.density01(z) * hermitePoly(getOrder(), z);
      }
      density *= norma;
      return density;
   }

   @Override
   public double[] evalDensity(double[] evalPoints) {
      double z = 0.0;
      double sign = getOrder() % 2 == 0 ? 1.0 : -1.0;
      int k = evalPoints.length;
      double[] density = new double[k];
      int n = data.length;
      double nInv = 1.0 / (double) n;
      double hInv = 1.0 / getH();
 
      double norma = sign * nInv * Math.pow(hInv, getOrder() + 1.0);
 
      for (int j = 0; j < k; j++) { // evalPoints indexed by j
         for (int i = 0; i < n; i++) { // data points indexed by i
            z = (evalPoints[j] - data[i]) * hInv;
            density[j] += NormalDist.density01(z) * hermitePoly(getOrder(), z);
         }
 
         density[j] *= norma;
      }
 
      return density;
   }

   @Override
   public String toString() {
      return "DDE [Gaussian kernel with h = " + getH() + "]";
   }
 
   // STATIC METHODS
   public static double evalDensity(double x, int order, double h, double[] data) {
      double z;
      double sign = order % 2 == 0 ? 1.0 : -1.0;
      double density = 0.0;
      int n = data.length;
      double nInv = 1.0 / (double) n;
      double hInv = 1.0 / h;
 
      double norma = sign * nInv * Math.pow(hInv, order + 1.0);
 
      for (int i = 0; i < n; i++) {
         z = (x - data[i]) * hInv;
 
         density += NormalDist.density01(z) * hermitePoly(order, z);
      }
      density *= norma;
      return density;
   }

   public static double[] evalDensity(double[] evalPoints, int order, double h, double[] data) {
      double z = 0.0;
      double sign = order == 0 ? 1.0 : -1.0;
      int k = evalPoints.length;
      double[] density = new double[k];
      int n = data.length;
      double nInv = 1.0 / (double) n;
      double hInv = 1.0 / h;
 
      double norma = sign * nInv * Math.pow(hInv, order + 1.0);
 
      for (int j = 0; j < k; j++) { // evalPoints indexed by j
         for (int i = 0; i < n; i++) { // data points indexed by i
            z = (evalPoints[j] - data[i]) * hInv;
            density[j] += NormalDist.density01(z) * hermitePoly(order, z);
         }
 
         density[j] *= norma;
      }
 
      return density;
   }

   public static double[][] evalDensity(double[] evalPoints, int order, double h, double[][] data) {
      int m = data.length;
      double[][] density = new double[m][];
      for (int r = 0; r < m; r++)
         density[r] = evalDensity(evalPoints, order, h, data[r]);
      return density;
   }

   public static double hermitePoly(int r, double x) {
      if (r == 0)
         return 1.0;
      else if (r == 1)
         return x;
      else
         return hermitePoly(r - 1, x) * x - ((double) r - 1.0) * hermitePoly(r - 2, x);
   }
 
}