ScaledHistogram.java

package umontreal.ssj.stat;
 
/*
 * Class:        ScaledHistogram
 * Description:  A histogram whose bin heights are arbitrary real numbers
 * Environment:  Java
 * Software:     SSJ
 * Copyright (C) 2016  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author       Pierre L'Ecuyer and Mamadou Thiongane
 * @since        December 2016
 *
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 */
 
// import java.util.logging.Logger;
// import umontreal.ssj.util.PrintfFormat;

 
public class ScaledHistogram {
   protected int numBins; // number of bins
   protected double m_h; // width of each bin
   protected double m_a; // left boundary of first bin
   protected double m_b; // right boundary of last bin
   protected double[] height; // rescaled counters: height[j] is the height of bin j.
   protected double integral; // Total area under the histogram, = (b-a) x sum of heights.
 
   private ScaledHistogram() {
   }

   public ScaledHistogram(double a, double b, int numBins) {
      init(a, b, numBins);
   }

   public ScaledHistogram(TallyHistogram hist, double integral) {
      init(hist, integral);
   }

   public void init(double a, double b, int numBins) {
      if (b <= a)
         throw new IllegalArgumentException("   b <= a");
      this.numBins = numBins;
      m_h = (b - a) / numBins;
      m_a = a;
      m_b = b;
      height = new double[numBins];
      for (int i = 0; i < numBins; i++)
         height[i] = 0;
      integral = 0.0;
   }

   public void init(TallyHistogram hist, double integral) {
      m_a = hist.getA();
      m_b = hist.getB();
      m_h = hist.getH();
      numBins = hist.numBins;
      height = new double[numBins];
      this.integral = integral;
      int count[] = hist.getCounters();
      double scaleFactor = integral / (hist.numberObs() * m_h);
      for (int i = 0; i < numBins; i++)
         height[i] = count[i] * scaleFactor;
   }

   public void init() {
      for (int i = 0; i < numBins; i++)
         height[i] = 0.0;
   }

   public void rescale(double integral) {
      double scaleFactor = integral / this.integral;
      for (int i = 0; i < numBins; i++)
         height[i] *= scaleFactor;
      this.integral = integral;
   }

   public ScaledHistogram averageShiftedHistogram(int r) {
      ScaledHistogram image = clone();
      double[] heightNew = image.getHeights();
      double rscale = 1.0 / (r * r); // Rescaling to be made for each bin.
      double sum = 0.0;
      for (int k = 0; k < numBins; k++) {
         heightNew[k] = r * height[k];
         for (int ell = 1; ell < r; ell++) {
            if (k - ell >= 0)
               heightNew[k] += (r - ell) * height[k - ell];
            if (k + ell < numBins)
               heightNew[k] += (r - ell) * height[k + ell];
         }
         heightNew[k] *= rscale;
         sum += heightNew[k];
      }
      image.height = heightNew;
      image.integral = sum * m_h;
      return image;
   }

   public ScaledHistogram averageShiftedHistogramTrunc(int r) {
      ScaledHistogram image = clone();
      double[] heightNew = image.getHeights();
      // double rscale = 1.0 / (r * r); // Rescaling to be made for each bin.
      double sum = 0.0;
      for (int k = 0; k < numBins; k++) {
         heightNew[k] = r * height[k]; //
         int weight = r;
         for (int ell = 1; ell < r; ell++) {
            if (k - ell >= 0) {
               heightNew[k] += (r - ell) * height[k - ell];
               weight += r - ell;
            }
            if (k + ell < numBins) {
               heightNew[k] += (r - ell) * height[k + ell];
               weight += r - ell;
            }
         }
         heightNew[k] *= 1.0 / (double) weight;
         sum += heightNew[k];
      }
      image.height = heightNew;
      image.integral = sum * m_h;
      return image;
   }

   public ScaledHistogram averageShiftedHistogram(int r, double[] w) {
      ScaledHistogram image = clone();
      double[] heightNew = image.getHeights();
      double weight = w[0];
      for (int ell = 1; ell < r; ell++)
         weight += 2.0 * w[ell];
      double rscale = 1.0 / weight; // Rescaling factor for each bin.
      double sum = 0.0;
      for (int k = 0; k < numBins; k++) {
         heightNew[k] = w[0] * height[k]; //
         for (int ell = 1; ell < r; ell++) {
            if (k - ell >= 0) {
               heightNew[k] += w[ell] * height[k - ell];
            }
            if (k + ell < numBins) {
               heightNew[k] += w[ell] * height[k + ell];
            }
         }
         heightNew[k] *= rscale;
         sum += heightNew[k];
      }
      image.height = heightNew;
      image.integral = sum * m_h;
      return image;
   }

   public ScaledHistogram averageShiftedHistogramTrunc(int r, double[] w) {
      ScaledHistogram image = clone();
      double[] heightNew = image.getHeights();
      // double rscale = 1.0 / (r * r); // Rescaling to be made for each bin.
      double sum = 0.0;
      for (int k = 0; k < numBins; k++) {
         heightNew[k] = w[0] * height[k]; //
         double weight = w[0];
         for (int ell = 1; ell < r; ell++) {
            if (k - ell >= 0) {
               heightNew[k] += w[ell] * height[k - ell];
               weight += w[ell];
            }
            if (k + ell < numBins) {
               heightNew[k] += w[ell] * height[k + ell];
               weight += w[ell];
            }
         }
         heightNew[k] *= 1.0 / weight;
         sum += heightNew[k];
      }
      image.height = heightNew;
      image.integral = sum * m_h;
      return image;
   }
 
   /*
    * public ScaledHistogram averageShiftedHistogram1 (int r) { // if (numBins % r
    * != 0) // throw new
    * IllegalArgumentException("r must divide the number of bins.");
    * ScaledHistogram image = clone(); double[] heightNew = image.getHeights();
    * double rscale = 1.0 / (r * r); // Rescaling to be made for each bin. double
    * sum = 0.0; double S1[] = new double[numBins]; double S2[] = new
    * double[numBins]; S1[0] = height[0]; S2[0] = 0.0; heightNew[0] = r * S1[0];
    * for (int ell = 1; ell < Math.min(r, numBins); ell++) { S2[0] += height[ell];
    * heightNew[0] += (r-ell) * height[ell]; } for (int k = 1; k < numBins; k++) {
    * for (int ell = 1; ell < r; ell++) { S1[k] = S1[k-1] + height[k]; if (k >= r)
    * S1[k] -= height[k-r]; S2[k] = S2[k-1] - height[k]; if (k+r-1 < numBins) S2[k]
    * += height[k+r-1]; heightNew[k] = heightNew[k-1] + S2[k-1] - S1[k-1]; }
    * heightNew[k] *= rscale; sum += heightNew[k]; } image.height = heightNew;
    * image.integral = sum*m_h; return image; }
    */

   public ScaledHistogram averageShiftedHistogram1(int r) {
      ScaledHistogram image = clone();
      double[] heightNew = image.getHeights();
      double rscale = 1.0 / (r * r); // Rescaling to be made for each bin.
      double sum = 0.0;
      double S1[] = new double[numBins];
      double S2[] = new double[numBins];
      S1[0] = height[0];
      S2[0] = 0.0;
      heightNew[0] = r * S1[0];
      for (int ell = 1; ell <= Math.min(r, numBins); ell++) {
         S2[0] += height[ell];
         heightNew[0] += (r - ell) * height[ell];
      }
      for (int k = 2; k <= numBins; k++) {
         S1[k - 1] = S1[k - 2] + height[k - 1];
         if (k >= r)
            S1[k - 1] -= height[k - r];
         S2[k - 1] = S2[k - 2] - height[k - 1];
         if (k + r < numBins)
            S2[k - 1] += height[k + r - 1];
         heightNew[k - 1] = heightNew[k - 2] + S2[k - 2] - S1[k - 2];
      }
      for (int k = 0; k < numBins; k++) {
         heightNew[k] *= rscale;
         sum += heightNew[k];
      }
      image.height = heightNew;
      image.integral = sum * m_h;
      return image;
   }

   public int getNumBins() {
      return numBins;
   }

   public double getA() {
      return m_a;
   }

   public double getB() {
      return m_b;
   }

   public double[] getHeights() {
      return height;
   }

   public double getIntegral() {
      return integral;
   }

   public double ISEvsU01() {
      double sum = 0.0;
      for (int j = 0; j < numBins; ++j)
         sum += (height[j] - 1.0) * (height[j] - 1.0);
      return sum / numBins;
   }

   public double ISEvsU01polygonal() {
      double w0[] = new double[numBins]; // histogram shifted down to zero mean.
      for (int j = 0; j < numBins; ++j) {
         w0[j] = height[j] - 1.0;
      }
      double a, b;
      double sum = 0.5 * (w0[0] * w0[0] + w0[numBins - 1] * w0[numBins - 1]);
      for (int j = 0; j < numBins - 2; ++j) {
         a = w0[j];
         b = w0[j + 1];
         sum += 0.3333333333333333333 * (b * b + a * a + a * b);
      }
      return sum / (double) numBins;
   }

   public ScaledHistogram clone() {
      ScaledHistogram image = new ScaledHistogram();
      image.numBins = numBins;
      image.m_h = m_h;
      image.m_a = m_a;
      image.m_b = m_b;
      image.height = new double[numBins];
      image.integral = integral;
      for (int j = 1; j < numBins; ++j)
         image.height[j] = height[j];
      return image;
   }
 
}