Tally.java

/*
 * Class:        Tally
 * Description:  statistical collector
 * Environment:  Java
 * Software:     SSJ
 * Copyright (C) 2001--2018  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author
 * @since
 *
 *
 * 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.
 *
 */
package umontreal.ssj.stat;
 
import umontreal.ssj.util.PrintfFormat;
import umontreal.ssj.probdist.StudentDist;
import umontreal.ssj.probdist.NormalDist;
import umontreal.ssj.probdist.ChiSquareDist;
import java.util.logging.Level;
import java.util.logging.Logger;

public class Tally extends StatProbe implements Cloneable {
   protected int numObs;
   // private double sumSquares;
   private double curAverage; // The average of the first numObs observations
   private double curSum2; // The sum (xi - average)^2 of the first numObs
                           // observations.
   private Logger log = Logger.getLogger("umontreal.ssj.stat");
 
   protected static enum CIType {
      CI_NONE, CI_NORMAL, CI_STUDENT
   };
 
   protected CIType confidenceInterval = CIType.CI_NONE;
   protected double level = 0.95;
   protected int digits = 3;

   public Tally() {
      super();
      init();
   }

   public Tally(String name) {
      super();
      this.name = name;
      init();
   }

   public void setName(String name) {
      this.name = name;
   }
 
   public void init() {
      maxValue = Double.NEGATIVE_INFINITY;
      minValue = Double.POSITIVE_INFINITY;
      sumValue = 0.0;
      // sumSquares = 0.0;
      curAverage = 0.0;
      curSum2 = 0.0;
      numObs = 0;
   }

   public void add(double x) {
      if (collect) {
         if (x < minValue)
            minValue = x;
         if (x > maxValue)
            maxValue = x;
         numObs++;
         // Algorithme dans Knuth ed. 3, p. 232; voir Wikipedia
         // http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#cite_note-1
         double y = x - curAverage;
         curAverage += y / numObs;
         curSum2 += y * (x - curAverage);
         // On pourrait utiliser l'algorithme correcteur de Kahan pour une
         // meilleure precision.
         // (voir http://en.wikipedia.org/wiki/Kahan_summation_algorithm)
      }
      notifyListeners(x);
   }

   public void add(double[] x, int number) {
      if (collect)
         for (int i = 0; i < number; i++)
            add(x[i]);
   }

   public int numberObs() {
      return numObs;
   }
 
   @Override
   public double sum() {
      return numObs * curAverage;
   }

   public double average() {
      if (numObs < 1) {
         // System.err.println (
         // "******* Tally " + name + ": calling average() with " + numObs +
         // " Observation");
         log.logp(Level.WARNING, "Tally", "average",
               "Tally " + name + ":   calling average() with " + numObs + " observation");
         return Double.NaN;
      }
      return curAverage;
   }

   public double variance() {
      // throws NumberObservationException {
      // if (numObs < 2) throw NumberObservationException;
      if (numObs < 2) {
         // System.err.println (
         // "******* Tally " + name + ": calling variance() with " + numObs +
         // " Observation");
         log.logp(Level.WARNING, "Tally", "variance",
               "Tally " + name + ":   calling variance() with " + numObs + " observation");
         return Double.NaN;
      }
      return curSum2 / (numObs - 1);
   }

   public double standardDeviation() {
      return Math.sqrt(variance());
   }

   public void confidenceIntervalNormal(double level, double[] centerAndRadius) {
      // Must return an array object, cannot return 2 doubles directly
      double z;
      if (numObs < 2)
         throw new RuntimeException("Tally " + name + ": Calling confidenceIntervalStudent with < 2 Observations");
      centerAndRadius[0] = average();
      z = NormalDist.inverseF01(0.5 * (level + 1.0));
      centerAndRadius[1] = z * Math.sqrt(variance() / (double) numObs);
   }

   public void confidenceIntervalStudent(double level, double[] centerAndRadius) {
      // Must return an array object, cannot return 2 doubles directly
      double t;
      if (numObs < 2)
         throw new RuntimeException("Tally " + name + ": Calling confidenceIntervalStudent with < 2 Observations");
      centerAndRadius[0] = average();
      t = StudentDist.inverseF(numObs - 1, 0.5 * (level + 1.0));
      centerAndRadius[1] = t * Math.sqrt(variance() / (double) numObs);
   }

   public String formatCINormal(double level, int d) {
      PrintfFormat str = new PrintfFormat();
      double ci[] = new double[2];
      confidenceIntervalNormal(level, ci);
      str.append("  " + (100 * level) + "%");
      str.append(" conf. interval for the mean (normal approx.): (");
      str.append(7 + d, d - 1, d, ci[0] - ci[1]).append(',');
      str.append(7 + d, d - 1, d, ci[0] + ci[1]).append(" )" + PrintfFormat.NEWLINE);
      return str.toString();
   }

   public String formatCINormal(double level) {
      return formatCINormal(level, 3);
   }

   public String formatCIStudent(double level, int d) {
      PrintfFormat str = new PrintfFormat();
      double ci[] = new double[2];
      confidenceIntervalStudent(level, ci);
      str.append("  " + (100 * level) + "%");
      str.append(" conf. interval for the mean (Student approx.): (");
      str.append(7 + d, d, d - 1, ci[0] - ci[1]).append(',');
      str.append(7 + d, d, d - 1, ci[0] + ci[1]).append(" )" + PrintfFormat.NEWLINE);
      return str.toString();
   }

   public String formatCIStudent(double level) {
      return formatCIStudent(level, 3);
   }

   public void confidenceIntervalVarianceChi2(double level, double[] interval) {
      // Must return an array object, cannot return 2 doubles directly
      if (numObs < 2)
         throw new RuntimeException(
               "Tally " + name + ":   calling confidenceIntervalVarianceChi2 with < 2 observations");
      double w = (numObs - 1) * variance();
      double x2 = ChiSquareDist.inverseF(numObs - 1, 0.5 * (1.0 + level));
      double x1 = ChiSquareDist.inverseF(numObs - 1, 0.5 * (1.0 - level));
      interval[0] = w / x2;
      interval[1] = w / x1;
   }

   public String formatCIVarianceChi2(double level, int d) {
      PrintfFormat str = new PrintfFormat();
      double ci[] = new double[2];
      confidenceIntervalVarianceChi2(level, ci);
      str.append("  " + (100 * level) + "%");
      str.append(" conf. interval for the variance (chi2 approx.): (");
      str.append(7 + d, d, d - 1, ci[0]).append(',');
      str.append(7 + d, d, d - 1, ci[1]).append(" )" + PrintfFormat.NEWLINE);
      return str.toString();
   }

   public String report() {
      return report(level, digits);
   }

   public String report(double level, int d) {
      PrintfFormat str = new PrintfFormat();
      str.append("REPORT on Tally stat. collector ==> " + name);
      str.append(
            PrintfFormat.NEWLINE + "    num. obs.      min          max        average     variance    standard dev."
                  + PrintfFormat.NEWLINE);
      str.append(7 + d, (int) numObs);
      str.append(" ");
      str.append(9 + d, d, d - 1, (double) minValue);
      str.append(" ");
      str.append(9 + d, d, d - 1, (double) maxValue);
      str.append(" ");
      str.append(9 + d, d, d - 1, (double) average());
      str.append(" ");
      str.append(9 + d, d, d - 1, (double) variance());
      str.append(" ");
      str.append(9 + d, d, d - 1, standardDeviation());
      str.append(PrintfFormat.NEWLINE);
 
      switch (confidenceInterval) {
      case CI_NORMAL:
         str.append(formatCINormal(level, d));
         break;
      case CI_STUDENT:
         str.append(formatCIStudent(level, d));
         break;
      case CI_NONE:
         break;
      }
 
      return str.toString();
   }
 
   public String shortReportHeader() {
      PrintfFormat pf = new PrintfFormat();
      if (showNobs)
         pf.append(-8, "num obs.").append("  ");
      pf.append(-8, "   min").append("   ");
      pf.append(-8, "   max").append("   ");
      pf.append(-8, "   average").append("   ");
      pf.append(-8, "   variance").append("   ");
      pf.append(-8, "   std. dev.");
      if (confidenceInterval != CIType.CI_NONE)
         pf.append("   ").append(-12, "conf. int.");
 
      return pf.toString();
   }

   public String shortReport() {
      PrintfFormat pf = new PrintfFormat();
      if (showNobs)
         pf.append(-8, numberObs());
      pf.append(9, 3, 2, min()).append("   ");
      pf.append(9, 3, 2, max()).append("   ");
      pf.append(10, 3, 2, average()).append("   ");
      if (numberObs() >= 2) {
         pf.append(10, 3, 2, variance()).append("   ");
         pf.append(11, 3, 2, standardDeviation());
      } else
         pf.append(21, "---");
 
      if (confidenceInterval != CIType.CI_NONE) {
         double[] ci = new double[2];
         switch (confidenceInterval) {
         case CI_NORMAL:
            confidenceIntervalNormal(level, ci);
            break;
         case CI_STUDENT:
            confidenceIntervalStudent(level, ci);
            break;
         case CI_NONE:
            break;
         }
         pf.append("   ").append((100 * level) + "% (");
         pf.append(9, 3, 2, ci[0] - ci[1]).append(',');
         pf.append(9, 3, 2, ci[0] + ci[1]).append(")");
      }
      return pf.toString();
   }

   public String reportAndCIStudent(double level, int d) {
      CIType oldCIType = confidenceInterval;
 
      try {
         confidenceInterval = CIType.CI_STUDENT;
         return report(level, d);
      } finally {
         confidenceInterval = oldCIType;
      }
   }

   public String reportAndCIStudent(double level) {
      return reportAndCIStudent(level, 3);
   }

   public double getConfidenceLevel() {
      return level;
   }

   public void setConfidenceLevel(double level) {
      if (level < 0.0)
         throw new IllegalArgumentException("level < 0");
      if (level >= 1.0)
         throw new IllegalArgumentException("level >= 1");
      this.level = level;
   }

   public void setConfidenceIntervalNone() {
      confidenceInterval = CIType.CI_NONE;
   }

   public void setConfidenceIntervalNormal() {
      confidenceInterval = CIType.CI_NORMAL;
   }

   public void setConfidenceIntervalStudent() {
      confidenceInterval = CIType.CI_STUDENT;
   }

   public void setShowNumberObs(boolean showNumObs) {
      showNobs = showNumObs;
   }

   public Tally clone() {
      try {
         return (Tally) super.clone();
      } catch (CloneNotSupportedException e) {
         throw new IllegalStateException("This Tally cannot be cloned");
      }
   }
 
}