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 * https://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
18 /*
19 * This is not the original file distributed by the Apache Software Foundation
20 * It has been modified by the Hipparchus project
21 */
22 package org.hipparchus.distribution.continuous;
23
24 import org.hipparchus.exception.LocalizedCoreFormats;
25 import org.hipparchus.exception.MathIllegalArgumentException;
26 import org.hipparchus.special.Beta;
27 import org.hipparchus.special.Gamma;
28 import org.hipparchus.util.FastMath;
29
30 /**
31 * Implementation of Student's t-distribution.
32 */
33 public class TDistribution extends AbstractRealDistribution {
34 /** Serializable version identifier */
35 private static final long serialVersionUID = 20160320L;
36 /** The degrees of freedom. */
37 private final double degreesOfFreedom;
38 /** Static computation factor based on degreesOfFreedom. */
39 private final double factor;
40
41 /**
42 * Create a t distribution using the given degrees of freedom.
43 *
44 * @param degreesOfFreedom Degrees of freedom.
45 * @throws MathIllegalArgumentException if {@code degreesOfFreedom <= 0}
46 */
47 public TDistribution(double degreesOfFreedom)
48 throws MathIllegalArgumentException {
49 this(degreesOfFreedom, DEFAULT_SOLVER_ABSOLUTE_ACCURACY);
50 }
51
52 /**
53 * Create a t distribution using the given degrees of freedom and the
54 * specified inverse cumulative probability absolute accuracy.
55 *
56 * @param degreesOfFreedom Degrees of freedom.
57 * @param inverseCumAccuracy the maximum absolute error in inverse
58 * cumulative probability estimates
59 * (defaults to {@link #DEFAULT_SOLVER_ABSOLUTE_ACCURACY}).
60 * @throws MathIllegalArgumentException if {@code degreesOfFreedom <= 0}
61 */
62 public TDistribution(double degreesOfFreedom, double inverseCumAccuracy)
63 throws MathIllegalArgumentException {
64 super(inverseCumAccuracy);
65
66 if (degreesOfFreedom <= 0) {
67 throw new MathIllegalArgumentException(LocalizedCoreFormats.DEGREES_OF_FREEDOM,
68 degreesOfFreedom);
69 }
70 this.degreesOfFreedom = degreesOfFreedom;
71
72 final double n = degreesOfFreedom;
73 final double nPlus1Over2 = (n + 1) / 2;
74 factor = Gamma.logGamma(nPlus1Over2) -
75 0.5 * (FastMath.log(FastMath.PI) + FastMath.log(n)) -
76 Gamma.logGamma(n / 2);
77 }
78
79 /**
80 * Access the degrees of freedom.
81 *
82 * @return the degrees of freedom.
83 */
84 public double getDegreesOfFreedom() {
85 return degreesOfFreedom;
86 }
87
88 /** {@inheritDoc} */
89 @Override
90 public double density(double x) {
91 return FastMath.exp(logDensity(x));
92 }
93
94 /** {@inheritDoc} */
95 @Override
96 public double logDensity(double x) {
97 final double n = degreesOfFreedom;
98 final double nPlus1Over2 = (n + 1) / 2;
99 return factor - nPlus1Over2 * FastMath.log(1 + x * x / n);
100 }
101
102 /** {@inheritDoc} */
103 @Override
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 * {@inheritDoc}
126 *
127 * For degrees of freedom parameter {@code df}, the mean is
128 * <ul>
129 * <li>if {@code df > 1} then {@code 0},</li>
130 * <li>else undefined ({@code Double.NaN}).</li>
131 * </ul>
132 */
133 @Override
134 public double getNumericalMean() {
135 final double df = getDegreesOfFreedom();
136
137 if (df > 1) {
138 return 0;
139 }
140
141 return Double.NaN;
142 }
143
144 /**
145 * {@inheritDoc}
146 *
147 * For degrees of freedom parameter {@code df}, the variance is
148 * <ul>
149 * <li>if {@code df > 2} then {@code df / (df - 2)},</li>
150 * <li>if {@code 1 < df <= 2} then positive infinity
151 * ({@code Double.POSITIVE_INFINITY}),</li>
152 * <li>else undefined ({@code Double.NaN}).</li>
153 * </ul>
154 */
155 @Override
156 public double getNumericalVariance() {
157 final double df = getDegreesOfFreedom();
158
159 if (df > 2) {
160 return df / (df - 2);
161 }
162
163 if (df > 1 && df <= 2) {
164 return Double.POSITIVE_INFINITY;
165 }
166
167 return Double.NaN;
168 }
169
170 /**
171 * {@inheritDoc}
172 *
173 * The lower bound of the support is always negative infinity no matter the
174 * parameters.
175 *
176 * @return lower bound of the support (always
177 * {@code Double.NEGATIVE_INFINITY})
178 */
179 @Override
180 public double getSupportLowerBound() {
181 return Double.NEGATIVE_INFINITY;
182 }
183
184 /**
185 * {@inheritDoc}
186 *
187 * The upper bound of the support is always positive infinity no matter the
188 * parameters.
189 *
190 * @return upper bound of the support (always
191 * {@code Double.POSITIVE_INFINITY})
192 */
193 @Override
194 public double getSupportUpperBound() {
195 return Double.POSITIVE_INFINITY;
196 }
197
198 /**
199 * {@inheritDoc}
200 *
201 * The support of this distribution is connected.
202 *
203 * @return {@code true}
204 */
205 @Override
206 public boolean isSupportConnected() {
207 return true;
208 }
209 }