/* *
* Title:
* Description:
* Copyright: Copyright (c)
* Company:
* @author
* @version 1.0
*/
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;
}
}
Ads help maintain this website.