DirichletGen.java

/*
 * Class:        DirichletGen
 * Description:  multivariate generator for a Dirichlet 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.randvarmulti;
 
import umontreal.ssj.probdist.GammaDist;
import umontreal.ssj.randvar.RandomVariateGen;
import umontreal.ssj.randvar.GammaAcceptanceRejectionGen;
import umontreal.ssj.rng.RandomStream;

public class DirichletGen extends RandomMultivariateGen {
   private GammaAcceptanceRejectionGen[] ggens;
   private double[] alphas;
 
   // use log gamma if all alphas are smaller than this threshold
   private static final double ALPHA_THRESHOLD = 0.1;
 
   // determine if the Dirichlet multivariate must be generated using the
   // log gamma, when the shape parameter alphas are all small.
   private boolean useLogGamma = false;

   public DirichletGen(RandomStream stream, double[] alphas) {
      if (stream == null)
         throw new NullPointerException("stream is null");
      this.stream = stream;
      dimension = alphas.length;
      ggens = new GammaAcceptanceRejectionGen[alphas.length];
      for (int k = 0; k < alphas.length; k++)
         ggens[k] = new GammaAcceptanceRejectionGen(stream, new GammaDist(alphas[k], 1.0 / 2.0));
 
      this.alphas = new double[alphas.length];
      System.arraycopy(alphas, 0, this.alphas, 0, alphas.length);
      useLogGamma = canUseLogGamma(alphas);
   }

   public double getAlpha(int i) {
      return alphas[i];
   }

   public static void nextPoint(RandomStream stream, double[] alphas, double[] p) {
      if (canUseLogGamma(alphas)) {
         nextPointUsingLog(stream, alphas, p);
         return;
      }
 
      double total = 0;
      for (int i = 0; i < alphas.length; i++) {
         p[i] = GammaAcceptanceRejectionGen.nextDouble(stream, stream, alphas[i], 1.0 / 2.0);
         total += p[i];
      }
      for (int i = 0; i < alphas.length; i++)
         p[i] /= total;
   }

   public void nextPoint(double[] p) {
      if (useLogGamma) {
         nextPointUsingLog(stream, alphas, p);
         return;
      }
 
      int n = ggens.length;
      double total = 0;
      for (int i = 0; i < n; i++) {
         p[i] = ggens[i].nextDouble();
         total += p[i];
      }
      for (int i = 0; i < n; i++)
         p[i] /= total;
   }

   private static void nextPointUsingLog(RandomStream stream, double[] alphas, double[] p) {
      double total = 0;
      double[] log = new double[p.length]; // contains the log value
 
      // get the log value of the gammas
      for (int i = 0; i < alphas.length; i++) {
         log[i] = GammaAcceptanceRejectionGen.nextDoubleLog(stream, stream, alphas[i], 1.0 / 2.0);
      }
 
      // find the ratio for each p
      for (int i = 0; i < alphas.length; i++) {
         total = 0;
         for (int j = 0; j < alphas.length; j++)
            total += Math.exp(log[j] - log[i]);
         p[i] = 1.0 / total;
      }
   }

   private static boolean canUseLogGamma(double[] alphas) {
      for (int i = 0; i < alphas.length; i++)
         if (alphas[i] >= ALPHA_THRESHOLD)
            return false;
      return true;
   }
 
}