umontreal.iro.lecuyer.probdist.TriangularDist Maven / Gradle / Ivy

Go to download
Show more of this group Show more artifacts with this name
Show all versions of ssj Show documentation
SSJ is a Java library for stochastic simulation, developed under the direction of Pierre L'Ecuyer, in the Département d'Informatique et de Recherche Opérationnelle (DIRO), at the Université de Montréal. It provides facilities for generating uniform and nonuniform random variates, computing different measures related to probability distributions, performing goodness-of-fit tests, applying quasi-Monte Carlo methods, collecting (elementary) statistics, and programming discrete-event simulations with both events and processes.
The newest version!


/*
 * Class:        TriangularDist
 * Description:  triangular distribution
 * Environment:  Java
 * Software:     SSJ 
 * Copyright (C) 2001  Pierre L'Ecuyer and Université de Montréal
 * Organization: DIRO, Université de Montréal
 * @author       
 * @since

 * 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
   GPL licence site.
 */

package umontreal.iro.lecuyer.probdist;

import java.util.Arrays;



/**
 * Extends the class {@link ContinuousDistribution} for
 * the triangular distribution with domain [a, b] and mode 
 * (or shape parameter) m, where  
 * a <= m <= b. 
 * The density function is
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * f (x) = 2(x - a)/[(b - a)(m - a)]          for a <= x <= m,
f (x) = 2(b - x)/[(b - a)(b - m)]          for m <= x <= b,
f (x) = 0          elsewhere,
 * 

 * the distribution function is
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * F(x) = 0          for x < a,
F(x) = (x - a)²/[(b - a)(m - a)]          if a <= x <= m,
F(x) = 1 - (b - x)²/[(b - a)(b - m)]          if m <= x <= b,
F(x) = 1          for x > b,
 * 

 * and the inverse distribution function is given by
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * F^-1(u) = a + ((b - a)(m - a)u)^1/2          if 0 <= u <= (m - a)/(b - a),
F^-1(u) = b - ((b - a)(b - m)(1 - u))^1/2          if (m - a)/(b - a <= u <= 1.
 * 

 * 
 */
public class TriangularDist extends ContinuousDistribution {
   private double a;
   private double b;
   private double m;



   /**
    * Constructs a TriangularDist object with default parameters
    *    a = 0, b = 1, and m = 0.5.
    * 
    */
   public TriangularDist() {
      setParams (0.0, 1.0, 0.5);
   }


   /**
    * Constructs a TriangularDist object with parameters a = 0 ,
    *  b = 1 and m = m.
    * 
    */
   public TriangularDist (double m) {
      setParams (0.0, 1.0, m);
   }


   /**
    * Constructs a TriangularDist object with
    *    parameters a, b and m.
    * 
    */
   public TriangularDist (double a, double b, double m) {
      setParams (a, b, m);
   }


   public double density (double x) {
      return density (a, b, m, x);
   }
       
   public double cdf (double x) {
      return cdf (a, b, m, x);
   }

   public double barF (double x) {
      return barF (a, b, m, x);
   }
     
   public double inverseF (double u){
      return inverseF (a, b, m, u);
   }

   public double getMean() {
      return TriangularDist.getMean (a, b, m);
   }

   public double getVariance() {
      return TriangularDist.getVariance (a, b, m);
   }

   public double getStandardDeviation() {
      return TriangularDist.getStandardDeviation (a, b, m);
   }

   /**
    * Computes the density function.
    * 
    */
   public static double density (double a, double b, double m, double x) {
      if (m < a || m > b)
         throw new IllegalArgumentException ("m is not in [a,b]");
      if (x < a || x > b)
         return 0.0;
      else if (x <= m && m != a)
         return 2.0*(x - a)/((b - a)*(m - a));
      else
         return 2.0*(b - x)/((b - a)*(b - m));
   }
      

   /**
    * Computes the  distribution function.
    * 
    */
   public static double cdf (double a, double b, double m, double x) {
      if (m < a || m > b)
         throw new IllegalArgumentException ("m is not in [a,b]");
      if (x <= a)
         return 0.0;
      else if (x <= m && m != a)
         return (x - a)*(x - a)/((b - a)*(m - a));
      else if (x < b)
         return 1.0 - (b - x)*(b - x)/((b - a)*(b - m));
      else
         return 1.0;
   }


   /**
    * Computes the complementary distribution function.
    * 
    */
   public static double barF (double a, double b, double m, double x) {
      if (m < a || m > b)
         throw new IllegalArgumentException ("m is not in [a,b]");
      if (x <= a)
         return 1.0;
      else if (x <= m && m != a)
         return 1.0 - (x - a)*(x - a)/((b - a)*(m - a));
      else if (x < b)
         return (b - x)*(b - x)/((b - a)*(b - m));
      else
         return 0.0;
   }


   /**
    * Computes the inverse distribution function.
    * 
    */
   public static double inverseF (double a, double b, double m, double u) {
      if (m < a || m > b)
         throw new IllegalArgumentException ("m is not in [a,b]");
      if (u < 0.0 || u > 1.0)
         throw new IllegalArgumentException ("u is not in [0,1]");
      if (u <= 0.0)
         return a;
      if (u >= 1.0)
         return b;
       // the code is taken and adapted from unuran
       // file /distributions/c_triangular_gen.c
       double h = (m - a)/(b - a);
       return u <= h && m != a ? a + Math.sqrt ((b - a)*(m - a)*u) 
                : b - Math.sqrt ((b - a)*(b - m)*(1 - u));
   }


   /**
    * Estimates the parameter m of the triangular distribution using the
    *    maximum likelihood method, from the n observations x[i],
    *    
    * i = 0, 1,…, n - 1. The estimated parameter is returned in a one-element
    *    array: [hat(m)].  See.
    * 
    * @param x the list of observations used to evaluate parameters
    * 
    *    @param n the number of observations used to evaluate parameters
    * 
    *    @param a lower limit of range
    * 
    *    @param b upper limit of range
    * 
    *    @return returns the parameter [m]
    * 
    */
   public static double[] getMLE (double[] x, int n, double a, double b) {
      if (n <= 0)
         throw new IllegalArgumentException ("n <= 0");
      double[] Y = new double[n];   // sorted x[i]
      System.arraycopy (x, 0, Y, 0, n);
      Arrays.sort (Y);

      int rmax = -1;
      double prodmax = -1.0e300;
      final double ba = b - a;
      double z;
      int i;
      for (int r = 0; r < n; r++) {
         z = (Y[r] - a) / ba;
         if ((z <= (double)r/n) || (z >= (double)(r + 1)/n))
            continue;    // MLE cannot be there
         double prod = 1.0;
         double d = Y[r] - a;
         for (i = 0; i < r; i++)
            prod *= (Y[i] - a)/d;
         
         d = b - Y[r];
         for (i = r+1; i < n; i++)
            prod *= (b - Y[i])/d;

         if (prod > prodmax) {
            prodmax = prod;
            rmax = r;
         }
      }
      
      if (rmax < 0)
         throw new UnsupportedOperationException (
            "   data cannot fit a triangular distribution");

      double[] param = new double[1];
      param[0] = Y[rmax];
      return param;
   }


   /**
    * Creates a new instance of a triangular distribution with parameters
    *   a and b.  m is estimated using the maximum
    *   likelihood method based on the n observations
    *   x[i], 
    * i = 0, 1,…, n - 1.
    * 
    * @param x the list of observations used to evaluate parameters
    * 
    *    @param n the number of observations used to evaluate parameters
    * 
    *    @param a lower limit of range
    * 
    *    @param b upper limit of range
    * 
    * 
    */
   public static TriangularDist getInstanceFromMLE (double[] x, int n,
                                                    double a, double b) {
      double param[] = getMLE (x, n, a, b);
      return new TriangularDist (a, b, param[0]);
   }


   /**
    * Computes and returns the mean 
    * E[X] = (a + b + m)/3
    *    of the triangular distribution with parameters a, b, m.
    * 
    * @return the mean of the triangular distribution
    * 
    */
   public static double getMean (double a, double b, double m) {
      if ((a == 0.0 && b == 1.0) && (m < 0 || m > 1))
         throw new IllegalArgumentException ("m is not in [0,1]");
      else if (m < a || m > b) 
         throw new IllegalArgumentException ("m is not in [a,b]");

      return ((a + b + m) / 3.0);
   }


   /**
    * Computes and returns the variance 
    *    
    * Var[X] = (a² + b² + m² - ab - am - bm)/18
    *    of the triangular distribution with parameters a, b, m.
    * 
    * @return the variance of the triangular distribution
    * 
    */
   public static double getVariance (double a, double b, double m) {
      if ((a == 0.0 && b == 1.0) && (m < 0 || m > 1))
         throw new IllegalArgumentException ("m is not in [0,1]");
      else if (m < a || m > b) 
         throw new IllegalArgumentException ("m is not in [a,b]");

      return ((a * a + b * b + m * m - a * b - a * m - b * m) / 18.0);
   }


   /**
    * Computes and returns the standard deviation
    *    of the triangular distribution with parameters a, b, m.
    * 
    * @return the standard deviation of the triangular distribution
    * 
    */
   public static double getStandardDeviation (double a, double b, double m) {
      return Math.sqrt (TriangularDist.getVariance (a, b, m));
   }
      

   /**
    * Returns the value of a for this object.
    * 
    */
   public double getA() {
      return a;
   }
      

   /**
    * Returns the value of b for this object.
    * 
    */
   public double getB() {
      return b;
   }
      

   /**
    * Returns the value of m for this object.
    * 
    */
   public double getM() {
      return m;
   }
      

   /**
    * Sets the value of the parameters a, b and m for this object.
    * 
    */
   public void setParams (double a, double b, double m) {
      if ((a == 0.0 && b == 1.0) && (m < 0 || m > 1))
         throw new IllegalArgumentException ("m is not in [0,1]");
      else if (a >= b)
         throw new IllegalArgumentException ("a >= b");
      else if (m < a || m > b) 
         throw new IllegalArgumentException ("m is not in [a,b]");
      this.a = a;
      this.b = b;
      this.m = m;
      supportA = a;
      supportB = b;
   }


   /**
    * Return a table containing the parameters of the current distribution.
    *    This table is put in regular order: [a, b, m].
    * 
    */
   public double[] getParams () {
      double[] retour = {a, b, m};
      return retour;
   }


   /**
    * Returns a String containing information about the current distribution.
    * 
    */
   public String toString () {
      return getClass().getSimpleName() + " : a = " + a + ", b = " + b + ", m = " + m;
   }

}
f (x) =	2(x - a)/[(b - a)(m - a)]	for a <= x <= m,
f (x) =	2(b - x)/[(b - a)(b - m)]	for m <= x <= b,
f (x) =	0	elsewhere,
F(x) =	0	for x < a,
F(x) =	(x - a)²/[(b - a)(m - a)]	if a <= x <= m,
F(x) =	1 - (b - x)²/[(b - a)(b - m)]	if m <= x <= b,
F(x) =	1	for x > b,
F^-1(u) =	a + ((b - a)(m - a)u)^1/2	if 0 <= u <= (m - a)/(b - a),
F^-1(u) =	b - ((b - a)(b - m)(1 - u))^1/2	if (m - a)/(b - a <= u <= 1.