001 /* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017 018 package org.apache.commons.math3.distribution; 019 020 import org.apache.commons.math3.exception.NotStrictlyPositiveException; 021 import org.apache.commons.math3.exception.util.LocalizedFormats; 022 import org.apache.commons.math3.special.Beta; 023 import org.apache.commons.math3.util.FastMath; 024 import org.apache.commons.math3.random.RandomGenerator; 025 import org.apache.commons.math3.random.Well19937c; 026 027 /** 028 * Implementation of the F-distribution. 029 * 030 * @see <a href="http://en.wikipedia.org/wiki/F-distribution">F-distribution (Wikipedia)</a> 031 * @see <a href="http://mathworld.wolfram.com/F-Distribution.html">F-distribution (MathWorld)</a> 032 * @version $Id: FDistribution.java 1416643 2012-12-03 19:37:14Z tn $ 033 */ 034 public class FDistribution extends AbstractRealDistribution { 035 /** 036 * Default inverse cumulative probability accuracy. 037 * @since 2.1 038 */ 039 public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9; 040 /** Serializable version identifier. */ 041 private static final long serialVersionUID = -8516354193418641566L; 042 /** The numerator degrees of freedom. */ 043 private final double numeratorDegreesOfFreedom; 044 /** The numerator degrees of freedom. */ 045 private final double denominatorDegreesOfFreedom; 046 /** Inverse cumulative probability accuracy. */ 047 private final double solverAbsoluteAccuracy; 048 /** Cached numerical variance */ 049 private double numericalVariance = Double.NaN; 050 /** Whether or not the numerical variance has been calculated */ 051 private boolean numericalVarianceIsCalculated = false; 052 053 /** 054 * Creates an F distribution using the given degrees of freedom. 055 * 056 * @param numeratorDegreesOfFreedom Numerator degrees of freedom. 057 * @param denominatorDegreesOfFreedom Denominator degrees of freedom. 058 * @throws NotStrictlyPositiveException if 059 * {@code numeratorDegreesOfFreedom <= 0} or 060 * {@code denominatorDegreesOfFreedom <= 0}. 061 */ 062 public FDistribution(double numeratorDegreesOfFreedom, 063 double denominatorDegreesOfFreedom) 064 throws NotStrictlyPositiveException { 065 this(numeratorDegreesOfFreedom, denominatorDegreesOfFreedom, 066 DEFAULT_INVERSE_ABSOLUTE_ACCURACY); 067 } 068 069 /** 070 * Creates an F distribution using the given degrees of freedom 071 * and inverse cumulative probability accuracy. 072 * 073 * @param numeratorDegreesOfFreedom Numerator degrees of freedom. 074 * @param denominatorDegreesOfFreedom Denominator degrees of freedom. 075 * @param inverseCumAccuracy the maximum absolute error in inverse 076 * cumulative probability estimates. 077 * @throws NotStrictlyPositiveException if 078 * {@code numeratorDegreesOfFreedom <= 0} or 079 * {@code denominatorDegreesOfFreedom <= 0}. 080 * @since 2.1 081 */ 082 public FDistribution(double numeratorDegreesOfFreedom, 083 double denominatorDegreesOfFreedom, 084 double inverseCumAccuracy) 085 throws NotStrictlyPositiveException { 086 this(new Well19937c(), numeratorDegreesOfFreedom, 087 denominatorDegreesOfFreedom, inverseCumAccuracy); 088 } 089 090 /** 091 * Creates an F distribution. 092 * 093 * @param rng Random number generator. 094 * @param numeratorDegreesOfFreedom Numerator degrees of freedom. 095 * @param denominatorDegreesOfFreedom Denominator degrees of freedom. 096 * @param inverseCumAccuracy the maximum absolute error in inverse 097 * cumulative probability estimates. 098 * @throws NotStrictlyPositiveException if 099 * {@code numeratorDegreesOfFreedom <= 0} or 100 * {@code denominatorDegreesOfFreedom <= 0}. 101 * @since 3.1 102 */ 103 public FDistribution(RandomGenerator rng, 104 double numeratorDegreesOfFreedom, 105 double denominatorDegreesOfFreedom, 106 double inverseCumAccuracy) 107 throws NotStrictlyPositiveException { 108 super(rng); 109 110 if (numeratorDegreesOfFreedom <= 0) { 111 throw new NotStrictlyPositiveException(LocalizedFormats.DEGREES_OF_FREEDOM, 112 numeratorDegreesOfFreedom); 113 } 114 if (denominatorDegreesOfFreedom <= 0) { 115 throw new NotStrictlyPositiveException(LocalizedFormats.DEGREES_OF_FREEDOM, 116 denominatorDegreesOfFreedom); 117 } 118 this.numeratorDegreesOfFreedom = numeratorDegreesOfFreedom; 119 this.denominatorDegreesOfFreedom = denominatorDegreesOfFreedom; 120 solverAbsoluteAccuracy = inverseCumAccuracy; 121 } 122 123 /** 124 * {@inheritDoc} 125 * 126 * @since 2.1 127 */ 128 public double density(double x) { 129 final double nhalf = numeratorDegreesOfFreedom / 2; 130 final double mhalf = denominatorDegreesOfFreedom / 2; 131 final double logx = FastMath.log(x); 132 final double logn = FastMath.log(numeratorDegreesOfFreedom); 133 final double logm = FastMath.log(denominatorDegreesOfFreedom); 134 final double lognxm = FastMath.log(numeratorDegreesOfFreedom * x + 135 denominatorDegreesOfFreedom); 136 return FastMath.exp(nhalf * logn + nhalf * logx - logx + 137 mhalf * logm - nhalf * lognxm - mhalf * lognxm - 138 Beta.logBeta(nhalf, mhalf)); 139 } 140 141 /** 142 * {@inheritDoc} 143 * 144 * The implementation of this method is based on 145 * <ul> 146 * <li> 147 * <a href="http://mathworld.wolfram.com/F-Distribution.html"> 148 * F-Distribution</a>, equation (4). 149 * </li> 150 * </ul> 151 */ 152 public double cumulativeProbability(double x) { 153 double ret; 154 if (x <= 0) { 155 ret = 0; 156 } else { 157 double n = numeratorDegreesOfFreedom; 158 double m = denominatorDegreesOfFreedom; 159 160 ret = Beta.regularizedBeta((n * x) / (m + n * x), 161 0.5 * n, 162 0.5 * m); 163 } 164 return ret; 165 } 166 167 /** 168 * Access the numerator degrees of freedom. 169 * 170 * @return the numerator degrees of freedom. 171 */ 172 public double getNumeratorDegreesOfFreedom() { 173 return numeratorDegreesOfFreedom; 174 } 175 176 /** 177 * Access the denominator degrees of freedom. 178 * 179 * @return the denominator degrees of freedom. 180 */ 181 public double getDenominatorDegreesOfFreedom() { 182 return denominatorDegreesOfFreedom; 183 } 184 185 /** {@inheritDoc} */ 186 @Override 187 protected double getSolverAbsoluteAccuracy() { 188 return solverAbsoluteAccuracy; 189 } 190 191 /** 192 * {@inheritDoc} 193 * 194 * For denominator degrees of freedom parameter {@code b}, the mean is 195 * <ul> 196 * <li>if {@code b > 2} then {@code b / (b - 2)},</li> 197 * <li>else undefined ({@code Double.NaN}). 198 * </ul> 199 */ 200 public double getNumericalMean() { 201 final double denominatorDF = getDenominatorDegreesOfFreedom(); 202 203 if (denominatorDF > 2) { 204 return denominatorDF / (denominatorDF - 2); 205 } 206 207 return Double.NaN; 208 } 209 210 /** 211 * {@inheritDoc} 212 * 213 * For numerator degrees of freedom parameter {@code a} and denominator 214 * degrees of freedom parameter {@code b}, the variance is 215 * <ul> 216 * <li> 217 * if {@code b > 4} then 218 * {@code [2 * b^2 * (a + b - 2)] / [a * (b - 2)^2 * (b - 4)]}, 219 * </li> 220 * <li>else undefined ({@code Double.NaN}). 221 * </ul> 222 */ 223 public double getNumericalVariance() { 224 if (!numericalVarianceIsCalculated) { 225 numericalVariance = calculateNumericalVariance(); 226 numericalVarianceIsCalculated = true; 227 } 228 return numericalVariance; 229 } 230 231 /** 232 * used by {@link #getNumericalVariance()} 233 * 234 * @return the variance of this distribution 235 */ 236 protected double calculateNumericalVariance() { 237 final double denominatorDF = getDenominatorDegreesOfFreedom(); 238 239 if (denominatorDF > 4) { 240 final double numeratorDF = getNumeratorDegreesOfFreedom(); 241 final double denomDFMinusTwo = denominatorDF - 2; 242 243 return ( 2 * (denominatorDF * denominatorDF) * (numeratorDF + denominatorDF - 2) ) / 244 ( (numeratorDF * (denomDFMinusTwo * denomDFMinusTwo) * (denominatorDF - 4)) ); 245 } 246 247 return Double.NaN; 248 } 249 250 /** 251 * {@inheritDoc} 252 * 253 * The lower bound of the support is always 0 no matter the parameters. 254 * 255 * @return lower bound of the support (always 0) 256 */ 257 public double getSupportLowerBound() { 258 return 0; 259 } 260 261 /** 262 * {@inheritDoc} 263 * 264 * The upper bound of the support is always positive infinity 265 * no matter the parameters. 266 * 267 * @return upper bound of the support (always Double.POSITIVE_INFINITY) 268 */ 269 public double getSupportUpperBound() { 270 return Double.POSITIVE_INFINITY; 271 } 272 273 /** {@inheritDoc} */ 274 public boolean isSupportLowerBoundInclusive() { 275 return false; 276 } 277 278 /** {@inheritDoc} */ 279 public boolean isSupportUpperBoundInclusive() { 280 return false; 281 } 282 283 /** 284 * {@inheritDoc} 285 * 286 * The support of this distribution is connected. 287 * 288 * @return {@code true} 289 */ 290 public boolean isSupportConnected() { 291 return true; 292 } 293 }