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&#39;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}