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 }