FunctionOfMultipleMeansTally.java

/*
 * Class:        FunctionOfMultipleMeansTally
 * Description:  statistical collector for a function of multiple means
 * Environment:  Java
 * Software:     SSJ
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author
 * @since
 
 * SSJ is free software: you can redistribute it and/or modify it under
 * the terms of the GNU General Public License (GPL) as published by the
 * Free Software Foundation, either version 3 of the License, or
 * any later version.
 
 * SSJ is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 
 * A copy of the GNU General Public License is available at
   <a href="http://www.gnu.org/licenses">GPL licence site</a>.
 */
package umontreal.ssj.stat;
 
import umontreal.ssj.util.MultivariateFunction;
import umontreal.ssj.stat.list.ListOfTalliesWithCovariance;
import umontreal.ssj.util.PrintfFormat;
import umontreal.ssj.probdist.NormalDist;
import java.util.logging.Level;
import java.util.logging.Logger;

public class FunctionOfMultipleMeansTally extends StatProbe implements Cloneable {
   private ListOfTalliesWithCovariance<Tally> ta;
   protected MultivariateFunction func;
   private double[] temp;
   private double[] delta;
 
   private Logger log = Logger.getLogger("umontreal.ssj.stat");
 
   private static enum CIType {
      CI_NONE, CI_DELTA
   };
 
   protected CIType confidenceInterval = CIType.CI_NONE;
   protected double level = 0.95;

   public FunctionOfMultipleMeansTally(MultivariateFunction func, int d) {
      ta = ListOfTalliesWithCovariance.createWithTally(d);
      ta.setUnmodifiable();
      this.func = func;
      internalInit();
   }

   public FunctionOfMultipleMeansTally(MultivariateFunction func, String name, int d) {
      ta = ListOfTalliesWithCovariance.createWithTally(d);
      ta.setName(name);
      ta.setUnmodifiable();
      this.func = func;
      this.name = name;
      internalInit();
   }

   public FunctionOfMultipleMeansTally(MultivariateFunction func, ListOfTalliesWithCovariance<Tally> ta) {
      if (ta == null)
         throw new NullPointerException("The list of tallies cannot be null");
      this.ta = ta;
      this.func = func;
      ta.setUnmodifiable();
      name = ta.getName();
      internalInit();
   }
 
   public void setName(String name) {
      super.setName(name);
      ta.setName(name);
   }

   public ListOfTalliesWithCovariance<Tally> getListOfTallies() {
      return ta;
   }

   public MultivariateFunction getFunction() {
      return func;
   }

   public int getDimension() {
      return ta.size();
   }

   public void add(double... x) {
      if (x.length != ta.size())
         throw new IllegalArgumentException("Incompatible length of vectors of observations: given length is "
               + x.length + ", but required length is " + ta.size());
      if (collect)
         ta.add(x);
   }

   public int numberObs() {
      return ta.numberObs();
   }

   public double average() {
      if (temp.length != ta.size())
         // The dimension of the internal tally could (rarely) change.
         temp = new double[ta.size()];
      ta.average(temp);
      return func.evaluate(temp);
   }

   public double variance() {
      if (ta.numberObs() < 2) {
         // System.out.println
         // ("******* FunctionOfMultipleMeansTally with name " + getName()
         // + ": Calling variance() with "
         // + ta.numberObs() + " observation");
         log.logp(Level.WARNING, "FunctionOfMultipleMeansTally", "variance", "FunctionOfMultipleMeansTally " + name
               + ":   calling variance() with " + ta.numberObs() + " observation");
         return Double.NaN;
      }
      if (temp.length != ta.size()) {
         // The dimension of the internal tally could (rarely) change.
         temp = new double[ta.size()];
         delta = new double[ta.size()];
      }
      ta.average(temp);
      for (int i = 0; i < delta.length; i++)
         delta[i] = func.evaluateGradient(i, temp);
      for (int i = 0; i < temp.length; i++) {
         temp[i] = 0;
         for (int j = 0; j < delta.length; j++)
            temp[i] += delta[j] * ta.covariance(j, i);
      }
      double sigma = 0;
      for (int i = 0; i < temp.length; i++)
         sigma += temp[i] * delta[i];
      // sigma could be negative but near 0, e.g., -1e-13, so round it
      // to 0 if it is negative.
      return sigma < 0 ? 0 : sigma;
   }

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

   public void confidenceIntervalDelta(double level, double[] centerAndRadius) {
      // Must return an array object, cannot return 2 doubles directly.
      centerAndRadius[0] = average();
      double z = NormalDist.inverseF01(0.5 * (level + 1.0));
      centerAndRadius[1] = z * Math.sqrt(variance() / (double) numberObs());
   }

   public String formatCIDelta(double level, int d) {
      PrintfFormat str = new PrintfFormat();
      double ci[] = new double[2];
      confidenceIntervalDelta(level, ci);
 
      str.append("  " + (100 * level) + "%");
      str.append(" confidence interval for function of means: (");
      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 formatCIDelta(double level) {
      return formatCIDelta(level, 3);
   }

   public String report() {
      PrintfFormat str = new PrintfFormat();
      str.append("REPORT on Tally stat. collector ==> " + name);
      str.append(PrintfFormat.NEWLINE + "    func. of averages    standard dev.  ");
      str.append("num. obs." + PrintfFormat.NEWLINE);
      str.append(20, 3, 2, (double) average());
      str.append(13, 3, 2, standardDeviation());
      str.append(13, (int) numberObs()).append(PrintfFormat.NEWLINE);
 
      if (confidenceInterval == CIType.CI_DELTA) {
         str.append(formatCIDelta(level));
      }
 
      return str.toString();
   }

   public String reportAndCIDelta(double level, int d) {
      CIType oldCIType = confidenceInterval;
      double oldLevel = this.level;
 
      try {
         confidenceInterval = CIType.CI_DELTA;
         this.level = level;
         return report();
      } finally {
         confidenceInterval = oldCIType;
         this.level = oldLevel;
      }
   }

   public String reportAndCIDelta(double level) {
      return reportAndCIDelta(level, 3);
   }
 
   public String shortReportHeader() {
      PrintfFormat pf = new PrintfFormat();
      pf.append(-20, " func. of averages").append("   ");
      pf.append(-9, "std. dev.");
      if (showNobs)
         pf.append("   ").append(-5, "nobs.");
      if (confidenceInterval != CIType.CI_NONE)
         pf.append("   ").append(-25, "conf. int.");
 
      return pf.toString();
   }

   public String shortReport() {
      PrintfFormat pf = new PrintfFormat();
      pf.append(20, 3, 2, average()).append("   ");
      if (numberObs() >= 2)
         pf.append(9, 3, 2, standardDeviation());
      else
         pf.append(9, "---");
      if (showNobs)
         pf.append("   ").append(5, numberObs());
 
      if (confidenceInterval == CIType.CI_DELTA) {
         double[] ci = new double[2];
         confidenceIntervalDelta(level, ci);
         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 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 setConfidenceIntervalDelta() {
      confidenceInterval = CIType.CI_DELTA;
   }

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

   public double min() {
      return Double.NaN;
   }

   public double max() {
      return Double.NaN;
   }

   public double sum() {
      return Double.NaN;
   }
 
   private void internalInit() {
      if ((func.getDimension() != ta.size()) && (func.getDimension() != -1))
         throw new IllegalArgumentException("The dimension of the function must be equal to d, or equal to -1");
      temp = new double[ta.size()];
      delta = new double[ta.size()];
   }
 
   public void init() {
      ta.init();
      internalInit();
   }

   public FunctionOfMultipleMeansTally clone() {
      FunctionOfMultipleMeansTally mta;
      try {
         mta = (FunctionOfMultipleMeansTally) super.clone();
      } catch (CloneNotSupportedException cne) {
         throw new IllegalStateException("CloneNotSupportedException for a class implementing Cloneable");
      }
 
      mta.ta = new ListOfTalliesWithCovariance<Tally>();
      for (int i = 0; i < ta.size(); i++) {
         mta.ta.add(ta.get(i).clone());
      }
 
      mta.temp = new double[ta.size()];
      mta.delta = new double[ta.size()];
      return mta;
   }
}