MultinomialDist.java

/*
 * Class:        MultinomialDist
 * Description:  multinomial distribution
 * Environment:  Java
 * Software:     SSJ
 * Copyright (C) 2001  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.probdistmulti;
 
import umontreal.ssj.util.Num;

public class MultinomialDist extends DiscreteDistributionIntMulti {
   protected int n;
   protected double p[];

   public MultinomialDist(int n, double p[]) {
      setParams(n, p);
   }
 
   public double prob(int x[]) {
      return prob_(n, p, x);
   }
 
   public double cdf(int x[]) {
      return cdf_(n, p, x);
   }
 
   public double[] getMean() {
      return getMean_(n, p);
   }
 
   public double[][] getCovariance() {
      return getCovariance_(n, p);
   }
 
   public double[][] getCorrelation() {
      return getCorrelation_(n, p);
   }
 
   private static void verifParam(int n, double p[]) {
      if (n <= 0)
         throw new IllegalArgumentException("n <= 0");
 
      double sumPi = 0.0;
      for (int i = 0; i < p.length; i++) {
         if ((p[i] < 0) || (p[i] > 1))
            throw new IllegalArgumentException("p is not a probability vector");
         sumPi += p[i];
      }
 
      if (Math.abs(sumPi - 1.0) > 1.0e-15)
         throw new IllegalArgumentException("p is not a probability vector");
   }
 
   private static double prob_(int n, double p[], int x[]) {
      if (x.length != p.length)
         throw new IllegalArgumentException("x and p must have the same dimension");
 
      double sumXFact = 0.0;
      int sumX = 0;
      double sumPX = 0.0;
 
      for (int i = 0; i < p.length; i++) {
         sumX += x[i];
         sumXFact += Num.lnFactorial(x[i]);
         sumPX += (x[i] * Math.log(p[i]));
      }
 
      if (sumX != n)
         return 0.0;
      else {
         return Math.exp(Num.lnFactorial(n) - sumXFact + sumPX);
      }
   }

   public static double prob(int n, double p[], int x[]) {
      verifParam(n, p);
      return prob_(n, p, x);
   }
 
   private static double cdf_(int n, double p[], int x[]) {
      boolean end = false;
      double sum = 0.0;
      int j;
 
      if (x.length != p.length)
         throw new IllegalArgumentException("x and p must have the same dimension");
 
      int is[] = new int[x.length];
      for (int i = 0; i < is.length; i++)
         is[i] = 0;
 
      sum = 0.0;
      while (!end) {
         sum += prob(n, p, is);
         is[0]++;
 
         if (is[0] > x[0]) {
            is[0] = 0;
            j = 1;
            while (j < x.length && is[j] == x[j])
               is[j++] = 0;
 
            if (j == x.length)
               end = true;
            else
               is[j]++;
         }
      }
 
      return sum;
   }

   public static double cdf(int n, double p[], int x[]) {
      verifParam(n, p);
      return cdf_(n, p, x);
   }
 
   private static double[] getMean_(int n, double[] p) {
      double mean[] = new double[p.length];
 
      for (int i = 0; i < p.length; i++)
         mean[i] = n * p[i];
 
      return mean;
   }

   public static double[] getMean(int n, double[] p) {
      verifParam(n, p);
 
      return getMean_(n, p);
   }
 
   private static double[][] getCovariance_(int n, double[] p) {
      double cov[][] = new double[p.length][p.length];
 
      for (int i = 0; i < p.length; i++) {
         for (int j = 0; j < p.length; j++)
            cov[i][j] = -n * p[i] * p[j];
 
         cov[i][i] = n * p[i] * (1.0 - p[i]);
      }
      return cov;
   }

   public static double[][] getCovariance(int n, double[] p) {
      verifParam(n, p);
      return getCovariance_(n, p);
   }
 
   private static double[][] getCorrelation_(int n, double[] p) {
      double corr[][] = new double[p.length][p.length];
 
      for (int i = 0; i < p.length; i++) {
         for (int j = 0; j < p.length; j++)
            corr[i][j] = -Math.sqrt(p[i] * p[j] / ((1.0 - p[i]) * (1.0 - p[j])));
         corr[i][i] = 1.0;
      }
      return corr;
   }

   public static double[][] getCorrelation(int n, double[] p) {
      verifParam(n, p);
      return getCorrelation_(n, p);
   }

   public static double[] getMLE(int x[][], int m, int d, int n) {
      double parameters[] = new double[d];
      double xBar[] = new double[d];
      double N = 0.0;
 
      if (m <= 0)
         throw new IllegalArgumentException("m <= 0");
      if (d <= 0)
         throw new IllegalArgumentException("d <= 0");
 
      for (int i = 0; i < d; i++)
         xBar[i] = 0;
 
      for (int v = 0; v < m; v++)
         for (int c = 0; c < d; c++)
            xBar[c] += x[v][c];
 
      for (int i = 0; i < d; i++) {
         xBar[i] = xBar[i] / (double) n;
         N += xBar[i];
      }
      if (N != (double) n)
         throw new IllegalArgumentException("n is not correct");
 
      for (int i = 0; i < d; i++)
         parameters[i] = xBar[i] / (double) n;
 
      return parameters;
   }

   public int getN() {
      return n;
   }

   public double[] getP() {
      return p;
   }

   public void setParams(int n, double p[]) {
      double sumP = 0.0;
 
      if (n <= 0)
         throw new IllegalArgumentException("n <= 0");
      if (p.length < 2)
         throw new IllegalArgumentException("p.length < 2");
 
      this.n = n;
      this.dimension = p.length;
      this.p = new double[dimension];
      for (int i = 0; i < dimension; i++) {
         if ((p[i] < 0) || (p[i] > 1))
            throw new IllegalArgumentException("p is not a probability vector");
 
         this.p[i] = p[i];
         sumP += p[i];
      }
 
      if (Math.abs(sumP - 1.0) > 1.0e-15)
         throw new IllegalArgumentException("p is not a probability vector");
   }
 
}