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