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 package org.apache.commons.math.distribution; 018 019 import java.io.Serializable; 020 021 import org.apache.commons.math.MathException; 022 import org.apache.commons.math.MathRuntimeException; 023 import org.apache.commons.math.special.Beta; 024 025 /** 026 * Default implementation of 027 * {@link org.apache.commons.math.distribution.FDistribution}. 028 * 029 * @version $Revision: 925897 $ $Date: 2010-03-21 17:06:46 -0400 (Sun, 21 Mar 2010) $ 030 */ 031 public class FDistributionImpl 032 extends AbstractContinuousDistribution 033 implements FDistribution, Serializable { 034 035 /** 036 * Default inverse cumulative probability accuracy 037 * @since 2.1 038 */ 039 public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9; 040 041 /** Message for non positive degrees of freddom. */ 042 private static final String NON_POSITIVE_DEGREES_OF_FREEDOM_MESSAGE = 043 "degrees of freedom must be positive ({0})"; 044 045 /** Serializable version identifier */ 046 private static final long serialVersionUID = -8516354193418641566L; 047 048 /** The numerator degrees of freedom*/ 049 private double numeratorDegreesOfFreedom; 050 051 /** The numerator degrees of freedom*/ 052 private double denominatorDegreesOfFreedom; 053 054 /** Inverse cumulative probability accuracy */ 055 private final double solverAbsoluteAccuracy; 056 057 /** 058 * Create a F distribution using the given degrees of freedom. 059 * @param numeratorDegreesOfFreedom the numerator degrees of freedom. 060 * @param denominatorDegreesOfFreedom the denominator degrees of freedom. 061 */ 062 public FDistributionImpl(double numeratorDegreesOfFreedom, 063 double denominatorDegreesOfFreedom) { 064 this(numeratorDegreesOfFreedom, denominatorDegreesOfFreedom, DEFAULT_INVERSE_ABSOLUTE_ACCURACY); 065 } 066 067 /** 068 * Create a F distribution using the given degrees of freedom and inverse cumulative probability accuracy. 069 * @param numeratorDegreesOfFreedom the numerator degrees of freedom. 070 * @param denominatorDegreesOfFreedom the denominator degrees of freedom. 071 * @param inverseCumAccuracy the maximum absolute error in inverse cumulative probability estimates 072 * (defaults to {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}) 073 * @since 2.1 074 */ 075 public FDistributionImpl(double numeratorDegreesOfFreedom, double denominatorDegreesOfFreedom, 076 double inverseCumAccuracy) { 077 super(); 078 setNumeratorDegreesOfFreedomInternal(numeratorDegreesOfFreedom); 079 setDenominatorDegreesOfFreedomInternal(denominatorDegreesOfFreedom); 080 solverAbsoluteAccuracy = inverseCumAccuracy; 081 } 082 083 /** 084 * Returns the probability density for a particular point. 085 * 086 * @param x The point at which the density should be computed. 087 * @return The pdf at point x. 088 * @since 2.1 089 */ 090 @Override 091 public double density(double x) { 092 final double nhalf = numeratorDegreesOfFreedom / 2; 093 final double mhalf = denominatorDegreesOfFreedom / 2; 094 final double logx = Math.log(x); 095 final double logn = Math.log(numeratorDegreesOfFreedom); 096 final double logm = Math.log(denominatorDegreesOfFreedom); 097 final double lognxm = Math.log(numeratorDegreesOfFreedom * x + denominatorDegreesOfFreedom); 098 return Math.exp(nhalf*logn + nhalf*logx - logx + mhalf*logm - nhalf*lognxm - 099 mhalf*lognxm - Beta.logBeta(nhalf, mhalf)); 100 } 101 102 /** 103 * For this distribution, X, this method returns P(X < x). 104 * 105 * The implementation of this method is based on: 106 * <ul> 107 * <li> 108 * <a href="http://mathworld.wolfram.com/F-Distribution.html"> 109 * F-Distribution</a>, equation (4).</li> 110 * </ul> 111 * 112 * @param x the value at which the CDF is evaluated. 113 * @return CDF for this distribution. 114 * @throws MathException if the cumulative probability can not be 115 * computed due to convergence or other numerical errors. 116 */ 117 public double cumulativeProbability(double x) throws MathException { 118 double ret; 119 if (x <= 0.0) { 120 ret = 0.0; 121 } else { 122 double n = numeratorDegreesOfFreedom; 123 double m = denominatorDegreesOfFreedom; 124 125 ret = Beta.regularizedBeta((n * x) / (m + n * x), 126 0.5 * n, 127 0.5 * m); 128 } 129 return ret; 130 } 131 132 /** 133 * For this distribution, X, this method returns the critical point x, such 134 * that P(X < x) = <code>p</code>. 135 * <p> 136 * Returns 0 for p=0 and <code>Double.POSITIVE_INFINITY</code> for p=1.</p> 137 * 138 * @param p the desired probability 139 * @return x, such that P(X < x) = <code>p</code> 140 * @throws MathException if the inverse cumulative probability can not be 141 * computed due to convergence or other numerical errors. 142 * @throws IllegalArgumentException if <code>p</code> is not a valid 143 * probability. 144 */ 145 @Override 146 public double inverseCumulativeProbability(final double p) 147 throws MathException { 148 if (p == 0) { 149 return 0d; 150 } 151 if (p == 1) { 152 return Double.POSITIVE_INFINITY; 153 } 154 return super.inverseCumulativeProbability(p); 155 } 156 157 /** 158 * Access the domain value lower bound, based on <code>p</code>, used to 159 * bracket a CDF root. This method is used by 160 * {@link #inverseCumulativeProbability(double)} to find critical values. 161 * 162 * @param p the desired probability for the critical value 163 * @return domain value lower bound, i.e. 164 * P(X < <i>lower bound</i>) < <code>p</code> 165 */ 166 @Override 167 protected double getDomainLowerBound(double p) { 168 return 0.0; 169 } 170 171 /** 172 * Access the domain value upper bound, based on <code>p</code>, used to 173 * bracket a CDF root. This method is used by 174 * {@link #inverseCumulativeProbability(double)} to find critical values. 175 * 176 * @param p the desired probability for the critical value 177 * @return domain value upper bound, i.e. 178 * P(X < <i>upper bound</i>) > <code>p</code> 179 */ 180 @Override 181 protected double getDomainUpperBound(double p) { 182 return Double.MAX_VALUE; 183 } 184 185 /** 186 * Access the initial domain value, based on <code>p</code>, used to 187 * bracket a CDF root. This method is used by 188 * {@link #inverseCumulativeProbability(double)} to find critical values. 189 * 190 * @param p the desired probability for the critical value 191 * @return initial domain value 192 */ 193 @Override 194 protected double getInitialDomain(double p) { 195 double ret = 1.0; 196 double d = denominatorDegreesOfFreedom; 197 if (d > 2.0) { 198 // use mean 199 ret = d / (d - 2.0); 200 } 201 return ret; 202 } 203 204 /** 205 * Modify the numerator degrees of freedom. 206 * @param degreesOfFreedom the new numerator degrees of freedom. 207 * @throws IllegalArgumentException if <code>degreesOfFreedom</code> is not 208 * positive. 209 * @deprecated as of 2.1 (class will become immutable in 3.0) 210 */ 211 @Deprecated 212 public void setNumeratorDegreesOfFreedom(double degreesOfFreedom) { 213 setNumeratorDegreesOfFreedomInternal(degreesOfFreedom); 214 } 215 216 /** 217 * Modify the numerator degrees of freedom. 218 * @param degreesOfFreedom the new numerator degrees of freedom. 219 * @throws IllegalArgumentException if <code>degreesOfFreedom</code> is not 220 * positive. 221 */ 222 private void setNumeratorDegreesOfFreedomInternal(double degreesOfFreedom) { 223 if (degreesOfFreedom <= 0.0) { 224 throw MathRuntimeException.createIllegalArgumentException( 225 NON_POSITIVE_DEGREES_OF_FREEDOM_MESSAGE, degreesOfFreedom); 226 } 227 this.numeratorDegreesOfFreedom = degreesOfFreedom; 228 } 229 230 /** 231 * Access the numerator degrees of freedom. 232 * @return the numerator degrees of freedom. 233 */ 234 public double getNumeratorDegreesOfFreedom() { 235 return numeratorDegreesOfFreedom; 236 } 237 238 /** 239 * Modify the denominator degrees of freedom. 240 * @param degreesOfFreedom the new denominator degrees of freedom. 241 * @throws IllegalArgumentException if <code>degreesOfFreedom</code> is not 242 * positive. 243 * @deprecated as of 2.1 (class will become immutable in 3.0) 244 */ 245 @Deprecated 246 public void setDenominatorDegreesOfFreedom(double degreesOfFreedom) { 247 setDenominatorDegreesOfFreedomInternal(degreesOfFreedom); 248 } 249 250 /** 251 * Modify the denominator degrees of freedom. 252 * @param degreesOfFreedom the new denominator degrees of freedom. 253 * @throws IllegalArgumentException if <code>degreesOfFreedom</code> is not 254 * positive. 255 */ 256 private void setDenominatorDegreesOfFreedomInternal(double degreesOfFreedom) { 257 if (degreesOfFreedom <= 0.0) { 258 throw MathRuntimeException.createIllegalArgumentException( 259 NON_POSITIVE_DEGREES_OF_FREEDOM_MESSAGE, degreesOfFreedom); 260 } 261 this.denominatorDegreesOfFreedom = degreesOfFreedom; 262 } 263 264 /** 265 * Access the denominator degrees of freedom. 266 * @return the denominator degrees of freedom. 267 */ 268 public double getDenominatorDegreesOfFreedom() { 269 return denominatorDegreesOfFreedom; 270 } 271 272 /** 273 * Return the absolute accuracy setting of the solver used to estimate 274 * inverse cumulative probabilities. 275 * 276 * @return the solver absolute accuracy 277 * @since 2.1 278 */ 279 @Override 280 protected double getSolverAbsoluteAccuracy() { 281 return solverAbsoluteAccuracy; 282 } 283 }