umontreal.iro.lecuyer.probdist.WeibullDist 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: WeibullDist
* Description: Weibull 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;
import umontreal.iro.lecuyer.util.RootFinder;
import umontreal.iro.lecuyer.functions.MathFunction;
/**
* This class extends the class {@link ContinuousDistribution} for
* the Weibull distribution with shape parameter
*
* α > 0, location parameter δ, and scale parameter
*
* λ > 0.
* The density function is
*
*
*
* f (x) = αλα(x - δ)α-1e-(λ(x-δ))α for x > δ.
*
* the distribution function is
*
*
*
* F(x) = 1 - e-(λ(x-δ))α for x > δ,
*
* and the inverse distribution function is
*
*
*
* F-1(u) = (- ln(1 - u))1/α/λ + δ for 0 <= u < 1.
*
*
*/
public class WeibullDist extends ContinuousDistribution {
private double alpha;
private double lambda;
private double delta;
private static class Function implements MathFunction {
private int n;
private double xi[];
private double lnXi[];
private double sumLnXi = 0.0;
private final double LN_EPS = Num.LN_DBL_MIN - Num.LN2;
public Function (double x[], int n)
{
this.n = n;
this.xi = new double[n];
this.lnXi = new double[n];
for (int i = 0; i < n; i++)
{
this.xi[i] = x[i];
if (x[i] > 0.0)
this.lnXi[i] = Math.log (x[i]);
else
this.lnXi[i] = LN_EPS;
sumLnXi += this.lnXi[i];
}
}
public double evaluate (double x)
{
if (x <= 0.0) return 1.0e200;
double sumXiLnXi = 0.0;
double sumXi = 0.0;
double xalpha;
for (int i = 0; i < n; i++)
{
xalpha = Math.pow (this.xi[i], x);
sumXiLnXi += xalpha * lnXi[i];
sumXi += xalpha;
}
return (x * (n * sumXiLnXi - sumLnXi * sumXi) - n * sumXi);
}
}
/**
* Constructs a WeibullDist object with parameters
* α = alpha, λ = 1, and δ = 0.
*
*/
public WeibullDist (double alpha) {
setParams (alpha, 1.0, 0.0);
}
/**
* Constructs a WeibullDist object with parameters
* α = alpha,
* λ = lambda, and δ = delta.
*
*/
public WeibullDist (double alpha, double lambda, double delta) {
setParams (alpha, lambda, delta);
}
public double density (double x) {
return density (alpha, lambda, delta, x);
}
public double cdf (double x) {
return cdf (alpha, lambda, delta, x);
}
public double barF (double x) {
return barF (alpha, lambda, delta, x);
}
public double inverseF (double u) {
return inverseF (alpha, lambda, delta, u);
}
public double getMean() {
return WeibullDist.getMean (alpha, lambda, delta);
}
public double getVariance() {
return WeibullDist.getVariance (alpha, lambda, delta);
}
public double getStandardDeviation() {
return WeibullDist.getStandardDeviation (alpha, lambda, delta);
}
/**
* Computes the density function.
*
*/
public static double density (double alpha, double lambda,
double delta, double x) {
if (alpha <= 0.0)
throw new IllegalArgumentException ("alpha <= 0");
if (lambda <= 0.0)
throw new IllegalArgumentException ("lambda <= 0");
if (x <= delta)
return 0.0;
double y = Math.log(lambda*(x - delta)) * alpha;
if (y >= 7.0)
return 0.0;
y = Math.exp(y);
return alpha * (y / (x - delta)) * Math.exp(-y);
}
/**
* Same as density (alpha, 1, 0, x).
*
*/
public static double density (double alpha, double x) {
return density (alpha, 1.0, 0.0, x);
}
/**
* Computes the distribution function.
*
*/
public static double cdf (double alpha, double lambda,
double delta, double x) {
if (alpha <= 0.0)
throw new IllegalArgumentException ("alpha <= 0");
if (lambda <= 0.0)
throw new IllegalArgumentException ("lambda <= 0");
if (x <= delta)
return 0.0;
if ((lambda*(x - delta) >= XBIG) && (alpha >= 1.0))
return 1.0;
double y = Math.log(lambda*(x - delta)) * alpha;
if (y >= 3.65)
return 1.0;
y = Math.exp(y);
return -Math.expm1 (-y); // in JDK-1.5
}
/**
* Same as cdf (alpha, 1, 0, x).
*
*/
public static double cdf (double alpha, double x) {
return cdf (alpha, 1.0, 0.0, x);
}
/**
* Computes the complementary distribution function.
*
*/
public static double barF (double alpha, double lambda,
double delta, double x) {
if (alpha <= 0)
throw new IllegalArgumentException ("alpha <= 0");
if (lambda <= 0)
throw new IllegalArgumentException ("lambda <= 0");
if (x <= delta)
return 1.0;
if (alpha >= 1.0 && x >= Num.DBL_MAX_EXP*2)
return 0.0;
double temp = Math.log (lambda*(x - delta)) * alpha;
if (temp >= Num.DBL_MAX_EXP * Num.LN2)
return 0.0;
temp = Math.exp(temp);
return Math.exp (-temp);
}
/**
* Same as barF (alpha, 1, 0, x).
*
*/
public static double barF (double alpha, double x) {
return barF (alpha, 1.0, 0.0, x);
}
/**
* Computes the inverse of the distribution function.
*
*/
public static double inverseF (double alpha, double lambda,
double delta, double u) {
double t;
if (alpha <= 0.0)
throw new IllegalArgumentException ("alpha <= 0");
if (lambda <= 0.0)
throw new IllegalArgumentException ("lambda <= 0");
if (u < 0.0 || u > 1.0)
throw new IllegalArgumentException ("u not in [0, 1]");
if (u <= 0.0)
return 0.0;
if (u >= 1.0)
return Double.POSITIVE_INFINITY;
t = -Math.log1p (-u);
if (Math.log (t)/Math.log (10) >= alpha*Num.DBL_MAX_10_EXP)
throw new ArithmeticException
("inverse function cannot be positive infinity");
return Math.pow (t, 1.0/alpha)/lambda + delta;
}
/**
* Same as inverseF (alpha, 1, 0, x).
*
*/
public static double inverseF (double alpha, double x) {
return inverseF (alpha, 1.0, 0.0, x);
}
private static double[] getMaximumLikelihoodEstimate (double[] x, int n,
double delta) {
if (n <= 0)
throw new IllegalArgumentException ("n <= 0");
if (delta != 0.0)
throw new IllegalArgumentException ("delta must be equal to 0");
// Verifier cette fonction si delta != 0.
final double LN_EPS = Num.LN_DBL_MIN - Num.LN2;
double sumLn = 0.0;
double sumLn2 = 0.0;
double lnxi;
for (int i = 0; i < n; i++) {
if (x[i] <= delta)
lnxi = LN_EPS;
else
lnxi = Math.log (x[i]);
sumLn += lnxi;
sumLn2 += lnxi * lnxi;
}
double alpha0 = Math.sqrt ((double) n / ((6.0 / (Math.PI * Math.PI)) *
(sumLn2 - sumLn * sumLn / (double) n)));
double a = alpha0 - 20.0;
if (a <= delta)
a = delta + 1.0e-5;
double param[] = new double[3];
param[2] = 0.0;
Function f = new Function (x, n);
param[0] = RootFinder.brentDekker (a, alpha0 + 20.0, f, 1e-5);
double sumXalpha = 0.0;
for (int i = 0; i < n; i++)
sumXalpha += Math.pow (x[i], param[0]);
param[1] = Math.pow ((double) n / sumXalpha, 1.0 / param[0]);
return param;
}
/**
* Estimates the parameters
* (α, λ) of the Weibull distribution,
* assuming that
* δ = 0,
* 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 to use to evaluate parameters
*
* @param n the number of observations to use to evaluate parameters
*
* @return returns the parameter [
* hat(α),
* hat(λ),
* hat(δ) = 0]
*
*/
public static double[] getMLE (double[] x, int n)
{
return getMaximumLikelihoodEstimate (x, n, 0.0);
}
/**
* Creates a new instance of a Weibull distribution with parameters α,
* λ and
* δ = 0
* 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 WeibullDist getInstanceFromMLE (double[] x, int n) {
double param[] = getMLE (x, n);
return new WeibullDist (param[0], param[1], param[2]);
}
/**
* Computes and returns the mean
* of the Weibull distribution with parameters α, λ and δ.
*
* @return the mean of the Weibull distribution
*
* E[X] = δ + Γ(1 + 1/α)/λ
*
*/
public static double getMean (double alpha, double lambda, double delta) {
if (alpha <= 0.0)
throw new IllegalArgumentException ("alpha <= 0");
if (lambda <= 0.0)
throw new IllegalArgumentException ("lambda <= 0");
return (delta + Math.exp (Num.lnGamma(1.0 + 1.0 / alpha)) / lambda);
}
/**
* Computes and returns the variance
* of the Weibull distribution with parameters α, λ and δ.
*
* @return the variance of the Weibull distribution
*
* Var[X] = 1/λ2| Γ(2/α +1) - Γ2(1/α + 1)|
*
*/
public static double getVariance (double alpha, double lambda,
double delta) {
double gAlpha;
if (alpha <= 0.0)
throw new IllegalArgumentException ("alpha <= 0");
if (lambda <= 0.0)
throw new IllegalArgumentException ("lambda <= 0");
gAlpha = Math.exp (Num.lnGamma (1.0 / alpha + 1.0));
return (Math.abs (Math.exp (Num.lnGamma(2 / alpha + 1)) - gAlpha * gAlpha) / (lambda * lambda));
}
/**
* Computes and returns the standard deviation
* of the Weibull distribution with parameters α, λ and δ.
*
* @return the standard deviation of the Weibull distribution
*
*/
public static double getStandardDeviation (double alpha, double lambda,
double delta) {
return Math.sqrt (WeibullDist.getVariance (alpha, lambda, delta));
}
/**
* Returns the parameter α.
*
*/
public double getAlpha() {
return alpha;
}
/**
* Returns the parameter λ.
*
*/
public double getLambda() {
return lambda;
}
/**
* Returns the parameter δ.
*
*/
public double getDelta() {
return delta;
}
/**
* Sets the parameters α, λ and δ for this
* object.
*
*/
public void setParams (double alpha, double lambda, double delta) {
if (alpha <= 0.0)
throw new IllegalArgumentException ("alpha <= 0");
if (lambda <= 0.0)
throw new IllegalArgumentException ("lambda <= 0");
this.alpha = alpha;
this.lambda = lambda;
this.delta = delta;
supportA = delta;
}
/**
* Return a table containing the parameters of the current distribution.
* This table is put in regular order: [α, λ, δ].
*
*
*/
public double[] getParams () {
double[] retour = {alpha, lambda, delta};
return retour;
}
public String toString () {
return getClass().getSimpleName() + " : alpha = " + alpha + ", lambda = " + lambda + ", delta = " + delta;
}
}
© 2015 - 2025 Weber Informatics LLC | Privacy Policy