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 */ 17 package org.apache.commons.math.distribution; 18 19 import java.io.Serializable; 20 21 import org.apache.commons.math.MathException; 22 import org.apache.commons.math.MathRuntimeException; 23 import org.apache.commons.math.exception.util.LocalizedFormats; 24 import org.apache.commons.math.special.Beta; 25 import org.apache.commons.math.special.Gamma; 26 import org.apache.commons.math.util.FastMath; 27 28 /** 29 * Default implementation of 30 * {@link org.apache.commons.math.distribution.TDistribution}. 31 * 32 * @version $Revision: 1054524 $ $Date: 2011-01-03 05:59:18 +0100 (lun. 03 janv. 2011) $ 33 */ 34 public class TDistributionImpl 35 extends AbstractContinuousDistribution 36 implements TDistribution, Serializable { 37 38 /** 39 * Default inverse cumulative probability accuracy 40 * @since 2.1 41 */ 42 public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9; 43 44 /** Serializable version identifier */ 45 private static final long serialVersionUID = -5852615386664158222L; 46 47 /** The degrees of freedom*/ 48 private double degreesOfFreedom; 49 50 /** Inverse cumulative probability accuracy */ 51 private final double solverAbsoluteAccuracy; 52 53 /** 54 * Create a t distribution using the given degrees of freedom and the 55 * specified inverse cumulative probability absolute accuracy. 56 * 57 * @param degreesOfFreedom the degrees of freedom. 58 * @param inverseCumAccuracy the maximum absolute error in inverse cumulative probability estimates 59 * (defaults to {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}) 60 * @since 2.1 61 */ 62 public TDistributionImpl(double degreesOfFreedom, double inverseCumAccuracy) { 63 super(); 64 setDegreesOfFreedomInternal(degreesOfFreedom); 65 solverAbsoluteAccuracy = inverseCumAccuracy; 66 } 67 68 /** 69 * Create a t distribution using the given degrees of freedom. 70 * @param degreesOfFreedom the degrees of freedom. 71 */ 72 public TDistributionImpl(double degreesOfFreedom) { 73 this(degreesOfFreedom, DEFAULT_INVERSE_ABSOLUTE_ACCURACY); 74 } 75 76 /** 77 * Modify the degrees of freedom. 78 * @param degreesOfFreedom the new degrees of freedom. 79 * @deprecated as of 2.1 (class will become immutable in 3.0) 80 */ 81 @Deprecated 82 public void setDegreesOfFreedom(double degreesOfFreedom) { 83 setDegreesOfFreedomInternal(degreesOfFreedom); 84 } 85 86 /** 87 * Modify the degrees of freedom. 88 * @param newDegreesOfFreedom the new degrees of freedom. 89 */ 90 private void setDegreesOfFreedomInternal(double newDegreesOfFreedom) { 91 if (newDegreesOfFreedom <= 0.0) { 92 throw MathRuntimeException.createIllegalArgumentException( 93 LocalizedFormats.NOT_POSITIVE_DEGREES_OF_FREEDOM, 94 newDegreesOfFreedom); 95 } 96 this.degreesOfFreedom = newDegreesOfFreedom; 97 } 98 99 /** 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 } 304