Pearson6Dist.java

/*
 * Class:        Pearson6Dist
 * Description:  Pearson type VI 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.Num;
import optimization.*;

public class Pearson6Dist extends ContinuousDistribution {
   protected double alpha1;
   protected double alpha2;
   protected double beta;
   protected double logBeta; // Ln (Beta (alpha1, alpha2))
 
   private static class Optim implements Uncmin_methods {
      private int n;
      private double[] x;
 
      public Optim(double[] x, int n) {
         this.n = n;
         this.x = new double[n];
         System.arraycopy(x, 0, this.x, 0, n);
      }
 
      public double f_to_minimize(double[] param) {
 
         if ((param[1] <= 0.0) || (param[2] <= 0.0) || (param[3] <= 0.0))
            return 1e200;
 
         double sumLogY = 0.0;
         double sumLog1_Y = 0.0;
         for (int i = 0; i < n; i++) {
            if (x[i] > 0.0)
               sumLogY += Math.log(x[i] / param[3]);
            else
               sumLogY -= 709.0;
            sumLog1_Y += Math.log1p(x[i] / param[3]);
         }
 
         return (n * (Math.log(param[3]) + Num.lnBeta(param[1], param[2])) - (param[1] - 1.0) * sumLogY
               + (param[1] + param[2]) * sumLog1_Y);
      }
 
      public void gradient(double[] x, double[] g) {
      }
 
      public void hessian(double[] x, double[][] h) {
      }
   }

   public Pearson6Dist(double alpha1, double alpha2, double beta) {
      setParam(alpha1, alpha2, beta);
   }
 
   public double density(double x) {
      if (x <= 0.0)
         return 0.0;
      return Math.exp((alpha1 - 1.0) * Math.log(x / beta) - (logBeta + (alpha1 + alpha2) * Math.log1p(x / beta)))
            / beta;
   }
 
   public double cdf(double x) {
      return cdf(alpha1, alpha2, beta, x);
   }
 
   public double barF(double x) {
      return barF(alpha1, alpha2, beta, x);
   }
 
   public double inverseF(double u) {
      return inverseF(alpha1, alpha2, beta, u);
   }
 
   public double getMean() {
      return getMean(alpha1, alpha2, beta);
   }
 
   public double getVariance() {
      return getVariance(alpha1, alpha2, beta);
   }
 
   public double getStandardDeviation() {
      return getStandardDeviation(alpha1, alpha2, beta);
   }

   public static double density(double alpha1, double alpha2, double beta, double x) {
      if (alpha1 <= 0.0)
         throw new IllegalArgumentException("alpha1 <= 0");
      if (alpha2 <= 0.0)
         throw new IllegalArgumentException("alpha2 <= 0");
      if (beta <= 0.0)
         throw new IllegalArgumentException("beta <= 0");
      if (x <= 0.0)
         return 0.0;
 
      return Math.exp((alpha1 - 1.0) * Math.log(x / beta)
            - (Num.lnBeta(alpha1, alpha2) + (alpha1 + alpha2) * Math.log1p(x / beta))) / beta;
   }

   public static double cdf(double alpha1, double alpha2, double beta, double x) {
      if (alpha1 <= 0.0)
         throw new IllegalArgumentException("alpha1 <= 0");
      if (alpha2 <= 0.0)
         throw new IllegalArgumentException("alpha2 <= 0");
      if (beta <= 0.0)
         throw new IllegalArgumentException("beta <= 0");
      if (x <= 0.0)
         return 0.0;
 
      return BetaDist.cdf(alpha1, alpha2, x / (x + beta));
   }

   public static double barF(double alpha1, double alpha2, double beta, double x) {
      if (alpha1 <= 0.0)
         throw new IllegalArgumentException("alpha1 <= 0");
      if (alpha2 <= 0.0)
         throw new IllegalArgumentException("alpha2 <= 0");
      if (beta <= 0.0)
         throw new IllegalArgumentException("beta <= 0");
      if (x <= 0.0)
         return 1.0;
 
      return BetaDist.barF(alpha1, alpha2, x / (x + beta));
   }

   public static double inverseF(double alpha1, double alpha2, double beta, double u) {
      if (alpha1 <= 0.0)
         throw new IllegalArgumentException("alpha1 <= 0");
      if (alpha2 <= 0.0)
         throw new IllegalArgumentException("alpha2 <= 0");
      if (beta <= 0.0)
         throw new IllegalArgumentException("beta <= 0");
 
      double y = BetaDist.inverseF(alpha1, alpha2, u);
 
      return ((y * beta) / (1.0 - y));
   }

   public static double[] getMLE(double[] x, int n) {
      if (n <= 0)
         throw new IllegalArgumentException("n <= 0");
 
      double[] parameters = new double[3];
      double[] xpls = new double[4];
      double[] param = new double[4];
      double[] fpls = new double[4];
      double[] gpls = new double[4];
      int[] itrcmd = new int[2];
      double[][] h = new double[4][4];
      double[] udiag = new double[4];
 
      Optim system = new Optim(x, n);
 
      double mean = 0.0;
      double mean2 = 0.0;
      double mean3 = 0.0;
      for (int i = 0; i < n; i++) {
         mean += x[i];
         mean2 += x[i] * x[i];
         mean3 += x[i] * x[i] * x[i];
      }
      mean /= (double) n;
      mean2 /= (double) n;
      mean3 /= (double) n;
 
      double r1 = mean2 / (mean * mean);
      double r2 = mean2 * mean / mean3;
 
      param[1] = -(2.0 * (-1.0 + r1 * r2)) / (-2.0 + r1 + r1 * r2);
      if (param[1] <= 0)
         param[1] = 1;
      param[2] = (-3.0 - r2 + 4.0 * r1 * r2) / (-1.0 - r2 + 2.0 * r1 * r2);
      if (param[2] <= 0)
         param[2] = 1;
      param[3] = (param[2] - 1.0) * mean / param[1];
      if (param[3] <= 0)
         param[3] = 1;
 
      Uncmin_f77.optif0_f77(3, param, system, xpls, fpls, gpls, itrcmd, h, udiag);
 
      for (int i = 0; i < 3; i++)
         parameters[i] = xpls[i + 1];
 
      return parameters;
   }

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

   public static double getMean(double alpha1, double alpha2, double beta) {
      if (alpha1 <= 0.0)
         throw new IllegalArgumentException("alpha1 <= 0");
      if (alpha2 <= 1.0)
         throw new IllegalArgumentException("alpha2 <= 1");
      if (beta <= 0.0)
         throw new IllegalArgumentException("beta <= 0");
 
      return ((beta * alpha1) / (alpha2 - 1.0));
   }

   public static double getVariance(double alpha1, double alpha2, double beta) {
      if (alpha1 <= 0.0)
         throw new IllegalArgumentException("alpha1 <= 0");
      if (alpha2 <= 0.0)
         throw new IllegalArgumentException("alpha2 <= 2");
      if (beta <= 0.0)
         throw new IllegalArgumentException("beta <= 0");
 
      return (((beta * beta) * alpha1 * (alpha1 + alpha2 - 1.0)) / ((alpha2 - 1.0) * (alpha2 - 1.0) * (alpha2 - 2.0)));
   }

   public static double getStandardDeviation(double alpha1, double alpha2, double beta) {
      return Math.sqrt(getVariance(alpha1, alpha2, beta));
   }

   public double getAlpha1() {
      return alpha1;
   }

   public double getAlpha2() {
      return alpha2;
   }

   public double getBeta() {
      return beta;
   }

   public void setParam(double alpha1, double alpha2, double beta) {
      if (alpha1 <= 0.0)
         throw new IllegalArgumentException("alpha1 <= 0");
      if (alpha2 <= 0.0)
         throw new IllegalArgumentException("alpha2 <= 0");
      if (beta <= 0.0)
         throw new IllegalArgumentException("beta <= 0");
      supportA = 0.0;
      this.alpha1 = alpha1;
      this.alpha2 = alpha2;
      this.beta = beta;
      logBeta = Num.lnBeta(alpha1, alpha2);
   }

   public double[] getParams() {
      double[] retour = { alpha1, alpha2, beta };
      return retour;
   }

   public String toString() {
      return getClass().getSimpleName() + " : alpha1 = " + alpha1 + ", alpha2 = " + alpha2 + ", beta = " + beta;
   }
 
}