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 * http://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 package org.apache.commons.rng.sampling.distribution;
18
19 import java.util.function.DoubleUnaryOperator;
20 import org.apache.commons.rng.UniformRandomProvider;
21
22 /**
23 * Sampling from a T distribution.
24 *
25 * <p>Uses Bailey's algorithm for t-distribution sampling:</p>
26 *
27 * <blockquote>
28 * <pre>
29 * Bailey, R. W. (1994)
30 * "Polar Generation of Random Variates with the t-Distribution."
31 * Mathematics of Computation 62, 779-781.
32 * </pre>
33 * </blockquote>
34 *
35 * <p>Sampling uses {@link UniformRandomProvider#nextLong()}.</p>
36 *
37 * @see <a href="https://en.wikipedia.org/wiki/Student%27s_t-distribution">Student's T distribution (wikipedia)</a>
38 * @see <a href="https://doi.org/10.2307/2153537">Mathematics of Computation, 62, 779-781</a>
39 * @since 1.5
40 */
41 public abstract class TSampler implements SharedStateContinuousSampler {
42 /** Threshold for huge degrees of freedom. Above this value the CDF of the t distribution
43 * matches the normal distribution. Value is 2/eps (where eps is the machine epsilon)
44 * or approximately 9.0e15. */
45 private static final double HUGE_DF = 0x1.0p53;
46
47 /** Source of randomness. */
48 private final UniformRandomProvider rng;
49
50 /**
51 * Sample from a t-distribution using Bailey's algorithm.
52 */
53 private static final class StudentsTSampler extends TSampler {
54 /** Threshold for large degrees of freedom. */
55 private static final double LARGE_DF = 25;
56 /** The multiplier to convert the least significant 53-bits of a {@code long} to a
57 * uniform {@code double}. */
58 private static final double DOUBLE_MULTIPLIER = 0x1.0p-53;
59
60 /** Degrees of freedom. */
61 private final double df;
62 /** Function to compute pow(x, -2/v) - 1, where v = degrees of freedom. */
63 private final DoubleUnaryOperator powm1;
64
65 /**
66 * @param rng Generator of uniformly distributed random numbers.
67 * @param v Degrees of freedom.
68 */
69 StudentsTSampler(UniformRandomProvider rng,
70 double v) {
71 super(rng);
72 df = v;
73
74 // The sampler requires pow(w, -2/v) - 1 with
75 // 0 <= w <= 1; Expected(w) = sqrt(0.5).
76 // When the exponent is small then pow(x, y) -> 1.
77 // This affects large degrees of freedom.
78 final double exponent = -2 / v;
79 powm1 = v > LARGE_DF ?
80 x -> Math.expm1(Math.log(x) * exponent) :
81 x -> Math.pow(x, exponent) - 1;
82 }
83
84 /**
85 * @param rng Generator of uniformly distributed random numbers.
86 * @param source Source to copy.
87 */
88 private StudentsTSampler(UniformRandomProvider rng,
89 StudentsTSampler source) {
90 super(rng);
91 df = source.df;
92 powm1 = source.powm1;
93 }
94
95 /** {@inheritDoc} */
96 @Override
97 public double sample() {
98 // Require u and v in [0, 1] and a random sign.
99 // 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 }