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.rng.sampling.distribution; 018 019import java.util.function.DoubleUnaryOperator; 020import org.apache.commons.rng.UniformRandomProvider; 021 022/** 023 * Sampling from a T distribution. 024 * 025 * <p>Uses Bailey's algorithm for t-distribution sampling:</p> 026 * 027 * <blockquote> 028 * <pre> 029 * Bailey, R. W. (1994) 030 * "Polar Generation of Random Variates with the t-Distribution." 031 * Mathematics of Computation 62, 779-781. 032 * </pre> 033 * </blockquote> 034 * 035 * <p>Sampling uses {@link UniformRandomProvider#nextLong()}.</p> 036 * 037 * @see <a href="https://en.wikipedia.org/wiki/Student%27s_t-distribution">Student's T distribution (wikipedia)</a> 038 * @see <a href="https://doi.org/10.2307/2153537">Mathematics of Computation, 62, 779-781</a> 039 * @since 1.5 040 */ 041public abstract class TSampler implements SharedStateContinuousSampler { 042 /** Threshold for huge degrees of freedom. Above this value the CDF of the t distribution 043 * matches the normal distribution. Value is 2/eps (where eps is the machine epsilon) 044 * or approximately 9.0e15. */ 045 private static final double HUGE_DF = 0x1.0p53; 046 047 /** Source of randomness. */ 048 private final UniformRandomProvider rng; 049 050 /** 051 * Sample from a t-distribution using Bailey's algorithm. 052 */ 053 private static final class StudentsTSampler extends TSampler { 054 /** Threshold for large degrees of freedom. */ 055 private static final double LARGE_DF = 25; 056 /** The multiplier to convert the least significant 53-bits of a {@code long} to a 057 * uniform {@code double}. */ 058 private static final double DOUBLE_MULTIPLIER = 0x1.0p-53; 059 060 /** Degrees of freedom. */ 061 private final double df; 062 /** Function to compute pow(x, -2/v) - 1, where v = degrees of freedom. */ 063 private final DoubleUnaryOperator powm1; 064 065 /** 066 * @param rng Generator of uniformly distributed random numbers. 067 * @param v Degrees of freedom. 068 */ 069 StudentsTSampler(UniformRandomProvider rng, 070 double v) { 071 super(rng); 072 df = v; 073 074 // The sampler requires pow(w, -2/v) - 1 with 075 // 0 <= w <= 1; Expected(w) = sqrt(0.5). 076 // When the exponent is small then pow(x, y) -> 1. 077 // This affects large degrees of freedom. 078 final double exponent = -2 / v; 079 powm1 = v > LARGE_DF ? 080 x -> Math.expm1(Math.log(x) * exponent) : 081 x -> Math.pow(x, exponent) - 1; 082 } 083 084 /** 085 * @param rng Generator of uniformly distributed random numbers. 086 * @param source Source to copy. 087 */ 088 private StudentsTSampler(UniformRandomProvider rng, 089 StudentsTSampler source) { 090 super(rng); 091 df = source.df; 092 powm1 = source.powm1; 093 } 094 095 /** {@inheritDoc} */ 096 @Override 097 public double sample() { 098 // Require u and v in [0, 1] and a random sign. 099 // Create u in (0, 1] to avoid generating nan 100 // from u*u/w (0/0) or r2*c2 (inf*0). 101 final double u = InternalUtils.makeNonZeroDouble(nextLong()); 102 final double v = makeSignedDouble(nextLong()); 103 final double w = u * u + v * v; 104 if (w > 1) { 105 // Rejection frequency = 1 - pi/4 = 0.215. 106 // Recursion will generate stack overflow given a broken RNG 107 // and avoids an infinite loop. 108 return sample(); 109 } 110 // Sidestep a square-root calculation. 111 final double c2 = u * u / w; 112 final double r2 = df * powm1.applyAsDouble(w); 113 // Choose sign at random from the sign of v. 114 return Math.copySign(Math.sqrt(r2 * c2), v); 115 } 116 117 /** {@inheritDoc} */ 118 @Override 119 public StudentsTSampler withUniformRandomProvider(UniformRandomProvider rng) { 120 return new StudentsTSampler(rng, this); 121 } 122 123 /** 124 * Creates a signed double in the range {@code [-1, 1)}. The magnitude is sampled evenly 125 * from the 2<sup>54</sup> dyadic rationals in the range. 126 * 127 * <p>Note: This method will not return samples for both -0.0 and 0.0. 128 * 129 * @param bits the bits 130 * @return the double 131 */ 132 private static double makeSignedDouble(long bits) { 133 // As per o.a.c.rng.core.utils.NumberFactory.makeDouble(long) but using a signed 134 // shift of 10 in place of an unsigned shift of 11. 135 // Use the upper 54 bits on the assumption they are more random. 136 // The sign bit is maintained by the signed shift. 137 // The next 53 bits generates a magnitude in the range [0, 2^53) or [-2^53, 0). 138 return (bits >> 10) * DOUBLE_MULTIPLIER; 139 } 140 } 141 142 /** 143 * Sample from a t-distribution using a normal distribution. 144 * This is used when the degrees of freedom is extremely large (e.g. {@code > 1e16}). 145 */ 146 private static final class NormalTSampler extends TSampler { 147 /** Underlying normalized Gaussian sampler. */ 148 private final NormalizedGaussianSampler sampler; 149 150 /** 151 * @param rng Generator of uniformly distributed random numbers. 152 */ 153 NormalTSampler(UniformRandomProvider rng) { 154 super(rng); 155 this.sampler = ZigguratSampler.NormalizedGaussian.of(rng); 156 } 157 158 /** {@inheritDoc} */ 159 @Override 160 public double sample() { 161 return sampler.sample(); 162 163 } 164 165 /** {@inheritDoc} */ 166 @Override 167 public NormalTSampler withUniformRandomProvider(UniformRandomProvider rng) { 168 return new NormalTSampler(rng); 169 } 170 } 171 172 /** 173 * @param rng Generator of uniformly distributed random numbers. 174 */ 175 TSampler(UniformRandomProvider rng) { 176 this.rng = rng; 177 } 178 179 /** {@inheritDoc} */ 180 // Redeclare the signature to return a TSampler not a SharedStateContinuousSampler 181 @Override 182 public abstract TSampler withUniformRandomProvider(UniformRandomProvider rng); 183 184 /** 185 * Generates a {@code long} value. 186 * Used by algorithm implementations without exposing access to the RNG. 187 * 188 * @return the next random value 189 */ 190 long nextLong() { 191 return rng.nextLong(); 192 } 193 194 /** {@inheritDoc} */ 195 @Override 196 public String toString() { 197 return "Student's t deviate [" + rng.toString() + "]"; 198 } 199 200 /** 201 * Create a new t distribution sampler. 202 * 203 * @param rng Generator of uniformly distributed random numbers. 204 * @param degreesOfFreedom Degrees of freedom. 205 * @return the sampler 206 * @throws IllegalArgumentException if {@code degreesOfFreedom <= 0} 207 */ 208 public static TSampler of(UniformRandomProvider rng, 209 double degreesOfFreedom) { 210 if (degreesOfFreedom > HUGE_DF) { 211 return new NormalTSampler(rng); 212 } else if (degreesOfFreedom > 0) { 213 return new StudentsTSampler(rng, degreesOfFreedom); 214 } else { 215 // df <= 0 or nan 216 throw new IllegalArgumentException( 217 "degrees of freedom is not strictly positive: " + degreesOfFreedom); 218 } 219 } 220}