umontreal.iro.lecuyer.probdist.ParetoDist Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of ssj Show documentation
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: ParetoDist
* Description: Pareto 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 umontreal.iro.lecuyer.util.Num;
/**
* Extends the class {@link ContinuousDistribution} for a distribution
* from the Pareto family, with
* shape parameter
* α > 0 and location parameter β > 0.
* The density for this type of Pareto distribution is
*
*
*
* f (x) = αβα/xα+1 for x >= β,
*
* and 0 otherwise. The distribution function is
*
*
*
* F(x) = 1 - (β/x)α for x >= β,
*
* and the inverse distribution function is
*
*
*
* F-1(u) = β(1 - u)-1/α for 0 <= u < 1.
*
*
*/
public class ParetoDist extends ContinuousDistribution {
private double alpha;
private double beta;
/**
* Constructs a ParetoDist object with parameters α =
* alpha and β = 1.
*
*/
public ParetoDist (double alpha) {
setParams (alpha, 1.0);
}
/**
* Constructs a ParetoDist object with parameters α =
* alpha and β = beta.
*
*/
public ParetoDist (double alpha, double beta) {
setParams (alpha, beta);
}
public double density (double x) {
return density (alpha, beta, x);
}
public double cdf (double x) {
return cdf (alpha, beta, x);
}
public double barF (double x) {
return barF (alpha, beta, x);
}
public double inverseF (double u) {
return inverseF (alpha, beta, u);
}
public double getMean() {
return ParetoDist.getMean (alpha, beta);
}
public double getVariance() {
return ParetoDist.getVariance (alpha, beta);
}
public double getStandardDeviation() {
return ParetoDist.getStandardDeviation (alpha, beta);
}
/**
* Computes the density function.
*
*/
public static double density (double alpha, double beta, double x) {
if (alpha <= 0.0)
throw new IllegalArgumentException ("alpha <= 0");
if (beta <= 0.0)
throw new IllegalArgumentException ("beta <= 0");
return x < beta ? 0 : alpha*Math.pow (beta/x, alpha)/x;
}
/**
* Computes the distribution function.
*
*/
public static double cdf (double alpha, double beta, double x) {
if (alpha <= 0.0)
throw new IllegalArgumentException ("alpha <= 0");
if (beta <= 0.0)
throw new IllegalArgumentException ("beta <= 0");
if (x <= beta)
return 0.0;
return 1.0 - Math.pow (beta/x, alpha);
}
/**
* Computes the complementary distribution function.
*
*/
public static double barF (double alpha, double beta, double x) {
if (alpha <= 0)
throw new IllegalArgumentException ("c <= 0");
if (beta <= 0.0)
throw new IllegalArgumentException ("beta <= 0");
if (x <= beta)
return 1.0;
return Math.pow (beta/x, alpha);
}
/**
* Computes the inverse of the distribution function.
*
*/
public static double inverseF (double alpha, double beta, double u) {
if (alpha <= 0)
throw new IllegalArgumentException ("c <= 0");
if (beta <= 0.0)
throw new IllegalArgumentException ("beta <= 0");
if (u < 0.0 || u > 1.0)
throw new IllegalArgumentException ("u not in [0,1]");
if (u <= 0.0)
return beta;
double t;
t = -Math.log1p (-u);
if ((u >= 1.0) || t/Math.log(10) >= alpha * Num.DBL_MAX_10_EXP)
return Double.POSITIVE_INFINITY;
return beta / Math.pow (1 - u, 1.0/alpha);
}
/**
* Estimates the parameters
* (α, β) of the Pareto distribution
* using the maximum likelihood method, from the n observations
* x[i],
* i = 0, 1,…, n - 1. The estimates are returned in a two-element
* array, in regular order: [α, β].
*
* @param x the list of observations used to evaluate parameters
*
* @param n the number of observations used to evaluate parameters
*
* @return returns the parameters [
* hat(α),
* hat(β)]
*
*/
public static double[] getMLE (double[] x, int n) {
if (n <= 0)
throw new IllegalArgumentException ("n <= 0");
double [] parameters = new double[2];
parameters[1] = Double.POSITIVE_INFINITY;
for (int i = 0; i < n; i++) {
if (x[i] < parameters[1])
parameters[1] = x[i];
}
double sum = 0.0;
for (int i = 0; i < n; i++) {
if (x[i] > 0.0)
sum += Math.log (x[i] / parameters[1]);
else
sum -= 709.0;
}
parameters[0] = n / sum;
return parameters;
}
/**
* Creates a new instance of a Pareto distribution with parameters α and β
* estimated using the maximum likelihood method based on the n observations
* x[i],
* i = 0, 1,…, n - 1.
*
* @param x the list of observations to use to evaluate parameters
*
* @param n the number of observations to use to evaluate parameters
*
*
*/
public static ParetoDist getInstanceFromMLE (double[] x, int n) {
double parameters[] = getMLE (x, n);
return new ParetoDist (parameters[0], parameters[1]);
}
/**
* Computes and returns the mean
* E[X] = αβ/(α - 1)
* of the Pareto distribution with parameters α and β.
*
* @return the mean of the Pareto distribution
*
*/
public static double getMean (double alpha, double beta) {
if (alpha <= 1.0)
throw new IllegalArgumentException("alpha <= 1");
if (beta <= 0.0)
throw new IllegalArgumentException ("beta <= 0");
return ((alpha * beta) / (alpha - 1.0));
}
/**
* Computes and returns the variance
* of the Pareto distribution with parameters α and β.
*
* @return the variance of the Pareto distribution
*
* Var[X] = αβ2/[(α -2)(α - 1)]
*
*/
public static double getVariance (double alpha, double beta) {
if (alpha <= 2)
throw new IllegalArgumentException("alpha <= 2");
if (beta <= 0.0)
throw new IllegalArgumentException ("beta <= 0");
return ((alpha * beta * beta) / ((alpha - 2.0) * (alpha - 1.0) * (alpha - 1.0)));
}
/**
* Computes and returns the standard deviation of the Pareto
* distribution with parameters α and β.
*
* @return the standard deviation of the Pareto distribution
*
*/
public static double getStandardDeviation (double alpha, double beta) {
return Math.sqrt (ParetoDist.getVariance (alpha, beta));
}
/**
* Returns the parameter α.
*
*/
public double getAlpha() {
return alpha;
}
/**
* Returns the parameter β.
*
*/
public double getBeta() {
return beta;
}
/**
* Sets the parameter α and β for this object.
*
*/
public void setParams (double alpha, double beta) {
if (alpha <= 0.0)
throw new IllegalArgumentException ("alpha <= 0");
if (beta <= 0.0)
throw new IllegalArgumentException ("beta <= 0");
this.alpha = alpha;
this.beta = beta;
supportA = beta;
}
/**
* Return a table containing the parameters of the current distribution.
* This table is put in regular order: [α, β].
*
*
*/
public double[] getParams () {
double[] retour = {alpha, beta};
return retour;
}
public String toString () {
return getClass().getSimpleName() + " : alpha = " + alpha + ", beta = " + beta;
}
}