DensityEstimator.java

 
package umontreal.ssj.stat.density;
 
import java.util.ArrayList;
import umontreal.ssj.probdist.ContinuousDistribution;
import umontreal.ssj.stat.PgfDataTable;

 
public abstract class DensityEstimator {

   protected double[] data;

   public abstract void setData(double[] data);

   public double[] getData() {
      return data;
   }

 
   public abstract double evalDensity(double x);

   public double[] evalDensity(double[] evalPoints) {
      int k = evalPoints.length;
      double[] dens = new double[k];
      for (int j = 0; j < k; j++)
         dens[j] = evalDensity(evalPoints[j]);
      return dens;
   }

   public double[] evalDensity(double[] evalPoints, double[] data) {
      setData(data);
      return evalDensity(evalPoints);
   }

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

   public static void evalDensity(ArrayList<DensityEstimator> listDE, double[] evalPoints, double[][] data,
         ArrayList<double[][]> listDensity) {
      for (DensityEstimator de : listDE)
         listDensity.add(de.evalDensity(evalPoints, data));
   }

   public static double[] computeVariance(double[][] density) {
      int m = density.length; // number of indep. replications
      int k = density[0].length; // number of evaluation points
 
      double x, y;
 
      double meanDens[] = new double[k]; // Average value over the rep replicates
      double varDens[] = new double[k]; // Variance at each evaluation point
      // Arrays.fill(meanDens, 0.0);
      // Arrays.fill(varDens, 0.0);
      for (int r = 0; r < m; r++) {
         // Update the empirical mean and variance at each evaluation point.
         for (int j = 0; j < k; j++) {
            x = density[r][j];
            y = x - meanDens[j];
            meanDens[j] += y / (double) (r + 1);
            varDens[j] += y * (x - meanDens[j]);
         }
      }
 
      for (int j = 0; j < k; j++) // normalize
         varDens[j] /= (double) (m - 1);
 
      return varDens;
   }

   public static double computeIV(double[][] density, double a, double b, double[] variance) {
      variance = computeVariance(density);
      int k = density[0].length;
      double iv = 0.0;
      for (double var : variance)
         iv += var;
      return iv * (b - a) / (double) k;
   }

   public static void computeIV(ArrayList<double[][]> listDensity, double a, double b, ArrayList<Double> listIV) {
 
      // ArrayList<Double> returnList = new ArrayList<Double>();
      int k = (listDensity.get(0))[0].length;
      double[] variance = new double[k];
      for (double[][] density : listDensity) {
         // returnList.add(computeIV(density, a, b, variance));
         listIV.add(computeIV(density, a, b, variance));
      }
   }

   public static double[] computeMISE(ContinuousDistribution dist, double[] evalPoints, double[][] density, double a,
         double b, double[] variance, double[] sqBias, double[] mse) {
      int m = density.length;
      int k = evalPoints.length;
 
      double x, y, z;
      double trueDensity;
 
      double meanDens[] = new double[k]; // Average value over m
      // replicates
 
      for (int r = 0; r < m; r++) {
         // Update the empirical mean, sum of squares, and mse of
         // observations at each evaluation point.
         for (int j = 0; j < k; j++) {
            x = density[r][j];
            y = x - meanDens[j];
            trueDensity = dist.density(evalPoints[j]);
            z = x - trueDensity;
 
            meanDens[j] += y / (double) (r + 1);
            variance[j] += y * (x - meanDens[j]);
            mse[j] += z * z;
         }
 
      }
 
      double iv = 0.0;
      double mise = 0.0;
      for (int j = 0; j < k; j++) {
         variance[j] /= (double) (m - 1.0);
         mse[j] /= (double) m;
         sqBias[j] = mse[j] - variance[j];
 
         iv += variance[j];
         mise += mse[j];
      }
 
      double fact = (b - a) / (double) k;
      iv *= fact;
      mise *= (b - a) / ((double) k);
 
      double[] res = { iv, mise - iv, mise };
      return res;
 
   }

 
   public static void computeMISE(ContinuousDistribution dist, double[] evalPoints, ArrayList<double[][]> listDensity,
         double a, double b, ArrayList<double[]> listMISE) {
      // ArrayList<double[]> returnList = new ArrayList<double[]>();
      int k = evalPoints.length;
      // double[] variance = new double[k];
      // double[] sqBias = new double[k];
      // double[] mse = new double[k];
      double[] tmp = new double[k];
      for (double[][] density : listDensity) {
         // returnList.add(computeMISE(dist, evalPoints, density, a, b, tmp, tmp, tmp));
         listMISE.add(computeMISE(dist, evalPoints, density, a, b, tmp, tmp, tmp));
      }
      // return returnList;
   }

   public abstract String toString();

   public static String plotDensity(double[] evalPoints, double[] density, String plotTitle, String[] axisTitles) {
      double[][] plotData = new double[evalPoints.length][];
      for (int i = 0; i < evalPoints.length; i++) {
         plotData[i] = new double[2];
         plotData[i][0] = evalPoints[i];
         plotData[i][1] = density[i];
      }
      PgfDataTable table = new PgfDataTable(plotTitle, "", axisTitles, plotData);
      StringBuffer sb = new StringBuffer("");
      sb.append(PgfDataTable.pgfplotFileHeader());
      sb.append(table.drawPgfPlotSingleCurve(plotTitle, "axis", 0, 1, 2, "", ""));
      sb.append(PgfDataTable.pgfplotEndDocument());
      return sb.toString();
   }

   public static double roughnessFunctional(double[] density, double a, double b) {
      double fac = (b - a) / (double) density.length;
      double sum = 0.0;
      for (double d : density) {
         sum += d * d;
      }
      return fac * sum;
   }
 
}