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 */ 017package org.apache.commons.statistics.distribution; 018 019import org.apache.commons.numbers.gamma.RegularizedBeta; 020import org.apache.commons.rng.UniformRandomProvider; 021import org.apache.commons.rng.sampling.distribution.TSampler; 022import org.apache.commons.numbers.gamma.Beta; 023import org.apache.commons.numbers.gamma.LogBeta; 024 025/** 026 * Implementation of Student's t-distribution. 027 * 028 * <p>The probability density function of \( X \) is: 029 * 030 * <p>\[ f(x; v) = \frac{\Gamma(\frac{\nu+1}{2})} {\sqrt{\nu\pi}\,\Gamma(\frac{\nu}{2})} \left(1+\frac{t^2}{\nu} \right)^{\!-\frac{\nu+1}{2}} \] 031 * 032 * <p>for \( v > 0 \) the degrees of freedom, 033 * \( \Gamma \) is the gamma function, and 034 * \( x \in (-\infty, \infty) \). 035 * 036 * @see <a href="https://en.wikipedia.org/wiki/Student%27s_t-distribution">Student's t-distribution (Wikipedia)</a> 037 * @see <a href="https://mathworld.wolfram.com/Studentst-Distribution.html">Student's t-distribution (MathWorld)</a> 038 */ 039public abstract class TDistribution extends AbstractContinuousDistribution { 040 /** A standard normal distribution used for calculations. 041 * This is immutable and thread-safe and can be used across instances. */ 042 static final NormalDistribution STANDARD_NORMAL = NormalDistribution.of(0, 1); 043 044 /** The degrees of freedom. */ 045 private final double degreesOfFreedom; 046 047 /** 048 * Specialisation of the T-distribution used when there are infinite degrees of freedom. 049 * In this case the distribution matches a normal distribution. This is used when the 050 * variance is not different from 1.0. 051 * 052 * <p>This delegates all methods to the standard normal distribution. Instances are 053 * allowed to provide access to the degrees of freedom used during construction. 054 */ 055 private static class NormalTDistribution extends TDistribution { 056 057 /** 058 * @param degreesOfFreedom Degrees of freedom. 059 */ 060 NormalTDistribution(double degreesOfFreedom) { 061 super(degreesOfFreedom); 062 } 063 064 @Override 065 public double density(double x) { 066 return STANDARD_NORMAL.density(x); 067 } 068 069 @Override 070 public double probability(double x0, double x1) { 071 return STANDARD_NORMAL.probability(x0, x1); 072 } 073 074 @Override 075 public double logDensity(double x) { 076 return STANDARD_NORMAL.logDensity(x); 077 } 078 079 @Override 080 public double cumulativeProbability(double x) { 081 return STANDARD_NORMAL.cumulativeProbability(x); 082 } 083 084 @Override 085 public double inverseCumulativeProbability(double p) { 086 return STANDARD_NORMAL.inverseCumulativeProbability(p); 087 } 088 089 // Survival probability functions inherit the symmetry operations from the TDistribution 090 091 @Override 092 public double getMean() { 093 return 0; 094 } 095 096 @Override 097 public double getVariance() { 098 return 1.0; 099 } 100 101 @Override 102 public Sampler createSampler(UniformRandomProvider rng) { 103 return STANDARD_NORMAL.createSampler(rng); 104 } 105 } 106 107 /** 108 * Implementation of Student's T-distribution. 109 */ 110 private static class StudentsTDistribution extends TDistribution { 111 /** 2. */ 112 private static final double TWO = 2; 113 /** The threshold for the density function where the 114 * power function base minus 1 is close to zero. */ 115 private static final double CLOSE_TO_ZERO = 0.25; 116 117 /** -(v + 1) / 2, where v = degrees of freedom. */ 118 private final double mvp1Over2; 119 /** Density normalisation factor, sqrt(v) * beta(1/2, v/2), where v = degrees of freedom. */ 120 private final double densityNormalisation; 121 /** Log density normalisation term, 0.5 * log(v) + log(beta(1/2, v/2)), where v = degrees of freedom. */ 122 private final double logDensityNormalisation; 123 /** Cached value for inverse probability function. */ 124 private final double mean; 125 /** Cached value for inverse probability function. */ 126 private final double variance; 127 128 /** 129 * @param degreesOfFreedom Degrees of freedom. 130 * @param variance Precomputed variance 131 */ 132 StudentsTDistribution(double degreesOfFreedom, double variance) { 133 super(degreesOfFreedom); 134 135 mvp1Over2 = -0.5 * (degreesOfFreedom + 1); 136 densityNormalisation = Math.sqrt(degreesOfFreedom) * Beta.value(0.5, degreesOfFreedom / 2); 137 logDensityNormalisation = 0.5 * Math.log(degreesOfFreedom) + LogBeta.value(0.5, degreesOfFreedom / 2); 138 mean = degreesOfFreedom > 1 ? 0 : Double.NaN; 139 this.variance = variance; 140 } 141 142 /** 143 * @param degreesOfFreedom Degrees of freedom. 144 * @return the variance 145 */ 146 static double computeVariance(double degreesOfFreedom) { 147 if (degreesOfFreedom == Double.POSITIVE_INFINITY) { 148 return 1; 149 } else if (degreesOfFreedom > TWO) { 150 return degreesOfFreedom / (degreesOfFreedom - 2); 151 } else if (degreesOfFreedom > 1) { 152 return Double.POSITIVE_INFINITY; 153 } else { 154 return Double.NaN; 155 } 156 } 157 158 @Override 159 public double density(double x) { 160 // 1 -(v+1)/2 161 // ------------------- * (1 + t^2/v) 162 // sqrt(v) B(1/2, v/2) 163 164 final double t2OverV = x * x / getDegreesOfFreedom(); 165 if (t2OverV < CLOSE_TO_ZERO) { 166 // Avoid power function when the base is close to 1 167 return Math.exp(Math.log1p(t2OverV) * mvp1Over2) / densityNormalisation; 168 } 169 return Math.pow(1 + t2OverV, mvp1Over2) / densityNormalisation; 170 } 171 172 @Override 173 public double logDensity(double x) { 174 return Math.log1p(x * x / getDegreesOfFreedom()) * mvp1Over2 - logDensityNormalisation; 175 } 176 177 @Override 178 public double cumulativeProbability(double x) { 179 if (x == 0) { 180 return 0.5; 181 } 182 final double v = getDegreesOfFreedom(); 183 184 // cdf(t) = 1 - 0.5 * I_x(t)(v/2, 1/2) 185 // where x(t) = v / (v + t^2) 186 // 187 // When t^2 << v using the regularized beta results 188 // in loss of precision. Use the complement instead: 189 // I[x](a,b) = 1 - I[1-x](b,a) 190 // x = v / (v + t^2) 191 // 1-x = t^2 / (v + t^2) 192 // Use the threshold documented in the Boost t-distribution: 193 // https://www.boost.org/doc/libs/1_78_0/libs/math/doc/html/math_toolkit/dist_ref/dists/students_t_dist.html 194 195 final double t2 = x * x; 196 double z; 197 if (v < 2 * t2) { 198 z = RegularizedBeta.value(v / (v + t2), v / 2, 0.5) / 2; 199 } else { 200 z = RegularizedBeta.complement(t2 / (v + t2), 0.5, v / 2) / 2; 201 } 202 return x > 0 ? 1 - z : z; 203 } 204 205 @Override 206 public double getMean() { 207 return mean; 208 } 209 210 @Override 211 public double getVariance() { 212 return variance; 213 } 214 215 @Override 216 double getMedian() { 217 // Overridden for the probability(double, double) method. 218 // This is intentionally not a public method. 219 return 0; 220 } 221 222 @Override 223 public Sampler createSampler(UniformRandomProvider rng) { 224 // T distribution sampler. 225 return TSampler.of(rng, getDegreesOfFreedom())::sample; 226 } 227 } 228 229 /** 230 * @param degreesOfFreedom Degrees of freedom. 231 */ 232 private TDistribution(double degreesOfFreedom) { 233 this.degreesOfFreedom = degreesOfFreedom; 234 } 235 236 /** 237 * Creates a Student's t-distribution. 238 * 239 * @param degreesOfFreedom Degrees of freedom. 240 * @return the distribution 241 * @throws IllegalArgumentException if {@code degreesOfFreedom <= 0} 242 */ 243 public static TDistribution of(double degreesOfFreedom) { 244 if (degreesOfFreedom <= 0) { 245 throw new DistributionException(DistributionException.NOT_STRICTLY_POSITIVE, 246 degreesOfFreedom); 247 } 248 // If the variance converges to 1 use a NormalDistribution. 249 // Occurs at 2^55 = 3.60e16 250 final double variance = StudentsTDistribution.computeVariance(degreesOfFreedom); 251 if (variance == 1) { 252 return new NormalTDistribution(degreesOfFreedom); 253 } 254 return new StudentsTDistribution(degreesOfFreedom, variance); 255 } 256 257 /** 258 * Gets the degrees of freedom parameter of this distribution. 259 * 260 * @return the degrees of freedom. 261 */ 262 public double getDegreesOfFreedom() { 263 return degreesOfFreedom; 264 } 265 266 /** {@inheritDoc} */ 267 @Override 268 public double survivalProbability(double x) { 269 // Exploit symmetry 270 return cumulativeProbability(-x); 271 } 272 273 /** {@inheritDoc} */ 274 @Override 275 public double inverseSurvivalProbability(double p) { 276 // Exploit symmetry 277 return -inverseCumulativeProbability(p); 278 } 279 280 /** 281 * {@inheritDoc} 282 * 283 * <p>For degrees of freedom parameter \( v \), the mean is: 284 * 285 * <p>\[ \mathbb{E}[X] = \begin{cases} 286 * 0 & \text{for } v \gt 1 \\ 287 * \text{undefined} & \text{otherwise} 288 * \end{cases} \] 289 * 290 * @return the mean, or {@link Double#NaN NaN} if it is not defined. 291 */ 292 @Override 293 public abstract double getMean(); 294 295 /** 296 * {@inheritDoc} 297 * 298 * <p>For degrees of freedom parameter \( v \), the variance is: 299 * 300 * <p>\[ \operatorname{var}[X] = \begin{cases} 301 * \frac{v}{v - 2} & \text{for } v \gt 2 \\ 302 * \infty & \text{for } 1 \lt v \le 2 \\ 303 * \text{undefined} & \text{otherwise} 304 * \end{cases} \] 305 * 306 * @return the variance, or {@link Double#NaN NaN} if it is not defined. 307 */ 308 @Override 309 public abstract double getVariance(); 310 311 /** 312 * {@inheritDoc} 313 * 314 * <p>The lower bound of the support is always negative infinity. 315 * 316 * @return {@link Double#NEGATIVE_INFINITY negative infinity}. 317 */ 318 @Override 319 public double getSupportLowerBound() { 320 return Double.NEGATIVE_INFINITY; 321 } 322 323 /** 324 * {@inheritDoc} 325 * 326 * <p>The upper bound of the support is always positive infinity. 327 * 328 * @return {@link Double#POSITIVE_INFINITY positive infinity}. 329 */ 330 @Override 331 public double getSupportUpperBound() { 332 return Double.POSITIVE_INFINITY; 333 } 334}