umontreal.iro.lecuyer.probdist.FrechetDist 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: FrechetDist
* Description: Fréchet 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.*;
import optimization.*;
import umontreal.iro.lecuyer.functions.MathFunction;
/**
* Extends the class {@link ContinuousDistribution} for the Fréchet
* distribution, with location parameter δ, scale
* parameter β > 0, and shape parameter
* α > 0, where we use
* the notation
* z = (x - δ)/β. It has density
*
*
*
* f (x) = αe-z-α/(βzα+1), for x > δ
*
* and distribution function
*
*
*
* F(x) = e-z-α, for x > δ.
*
* Both the density and the distribution are 0 for
* x <= δ.
*
*
* The mean is given by
*
*
*
* E[X] = δ + βΓ(1 - 1/α),
*
* where Γ(x) is the gamma function. The variance is
*
*
*
* Var[X] = β2[Γ(1 - 2/α) - (Γ(1 - 1/α))2].
*
*
*/
public class FrechetDist extends ContinuousDistribution {
private double delta;
private double beta;
private double alpha;
private static class Optim implements Lmder_fcn {
protected double[] x;
protected int n;
public Optim (double[] x, int n) {
this.n = n;
this.x = x;
}
public void fcn (int m, int n, double[] par, double[] fvec, double[][] fjac, int iflag[])
{
if (par[1] <= 0.0 || par[2] <= 0.0) {
final double BIG = 1.0e100;
fvec[1] = BIG;
fvec[2] = BIG;
fvec[3] = BIG;
return;
}
double sum1, sum2, sumb, sum4, sum5;
double z, w, v;
double alpha = par[1];
double beta = par[2];
double mu = par[3];
if (iflag[1] == 1) {
sum1 = sum2 = sumb = sum4 = sum5 = 0;
for (int i = 0; i < n; i++) {
z = (x[i] - mu) / beta;
sum1 += 1.0 / z;
v = Math.pow(z, -alpha);
sum2 += v / z;
sumb += v;
w = Math.log(z);
sum4 += w;
sum5 += v * w;
}
fvec[2] = sumb - n; // eq. for beta
fvec[3] = (alpha + 1) * sum1 - alpha * sum2; // eq. for mu
fvec[1] = n / alpha + sum5 - sum4; // eq. for alpha
} else if (iflag[1] == 2) {
throw new IllegalArgumentException ("iflag = 2");
// The 3 X 3 Jacobian must be calculated and put in fjac
}
}
}
private static class Function implements MathFunction {
private int n;
private double[] x;
private double delta;
public double sumxi;
public double dif;
public Function (double[] y, int n, double delta) {
this.n = n;
this.x = y;
this.delta = delta;
double xmin = Double.MAX_VALUE;
for (int i = 0; i < n; i++) {
if ((y[i] < xmin) && (y[i] > delta))
xmin = y[i];
}
dif = xmin - delta;
}
public double evaluate (double alpha) {
if (alpha <= 0.0) return 1.0e100;
double v, w;
double sum1 = 0, sum2 = 0, sum3 = 0;
for (int i = 0; i < n; i++) {
if (x[i] <= delta)
continue;
v = Math.log(x[i] - delta);
w = Math.pow(dif / (x[i] - delta), alpha);
sum1 += v;
sum2 += w;
sum3 += v * w;
}
sum1 /= n;
sumxi = sum2 / n;
return 1 / alpha + sum3 / sum2 - sum1;
}
}
/**
* Constructor for the standard Fréchet
* distribution with parameters β = 1 and δ = 0.
*
*/
public FrechetDist (double alpha) {
setParams (alpha, 1.0, 0.0);
}
/**
* Constructs a FrechetDist object with parameters
* α = alpha, β = beta and δ = delta.
*
*/
public FrechetDist (double alpha, double beta, double delta) {
setParams (alpha, beta, delta);
}
public double density (double x) {
return density (alpha, beta, delta, x);
}
public double cdf (double x) {
return cdf (alpha, beta, delta, x);
}
public double barF (double x) {
return barF (alpha, beta, delta, x);
}
public double inverseF (double u) {
return inverseF (alpha, beta, delta, u);
}
public double getMean() {
return getMean (alpha, beta, delta);
}
public double getVariance() {
return getVariance (alpha, beta, delta);
}
public double getStandardDeviation() {
return getStandardDeviation (alpha, beta, delta);
}
/**
* Computes and returns the density function.
*
*/
public static double density (double alpha, double beta, double delta,
double x) {
if (beta <= 0)
throw new IllegalArgumentException ("beta <= 0");
if (alpha <= 0)
throw new IllegalArgumentException ("alpha <= 0");
final double z = (x - delta)/beta;
if (z <= 0.0)
return 0.0;
double t = Math.pow (z, -alpha);
return alpha * t * Math.exp (-t) / (z * beta);
}
/**
* Computes and returns the distribution function.
*
*/
public static double cdf (double alpha, double beta, double delta,
double x) {
if (beta <= 0)
throw new IllegalArgumentException ("beta <= 0");
if (alpha <= 0)
throw new IllegalArgumentException ("alpha <= 0");
final double z = (x - delta)/beta;
if (z <= 0.0)
return 0.0;
double t = Math.pow (z, -alpha);
return Math.exp (-t);
}
/**
* Computes and returns the complementary distribution function 1 - F(x).
*
*/
public static double barF (double alpha, double beta, double delta,
double x) {
if (beta <= 0)
throw new IllegalArgumentException ("beta <= 0");
if (alpha <= 0)
throw new IllegalArgumentException ("alpha <= 0");
final double z = (x - delta)/beta;
if (z <= 0.0)
return 1.0;
double t = Math.pow (z, -alpha);
return -Math.expm1 (-t);
}
/**
* Computes and returns the inverse distribution function.
*
*/
public static double inverseF (double alpha, double beta, double delta,
double u) {
if (u < 0.0 || u > 1.0)
throw new IllegalArgumentException ("u not in [0, 1]");
if (beta <= 0)
throw new IllegalArgumentException ("beta <= 0");
if (alpha <= 0)
throw new IllegalArgumentException ("alpha <= 0");
if (u >= 1.0)
return Double.POSITIVE_INFINITY;
if (u <= 0.0)
return delta;
double t = Math.pow (-Math.log (u), 1.0/alpha);
if (t <= Double.MIN_NORMAL)
return Double.MAX_VALUE;
return delta + beta / t;
}
/**
* Given δ = delta, estimates the parameters
* (α, β)
* of the Fréchet distribution
* using the maximum likelihood method with 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
*
* @param delta location parameter
*
* @return returns the parameters [
* hat(α),
* hat(β)]
*
*/
public static double[] getMLE (double[] x, int n, double delta) {
if (n <= 1)
throw new IllegalArgumentException ("n <= 1");
Function func = new Function (x, n, delta);
double a = 1e-4;
double b = 1.0e12;
double alpha = RootFinder.brentDekker (a, b, func, 1e-12);
double par[] = new double[2];
par[0] = alpha;
par[1] = func.dif * Math.pow (func.sumxi, -1.0/alpha);
return par;
}
/**
* Given δ = delta, creates a new instance of a Fréchet
* 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
*
* @param delta location parameter
*
*/
public static FrechetDist getInstanceFromMLE (double[] x, int n,
double delta) {
double par[] = getMLE (x, n, delta);
return new FrechetDist (par[0], par[1], delta);
}
/**
* Returns the mean of the Fréchet distribution with
* parameters α, β and δ.
*
* @return the mean
*
*/
public static double getMean (double alpha, double beta, double delta) {
if (beta <= 0)
throw new IllegalArgumentException ("beta <= 0");
if (alpha <= 1)
throw new IllegalArgumentException ("alpha <= 1");
double t = Num.lnGamma(1.0 - 1.0/alpha);
return delta + beta * Math.exp(t);
}
/**
* Returns the variance of the Fréchet distribution with parameters
* α, β and δ.
*
* @return the variance
*
*/
public static double getVariance (double alpha, double beta,
double delta) {
if (beta <= 0)
throw new IllegalArgumentException ("beta <= 0");
if (alpha <= 2)
throw new IllegalArgumentException ("alpha <= 2");
double t = Num.lnGamma(1.0 - 1.0/alpha);
double mu = Math.exp(t);
double v = Math.exp(Num.lnGamma(1.0 - 2.0/alpha));
return beta * beta * (v - mu * mu);
}
/**
* Returns the standard deviation of the Fréchet distribution
* with parameters α, β and δ.
*
* @return the standard deviation
*
*/
public static double getStandardDeviation (double alpha, double beta,
double delta) {
return Math.sqrt(getVariance (alpha, beta, delta));
}
/**
* Returns the parameter α of this object.
*
*/
public double getAlpha() {
return alpha;
}
/**
* Returns the parameter β of this object.
*
*/
public double getBeta() {
return beta;
}
/**
* Returns the parameter δ of this object.
*
*/
public double getDelta() {
return delta;
}
/**
* Sets the parameters α, β and δ of this object.
*
*/
public void setParams (double alpha, double beta, double delta) {
if (beta <= 0)
throw new IllegalArgumentException ("beta <= 0");
if (alpha <= 0)
throw new IllegalArgumentException ("alpha <= 0");
this.delta = delta;
this.beta = beta;
this.alpha = alpha;
}
/**
* Return an array containing the parameters of the current object
* in regular order: [α, β, δ].
*
*
*/
public double[] getParams() {
double[] retour = {alpha, beta, delta};
return retour;
}
public String toString () {
return getClass().getSimpleName() + " : alpha = " + alpha + ", beta = " + beta + ", delta = " + delta;
}
}