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.math3.distribution;
018
019 import org.apache.commons.math3.exception.NotStrictlyPositiveException;
020 import org.apache.commons.math3.exception.util.LocalizedFormats;
021 import org.apache.commons.math3.special.Beta;
022 import org.apache.commons.math3.special.Gamma;
023 import org.apache.commons.math3.util.FastMath;
024 import org.apache.commons.math3.random.RandomGenerator;
025 import org.apache.commons.math3.random.Well19937c;
026
027 /**
028 * Implementation of Student's t-distribution.
029 *
030 * @see "<a href='http://en.wikipedia.org/wiki/Student's_t-distribution'>Student's t-distribution (Wikipedia)</a>"
031 * @see "<a href='http://mathworld.wolfram.com/Studentst-Distribution.html'>Student's t-distribution (MathWorld)</a>"
032 * @version $Id: TDistribution.java 1416643 2012-12-03 19:37:14Z tn $
033 */
034 public class TDistribution extends AbstractRealDistribution {
035 /**
036 * Default inverse cumulative probability accuracy.
037 * @since 2.1
038 */
039 public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9;
040 /** Serializable version identifier */
041 private static final long serialVersionUID = -5852615386664158222L;
042 /** The degrees of freedom. */
043 private final double degreesOfFreedom;
044 /** Inverse cumulative probability accuracy. */
045 private final double solverAbsoluteAccuracy;
046
047 /**
048 * Create a t distribution using the given degrees of freedom.
049 *
050 * @param degreesOfFreedom Degrees of freedom.
051 * @throws NotStrictlyPositiveException if {@code degreesOfFreedom <= 0}
052 */
053 public TDistribution(double degreesOfFreedom)
054 throws NotStrictlyPositiveException {
055 this(degreesOfFreedom, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
056 }
057
058 /**
059 * Create a t distribution using the given degrees of freedom and the
060 * specified inverse cumulative probability absolute accuracy.
061 *
062 * @param degreesOfFreedom Degrees of freedom.
063 * @param inverseCumAccuracy the maximum absolute error in inverse
064 * cumulative probability estimates
065 * (defaults to {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}).
066 * @throws NotStrictlyPositiveException if {@code degreesOfFreedom <= 0}
067 * @since 2.1
068 */
069 public TDistribution(double degreesOfFreedom, double inverseCumAccuracy)
070 throws NotStrictlyPositiveException {
071 this(new Well19937c(), degreesOfFreedom, inverseCumAccuracy);
072 }
073
074 /**
075 * Creates a t distribution.
076 *
077 * @param rng Random number generator.
078 * @param degreesOfFreedom Degrees of freedom.
079 * @param inverseCumAccuracy the maximum absolute error in inverse
080 * cumulative probability estimates
081 * (defaults to {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}).
082 * @throws NotStrictlyPositiveException if {@code degreesOfFreedom <= 0}
083 * @since 3.1
084 */
085 public TDistribution(RandomGenerator rng,
086 double degreesOfFreedom,
087 double inverseCumAccuracy)
088 throws NotStrictlyPositiveException {
089 super(rng);
090
091 if (degreesOfFreedom <= 0) {
092 throw new NotStrictlyPositiveException(LocalizedFormats.DEGREES_OF_FREEDOM,
093 degreesOfFreedom);
094 }
095 this.degreesOfFreedom = degreesOfFreedom;
096 solverAbsoluteAccuracy = inverseCumAccuracy;
097 }
098
099 /**
100 * Access the degrees of freedom.
101 *
102 * @return the degrees of freedom.
103 */
104 public double getDegreesOfFreedom() {
105 return degreesOfFreedom;
106 }
107
108 /** {@inheritDoc} */
109 public double density(double x) {
110 final double n = degreesOfFreedom;
111 final double nPlus1Over2 = (n + 1) / 2;
112 return FastMath.exp(Gamma.logGamma(nPlus1Over2) -
113 0.5 * (FastMath.log(FastMath.PI) +
114 FastMath.log(n)) -
115 Gamma.logGamma(n / 2) -
116 nPlus1Over2 * FastMath.log(1 + x * x / n));
117 }
118
119 /** {@inheritDoc} */
120 public double cumulativeProbability(double x) {
121 double ret;
122 if (x == 0) {
123 ret = 0.5;
124 } else {
125 double t =
126 Beta.regularizedBeta(
127 degreesOfFreedom / (degreesOfFreedom + (x * x)),
128 0.5 * degreesOfFreedom,
129 0.5);
130 if (x < 0.0) {
131 ret = 0.5 * t;
132 } else {
133 ret = 1.0 - 0.5 * t;
134 }
135 }
136
137 return ret;
138 }
139
140 /** {@inheritDoc} */
141 @Override
142 protected double getSolverAbsoluteAccuracy() {
143 return solverAbsoluteAccuracy;
144 }
145
146 /**
147 * {@inheritDoc}
148 *
149 * For degrees of freedom parameter {@code df}, the mean is
150 * <ul>
151 * <li>if {@code df > 1} then {@code 0},</li>
152 * <li>else undefined ({@code Double.NaN}).</li>
153 * </ul>
154 */
155 public double getNumericalMean() {
156 final double df = getDegreesOfFreedom();
157
158 if (df > 1) {
159 return 0;
160 }
161
162 return Double.NaN;
163 }
164
165 /**
166 * {@inheritDoc}
167 *
168 * For degrees of freedom parameter {@code df}, the variance is
169 * <ul>
170 * <li>if {@code df > 2} then {@code df / (df - 2)},</li>
171 * <li>if {@code 1 < df <= 2} then positive infinity
172 * ({@code Double.POSITIVE_INFINITY}),</li>
173 * <li>else undefined ({@code Double.NaN}).</li>
174 * </ul>
175 */
176 public double getNumericalVariance() {
177 final double df = getDegreesOfFreedom();
178
179 if (df > 2) {
180 return df / (df - 2);
181 }
182
183 if (df > 1 && df <= 2) {
184 return Double.POSITIVE_INFINITY;
185 }
186
187 return Double.NaN;
188 }
189
190 /**
191 * {@inheritDoc}
192 *
193 * The lower bound of the support is always negative infinity no matter the
194 * parameters.
195 *
196 * @return lower bound of the support (always
197 * {@code Double.NEGATIVE_INFINITY})
198 */
199 public double getSupportLowerBound() {
200 return Double.NEGATIVE_INFINITY;
201 }
202
203 /**
204 * {@inheritDoc}
205 *
206 * The upper bound of the support is always positive infinity no matter the
207 * parameters.
208 *
209 * @return upper bound of the support (always
210 * {@code Double.POSITIVE_INFINITY})
211 */
212 public double getSupportUpperBound() {
213 return Double.POSITIVE_INFINITY;
214 }
215
216 /** {@inheritDoc} */
217 public boolean isSupportLowerBoundInclusive() {
218 return false;
219 }
220
221 /** {@inheritDoc} */
222 public boolean isSupportUpperBoundInclusive() {
223 return false;
224 }
225
226 /**
227 * {@inheritDoc}
228 *
229 * The support of this distribution is connected.
230 *
231 * @return {@code true}
232 */
233 public boolean isSupportConnected() {
234 return true;
235 }
236 }