/* *
* 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.functions.MathFunction;
/**
* Extends the class {@link ContinuousDistribution} for
* the Student t-distribution
* with n degrees of freedom, where n is a positive integer.
* Its density is
*
*
*
* f (x) = [Γ((n + 1)/2)/(Γ(n/2)(πn)1/2)][1 + x2/n]-(n+1)/2 for - ∞ < x < ∞,
*
* where Γ(x) is the gamma function defined in
* {@link GammaDist}.
*
*/
public class StudentDist extends ContinuousDistribution {
protected int n;
private double factor;
private static final int NLIM1 = 100000;
/*
private static double cdfPeizer (int n, double x) {
// Peizer-Pratt normal approximation for the cdf (n, u)
// \cite{tPEI68a}
double v = Math.log1p(x*x/n) / (n - 5.0/6.0);
double z = -(n - 2.0/3.0 + 0.1/n) * Math.sqrt(v);
double u = NormalDist.cdf01 (z);
if (x >= 0.0)
return 1.0 - u;
return u;
}
private static double invPeizer (int n, double u) {
// Peizer-Pratt normal approximation for the inverseF (n, u)
// \cite{tPEI68a}
double z = NormalDist.inverseF01 (u);
double q = z / (n - 2.0/3.0 + 0.1/n);
double v = q*q*(n - 5.0/6.0);
double t = Math.sqrt(n * Math.expm1(v));
if (u >= 0.5)
return t;
else
return -t;
}
*/
private static double cdfGaver (int n, double x) {
// Gaver-Kafadar normal approximation for the cdf
// \cite{tGAV84a}
double v = Math.log1p(x * x / n) / (n - 1.5);
double z = -(n - 1) * Math.sqrt(v);
double u = NormalDist.cdf01 (z);
if (x >= 0.0)
return 1.0 - u;
return u;
}
private static double invGaver (int n, double u) {
// Gaver-Kafadar normal approximation for the inverse
// \cite{tGAV84a}
double z = NormalDist.inverseF01 (u);
double q = z / (n - 1.0);
double v = q * q * (n - 1.5);
double t = Math.sqrt(n * Math.expm1(v));
if (u >= 0.5)
return t;
else
return -t;
}
private static class Function implements MathFunction {
private int n;
private double[] xi;
public Function (double[] x, int n) {
this.n = n;
this.xi = new double[n];
System.arraycopy (x, 0, this.xi, 0, n);
}
public double evaluate (double x) {
if (x <= 0.0)
return 1e200;
double sum = 0.0;
for (int i = 0; i < n; i++)
sum += Math.log (density ((int) Math.round (x), xi[i]));
return sum;
}
}
/**
* Constructs a StudentDist object with n degrees of freedom.
*
*/
public StudentDist (int n) {
setN (n);
}
public double density (double x) {
return factor*Math.pow (1.0 / (1.0 + x*x/n), (n + 1)/2.0);
}
public double cdf (double x) {
return cdf (n, x);
}
public double barF (double x) {
return barF (n, x);
}
public double inverseF (double u) {
return inverseF (n, u);
}
public double getMean() {
return StudentDist.getMean (n);
}
public double getVariance() {
return StudentDist.getVariance (n);
}
public double getStandardDeviation() {
return StudentDist.getStandardDeviation (n);
}
/**
* Computes the density function of a
* Student t-distribution with n degrees of freedom.
*
*/
public static double density (int n, double x) {
double factor = Num.gammaRatioHalf(n/2.0)/ Math.sqrt (n*Math.PI);
return factor*Math.pow (1.0 / (1.0 + x*x/n), (n + 1)/2.0);
}
/**
* Computes the Student t-distribution function u = F(x) with n degrees of freedom.
* Gives 13 decimal digits of precision for
* n <= 105.
* For n > 105, gives at least 6 decimal digits of precision everywhere, and
* at least 9 decimal digits of precision for all
* u > 10-15.
*
*/
public static double cdf (int n, double x) {
if (n <= 0)
throw new IllegalArgumentException ("n <= 0");
if (n == 1)
return CauchyDist.cdf(0, 1, x);
if (x > 1.0e10)
return 1.0;
if (n > NLIM1)
return cdfGaver(n, x);
double r = Math.abs(x);
if (r < 1.0e15)
r = Math.sqrt (n + x*x);
if (n == 2) {
if (x >= 0.) {
return 0.5*(1.0 + x/r);
} else
return 1.0 / (r * (r - x));
}
if (x > 0.0)
return 1.0 - cdf (n, -x);
double z = 0.5*n / (r*(r - x));
return BetaSymmetricalDist.cdf (0.5*n, 15, z);
}
/**
* Same as {@link #cdf(int,double) cdf}(n, x).
*
*/
@Deprecated
public static double cdf2 (int n, int d, double x) {
if (d <= 0)
throw new IllegalArgumentException ("student2: d <= 0");
return cdf (n, x);
}
/**
* Computes the complementary distribution function
* v = bar(F)(x)
* with n degrees of freedom. Gives 13 decimal digits of precision for
* n <= 105.
* For n > 105, gives at least 6 decimal digits of precision everywhere, and
* at least 9 decimal digits of precision for all
* v > 10-15.
*
*/
public static double barF (int n, double x) {
if (n < 1)
throw new IllegalArgumentException ("n < 1");
if (n == 1)
return CauchyDist.barF(0, 1, x);
if (n == 2) {
double z = Math.abs(x);
if (z < 1.0e10)
z = Math.sqrt(2.0 + x*x);
if (x <= 0.) {
if (x < -1.0e10)
return 1.0;
return 0.5* (1.0 - x / z);
} else
return 1.0 / (z * (z + x));
}
return cdf (n, -x);
}
/**
* Returns the inverse
* x = F-1(u) of
* Student t-distribution function with n degrees of freedom.
* Gives 13 decimal digits of precision for
* n <= 105,
* and at least 9 decimal digits of precision for n > 105.
*
*/
public static double inverseF (int n, double u) {
if (n < 1)
throw new IllegalArgumentException ("Student: n < 1");
if (u > 1.0 || u < 0.0)
throw new IllegalArgumentException ("Student: u not in [0, 1]");
if (u <= 0.0)
return Double.NEGATIVE_INFINITY;
if (u >= 1.0)
return Double.POSITIVE_INFINITY;
if (1 == n)
return CauchyDist.inverseF(0, 1, u);
if (2 == n)
return (2.0*u - 1.0) / Math.sqrt(2.0*u*(1.0 - u));
if (n > NLIM1)
return invGaver(n, u);
double z = BetaSymmetricalDist.inverseF (0.5*n, u);
return (z - 0.5) * Math.sqrt(n / (z*(1.0 - z)));
}
/**
* Estimates the parameter n of the Student t-distribution
* using the maximum likelihood method, from the m observations
* x[i],
* i = 0, 1,…, m - 1. The estimate is returned in a one-element
* array.
*
* @param x the list of observations to use to evaluate parameters
*
* @param m the number of observations to use to evaluate parameters
*
* @return returns the parameter [hat(n)]
*
*/
public static double[] getMLE (double[] x, int m) {
double sum = 0.0;
double[] parameters = new double[1];
if (m <= 0)
throw new IllegalArgumentException ("m <= 0");
double var = 0.0;
for (int i = 0; i < m; i++)
var += x[i] * x[i];
var /= (double) m;
Function f = new Function (x, m);
double n0 = Math.round ((2.0 * var) / (var - 1.0));
double fn0 = f.evaluate (n0);
double min = fn0;
double fn1 = f.evaluate (n0 + 1.0);
double fn_1 = f.evaluate (n0 - 1.0);
parameters[0] = n0;
if (fn_1 > fn0) {
double n = n0 - 1.0;
double y;
while (((y = f.evaluate (n)) > min) && (n >= 1.0)) {
min = y;
parameters[0] = n;
n -= 1.0;
}
} else if (fn1 > fn0) {
double n = n0 + 1.0;
double y;
while ((y = f.evaluate (n)) > min) {
min = y;
parameters[0] = n;
n += 1.0;
}
}
return parameters;
}
/**
* Creates a new instance of a Student t-distribution with parameter n
* estimated using the maximum likelihood method based on the m observations
* x[i],
* i = 0, 1,…, m - 1.
*
* @param x the list of observations to use to evaluate parameters
*
* @param m the number of observations to use to evaluate parameters
*
*
*/
public static StudentDist getInstanceFromMLE (double[] x, int m) {
double parameters[] = getMLE (x, m);
return new StudentDist ((int) parameters[0]);
}
/**
* Returns the mean E[X] = 0 of the Student t-distribution with
* parameter n.
*
* @return the mean of the Student t-distribution E[X] = 0
*
*/
public static double getMean (int n) {
if (n < 2)
throw new IllegalArgumentException ("n <= 1");
return 0;
}
/**
* Computes and returns the variance
* Var[X] = n/(n - 2)
* of the Student t-distribution with parameter n.
*
* @return the variance of the Student t-distribution
*
* Var[X] = n/(n - 2)
*
*/
public static double getVariance (int n) {
if (n < 3)
throw new IllegalArgumentException("n <= 2");
return (n / (n - 2.0));
}
/**
* Computes and returns the standard deviation
* of the Student t-distribution with parameter n.
*
* @return the standard deviation of the Student t-distribution
*
*/
public static double getStandardDeviation (int n) {
return Math.sqrt (StudentDist.getVariance (n));
}
/**
* Returns the parameter n associated with this object.
*
*/
public int getN() {
return n;
}
/**
* Sets the parameter n associated with this object.
*
*/
public void setN (int n) {
if (n <= 0)
throw new IllegalArgumentException ("n <= 0");
this.n = n;
factor = Num.gammaRatioHalf(n/2.0) / Math.sqrt (n*Math.PI);
}
/**
* Return a table containing the parameter of the current distribution.
*
*
*/
public double[] getParams () {
double[] retour = {n};
return retour;
}
public String toString () {
return getClass().getSimpleName() + " : n = " + n;
}
}
Ads help maintain this website.