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 import org.apache.commons.math.special.Gamma; 025 026 /** 027 * Default implementation of 028 * {@link org.apache.commons.math.distribution.TDistribution}. 029 * 030 * @version $Revision: 925812 $ $Date: 2010-03-21 11:49:31 -0400 (Sun, 21 Mar 2010) $ 031 */ 032 public class TDistributionImpl 033 extends AbstractContinuousDistribution 034 implements TDistribution, Serializable { 035 036 /** 037 * Default inverse cumulative probability accuracy 038 * @since 2.1 039 */ 040 public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9; 041 042 /** Serializable version identifier */ 043 private static final long serialVersionUID = -5852615386664158222L; 044 045 /** The degrees of freedom*/ 046 private double degreesOfFreedom; 047 048 /** Inverse cumulative probability accuracy */ 049 private final double solverAbsoluteAccuracy; 050 051 /** 052 * Create a t distribution using the given degrees of freedom and the 053 * specified inverse cumulative probability absolute accuracy. 054 * 055 * @param degreesOfFreedom the degrees of freedom. 056 * @param inverseCumAccuracy the maximum absolute error in inverse cumulative probability estimates 057 * (defaults to {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}) 058 * @since 2.1 059 */ 060 public TDistributionImpl(double degreesOfFreedom, double inverseCumAccuracy) { 061 super(); 062 setDegreesOfFreedomInternal(degreesOfFreedom); 063 solverAbsoluteAccuracy = inverseCumAccuracy; 064 } 065 066 /** 067 * Create a t distribution using the given degrees of freedom. 068 * @param degreesOfFreedom the degrees of freedom. 069 */ 070 public TDistributionImpl(double degreesOfFreedom) { 071 this(degreesOfFreedom, DEFAULT_INVERSE_ABSOLUTE_ACCURACY); 072 } 073 074 /** 075 * Modify the degrees of freedom. 076 * @param degreesOfFreedom the new degrees of freedom. 077 * @deprecated as of 2.1 (class will become immutable in 3.0) 078 */ 079 @Deprecated 080 public void setDegreesOfFreedom(double degreesOfFreedom) { 081 setDegreesOfFreedomInternal(degreesOfFreedom); 082 } 083 /** 084 * Modify the degrees of freedom. 085 * @param newDegreesOfFreedom the new degrees of freedom. 086 */ 087 private void setDegreesOfFreedomInternal(double newDegreesOfFreedom) { 088 if (newDegreesOfFreedom <= 0.0) { 089 throw MathRuntimeException.createIllegalArgumentException( 090 "degrees of freedom must be positive ({0})", 091 newDegreesOfFreedom); 092 } 093 this.degreesOfFreedom = newDegreesOfFreedom; 094 } 095 096 /** 097 * Access the degrees of freedom. 098 * @return the degrees of freedom. 099 */ 100 public double getDegreesOfFreedom() { 101 return degreesOfFreedom; 102 } 103 104 /** 105 * Returns the probability density for a particular point. 106 * 107 * @param x The point at which the density should be computed. 108 * @return The pdf at point x. 109 * @since 2.1 110 */ 111 @Override 112 public double density(double x) { 113 final double n = degreesOfFreedom; 114 final double nPlus1Over2 = (n + 1) / 2; 115 return Math.exp(Gamma.logGamma(nPlus1Over2) - 0.5 * (Math.log(Math.PI) + Math.log(n)) - 116 Gamma.logGamma(n/2) - nPlus1Over2 * Math.log(1 + x * x /n)); 117 } 118 119 /** 120 * For this distribution, X, this method returns P(X < <code>x</code>). 121 * @param x the value at which the CDF is evaluated. 122 * @return CDF evaluted at <code>x</code>. 123 * @throws MathException if the cumulative probability can not be 124 * computed due to convergence or other numerical errors. 125 */ 126 public double cumulativeProbability(double x) throws MathException{ 127 double ret; 128 if (x == 0.0) { 129 ret = 0.5; 130 } else { 131 double t = 132 Beta.regularizedBeta( 133 degreesOfFreedom / (degreesOfFreedom + (x * x)), 134 0.5 * degreesOfFreedom, 135 0.5); 136 if (x < 0.0) { 137 ret = 0.5 * t; 138 } else { 139 ret = 1.0 - 0.5 * t; 140 } 141 } 142 143 return ret; 144 } 145 146 /** 147 * For this distribution, X, this method returns the critical point x, such 148 * that P(X < x) = <code>p</code>. 149 * <p> 150 * Returns <code>Double.NEGATIVE_INFINITY</code> for p=0 and 151 * <code>Double.POSITIVE_INFINITY</code> for p=1.</p> 152 * 153 * @param p the desired probability 154 * @return x, such that P(X < x) = <code>p</code> 155 * @throws MathException if the inverse cumulative probability can not be 156 * computed due to convergence or other numerical errors. 157 * @throws IllegalArgumentException if <code>p</code> is not a valid 158 * probability. 159 */ 160 @Override 161 public double inverseCumulativeProbability(final double p) 162 throws MathException { 163 if (p == 0) { 164 return Double.NEGATIVE_INFINITY; 165 } 166 if (p == 1) { 167 return Double.POSITIVE_INFINITY; 168 } 169 return super.inverseCumulativeProbability(p); 170 } 171 172 /** 173 * Access the domain value lower bound, based on <code>p</code>, used to 174 * bracket a CDF root. This method is used by 175 * {@link #inverseCumulativeProbability(double)} to find critical values. 176 * 177 * @param p the desired probability for the critical value 178 * @return domain value lower bound, i.e. 179 * P(X < <i>lower bound</i>) < <code>p</code> 180 */ 181 @Override 182 protected double getDomainLowerBound(double p) { 183 return -Double.MAX_VALUE; 184 } 185 186 /** 187 * Access the domain value upper bound, based on <code>p</code>, used to 188 * bracket a CDF root. This method is used by 189 * {@link #inverseCumulativeProbability(double)} to find critical values. 190 * 191 * @param p the desired probability for the critical value 192 * @return domain value upper bound, i.e. 193 * P(X < <i>upper bound</i>) > <code>p</code> 194 */ 195 @Override 196 protected double getDomainUpperBound(double p) { 197 return Double.MAX_VALUE; 198 } 199 200 /** 201 * Access the initial domain value, based on <code>p</code>, used to 202 * bracket a CDF root. This method is used by 203 * {@link #inverseCumulativeProbability(double)} to find critical values. 204 * 205 * @param p the desired probability for the critical value 206 * @return initial domain value 207 */ 208 @Override 209 protected double getInitialDomain(double p) { 210 return 0.0; 211 } 212 213 /** 214 * Return the absolute accuracy setting of the solver used to estimate 215 * inverse cumulative probabilities. 216 * 217 * @return the solver absolute accuracy 218 * @since 2.1 219 */ 220 @Override 221 protected double getSolverAbsoluteAccuracy() { 222 return solverAbsoluteAccuracy; 223 } 224 }