net.finmath.functions.NonCentralChiSquaredDistribution Maven / Gradle / Ivy

Show more of this group Show more artifacts with this name
Show all versions of finmath-lib Show documentation
finmath lib is a Mathematical Finance Library in Java. It provides algorithms and methodologies related to mathematical finance.
The newest version!
package net.finmath.functions;

import org.apache.commons.lang3.Validate;
import org.apache.commons.math3.special.Gamma;

/**
 * Implementation of the cumulative distribution function of the non-central Χ² distribution.
 *
 * The implementation is that of Ding's algorithm, splitting the sum since the first term
 * of the sum is not the largest one (and may be zero).
 *
 * The algorithm uses Kahan summation for the sums and limits the calculation time to avoid infinity loops.
 *
 * The method is currently intended to be used as a benchmark in unit tests only.
 *
 * Note: The implementation currently does not use an alternative algorithm
 * in cases of high non-centrality, like Benton and Krishnamoorthy. Since the number
 * of summation is limited, the method may be inaccurate for high non-centrality.
 *
 * @author Christian Fries
 * @author Ralph Rudd
 */
public class NonCentralChiSquaredDistribution {

	private final double nonCentralityHalf;
	private final double degreesOfFreedomHalf;

	private final int summationSplitIndex;
	private final int maxSummation;
	private final double pInitial;

	private static final double PRECISION = 1e-16;

	/**
	 * Create non-central Χ² distribution (non-central chi-squared distribution).
	 *
	 * @param degreesOfFreedom The number of degrees of freedom, positive
	 * @param nonCentrality The non-centrality parameter, not negative
	 */
	public NonCentralChiSquaredDistribution(final double degreesOfFreedom, double nonCentrality) {
		Validate.isTrue(degreesOfFreedom > 0, "Parameter degreesOfFreedom must be > 0 (given %d).", degreesOfFreedom);
		Validate.isTrue(nonCentrality >= 0, "Parameter nonCentrality must be >= 0 (given %d).", nonCentrality);

		degreesOfFreedomHalf = degreesOfFreedom / 2.0;
		nonCentralityHalf = nonCentrality / 2.0;
		summationSplitIndex = (int) Math.round(nonCentralityHalf);
		maxSummation = Math.max(summationSplitIndex, 10000);

		if (nonCentralityHalf == 0) {
			pInitial = 0.0;
		} else {
			pInitial = Math.exp(-nonCentralityHalf + summationSplitIndex * Math.log(nonCentralityHalf) - Gamma.logGamma(summationSplitIndex + 1));
		}
	}

	/**
	 * Cumulative distribution function of the non-central Χ² distribution
	 *
	 * @param x A sample point
	 * @return The probability of X being below x, given that X is non-central Χ² distributed
	 */
	public double cumulativeDistribution(final double x) {
		if (x < 0) {
			return 0.0;
		}

		final double xHalf		= x / 2.0;
		final double xHalfLog	= Math.log(xHalf);

		final double gammaRegularizedInitial = Gamma.regularizedGammaP(degreesOfFreedomHalf + summationSplitIndex, xHalf);

		double p = pInitial;
		double gammaRegularized = gammaRegularizedInitial;

		/*
		 * Sum of terms below summationSplitIndex, backward, using Kahan summation.
		 */
		double sum = p * gammaRegularized;
		double error = 0.0;
		double newSum = 0.0;
		for(int i=summationSplitIndex-1; i >= 0; i--) {

			p *= (i + 1) / nonCentralityHalf;
			gammaRegularized += Math.exp((degreesOfFreedomHalf + i) * xHalfLog - xHalf - Gamma.logGamma(degreesOfFreedomHalf + i + 1));

			// Kahan summation
			final double value = p * gammaRegularized - error;

			newSum = sum + value;
			error = (newSum - sum) - value;

			if(Math.abs(newSum - sum) / newSum <= PRECISION) {
				break;
			}
			sum = newSum;
		}

		p = pInitial;
		gammaRegularized = gammaRegularizedInitial;
		sum = newSum;

		/*
		 * Sum of terms above summationSplitIndex, forward, using Kahan summation.
		 * (to prevent infinity loops the sum has a MAXSUMMATION).
		 */
		for(int i=summationSplitIndex; i