DensityDerivativeEstimator.java

package umontreal.ssj.stat.density;

public abstract class DensityDerivativeEstimator extends DensityEstimator {

   private int order;

   private double h;

   public int getOrder() {
      return order;
   }

   public void setOrder(int order) {
      this.order = order;
   }

   public double getH() {
      return h;
   }

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

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

   public static double hAmiseR(int r, double mu2, double mu2Derivative, double init, int n) {
      double invN = 1.0 / (double) n;
      double sign = r % 2 == 0 ? 1.0 : -1.0;
      double h = sign * mu2Derivative * (2.0 * (double) r + 1.0) * invN / (mu2 * init);
      double exp = 1.0 / (2.0 * (double) r + 5.0);
      return Math.pow(h, exp);
   }

 
   public static double hAmiseR(int r, int t, double mu2, double[] mu2Derivative, double init,
         DensityDerivativeEstimator dde, double[] evalPoints, double a, double b) {
      double h;
      int k = evalPoints.length;
      double[] estDensity = new double[k];
      int n = dde.data.length;
 
      for (int tau = t - 1; tau >= 1; tau--) {
         h = hAmiseR(r + 2 * tau, mu2, mu2Derivative[tau], init, n);
         dde.setH(h);
         dde.setOrder(r + 2 * tau);
         estDensity = dde.evalDensity(evalPoints);
         init = roughnessFunctional(estDensity, a, b);
      }
 
      return hAmiseR(r, mu2, mu2Derivative[0], init, n);
   }

   // TODO: look for formulas with odd r.
   public static double densityFunctionalGaussian(int r, double sigma) {
      double sign = (r % 2 == 0 ? 1.0 : -1.0);
      double facTerm = 1.0;
      for (int i = r + 1; i <= 2 * r; i++) {
         facTerm *= (double) i;
      }
      double denom = Math.pow(2.0 * sigma, (double) (2.0 * r + 1.0));
 
      denom *= Math.sqrt(Math.PI);
 
      return sign * facTerm / denom;
   }

   public String toString() {
      return "DDE [h = " + h + "]";
   }
}