PascalDist.java

/*
 * Class:        PascalDist
 * Description:  Pascal 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.probdist;
 
import umontreal.ssj.util.RootFinder;
import umontreal.ssj.functions.MathFunction;

public class PascalDist extends NegativeBinomialDist {
   private static final double EPSI = 1.0E-10;
 
   private static class Function implements MathFunction {
      protected int m;
      protected int max;
      protected double mean;
      protected int[] Fj;
 
      public Function(int m, int max, double mean, int[] Fj) {
         this.m = m;
         this.max = max;
         this.mean = mean;
         this.Fj = new int[Fj.length];
         System.arraycopy(Fj, 0, this.Fj, 0, Fj.length);
      }
 
      public double evaluate(double p) {
         double sum = 0.0;
         double s = (p * mean) / (1.0 - p);
 
         for (int j = 0; j < max; j++)
            sum += Fj[j] / (s + (double) j);
 
         return sum + m * Math.log(p);
      }
 
      public double evaluateN(int n, double p) {
         double sum = 0.0;
 
         for (int j = 0; j < max; j++)
            sum += Fj[j] / (n + j);
 
         return sum + m * Math.log(p);
      }
   }

   public PascalDist(int n, double p) {
      setParams(n, p);
   }

   public static double[] getMLE(int[] x, int m) {
      if (m <= 0)
         throw new IllegalArgumentException("m <= 0");
 
      double sum = 0.0;
      int max = Integer.MIN_VALUE;
      for (int i = 0; i < m; i++) {
         sum += x[i];
         if (x[i] > max)
            max = x[i];
      }
      double mean = (double) sum / (double) m;
 
      double var = 0.0;
      for (int i = 0; i < m; i++)
         var += (x[i] - mean) * (x[i] - mean);
      var /= (double) m;
 
      if (mean >= var)
         throw new UnsupportedOperationException("mean >= variance");
 
      int[] Fj = new int[max];
      for (int j = 0; j < max; j++) {
         int prop = 0;
         for (int i = 0; i < m; i++)
            if (x[i] > j)
               prop++;
 
         Fj[j] = prop;
      }
 
      double[] parameters = new double[2];
      Function f = new Function(m, max, mean, Fj);
 
      parameters[1] = RootFinder.brentDekker(EPSI, 1 - EPSI, f, 1e-5);
      if (parameters[1] >= 1.0)
         parameters[1] = 1.0 - 1e-15;
 
      parameters[0] = Math.round((parameters[1] * mean) / (1.0 - parameters[1]));
      if (parameters[0] == 0)
         parameters[0] = 1;
 
      return parameters;
   }

   public static PascalDist getInstanceFromMLE(int[] x, int m) {
      double parameters[] = getMLE(x, m);
      return new PascalDist((int) parameters[0], parameters[1]);
   }

   public int getN1() {
      return (int) (n + 0.5);
   }

   public void setParams(int n, double p) {
      super.setParams((double) n, p);
   }

   public String toString() {
      return getClass().getSimpleName() + " : n = " + getN1() + ", p = " + getP();
   }
 
}