/* *
* Title:
* Description:
* Copyright: Copyright (c)
* Company:
* @author
* @version 1.0
*/
package umontreal.iro.lecuyer.probdist;
import umontreal.iro.lecuyer.util.*;
import umontreal.iro.lecuyer.functions.MathFunction;
/**
* Extends the class {@link ContinuousDistribution} for
* the chi-square distribution with n degrees of freedom,
* where n is a positive integer.
* Its density is
*
*
*
* f (x) = x(n/2)-1e-x/2/(2n/2Γ(n/2)), for x > 0.
*
* where Γ(x) is the gamma function defined in
* {@link GammaDist}.
* The chi-square distribution is a special case of the gamma
* distribution with shape parameter n/2 and scale parameter 1/2.
* Therefore, one can use the methods of {@link GammaDist} for this distribution.
*
*
* The non-static versions of the methods cdf, barF,
* and inverseF call the static version of the same name.
*
*/
public class ChiSquareDist extends ContinuousDistribution {
protected int n;
protected double C1;
private static class Function implements MathFunction {
protected int n;
protected double sumLog;
public Function (double s, int n)
{
this.n = n;
this.sumLog = s;
}
public double evaluate (double k)
{
if (k < 1.0) return 1.0e200;
return (sumLog + n * (Num.lnGamma (k / 2.0) - 0.5*Num.LN2 - Num.lnGamma ((k + 1.0) / 2.0)));
}
}
/**
* Constructs a chi-square distribution with n degrees of freedom.
*
*/
public ChiSquareDist (int n) {
setN (n);
}
public double density (double x) {
if (x <= 0)
return 0.0;
return Math.exp ((n/2.0 - 1)*Math.log (x) - x/2.0 - C1);
}
public double cdf (double x) {
return cdf (n, decPrec, x);
}
public double barF (double x) {
return barF (n, decPrec, x);
}
public double inverseF (double u) {
return inverseF (n, u);
}
public double getMean() {
return ChiSquareDist.getMean (n);
}
public double getVariance() {
return ChiSquareDist.getVariance (n);
}
public double getStandardDeviation() {
return ChiSquareDist.getStandardDeviation (n);
}
/**
* Computes the density function
* for a chi-square distribution with n degrees of freedom.
*
*/
public static double density (int n, double x) {
if (x <= 0)
return 0.0;
return Math.exp ((n/2.0 - 1)*Math.log (x) - x/2.0 -
(n/2.0)*Num.LN2 - Num.lnGamma(n/2.0));
}
/**
* Computes the chi-square distribution function with n degrees of freedom,
* evaluated at x. The method tries to return d decimals digits of precision,
* but there is no guarantee.
*
*/
public static double cdf (int n, int d, double x) {
if (n <= 0)
throw new IllegalArgumentException ("n <= 0");
if (x <= 0.0)
return 0.0;
if (x >= XBIG*n)
return 1.0;
return GammaDist.cdf (n/2.0, d, x/2.0);
}
/**
* Computes the complementary chi-square distribution function with n degrees
* of freedom, evaluated at x. The method tries to return d decimals digits
* of precision, but there is no guarantee.
*
*/
public static double barF (int n, int d, double x) {
if (n <= 0)
throw new IllegalArgumentException ("n <= 0");
if (x <= 0.0)
return 1.0;
return GammaDist.barF (n/2.0, d, x/2.0);
}
/**
* Computes an approximation of F-1(u), where F is the
* chi-square distribution with n degrees of freedom.
*
* It gives at least 6 decimal digits of precision, except far in the tails
* (that is, for
* u < 10-5 or
* u > 1 - 10-5) where the function
* calls the method GammaDist.inverseF (n/2, 7, u) and
* multiplies the result by 2.0.
* To get better precision, one may
* call GammaDist.inverseF, but this method is slower than the
* current method, especially for large n. For instance, for n = 16, 1024,
* and 65536, the GammaDist.inverseF method is 2, 5, and 8 times slower,
* respectively, than the current method.
*
*/
public static double inverseF (int n, double u) {
/*
* Returns an approximation of the inverse of Chi square cdf
* with n degrees of freedom.
* As in Figure L.23 of P.Bratley, B.L.Fox, and L.E.Schrage.
* A Guide to Simulation Springer-Verlag,
* New York, second edition, 1987.
*/
if (n <= 0)
throw new IllegalArgumentException ("n <= 0");
if (u < 0.0 || u > 1.0)
throw new IllegalArgumentException ("u must be in [0,1]");
if (u == 1.0)
return Double.POSITIVE_INFINITY;
if (u == 0.0)
return 0.0;
final double E = 0.5e-5; // Precision of this approximation
final double AA = 0.6931471805;
double A, XX, X, C, G, CH, Q, P1, P2, T, B, S1, S2, S3, S4, S5, S6;
if (u < 0.00001 || u > 1.0 - 1.0e-5)
return 2.0 * GammaDist.inverseF (n / 2.0, 7, u);
if (u >= 1.0)
return n * XBIG;
if (u >= 0.999998)
return (n + 4.0 * Math.sqrt (2.0 * n));
G = Num.lnGamma (n / 2.0);
XX = 0.5 * n;
C = XX - 1.0;
if (n >= -1.24 * Math.log (u)) {
X = NormalDist.inverseF01 (u);
P1 = 0.222222 / n;
Q = X * Math.sqrt (P1) + 1.0 - P1;
CH = n * Q * Q * Q;
if (CH > 2.2 * n + 6.0)
CH = -2.0 * (Math.log1p (-u) - C * Math.log (0.5 * CH) + G);
} else {
CH = Math.pow (u * XX * Math.exp (G + XX * AA), 1.0 / XX);
if (CH - E < 0)
return CH;
}
Q = CH;
P1 = 0.5 * CH;
P2 = u - GammaDist.cdf (XX, 5, P1);
if (GammaDist.cdf (XX, 5, P1) == -1.0)
throw new IllegalArgumentException ("RESULT = -1");
T = P2 * Math.exp (XX * AA + G + P1 - C * Math.log (CH));
B = T / CH;
A = 0.5 * T - B * C;
S1 = (210.0 + A * (140.0 +
A * (105.0 + A * (84.0 + A * (70.0 + 60.0 * A))))) / 420.0;
S2 = (420.0 + A * (735.0 + A * (966.0 + A * (1141.0 + 1278.0 * A))))
/ 2520.0;
S3 = (210.0 + A * (462.0 + A * (707.0 + 932.0 * A))) / 2520.0;
S4 = (252.0 + A * (672.0 + 1182.0 * A) +
C * (294.0 + A * (889.0 + 1740.0 * A))) / 5040.0;
S5 = (84.0 + 264.0 * A + C * (175.0 + 606.0 * A)) / 2520.0;
S6 = (120.0 + C * (346.0 + 127.0 * C)) / 5040.0;
CH = CH + T * (1.0 + 0.5 * T * S1 - B * C * (S1 - B * (S2 -
B * (S3 - B * (S4 - B * (S5 - B * S6))))));
double temp;
while (Math.abs (Q / CH - 1.0) > E) {
Q = CH;
P1 = 0.5 * CH;
temp = GammaDist.cdf (XX, 6, P1);
P2 = u - temp;
if (temp == -1.0)
return -1.0;
T = P2 * Math.exp (XX * AA + G + P1 - C * Math.log (CH));
B = T / CH;
A = 0.5 * T - B * C;
S1 = (210.0 + A * (140.0 + A * (105.0 + A * (84.0 +
A * (70.0 + 60.0 * A))))) / 420.0;
S2 = (420.0 + A * (735.0 + A * (966.0 + A * (1141.0 +
1278.0 * A)))) / 2520.0;
S3 = (210.0 + A * (462.0 + A * (707.0 + 932.0 * A))) / 2520.0;
S4 = (252.0 + A * (672.0 + 1182.0 * A) +
C * (294.0 + A * (889.0 + 1740.0 * A))) / 5040.0;
S5 = (84.0 + 264.0 * A + C * (175.0 + 606.0 * A)) / 2520.0;
S6 = (120.0 + C * (346.0 + 127.0 * C)) / 5040.0;
CH = CH + T * (1.0 + 0.5 * T * S1 - B * C * (S1 - B * (S2 -
B * (S3 - B * (S4 - B * (S5 - B * S6))))));
}
return CH;
}
/**
* Estimates the parameter n of the chi-square distribution
* using the maximum likelihood method, from the m observations
* x[i],
* i = 0, 1,…, m - 1. The estimate is returned in element 0
* of the returned 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) {
if (m <= 0)
throw new IllegalArgumentException ("m <= 0");
double[] parameters;
parameters = getMomentsEstimate (x, m);
double k = Math.round (parameters[0]) - 5.0;
if (k < 1.0)
k = 1.0;
double sum = 0.0;
for (int i = 0; i < m; i++) {
if (x[i] > 0.0)
sum += 0.5*Math.log (x[i]);
else
sum -= 709.0;
}
Function f = new Function (sum, m);
while (f.evaluate(k) > 0.0)
k++;
parameters[0] = k;
return parameters;
}
/**
* Creates a new instance of a chi-square 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 ChiSquareDist getInstanceFromMLE (double[] x, int m) {
double parameters[] = getMLE (x, m);
return new ChiSquareDist ((int) parameters[0]);
}
/**
* Computes and returns the mean E[X] = n of the
* chi-square distribution with parameter n.
*
* @return the mean of the Chi-square distribution E[X] = n
*
*/
public static double getMean (int n) {
if (n <= 0)
throw new IllegalArgumentException ("degrees of freedom " +
"must be non-null and positive.");
return n;
}
/**
* Estimates and returns the parameter [hat(n)] of the chi-square
* distribution using the moments method based on the m observations
* in table 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
*
* @return returns the parameter [hat(n)]
*
*/
public static double[] getMomentsEstimate (double[] x, int m) {
double[] parameters = new double[1];
double sum = 0.0;
for (int i = 0; i < m; i++)
sum += x[i];
parameters[0] = sum / (double) m;
return parameters;
}
/**
* Returns the variance
* Var[X] = 2n
* of the chi-square distribution with parameter n.
*
* @return the variance of the chi-square distribution
* VarX] = 2n
*
*/
public static double getVariance (int n) {
if (n <= 0)
throw new IllegalArgumentException ("degrees of freedom " +
"must be non-null and positive.");
return (2 * n);
}
/**
* Returns the standard deviation
* of the chi-square distribution with parameter n.
*
* @return the standard deviation of the chi-square distribution
*
*/
public static double getStandardDeviation (int n) {
if (n <= 0)
throw new IllegalArgumentException ("degrees of freedom " +
"must be non-null and positive.");
return Math.sqrt(2 * n);
}
/**
* Returns the parameter n of this object.
*
*/
public int getN() {
return n;
}
/**
* Sets the parameter n of this object.
*
*/
public void setN (int n) {
if (n <= 0)
throw new IllegalArgumentException ("degrees of freedom " +
"must be non-null and positive.");
this.n = n;
supportA = 0.0;
C1 = 0.5 * n *Num.LN2 + Num.lnGamma(n/2.0);
}
/**
* Return a table containing the parameters 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.