/* *
* Title:
* Description:
* Copyright: Copyright (c)
* Company:
* @author
* @version 1.0
*/
package umontreal.iro.lecuyer.probdist;
import umontreal.iro.lecuyer.probdist.NormalDist;
import optimization.*;
/**
* Extends the class {@link ContinuousDistribution} for
* the fatigue life distribution with location
* parameter μ, scale parameter β and shape
* parameter γ.
* Its density is
*
*
*
* f (x) = [(((x - μ)/β)1/2 + (β/(x - μ))1/2)/(2γ(x - μ))]φ((((x - μ)/β)1/2 - (β/(x - μ))1/2)/γ), for x > μ,
*
* where φ is the probability density of the standard normal distribution.
* The distribution function is given by
*
*
*
* F(x) = Φ((((x - μ)/β)1/2 - (β/(x - μ))1/2)/γ), for x > μ,
*
* where Φ is the standard normal distribution function.
* Restrictions: β > 0,
* γ > 0.
*
*
* The non-static versions of the methods cdf, barF,
* and inverseF call the static version of the same name.
*
*/
public class FatigueLifeDist extends ContinuousDistribution {
protected double mu;
protected double beta;
protected double gamma;
private static class Optim implements Uncmin_methods
{
private int n;
private double[] xi;
private double mu;
public Optim (double[] x, int n, double min)
{
this.n = n;
this.mu = min;
this.xi = new double[n];
System.arraycopy (x, 0, this.xi, 0, n);
}
public double f_to_minimize (double[] p)
{
double sum = 0.0;
if ((p[1] <= 0.0) || (p[2] <= 0.0))
return 1e200;
for (int i = 0; i < n; i++)
sum -= Math.log (density (mu, p[1], p[2], xi[i]));
return sum;
}
public void gradient (double[] x, double[] g)
{
}
public void hessian (double[] x, double[][] h)
{
}
}
/**
* Constructs a fatigue life distribution with parameters
* μ, β and γ.
*
*/
public FatigueLifeDist (double mu, double beta, double gamma) {
setParams (mu, beta, gamma);
}
public double density (double x) {
return FatigueLifeDist.density (mu, beta, gamma, x);
}
public double cdf (double x) {
return FatigueLifeDist.cdf (mu, beta, gamma, x);
}
public double barF (double x) {
return FatigueLifeDist.barF (mu, beta, gamma, x);
}
public double inverseF (double u) {
return FatigueLifeDist.inverseF (mu, beta, gamma, u);
}
public double getMean() {
return FatigueLifeDist.getMean (mu, beta, gamma);
}
public double getVariance() {
return FatigueLifeDist.getVariance (mu, beta, gamma);
}
public double getStandardDeviation() {
return FatigueLifeDist.getStandardDeviation (mu, beta, gamma);
}
/**
* Computes the density for the
* fatigue life distribution with parameters μ, β and γ.
*
*/
public static double density (double mu, double beta, double gamma,
double x) {
if (beta <= 0.0)
throw new IllegalArgumentException ("beta <= 0");
if (gamma <= 0.0)
throw new IllegalArgumentException ("gamma <= 0");
if (x <= mu)
return 0.0;
double y;
y = (Math.sqrt ((x - mu) / beta) - Math.sqrt (beta / (x - mu))) / gamma;
return (((Math.sqrt ((x - mu) / beta) + Math.sqrt (beta / (x - mu))) /
(2 * gamma * (x - mu))) * NormalDist.density (0.0, 1.0, y));
}
/**
* Computes the fatigue life distribution
* function with parameters μ, β and γ.
*
*/
public static double cdf (double mu, double beta, double gamma, double x) {
if (beta <= 0.0)
throw new IllegalArgumentException ("beta <= 0");
if (gamma <= 0.0)
throw new IllegalArgumentException ("gamma <= 0");
if (x <= mu)
return 0.0;
return NormalDist.cdf01 ((Math.sqrt ((x - mu) / beta) - Math.sqrt (beta / (x - mu))) / gamma);
}
/**
* Computes the complementary distribution function of the
* fatigue life distribution with parameters μ, β and γ.
*
*/
public static double barF (double mu, double beta, double gamma,
double x) {
if (beta <= 0.0)
throw new IllegalArgumentException ("beta <= 0");
if (gamma <= 0.0)
throw new IllegalArgumentException ("gamma <= 0");
if (x <= mu)
return 1.0;
return NormalDist.barF01 ((Math.sqrt ((x - mu) / beta) - Math.sqrt (beta / (x - mu))) / gamma);
}
/**
* Computes the inverse of the fatigue life distribution
* with parameters μ, β and γ.
*
*/
public static double inverseF (double mu, double beta, double gamma,
double u) {
if (beta <= 0.0)
throw new IllegalArgumentException ("beta <= 0");
if (gamma <= 0.0)
throw new IllegalArgumentException ("gamma <= 0");
if (u > 1.0 || u < 0.0)
throw new IllegalArgumentException ("u not in [0,1]");
if (u <= 0.0) // if u == 0, in fact
return mu;
if (u >= 1.0) // if u == 1, in fact
return Double.POSITIVE_INFINITY;
double w = gamma * NormalDist.inverseF01 (u);
double sqrtZ = 0.5 * (w + Math.sqrt (w * w + 4.0));
return (mu + sqrtZ * sqrtZ * beta);
}
/**
* Estimates the parameters (μ, β, γ) of the fatigue life
* distribution using the maximum likelihood method, from the n observations
* x[i],
* i = 0, 1,…, n - 1. The estimates are returned in a three-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
*
* @param mu the location parameter
*
* @return returns the parameters [
* hat(β),
* hat(γ)]
*
*/
public static double[] getMLE (double[] x, int n, double mu) {
double sum = 0.0;
if (n <= 0)
throw new IllegalArgumentException ("n <= 0");
double[] parameters = new double[3];
double[] xpls = new double[3];
double[] param = new double[3];
double[] fpls = new double[3];
double[] gpls = new double[3];
int[] itrcmd = new int[2];
double[][] h = new double[3][3];
double[] udiag = new double[3];
Optim system = new Optim (x, n, mu);
double mean = 0.0;
for (int i = 0; i < n; i++)
mean += x[i];
mean /= (double) n;
double var = 0.0;
for (int i = 0; i < n; i++)
var += (x[i] - mean) * (x[i] - mean);
var /= (double) n;
double loc2 = (mean - mu) * (mean - mu);
double a = 0.25 * (var - 5 * loc2);
double b = (var - loc2);
double c = var;
double delta = b * b - 4.0 * a * c;
double gamma2 = (- b - Math.sqrt (delta)) / (2.0 * a);
param[2] = Math.sqrt (gamma2);
param[1] = (mean - mu) / (1.0 + gamma2 / 2.0);
Uncmin_f77.optif0_f77 (2, param, system, xpls, fpls, gpls, itrcmd, h, udiag);
for (int i = 1; i < 3; i++)
parameters[i] = xpls[i];
parameters[0] = mu;
return parameters;
}
/**
* Computes and returns the mean
*
* E[X] = μ + β(1 + γ2/2)
* of the fatigue life distribution with parameters μ, β and γ.
* @return the mean of the fatigue life distribution
*
*/
public static double getMean (double mu, double beta, double gamma) {
if (beta <= 0.0)
throw new IllegalArgumentException ("beta <= 0");
if (gamma <= 0.0)
throw new IllegalArgumentException ("gamma <= 0");
return (mu + beta * (1 + 0.5 * gamma * gamma));
}
/**
* Computes and returns the variance
*
* Var[X] = β2γ2(1 + 5γ2/4) of the fatigue life distribution
* with parameters μ, β and γ.
* @return the variance of the fatigue life distribution
*
*/
public static double getVariance (double mu, double beta, double gamma) {
if (beta <= 0.0)
throw new IllegalArgumentException ("beta <= 0");
if (gamma <= 0.0)
throw new IllegalArgumentException ("gamma <= 0");
return (beta * beta * gamma * gamma * (1.0 + 5.0/4.0 * gamma * gamma));
}
/**
* Computes and returns the standard deviation
* of the fatigue life distribution
* with parameters μ, β and γ.
*
* @return the standard deviation of the fatigue life distribution
*
*/
public static double getStandardDeviation (double mu, double beta,
double gamma) {
return Math.sqrt (FatigueLifeDist.getVariance (mu, beta, gamma));
}
/**
* Returns the parameter β of this object.
*
*/
public double getBeta() {
return beta;
}
/**
* Returns the parameter γ of this object.
*
*/
public double getGamma() {
return gamma;
}
/**
* Returns the parameter μ of this object.
*
*/
public double getMu() {
return mu;
}
/**
* Sets the parameters μ, β and γ of this object.
*
*/
public void setParams (double mu, double beta, double gamma) {
if (beta <= 0.0)
throw new IllegalArgumentException ("beta <= 0");
if (gamma <= 0.0)
throw new IllegalArgumentException ("gamma <= 0");
this.mu = mu;
this.beta = beta;
this.gamma = gamma;
supportA = mu;
}
/**
* Return a table containing the parameters of the current distribution.
* This table is put in regular order: [μ, β, γ].
*
*
*/
public double[] getParams () {
double[] retour = {mu, beta, gamma};
return retour;
}
public String toString () {
return getClass().getSimpleName() + " : mu = " + mu + ", beta = " + beta + ", gamma = " + gamma;
}
}
Ads help maintain this website.