/* *
* Title:
* Description:
* Copyright: Copyright (c)
* Company:
* @author
* @version 1.0
*/
package umontreal.iro.lecuyer.probdist;
import umontreal.iro.lecuyer.util.Misc;
import optimization.*;
/**
* Extends the class {@link ContinuousDistribution} for
* the Cauchy distribution
* with location parameter α
* and scale parameter β > 0.
* The density function is given by
*
*
*
* f (x) = β/(π[(x - α)2 + β2]) for - ∞ < x < ∞.
*
* The distribution function is
*
*
*
* F(x) = 1/2 + arctan((x - α)/β)/π, for - ∞ < x < ∞,
*
* and its inverse is
*
*
*
* F-1(u) = α + βtan(π(u - 1/2)). for 0 < u < 1.
*
*
*/
public class CauchyDist 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)
{
double sum = 0.0;
if (p[2] <= 0.0) // barrier at 0
return 1.0e200;
for (int i = 0; i < n; i++)
sum -= Math.log (density (p[1], p[2], xi[i]));
return sum;
}
public void gradient (double[] x, double[] g)
{
}
public void hessian (double[] x, double[][] h)
{
}
}
/**
* Constructs a CauchyDist object
* with parameters α = 0 and β = 1.
*
*/
public CauchyDist() {
setParams (0.0, 1.0);
}
/**
* Constructs a CauchyDist object with parameters
* α = alpha and β = beta.
*
*/
public CauchyDist (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 CauchyDist.getMean (alpha, beta);
}
public double getVariance() {
return CauchyDist.getVariance (alpha, beta);
}
public double getStandardDeviation() {
return CauchyDist.getStandardDeviation (alpha, beta);
}
/**
* Computes the density function.
*
*/
public static double density (double alpha, double beta, double x) {
if (beta <= 0.0)
throw new IllegalArgumentException ("beta <= 0");
double t = (x - alpha)/beta;
return 1.0/(beta * Math.PI*(1 + t*t));
}
/**
* Computes the distribution function.
*
*/
public static double cdf (double alpha, double beta, double x) {
// The integral was computed analytically using Mathematica
if (beta <= 0.0)
throw new IllegalArgumentException ("beta <= 0");
double z = (x - alpha)/beta;
if (z < -0.5)
return Math.atan(-1.0/z)/Math.PI;
return Math.atan(z)/Math.PI + 0.5;
}
/**
* Computes the complementary distribution.
*
*/
public static double barF (double alpha, double beta, double x) {
if (beta <= 0.0)
throw new IllegalArgumentException ("beta <= 0");
double z = (x - alpha)/beta;
if (z > 0.5)
return Math.atan(1./z)/Math.PI;
return 0.5 - Math.atan(z)/Math.PI;
}
/**
* Computes the inverse of the distribution.
*
*/
public static double inverseF (double alpha, double beta, double u) {
if (beta <= 0.0)
throw new IllegalArgumentException ("beta <= 0");
if (u < 0.0 || u > 1.0)
throw new IllegalArgumentException ("u must be in [0,1]");
if (u <= 0.0)
return Double.NEGATIVE_INFINITY;
if (u >= 1.0)
return Double.POSITIVE_INFINITY;
if (u < 0.5)
return alpha - 1.0/Math.tan (Math.PI*u) * beta;
return alpha + Math.tan (Math.PI*(u - 0.5)) * beta;
}
/**
* Estimates the parameters
* (α, β) of the Cauchy 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[1] = EmpiricalDist.getMedian (x, n);
int m = Math.round ((float) n / 4.0f);
double q3 = Misc.quickSelect (x, n, 3 * m);
double q1 = Misc.quickSelect (x, n, m);
param[2] = (q3 - q1) / 2.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 Cauchy 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 CauchyDist getInstanceFromMLE (double[] x, int n) {
double parameters[] = getMLE (x, n);
return new CauchyDist (parameters[0], parameters[1]);
}
/**
* Throws an exception since the mean does not exist.
*
* @exception UnsupportedOperationException the mean of the Cauchy distribution is undefined.
*
*
*/
public static double getMean (double alpha, double beta) {
if (beta <= 0.0)
throw new IllegalArgumentException ("beta <= 0");
throw new UnsupportedOperationException("Undefined mean");
}
/**
* Returns ∞ since the variance does not exist.
*
* @return ∞.
*
*/
public static double getVariance (double alpha, double beta) {
if (beta <= 0.0)
throw new IllegalArgumentException ("beta <= 0");
return Double.POSITIVE_INFINITY;
}
/**
* Returns ∞ since the standard deviation does not exist.
*
* @return ∞
*
*/
public static double getStandardDeviation (double alpha, double beta) {
return Double.POSITIVE_INFINITY;
}
/**
* Returns the value of α for this object.
*
*/
public double getAlpha() {
return alpha;
}
/**
* Returns the value of β for this object.
*
*/
public double getBeta() {
return beta;
}
/**
* Sets the value of the parameters α and β for this object.
*
*/
public void setParams (double alpha, double beta) {
if (beta <= 0.0)
throw new IllegalArgumentException ("beta <= 0");
this.alpha = alpha;
this.beta = beta;
}
/**
* Return a table containing 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.