ListOfTalliesWithCV.java

/*
 * Class:        ListOfTalliesWithCV
 * Description:  List of tallies with control variables
 * Environment:  Java
 * Software:     SSJ 
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author       Eric Buist
 * @since        August 2007
 
 * 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.list.lincv;
 
import cern.colt.list.DoubleArrayList;
import cern.colt.matrix.DoubleMatrix1D;
import cern.colt.matrix.DoubleMatrix2D;
import cern.colt.matrix.impl.DenseDoubleMatrix2D;
import cern.colt.matrix.linalg.Algebra;
 
import umontreal.ssj.probdist.StudentDist;
import umontreal.ssj.stat.Tally;
import umontreal.ssj.stat.TallyStore;
import umontreal.ssj.stat.list.ListOfTallies;
import umontreal.ssj.stat.list.ListOfTalliesWithCovariance;
import umontreal.ssj.util.PrintfFormat;

public class ListOfTalliesWithCV<E extends Tally> extends ListOfTalliesWithCovariance<E> {
   private static Algebra alg = new Algebra();
   private DoubleMatrix2D beta;
   private double[] exp;
   private double[] tmp;
   private int q;
 
   private DoubleMatrix2D tempPP;
   private DoubleMatrix2D tempQQ;
   private DoubleMatrix2D tempPQ;
   private DoubleMatrix2D tempQP;

   public ListOfTalliesWithCV() {
      super();
   }

   public ListOfTalliesWithCV(String name) {
      super(name);
   }

   public static ListOfTalliesWithCV<Tally> createWithTally(int p, int q) {
      ListOfTalliesWithCV<Tally> list = new ListOfTalliesWithCV<Tally>();
      int size = p + q;
      for (int i = 0; i < size; i++)
         list.add(new Tally());
      list.setNumControlVariables(q);
      list.init();
      return list;
   }

   public static ListOfTalliesWithCV<TallyStore> createWithTallyStore(int p, int q) {
      ListOfTalliesWithCV<TallyStore> list = new ListOfTalliesWithCV<TallyStore>();
      int size = p + q;
      for (int i = 0; i < size; i++)
         list.add(new TallyStore());
      list.setNumControlVariables(q);
      list.init();
      return list;
   }
 
   public void init() {
      super.init();
      internalInit();
   }
 
   private void internalInit() {
      int p = sizeWithoutCV();
      beta = new DenseDoubleMatrix2D(q, p);
      exp = new double[q];
      tmp = new double[p + q];
 
      tempPP = new DenseDoubleMatrix2D(p, p);
      tempQQ = new DenseDoubleMatrix2D(q, q);
      tempPQ = new DenseDoubleMatrix2D(p, q);
      tempQP = new DenseDoubleMatrix2D(q, p);
      setUnmodifiable();
   }

   public DoubleMatrix2D getBeta() {
      return beta;
   }

   public void setBeta(DoubleMatrix2D beta) {
      if (beta.rows() != getNumControlVariables())
         throw new IllegalArgumentException("The number of rows in beta must be equal to q");
      if (beta.columns() != sizeWithoutCV())
         throw new IllegalArgumentException("The number of columns in beta must be equal to p");
      this.beta = beta;
   }

   public double getExpectedValue(int i) {
      return exp[i];
   }

   public void setExpectedValue(int i, double e) {
      exp[i] = e;
   }

   public double[] getExpectedValues() {
      return exp;
   }

   public void setExpectedValues(double[] exp) {
      if (exp.length != this.exp.length)
         throw new IllegalArgumentException("Invalid length of exp");
      this.exp = exp;
   }

   public int sizeWithoutCV() {
      return size() - q;
   }

   public int getNumControlVariables() {
      return q;
   }

   public void setNumControlVariables(int q) {
      if (beta != null)
         throw new IllegalArgumentException("Cannot change the number of control variables");
      if (q < 0 || q >= size())
         throw new IllegalArgumentException("q is negative or greater than or equal to " + size());
      this.q = q;
   }

   public void correlationX(DoubleMatrix2D c) {
      int l = sizeWithoutCV();
      if (c.rows() != l)
         throw new IllegalArgumentException("Invalid number of rows in covariance matrix");
      if (c.columns() != l)
         throw new IllegalArgumentException("Invalid number of columns in covariance matrix");
      for (int i1 = 0; i1 < l; i1++)
         for (int i2 = 0; i2 < l; i2++)
            c.setQuick(i1, i2, correlation(i1, i2));
   }

   public void covarianceX(DoubleMatrix2D c) {
      int l = sizeWithoutCV();
      if (c.rows() != l)
         throw new IllegalArgumentException("Invalid number of rows in covariance matrix");
      if (c.columns() != l)
         throw new IllegalArgumentException("Invalid number of columns in covariance matrix");
      for (int i1 = 0; i1 < l; i1++)
         for (int i2 = 0; i2 < l; i2++)
            c.setQuick(i1, i2, covariance(i1, i2));
   }

   public void correlationC(DoubleMatrix2D c) {
      int p = sizeWithoutCV();
      int q = getNumControlVariables();
      if (c.rows() != q)
         throw new IllegalArgumentException("Invalid number of rows in covariance matrix");
      if (c.columns() != q)
         throw new IllegalArgumentException("Invalid number of columns in covariance matrix");
      for (int i1 = 0; i1 < q; i1++)
         for (int i2 = 0; i2 < q; i2++)
            c.setQuick(i1, i2, correlation(p + i1, p + i2));
   }

   public void covarianceC(DoubleMatrix2D c) {
      int p = sizeWithoutCV();
      int q = getNumControlVariables();
      if (c.rows() != q)
         throw new IllegalArgumentException("Invalid number of rows in covariance matrix");
      if (c.columns() != q)
         throw new IllegalArgumentException("Invalid number of columns in covariance matrix");
      for (int i1 = 0; i1 < q; i1++)
         for (int i2 = 0; i2 < q; i2++)
            c.setQuick(i1, i2, covariance(p + i1, p + i2));
   }

   public void correlationCX(DoubleMatrix2D c) {
      int p = sizeWithoutCV();
      int q = getNumControlVariables();
      if (c.rows() != q)
         throw new IllegalArgumentException("Invalid number of rows in covariance matrix");
      if (c.columns() != p)
         throw new IllegalArgumentException("Invalid number of columns in covariance matrix");
      for (int i1 = 0; i1 < q; i1++)
         for (int i2 = 0; i2 < p; i2++)
            c.setQuick(i1, i2, correlation(p + i1, i2));
   }

   public void covarianceCX(DoubleMatrix2D c) {
      int p = sizeWithoutCV();
      int q = getNumControlVariables();
      if (c.rows() != q)
         throw new IllegalArgumentException("Invalid number of rows in covariance matrix");
      if (c.columns() != p)
         throw new IllegalArgumentException("Invalid number of columns in covariance matrix");
      for (int i1 = 0; i1 < q; i1++)
         for (int i2 = 0; i2 < p; i2++)
            c.setQuick(i1, i2, covariance(p + i1, i2));
   }

   public void add(double[] x, double[] c) {
      if (x.length != sizeWithoutCV())
         throw new IllegalArgumentException("Invalid length of x");
      if (c.length != getNumControlVariables())
         throw new IllegalArgumentException("Invalid length of c");
      System.arraycopy(x, 0, tmp, 0, x.length);
      System.arraycopy(c, 0, tmp, x.length, c.length);
      add(tmp);
   }

   public void add(double x, double[] c) {
      if (sizeWithoutCV() != 1)
         throw new IllegalArgumentException("Cannot use this method if p != 1");
      if (c.length != getNumControlVariables())
         throw new IllegalArgumentException("Invalid length of c");
      tmp[0] = x;
      System.arraycopy(c, 0, tmp, 1, c.length);
      add(tmp);
   }

   public void add(double x, double c) {
      if (sizeWithoutCV() != 1)
         throw new IllegalArgumentException("Cannot use this method if p != 1");
      if (getNumControlVariables() != 1)
         throw new IllegalArgumentException("Cannot use this method if q != 1");
      tmp[0] = x;
      tmp[1] = c;
      add(tmp);
   }

   public double averageWithCV(int i) {
      int p = sizeWithoutCV();
      if (i >= p)
         throw new ArrayIndexOutOfBoundsException(i);
      Tally tally = get(i);
      if (tally == null || tally.numberObs() == 0)
         return Double.NaN;
      else {
         double avg = tally.average();
         for (int j = 0; j < q; j++)
            avg -= beta.getQuick(j, i) * (get(p + j).average() - exp[j]);
         return avg;
      }
   }

   public void covarianceWithCV(DoubleMatrix2D covCV) {
      final int p = sizeWithoutCV();
      final int q = getNumControlVariables();
      if (covCV.rows() != p || covCV.columns() != p)  // ???????????
         throw new IllegalArgumentException("Invalid dimensions of covCV, (p,q) = " + p + " " + q);
      DoubleMatrix2D covX = covCV;
      covarianceX(covX);
      DoubleMatrix2D covC = tempQQ;
      covarianceC(covC);
      DoubleMatrix2D covCX = tempQP;
      covarianceCX(covCX);
 
      // covCV contains covX
      // Add beta^t*sigmaC*beta
      beta.viewDice().zMult(covC, tempPQ).zMult(beta, covCV, 1, 1, false, false);
      // Subtract 2beta^t*sigmaCX
      beta.viewDice().zMult(covCX, covCV, -2, 1, false, false);
   }

   public double covarianceWithCV(int i, int j) {
      final int p = sizeWithoutCV();
      // final int q = getNumControlVariables();
      if (i >= p || j >= p)
         throw new IllegalArgumentException("i >= p or j >= p");
      double cov = covariance(i, j);
      DoubleMatrix2D covC = tempQQ;
      covarianceC(covC);
      DoubleMatrix2D covCX = tempQP;
      covarianceCX(covCX);
      // Add beta_.i^t*sigmaC*beta.j
      cov += beta.viewColumn(i).zDotProduct(covC.zMult(beta.viewColumn(j), null));
      // Subtract sigmaCX_.i*beta_.j
      cov -= covCX.viewColumn(i).zDotProduct(beta.viewColumn(j));
      // Subtract sigmaCX_.j*beta_.i
      cov -= covCX.viewColumn(j).zDotProduct(beta.viewColumn(i));
      return cov;
   }

   public void averageWithCV(double[] a) {
      int p = sizeWithoutCV();
      if (a.length != p)
         throw new IllegalArgumentException(
               "Invalid length of the given array: given length is " + a.length + ", required length is " + p);
      for (int i = 0; i < p; i++)
         a[i] = averageWithCV(i);
   }

   public void averageX(double[] a) {
      int l = sizeWithoutCV();
      if (a.length != l)
         throw new IllegalArgumentException(
               "Invalid length of the given array: given length is " + a.length + ", required length is " + l);
      for (int i = 0; i < a.length; i++) {
         Tally tally = get(i);
         a[i] = tally == null ? Double.NaN : tally.average();
         if (tally.numberObs() == 0)
            a[i] = Double.NaN;
      }
   }

   public void averageC(double[] a) {
      int p = sizeWithoutCV();
      int l = getNumControlVariables();
      if (a.length != l)
         throw new IllegalArgumentException(
               "Invalid length of the given array: given length is " + a.length + ", required length is " + l);
      for (int i = 0; i < a.length; i++) {
         Tally tally = get(i + p);
         a[i] = tally == null ? Double.NaN : tally.average();
         if (tally.numberObs() == 0)
            a[i] = Double.NaN;
      }
   }

   public void standardDeviationWithCV(double[] std) {
      final int l = sizeWithoutCV();
      if (l != std.length)
         throw new IllegalArgumentException("Invalid length of given array");
      DoubleMatrix2D covCV = tempPP;
      covarianceWithCV(covCV);
      for (int i = 0; i < std.length; i++)
         std[i] = Math.sqrt(covCV.getQuick(i, i));
   }

   public void varianceWithCV(double[] v) {
      final int l = sizeWithoutCV();
      if (l != v.length)
         throw new IllegalArgumentException("Invalid length of given array");
      DoubleMatrix2D covCV = tempPP;
      covarianceWithCV(covCV);
      for (int i = 0; i < v.length; i++)
         v[i] = covCV.getQuick(i, i);
   }

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

   public void estimateBeta() {
      DoubleMatrix2D covC = tempQQ;
      covarianceC(covC);
      DoubleMatrix2D covCX = tempQP;
      covarianceCX(covCX);
      beta = alg.solve(covC, covCX);
      assert beta.rows() == getNumControlVariables() && beta.columns() == sizeWithoutCV();
   }

   /*
    * public void addControlledObservations (ListOfTallies<? extends Tally> a) {
    * int p = sizeWithoutCV(); if (a.size() != p) throw new
    * IllegalArgumentException ("Invalid length of the given array of tallies"); if
    * (p == 0) return; int numObs = numberObs(); double[] tmp = new double[p]; for
    * (int obs = 0; obs < numObs; obs++) { for (int i = 0; i < p; i++) { TallyStore
    * tallyStore = (TallyStore)get (i); double v = tallyStore.getArray().getQuick
    * (obs); for (int j = 0; j < q; j++) { tallyStore = (TallyStore)get (j + p);
    * double c = tallyStore.getArray().getQuick (obs); v -= beta.getQuick (j, i)*(c
    * - exp[j]); } tmp[i] = v; } a.add (tmp); } }
    */

   /*
    * public void addControlledObservations (Tally ta) { int p = sizeWithoutCV();
    * int q = getNumControlVariables(); if (p != 1) throw new
    * IllegalArgumentException ("This method can only be called with p=1"); int
    * numObs = get (0).numberObs(); ta.init(); TallyStore tallyStore =
    * (TallyStore)get (0); DoubleArrayList array = tallyStore.getArray(); for (int
    * obs = 0; obs < numObs; obs++) { double v = array.getQuick (obs); for (int j =
    * 0; j < q; j++) { tallyStore = (TallyStore)get (j + p); double c =
    * tallyStore.getArray().getQuick (obs); v -= beta.getQuick (j, 0)*(c - exp[j]);
    * } ta.add (v); } }
    */

   public ListOfTalliesWithCV<E> clone() {
      ListOfTalliesWithCV<E> l = (ListOfTalliesWithCV<E>) super.clone();
      l.beta = (DoubleMatrix2D) beta.clone();
      l.exp = exp.clone();
      l.tmp = tmp.clone();
      l.q = q;
 
      l.tempPP = (DoubleMatrix2D) tempPP.clone();
      l.tempQQ = (DoubleMatrix2D) tempQQ.clone();
      l.tempPQ = (DoubleMatrix2D) tempPQ.clone();
      l.tempQP = (DoubleMatrix2D) tempQP.clone();
 
      return l;
   }
 
}