| /* |
| * Licensed to the Apache Software Foundation (ASF) under one or more |
| * contributor license agreements. See the NOTICE file distributed with |
| * this work for additional information regarding copyright ownership. |
| * The ASF licenses this file to You under the Apache License, Version 2.0 |
| * (the "License"); you may not use this file except in compliance with |
| * the License. You may obtain a copy of the License at |
| * |
| * http://www.apache.org/licenses/LICENSE-2.0 |
| * |
| * Unless required by applicable law or agreed to in writing, software |
| * distributed under the License is distributed on an "AS IS" BASIS, |
| * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
| * See the License for the specific language governing permissions and |
| * limitations under the License. |
| */ |
| |
| package org.apache.commons.math3.distribution; |
| |
| import java.io.Serializable; |
| import java.math.BigDecimal; |
| |
| import org.apache.commons.math3.exception.MathArithmeticException; |
| import org.apache.commons.math3.exception.NotStrictlyPositiveException; |
| import org.apache.commons.math3.exception.NumberIsTooLargeException; |
| import org.apache.commons.math3.exception.util.LocalizedFormats; |
| import org.apache.commons.math3.fraction.BigFraction; |
| import org.apache.commons.math3.fraction.BigFractionField; |
| import org.apache.commons.math3.fraction.FractionConversionException; |
| import org.apache.commons.math3.linear.Array2DRowFieldMatrix; |
| import org.apache.commons.math3.linear.Array2DRowRealMatrix; |
| import org.apache.commons.math3.linear.FieldMatrix; |
| import org.apache.commons.math3.linear.RealMatrix; |
| import org.apache.commons.math3.util.FastMath; |
| |
| /** |
| * Implementation of the Kolmogorov-Smirnov distribution. |
| * |
| * <p> |
| * Treats the distribution of the two-sided {@code P(D_n < d)} where |
| * {@code D_n = sup_x |G(x) - G_n (x)|} for the theoretical cdf {@code G} and |
| * the empirical cdf {@code G_n}. |
| * </p> |
| * <p> |
| * This implementation is based on [1] with certain quick decisions for extreme |
| * values given in [2]. |
| * </p> |
| * <p> |
| * In short, when wanting to evaluate {@code P(D_n < d)}, the method in [1] is |
| * to write {@code d = (k - h) / n} for positive integer {@code k} and |
| * {@code 0 <= h < 1}. Then {@code P(D_n < d) = (n! / n^n) * t_kk}, where |
| * {@code t_kk} is the {@code (k, k)}'th entry in the special matrix |
| * {@code H^n}, i.e. {@code H} to the {@code n}'th power. |
| * </p> |
| * <p> |
| * References: |
| * <ul> |
| * <li>[1] <a href="http://www.jstatsoft.org/v08/i18/"> |
| * Evaluating Kolmogorov's Distribution</a> by George Marsaglia, Wai |
| * Wan Tsang, and Jingbo Wang</li> |
| * <li>[2] <a href="http://www.jstatsoft.org/v39/i11/"> |
| * Computing the Two-Sided Kolmogorov-Smirnov Distribution</a> by Richard Simard |
| * and Pierre L'Ecuyer</li> |
| * </ul> |
| * Note that [1] contains an error in computing h, refer to |
| * <a href="https://issues.apache.org/jira/browse/MATH-437">MATH-437</a> for details. |
| * </p> |
| * |
| * @see <a href="http://en.wikipedia.org/wiki/Kolmogorov-Smirnov_test"> |
| * Kolmogorov-Smirnov test (Wikipedia)</a> |
| * @deprecated to be removed in version 4.0 - |
| * use {@link org.apache.commons.math3.stat.inference.KolmogorovSmirnovTest} |
| */ |
| public class KolmogorovSmirnovDistribution implements Serializable { |
| |
| /** Serializable version identifier. */ |
| private static final long serialVersionUID = -4670676796862967187L; |
| |
| /** Number of observations. */ |
| private int n; |
| |
| /** |
| * @param n Number of observations |
| * @throws NotStrictlyPositiveException if {@code n <= 0} |
| */ |
| public KolmogorovSmirnovDistribution(int n) |
| throws NotStrictlyPositiveException { |
| if (n <= 0) { |
| throw new NotStrictlyPositiveException(LocalizedFormats.NOT_POSITIVE_NUMBER_OF_SAMPLES, n); |
| } |
| |
| this.n = n; |
| } |
| |
| /** |
| * Calculates {@code P(D_n < d)} using method described in [1] with quick |
| * decisions for extreme values given in [2] (see above). The result is not |
| * exact as with |
| * {@link KolmogorovSmirnovDistribution#cdfExact(double)} because |
| * calculations are based on {@code double} rather than |
| * {@link org.apache.commons.math3.fraction.BigFraction}. |
| * |
| * @param d statistic |
| * @return the two-sided probability of {@code P(D_n < d)} |
| * @throws MathArithmeticException if algorithm fails to convert {@code h} |
| * to a {@link org.apache.commons.math3.fraction.BigFraction} in expressing |
| * {@code d} as {@code (k - h) / m} for integer {@code k, m} and |
| * {@code 0 <= h < 1}. |
| */ |
| public double cdf(double d) throws MathArithmeticException { |
| return this.cdf(d, false); |
| } |
| |
| /** |
| * Calculates {@code P(D_n < d)} using method described in [1] with quick |
| * decisions for extreme values given in [2] (see above). The result is |
| * exact in the sense that BigFraction/BigReal is used everywhere at the |
| * expense of very slow execution time. Almost never choose this in real |
| * applications unless you are very sure; this is almost solely for |
| * verification purposes. Normally, you would choose |
| * {@link KolmogorovSmirnovDistribution#cdf(double)} |
| * |
| * @param d statistic |
| * @return the two-sided probability of {@code P(D_n < d)} |
| * @throws MathArithmeticException if algorithm fails to convert {@code h} |
| * to a {@link org.apache.commons.math3.fraction.BigFraction} in expressing |
| * {@code d} as {@code (k - h) / m} for integer {@code k, m} and |
| * {@code 0 <= h < 1}. |
| */ |
| public double cdfExact(double d) throws MathArithmeticException { |
| return this.cdf(d, true); |
| } |
| |
| /** |
| * Calculates {@code P(D_n < d)} using method described in [1] with quick |
| * decisions for extreme values given in [2] (see above). |
| * |
| * @param d statistic |
| * @param exact whether the probability should be calculated exact using |
| * {@link org.apache.commons.math3.fraction.BigFraction} everywhere at the |
| * expense of very slow execution time, or if {@code double} should be used |
| * convenient places to gain speed. Almost never choose {@code true} in real |
| * applications unless you are very sure; {@code true} is almost solely for |
| * verification purposes. |
| * @return the two-sided probability of {@code P(D_n < d)} |
| * @throws MathArithmeticException if algorithm fails to convert {@code h} |
| * to a {@link org.apache.commons.math3.fraction.BigFraction} in expressing |
| * {@code d} as {@code (k - h) / m} for integer {@code k, m} and |
| * {@code 0 <= h < 1}. |
| */ |
| public double cdf(double d, boolean exact) throws MathArithmeticException { |
| |
| final double ninv = 1 / ((double) n); |
| final double ninvhalf = 0.5 * ninv; |
| |
| if (d <= ninvhalf) { |
| |
| return 0; |
| |
| } else if (ninvhalf < d && d <= ninv) { |
| |
| double res = 1; |
| double f = 2 * d - ninv; |
| |
| // n! f^n = n*f * (n-1)*f * ... * 1*x |
| for (int i = 1; i <= n; ++i) { |
| res *= i * f; |
| } |
| |
| return res; |
| |
| } else if (1 - ninv <= d && d < 1) { |
| |
| return 1 - 2 * FastMath.pow(1 - d, n); |
| |
| } else if (1 <= d) { |
| |
| return 1; |
| } |
| |
| return exact ? exactK(d) : roundedK(d); |
| } |
| |
| /** |
| * Calculates the exact value of {@code P(D_n < d)} using method described |
| * in [1] and {@link org.apache.commons.math3.fraction.BigFraction} (see |
| * above). |
| * |
| * @param d statistic |
| * @return the two-sided probability of {@code P(D_n < d)} |
| * @throws MathArithmeticException if algorithm fails to convert {@code h} |
| * to a {@link org.apache.commons.math3.fraction.BigFraction} in expressing |
| * {@code d} as {@code (k - h) / m} for integer {@code k, m} and |
| * {@code 0 <= h < 1}. |
| */ |
| private double exactK(double d) throws MathArithmeticException { |
| |
| final int k = (int) FastMath.ceil(n * d); |
| |
| final FieldMatrix<BigFraction> H = this.createH(d); |
| final FieldMatrix<BigFraction> Hpower = H.power(n); |
| |
| BigFraction pFrac = Hpower.getEntry(k - 1, k - 1); |
| |
| for (int i = 1; i <= n; ++i) { |
| pFrac = pFrac.multiply(i).divide(n); |
| } |
| |
| /* |
| * BigFraction.doubleValue converts numerator to double and the |
| * denominator to double and divides afterwards. That gives NaN quite |
| * easy. This does not (scale is the number of digits): |
| */ |
| return pFrac.bigDecimalValue(20, BigDecimal.ROUND_HALF_UP).doubleValue(); |
| } |
| |
| /** |
| * Calculates {@code P(D_n < d)} using method described in [1] and doubles |
| * (see above). |
| * |
| * @param d statistic |
| * @return the two-sided probability of {@code P(D_n < d)} |
| * @throws MathArithmeticException if algorithm fails to convert {@code h} |
| * to a {@link org.apache.commons.math3.fraction.BigFraction} in expressing |
| * {@code d} as {@code (k - h) / m} for integer {@code k, m} and |
| * {@code 0 <= h < 1}. |
| */ |
| private double roundedK(double d) throws MathArithmeticException { |
| |
| final int k = (int) FastMath.ceil(n * d); |
| final FieldMatrix<BigFraction> HBigFraction = this.createH(d); |
| final int m = HBigFraction.getRowDimension(); |
| |
| /* |
| * Here the rounding part comes into play: use |
| * RealMatrix instead of FieldMatrix<BigFraction> |
| */ |
| final RealMatrix H = new Array2DRowRealMatrix(m, m); |
| |
| for (int i = 0; i < m; ++i) { |
| for (int j = 0; j < m; ++j) { |
| H.setEntry(i, j, HBigFraction.getEntry(i, j).doubleValue()); |
| } |
| } |
| |
| final RealMatrix Hpower = H.power(n); |
| |
| double pFrac = Hpower.getEntry(k - 1, k - 1); |
| |
| for (int i = 1; i <= n; ++i) { |
| pFrac *= (double) i / (double) n; |
| } |
| |
| return pFrac; |
| } |
| |
| /*** |
| * Creates {@code H} of size {@code m x m} as described in [1] (see above). |
| * |
| * @param d statistic |
| * @return H matrix |
| * @throws NumberIsTooLargeException if fractional part is greater than 1 |
| * @throws FractionConversionException if algorithm fails to convert |
| * {@code h} to a {@link org.apache.commons.math3.fraction.BigFraction} in |
| * expressing {@code d} as {@code (k - h) / m} for integer {@code k, m} and |
| * {@code 0 <= h < 1}. |
| */ |
| private FieldMatrix<BigFraction> createH(double d) |
| throws NumberIsTooLargeException, FractionConversionException { |
| |
| int k = (int) FastMath.ceil(n * d); |
| |
| int m = 2 * k - 1; |
| double hDouble = k - n * d; |
| |
| if (hDouble >= 1) { |
| throw new NumberIsTooLargeException(hDouble, 1.0, false); |
| } |
| |
| BigFraction h = null; |
| |
| try { |
| h = new BigFraction(hDouble, 1.0e-20, 10000); |
| } catch (FractionConversionException e1) { |
| try { |
| h = new BigFraction(hDouble, 1.0e-10, 10000); |
| } catch (FractionConversionException e2) { |
| h = new BigFraction(hDouble, 1.0e-5, 10000); |
| } |
| } |
| |
| final BigFraction[][] Hdata = new BigFraction[m][m]; |
| |
| /* |
| * Start by filling everything with either 0 or 1. |
| */ |
| for (int i = 0; i < m; ++i) { |
| for (int j = 0; j < m; ++j) { |
| if (i - j + 1 < 0) { |
| Hdata[i][j] = BigFraction.ZERO; |
| } else { |
| Hdata[i][j] = BigFraction.ONE; |
| } |
| } |
| } |
| |
| /* |
| * Setting up power-array to avoid calculating the same value twice: |
| * hPowers[0] = h^1 ... hPowers[m-1] = h^m |
| */ |
| final BigFraction[] hPowers = new BigFraction[m]; |
| hPowers[0] = h; |
| for (int i = 1; i < m; ++i) { |
| hPowers[i] = h.multiply(hPowers[i - 1]); |
| } |
| |
| /* |
| * First column and last row has special values (each other reversed). |
| */ |
| for (int i = 0; i < m; ++i) { |
| Hdata[i][0] = Hdata[i][0].subtract(hPowers[i]); |
| Hdata[m - 1][i] = Hdata[m - 1][i].subtract(hPowers[m - i - 1]); |
| } |
| |
| /* |
| * [1] states: "For 1/2 < h < 1 the bottom left element of the matrix |
| * should be (1 - 2*h^m + (2h - 1)^m )/m!" Since 0 <= h < 1, then if h > |
| * 1/2 is sufficient to check: |
| */ |
| if (h.compareTo(BigFraction.ONE_HALF) == 1) { |
| Hdata[m - 1][0] = Hdata[m - 1][0].add(h.multiply(2).subtract(1).pow(m)); |
| } |
| |
| /* |
| * Aside from the first column and last row, the (i, j)-th element is |
| * 1/(i - j + 1)! if i - j + 1 >= 0, else 0. 1's and 0's are already |
| * put, so only division with (i - j + 1)! is needed in the elements |
| * that have 1's. There is no need to calculate (i - j + 1)! and then |
| * divide - small steps avoid overflows. |
| * |
| * Note that i - j + 1 > 0 <=> i + 1 > j instead of j'ing all the way to |
| * m. Also note that it is started at g = 2 because dividing by 1 isn't |
| * really necessary. |
| */ |
| for (int i = 0; i < m; ++i) { |
| for (int j = 0; j < i + 1; ++j) { |
| if (i - j + 1 > 0) { |
| for (int g = 2; g <= i - j + 1; ++g) { |
| Hdata[i][j] = Hdata[i][j].divide(g); |
| } |
| } |
| } |
| } |
| |
| return new Array2DRowFieldMatrix<BigFraction>(BigFractionField.getInstance(), Hdata); |
| } |
| } |
