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.exception.util.LocalizedFormats; 024 import org.apache.commons.math.special.Beta; 025 import org.apache.commons.math.special.Gamma; 026 import org.apache.commons.math.util.FastMath; 027 028 /** 029 * Default implementation of 030 * {@link org.apache.commons.math.distribution.TDistribution}. 031 * 032 * @version $Revision: 1054524 $ $Date: 2011-01-03 05:59:18 +0100 (lun. 03 janv. 2011) $ 033 */ 034 public class TDistributionImpl 035 extends AbstractContinuousDistribution 036 implements TDistribution, Serializable { 037 038 /** 039 * Default inverse cumulative probability accuracy 040 * @since 2.1 041 */ 042 public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9; 043 044 /** Serializable version identifier */ 045 private static final long serialVersionUID = -5852615386664158222L; 046 047 /** The degrees of freedom*/ 048 private double degreesOfFreedom; 049 050 /** Inverse cumulative probability accuracy */ 051 private final double solverAbsoluteAccuracy; 052 053 /** 054 * Create a t distribution using the given degrees of freedom and the 055 * specified inverse cumulative probability absolute accuracy. 056 * 057 * @param degreesOfFreedom the degrees of freedom. 058 * @param inverseCumAccuracy the maximum absolute error in inverse cumulative probability estimates 059 * (defaults to {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}) 060 * @since 2.1 061 */ 062 public TDistributionImpl(double degreesOfFreedom, double inverseCumAccuracy) { 063 super(); 064 setDegreesOfFreedomInternal(degreesOfFreedom); 065 solverAbsoluteAccuracy = inverseCumAccuracy; 066 } 067 068 /** 069 * Create a t distribution using the given degrees of freedom. 070 * @param degreesOfFreedom the degrees of freedom. 071 */ 072 public TDistributionImpl(double degreesOfFreedom) { 073 this(degreesOfFreedom, DEFAULT_INVERSE_ABSOLUTE_ACCURACY); 074 } 075 076 /** 077 * Modify the degrees of freedom. 078 * @param degreesOfFreedom the new degrees of freedom. 079 * @deprecated as of 2.1 (class will become immutable in 3.0) 080 */ 081 @Deprecated 082 public void setDegreesOfFreedom(double degreesOfFreedom) { 083 setDegreesOfFreedomInternal(degreesOfFreedom); 084 } 085 086 /** 087 * Modify the degrees of freedom. 088 * @param newDegreesOfFreedom the new degrees of freedom. 089 */ 090 private void setDegreesOfFreedomInternal(double newDegreesOfFreedom) { 091 if (newDegreesOfFreedom <= 0.0) { 092 throw MathRuntimeException.createIllegalArgumentException( 093 LocalizedFormats.NOT_POSITIVE_DEGREES_OF_FREEDOM, 094 newDegreesOfFreedom); 095 } 096 this.degreesOfFreedom = newDegreesOfFreedom; 097 } 098 099 /** 100 * Access the degrees of freedom. 101 * @return the degrees of freedom. 102 */ 103 public double getDegreesOfFreedom() { 104 return degreesOfFreedom; 105 } 106 107 /** 108 * Returns the probability density for a particular point. 109 * 110 * @param x The point at which the density should be computed. 111 * @return The pdf at point x. 112 * @since 2.1 113 */ 114 @Override 115 public double density(double x) { 116 final double n = degreesOfFreedom; 117 final double nPlus1Over2 = (n + 1) / 2; 118 return FastMath.exp(Gamma.logGamma(nPlus1Over2) - 0.5 * (FastMath.log(FastMath.PI) + FastMath.log(n)) - 119 Gamma.logGamma(n/2) - nPlus1Over2 * FastMath.log(1 + x * x /n)); 120 } 121 122 /** 123 * For this distribution, X, this method returns P(X < <code>x</code>). 124 * @param x the value at which the CDF is evaluated. 125 * @return CDF evaluated at <code>x</code>. 126 * @throws MathException if the cumulative probability can not be 127 * computed due to convergence or other numerical errors. 128 */ 129 public double cumulativeProbability(double x) throws MathException{ 130 double ret; 131 if (x == 0.0) { 132 ret = 0.5; 133 } else { 134 double t = 135 Beta.regularizedBeta( 136 degreesOfFreedom / (degreesOfFreedom + (x * x)), 137 0.5 * degreesOfFreedom, 138 0.5); 139 if (x < 0.0) { 140 ret = 0.5 * t; 141 } else { 142 ret = 1.0 - 0.5 * t; 143 } 144 } 145 146 return ret; 147 } 148 149 /** 150 * For this distribution, X, this method returns the critical point x, such 151 * that P(X < x) = <code>p</code>. 152 * <p> 153 * Returns <code>Double.NEGATIVE_INFINITY</code> for p=0 and 154 * <code>Double.POSITIVE_INFINITY</code> for p=1.</p> 155 * 156 * @param p the desired probability 157 * @return x, such that P(X < x) = <code>p</code> 158 * @throws MathException if the inverse cumulative probability can not be 159 * computed due to convergence or other numerical errors. 160 * @throws IllegalArgumentException if <code>p</code> is not a valid 161 * probability. 162 */ 163 @Override 164 public double inverseCumulativeProbability(final double p) 165 throws MathException { 166 if (p == 0) { 167 return Double.NEGATIVE_INFINITY; 168 } 169 if (p == 1) { 170 return Double.POSITIVE_INFINITY; 171 } 172 return super.inverseCumulativeProbability(p); 173 } 174 175 /** 176 * Access the domain value lower bound, based on <code>p</code>, used to 177 * bracket a CDF root. This method is used by 178 * {@link #inverseCumulativeProbability(double)} to find critical values. 179 * 180 * @param p the desired probability for the critical value 181 * @return domain value lower bound, i.e. 182 * P(X < <i>lower bound</i>) < <code>p</code> 183 */ 184 @Override 185 protected double getDomainLowerBound(double p) { 186 return -Double.MAX_VALUE; 187 } 188 189 /** 190 * Access the domain value upper bound, based on <code>p</code>, used to 191 * bracket a CDF root. This method is used by 192 * {@link #inverseCumulativeProbability(double)} to find critical values. 193 * 194 * @param p the desired probability for the critical value 195 * @return domain value upper bound, i.e. 196 * P(X < <i>upper bound</i>) > <code>p</code> 197 */ 198 @Override 199 protected double getDomainUpperBound(double p) { 200 return Double.MAX_VALUE; 201 } 202 203 /** 204 * Access the initial domain value, based on <code>p</code>, used to 205 * bracket a CDF root. This method is used by 206 * {@link #inverseCumulativeProbability(double)} to find critical values. 207 * 208 * @param p the desired probability for the critical value 209 * @return initial domain value 210 */ 211 @Override 212 protected double getInitialDomain(double p) { 213 return 0.0; 214 } 215 216 /** 217 * Return the absolute accuracy setting of the solver used to estimate 218 * inverse cumulative probabilities. 219 * 220 * @return the solver absolute accuracy 221 * @since 2.1 222 */ 223 @Override 224 protected double getSolverAbsoluteAccuracy() { 225 return solverAbsoluteAccuracy; 226 } 227 228 /** 229 * Returns the lower bound of the support for the distribution. 230 * 231 * The lower bound of the support is always negative infinity 232 * no matter the parameters. 233 * 234 * @return lower bound of the support (always Double.NEGATIVE_INFINITY) 235 * @since 2.2 236 */ 237 public double getSupportLowerBound() { 238 return Double.NEGATIVE_INFINITY; 239 } 240 241 /** 242 * Returns the upper bound of the support for the distribution. 243 * 244 * The upper bound of the support is always positive infinity 245 * no matter the parameters. 246 * 247 * @return upper bound of the support (always Double.POSITIVE_INFINITY) 248 * @since 2.2 249 */ 250 public double getSupportUpperBound() { 251 return Double.POSITIVE_INFINITY; 252 } 253 254 /** 255 * Returns the mean. 256 * 257 * For degrees of freedom parameter df, the mean is 258 * <ul> 259 * <li>if <code>df > 1</code> then <code>0</code></li> 260 * <li>else <code>undefined</code></li> 261 * </ul> 262 * 263 * @return the mean 264 * @since 2.2 265 */ 266 public double getNumericalMean() { 267 final double df = getDegreesOfFreedom(); 268 269 if (df > 1) { 270 return 0; 271 } 272 273 return Double.NaN; 274 } 275 276 /** 277 * Returns the variance. 278 * 279 * For degrees of freedom parameter df, the variance is 280 * <ul> 281 * <li>if <code>df > 2</code> then <code>df / (df - 2)</code> </li> 282 * <li>if <code>1 < df <= 2</code> then <code>positive infinity</code></li> 283 * <li>else <code>undefined</code></li> 284 * </ul> 285 * 286 * @return the variance 287 * @since 2.2 288 */ 289 public double getNumericalVariance() { 290 final double df = getDegreesOfFreedom(); 291 292 if (df > 2) { 293 return df / (df - 2); 294 } 295 296 if (df > 1 && df <= 2) { 297 return Double.POSITIVE_INFINITY; 298 } 299 300 return Double.NaN; 301 } 302 303 }