/* *
* 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.Misc;
import optimization.*;
/**
* Extends the class {@link ContinuousDistribution} for the
* Log-Logistic distribution with shape parameter
* α > 0
* and scale parameter β > 0.
* Its density is
*
*
*
* f (x) = (α(x/β)α-1)/(β[1 + (x/β)α]2) for x > 0
*
* and its distribution function is
*
*
*
* F(x) = 1/(1 + (x/β)-α) for x > 0.
*
* The complementary distribution is
*
*
*
* bar(F)(x) = 1/(1 + (x/β)α) for x > 0.
*
*
*/
public class LoglogisticDist extends ContinuousDistribution {
private double alpha;
private double beta;
private static class Optim implements Uncmin_methods
{
private int n;
private double[] xi;
public Optim (double[] x, int n)
{
this.n = n;
this.xi = new double[n];
System.arraycopy (x, 0, this.xi, 0, n);
}
public double f_to_minimize (double[] p)
{
if ((p[1] <= 0.0) || (p[2] <= 0.0))
return 1e200;
double sum = 0.0;
for (int i = 0; i < n; i++) {
double tmp = density (p[1], p[2], xi[i]);
if (tmp > 0.0)
sum -= Math.log (tmp);
else
sum += 709.0; // log (Double.MIN_VALUE)
}
return sum;
}
public void gradient (double[] x, double[] g)
{
}
public void hessian (double[] x, double[][] h)
{
}
}
/**
* Constructs a log-logistic distribution with parameters
* α and β.
*
*/
public LoglogisticDist (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 getMean (alpha, beta);
}
public double getVariance() {
return getVariance (alpha, beta);
}
public double getStandardDeviation() {
return getStandardDeviation (alpha, beta);
}
/**
* Computes the density function
* for a log-logisitic distribution with parameters α
* and β.
*
*/
public static double density (double alpha, double beta, double x) {
double denominateur;
if (alpha <= 0.0)
throw new IllegalArgumentException ("alpha <= 0");
if (beta <= 0.0)
throw new IllegalArgumentException ("beta <= 0");
if (x <= 0.0 || x >= Double.MAX_VALUE / 2.0)
return 0.0;
if (x <= beta) {
double v = Math.pow (x / beta, alpha);
denominateur = 1.0 + v;
denominateur *= denominateur * beta;
return alpha * v * beta / (x * denominateur);
} else {
double v = Math.pow (beta / x, alpha);
denominateur = 1.0 + v;
denominateur *= denominateur * beta;
return alpha * v * beta / (x * denominateur);
}
}
/**
* Computes the distribution function of the
* log-logistic distribution with parameters α and β.
*
*/
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 <= 0.0)
return 0.0;
if (x >= Double.MAX_VALUE / 2.0)
return 1.0;
double z = x/beta;
if (z >= 1.0)
return 1.0 / (1.0 + Math.pow (1.0/z, alpha));
double v = Math.pow (z, alpha);
return v/(v + 1.0);
}
/**
* Computes the complementary distribution function
* of the log-logistic distribution with parameters α and β.
*
*/
public static double barF (double alpha, double beta, double x) {
double power;
if (alpha <= 0.0)
throw new IllegalArgumentException ("alpha <= 0");
if (beta <= 0.0)
throw new IllegalArgumentException ("beta <= 0");
if (x <= 0.0)
return 1.0;
if (x >= Double.MAX_VALUE / 2.0)
return 0.0;
double z = x/beta;
if (z <= 1.0)
return 1.0 / (1.0 + Math.pow (z, alpha));
double v = Math.pow (1.0/z, alpha);
return v/(v + 1.0);
}
/**
* Computes the inverse of the log-logistic distribution
* with parameters α and β.
*
*/
public static double inverseF (double alpha, double beta, double u) {
if (alpha <= 0.0)
throw new IllegalArgumentException ("alpha <= 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 >= 1.0)
return Double.POSITIVE_INFINITY;
if (u <= 0.0)
return 0.0;
if (u <= 0.5)
return (beta * Math.pow (u / (1.0 - u), 1.0 / alpha));
else
return (beta / Math.pow ((1.0 - u)/ u, 1.0 / alpha));
}
/**
* Estimates the parameters
* (α, β) of the log-logistic 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 to use to evaluate parameters
*
* @param n the number of observations to use to evaluate parameters
*
* @return returns the parameters [
* hat(α),
* hat(β)]
*
*/
public static double[] getMLE (double[] x, int n) {
double sum = 0.0;
if (n <= 0)
throw new IllegalArgumentException ("n <= 0");
Optim system = new Optim (x, n);
double[] parameters = new double[2];
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[][] a = new double[3][3];
double[] udiag = new double[3];
param[2] = EmpiricalDist.getMedian (x, n);
if (param[2] < 0) throw new IllegalArgumentException ("median < 0");
if (param[2] <= 0) param[2] = 1.0;
int m = Math.round ((float) n / 4.0f);
double q1 = Misc.quickSelect (x, n, m);
if (q1 < 0) throw new IllegalArgumentException ("x[i] < 0");
if (q1 > 0)
param[1] = Math.log (3) / (Math.log(param[2]) - Math.log(q1));
else
param[1] = 1.0;
Uncmin_f77.optif0_f77 (2, param, system, xpls, fpls, gpls, itrcmd, a, udiag);
for (int i = 0; i < 2; i++)
parameters[i] = xpls[i+1];
return parameters;
}
/**
* Creates a new instance of a log-logistic 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 LoglogisticDist getInstanceFromMLE (double[] x, int n) {
double parameters[] = getMLE (x, n);
return new LoglogisticDist (parameters[0], parameters[1]);
}
/**
* Computes and returns the mean
* of the log-logistic distribution with parameters α and β.
*
* @return the mean of the log-logistic distribution
*
* E[X] = βθ cosec(θ), where θ = π/α
*
*/
public static double getMean (double alpha, double beta) {
double theta;
if (alpha <= 1.0)
throw new IllegalArgumentException ("alpha <= 1");
if (beta <= 0.0)
throw new IllegalArgumentException ("beta <= 0");
theta = Math.PI / alpha;
return (beta * theta / Math.sin (theta));
}
/**
* Computes and returns the variance
* of the log-logistic distribution with parameters α and β.
*
* @return the variance of the log-logistic distribution
*
* Var[X] = β2θ(2cosec(2θ) - θ[cosec(θ)]2), where θ = π/α
*
*/
public static double getVariance (double alpha, double beta) {
double theta;
if (alpha <= 2.0)
throw new IllegalArgumentException ("alpha <= 2");
if (beta <= 0.0)
throw new IllegalArgumentException ("beta <= 0");
theta = Math.PI / alpha;
return (beta * beta * theta * ((2.0 / Math.sin (2.0 * theta)) -
(theta / (Math.sin (theta) * Math.sin (theta)))));
}
/**
* Computes and returns the standard deviation of the log-logistic
* distribution with parameters α and β.
*
* @return the standard deviation of the log-logistic distribution
*
*/
public static double getStandardDeviation (double alpha, double beta) {
return Math.sqrt (getVariance (alpha, beta));
}
/**
* Return the parameter α of this object.
*
*/
public double getAlpha() {
return alpha;
}
/**
* Returns the parameter β of this object.
*
*/
public double getBeta() {
return beta;
}
/**
* Sets the parameters α and β of 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 = 0.0;
}
/**
* 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;
}
}
Ads help maintain this website.