FDistribution.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */

  18. package org.apache.commons.math3.distribution;

  20. import org.apache.commons.math3.exception.NotStrictlyPositiveException;
  21. import org.apache.commons.math3.exception.util.LocalizedFormats;
  22. import org.apache.commons.math3.random.RandomGenerator;
  23. import org.apache.commons.math3.random.Well19937c;
  24. import org.apache.commons.math3.special.Beta;
  25. import org.apache.commons.math3.util.FastMath;

  27. /**
  28.  * Implementation of the F-distribution.
  29.  *
  30.  * @see <a href="http://en.wikipedia.org/wiki/F-distribution">F-distribution (Wikipedia)</a>
  31.  * @see <a href="http://mathworld.wolfram.com/F-Distribution.html">F-distribution (MathWorld)</a>
  32.  */
  33. public class FDistribution extends AbstractRealDistribution {
  34.     /**
  35.      * Default inverse cumulative probability accuracy.
  36.      * @since 2.1
  37.      */
  38.     public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9;
  39.     /** Serializable version identifier. */
  40.     private static final long serialVersionUID = -8516354193418641566L;
  41.     /** The numerator degrees of freedom. */
  42.     private final double numeratorDegreesOfFreedom;
  43.     /** The numerator degrees of freedom. */
  44.     private final double denominatorDegreesOfFreedom;
  45.     /** Inverse cumulative probability accuracy. */
  46.     private final double solverAbsoluteAccuracy;
  47.     /** Cached numerical variance */
  48.     private double numericalVariance = Double.NaN;
  49.     /** Whether or not the numerical variance has been calculated */
  50.     private boolean numericalVarianceIsCalculated = false;

  52.     /**
  53.      * Creates an F distribution using the given degrees of freedom.
  54.      * <p>
  55.      * <b>Note:</b> this constructor will implicitly create an instance of
  56.      * {@link Well19937c} as random generator to be used for sampling only (see
  57.      * {@link #sample()} and {@link #sample(int)}). In case no sampling is
  58.      * needed for the created distribution, it is advised to pass {@code null}
  59.      * as random generator via the appropriate constructors to avoid the
  60.      * additional initialisation overhead.
  61.      *
  62.      * @param numeratorDegreesOfFreedom Numerator degrees of freedom.
  63.      * @param denominatorDegreesOfFreedom Denominator degrees of freedom.
  64.      * @throws NotStrictlyPositiveException if
  65.      * {@code numeratorDegreesOfFreedom <= 0} or
  66.      * {@code denominatorDegreesOfFreedom <= 0}.
  67.      */
  68.     public FDistribution(double numeratorDegreesOfFreedom,
  69.                          double denominatorDegreesOfFreedom)
  70.         throws NotStrictlyPositiveException {
  71.         this(numeratorDegreesOfFreedom, denominatorDegreesOfFreedom,
  72.              DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
  73.     }

  75.     /**
  76.      * Creates an F distribution using the given degrees of freedom
  77.      * and inverse cumulative probability accuracy.
  78.      * <p>
  79.      * <b>Note:</b> this constructor will implicitly create an instance of
  80.      * {@link Well19937c} as random generator to be used for sampling only (see
  81.      * {@link #sample()} and {@link #sample(int)}). In case no sampling is
  82.      * needed for the created distribution, it is advised to pass {@code null}
  83.      * as random generator via the appropriate constructors to avoid the
  84.      * additional initialisation overhead.
  85.      *
  86.      * @param numeratorDegreesOfFreedom Numerator degrees of freedom.
  87.      * @param denominatorDegreesOfFreedom Denominator degrees of freedom.
  88.      * @param inverseCumAccuracy the maximum absolute error in inverse
  89.      * cumulative probability estimates.
  90.      * @throws NotStrictlyPositiveException if
  91.      * {@code numeratorDegreesOfFreedom <= 0} or
  92.      * {@code denominatorDegreesOfFreedom <= 0}.
  93.      * @since 2.1
  94.      */
  95.     public FDistribution(double numeratorDegreesOfFreedom,
  96.                          double denominatorDegreesOfFreedom,
  97.                          double inverseCumAccuracy)
  98.         throws NotStrictlyPositiveException {
  99.         this(new Well19937c(), numeratorDegreesOfFreedom,
  100.              denominatorDegreesOfFreedom, inverseCumAccuracy);
  101.     }

  103.     /**
  104.      * Creates an F distribution.
  105.      *
  106.      * @param rng Random number generator.
  107.      * @param numeratorDegreesOfFreedom Numerator degrees of freedom.
  108.      * @param denominatorDegreesOfFreedom Denominator degrees of freedom.
  109.      * @throws NotStrictlyPositiveException if {@code numeratorDegreesOfFreedom <= 0} or
  110.      * {@code denominatorDegreesOfFreedom <= 0}.
  111.      * @since 3.3
  112.      */
  113.     public FDistribution(RandomGenerator rng,
  114.                          double numeratorDegreesOfFreedom,
  115.                          double denominatorDegreesOfFreedom)
  116.         throws NotStrictlyPositiveException {
  117.         this(rng, numeratorDegreesOfFreedom, denominatorDegreesOfFreedom, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
  118.     }

  120.     /**
  121.      * Creates an F distribution.
  122.      *
  123.      * @param rng Random number generator.
  124.      * @param numeratorDegreesOfFreedom Numerator degrees of freedom.
  125.      * @param denominatorDegreesOfFreedom Denominator degrees of freedom.
  126.      * @param inverseCumAccuracy the maximum absolute error in inverse
  127.      * cumulative probability estimates.
  128.      * @throws NotStrictlyPositiveException if {@code numeratorDegreesOfFreedom <= 0} or
  129.      * {@code denominatorDegreesOfFreedom <= 0}.
  130.      * @since 3.1
  131.      */
  132.     public FDistribution(RandomGenerator rng,
  133.                          double numeratorDegreesOfFreedom,
  134.                          double denominatorDegreesOfFreedom,
  135.                          double inverseCumAccuracy)
  136.         throws NotStrictlyPositiveException {
  137.         super(rng);

  139.         if (numeratorDegreesOfFreedom <= 0) {
  140.             throw new NotStrictlyPositiveException(LocalizedFormats.DEGREES_OF_FREEDOM,
  141.                                                    numeratorDegreesOfFreedom);
  142.         }
  143.         if (denominatorDegreesOfFreedom <= 0) {
  144.             throw new NotStrictlyPositiveException(LocalizedFormats.DEGREES_OF_FREEDOM,
  145.                                                    denominatorDegreesOfFreedom);
  146.         }
  147.         this.numeratorDegreesOfFreedom = numeratorDegreesOfFreedom;
  148.         this.denominatorDegreesOfFreedom = denominatorDegreesOfFreedom;
  149.         solverAbsoluteAccuracy = inverseCumAccuracy;
  150.     }

  152.     /**
  153.      * {@inheritDoc}
  154.      *
  155.      * @since 2.1
  156.      */
  157.     public double density(double x) {
  158.         return FastMath.exp(logDensity(x));
  159.     }

  161.     /** {@inheritDoc} **/
  162.     @Override
  163.     public double logDensity(double x) {
  164.         final double nhalf = numeratorDegreesOfFreedom / 2;
  165.         final double mhalf = denominatorDegreesOfFreedom / 2;
  166.         final double logx = FastMath.log(x);
  167.         final double logn = FastMath.log(numeratorDegreesOfFreedom);
  168.         final double logm = FastMath.log(denominatorDegreesOfFreedom);
  169.         final double lognxm = FastMath.log(numeratorDegreesOfFreedom * x +
  170.                 denominatorDegreesOfFreedom);
  171.         return nhalf * logn + nhalf * logx - logx +
  172.                mhalf * logm - nhalf * lognxm - mhalf * lognxm -
  173.                Beta.logBeta(nhalf, mhalf);
  174.     }

  176.     /**
  177.      * {@inheritDoc}
  178.      *
  179.      * The implementation of this method is based on
  180.      * <ul>
  181.      *  <li>
  182.      *   <a href="http://mathworld.wolfram.com/F-Distribution.html">
  183.      *   F-Distribution</a>, equation (4).
  184.      *  </li>
  185.      * </ul>
  186.      */
  187.     public double cumulativeProbability(double x)  {
  188.         double ret;
  189.         if (x <= 0) {
  190.             ret = 0;
  191.         } else {
  192.             double n = numeratorDegreesOfFreedom;
  193.             double m = denominatorDegreesOfFreedom;

  195.             ret = Beta.regularizedBeta((n * x) / (m + n * x),
  196.                 0.5 * n,
  197.                 0.5 * m);
  198.         }
  199.         return ret;
  200.     }

  202.     /**
  203.      * Access the numerator degrees of freedom.
  204.      *
  205.      * @return the numerator degrees of freedom.
  206.      */
  207.     public double getNumeratorDegreesOfFreedom() {
  208.         return numeratorDegreesOfFreedom;
  209.     }

  211.     /**
  212.      * Access the denominator degrees of freedom.
  213.      *
  214.      * @return the denominator degrees of freedom.
  215.      */
  216.     public double getDenominatorDegreesOfFreedom() {
  217.         return denominatorDegreesOfFreedom;
  218.     }

  220.     /** {@inheritDoc} */
  221.     @Override
  222.     protected double getSolverAbsoluteAccuracy() {
  223.         return solverAbsoluteAccuracy;
  224.     }

  226.     /**
  227.      * {@inheritDoc}
  228.      *
  229.      * For denominator degrees of freedom parameter {@code b}, the mean is
  230.      * <ul>
  231.      *  <li>if {@code b > 2} then {@code b / (b - 2)},</li>
  232.      *  <li>else undefined ({@code Double.NaN}).
  233.      * </ul>
  234.      */
  235.     public double getNumericalMean() {
  236.         final double denominatorDF = getDenominatorDegreesOfFreedom();

  238.         if (denominatorDF > 2) {
  239.             return denominatorDF / (denominatorDF - 2);
  240.         }

  242.         return Double.NaN;
  243.     }

  245.     /**
  246.      * {@inheritDoc}
  247.      *
  248.      * For numerator degrees of freedom parameter {@code a} and denominator
  249.      * degrees of freedom parameter {@code b}, the variance is
  250.      * <ul>
  251.      *  <li>
  252.      *    if {@code b > 4} then
  253.      *    {@code [2 * b^2 * (a + b - 2)] / [a * (b - 2)^2 * (b - 4)]},
  254.      *  </li>
  255.      *  <li>else undefined ({@code Double.NaN}).
  256.      * </ul>
  257.      */
  258.     public double getNumericalVariance() {
  259.         if (!numericalVarianceIsCalculated) {
  260.             numericalVariance = calculateNumericalVariance();
  261.             numericalVarianceIsCalculated = true;
  262.         }
  263.         return numericalVariance;
  264.     }

  266.     /**
  267.      * used by {@link #getNumericalVariance()}
  268.      *
  269.      * @return the variance of this distribution
  270.      */
  271.     protected double calculateNumericalVariance() {
  272.         final double denominatorDF = getDenominatorDegreesOfFreedom();

  274.         if (denominatorDF > 4) {
  275.             final double numeratorDF = getNumeratorDegreesOfFreedom();
  276.             final double denomDFMinusTwo = denominatorDF - 2;

  278.             return ( 2 * (denominatorDF * denominatorDF) * (numeratorDF + denominatorDF - 2) ) /
  279.                    ( (numeratorDF * (denomDFMinusTwo * denomDFMinusTwo) * (denominatorDF - 4)) );
  280.         }

  282.         return Double.NaN;
  283.     }

  285.     /**
  286.      * {@inheritDoc}
  287.      *
  288.      * The lower bound of the support is always 0 no matter the parameters.
  289.      *
  290.      * @return lower bound of the support (always 0)
  291.      */
  292.     public double getSupportLowerBound() {
  293.         return 0;
  294.     }

  296.     /**
  297.      * {@inheritDoc}
  298.      *
  299.      * The upper bound of the support is always positive infinity
  300.      * no matter the parameters.
  301.      *
  302.      * @return upper bound of the support (always Double.POSITIVE_INFINITY)
  303.      */
  304.     public double getSupportUpperBound() {
  305.         return Double.POSITIVE_INFINITY;
  306.     }

  308.     /** {@inheritDoc} */
  309.     public boolean isSupportLowerBoundInclusive() {
  310.         return false;
  311.     }

  313.     /** {@inheritDoc} */
  314.     public boolean isSupportUpperBoundInclusive() {
  315.         return false;
  316.     }

  318.     /**
  319.      * {@inheritDoc}
  320.      *
  321.      * The support of this distribution is connected.
  322.      *
  323.      * @return {@code true}
  324.      */
  325.     public boolean isSupportConnected() {
  326.         return true;
  327.     }
  328. }
Created with JaCoCo 0.7.2.201409121644