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 org.apache.commons.rng.UniformRandomProvider;
020
021/**
022 * Samples from a stable distribution.
023 *
024 * <p>Several different parameterizations exist for the stable distribution.
025 * This sampler uses the 0-parameterization distribution described in Nolan (2020) "Univariate Stable
026 * Distributions: Models for Heavy Tailed Data". Springer Series in Operations Research and
027 * Financial Engineering. Springer. Sections 1.7 and 3.3.3.
028 *
029 * <p>The random variable \( X \) has
030 * the stable distribution \( S(\alpha, \beta, \gamma, \delta; 0) \) if its characteristic
031 * function is given by:
032 *
033 * <p>\[ E(e^{iuX}) = \begin{cases} \exp \left (- \gamma^\alpha |u|^\alpha \left [1 - i \beta (\tan \frac{\pi \alpha}{2})(\text{sgn}(u)) \right ] + i \delta u \right ) &amp; \alpha \neq 1 \\
034 * \exp \left (- \gamma |u| \left [1 + i \beta \frac{2}{\pi} (\text{sgn}(u)) \log |u| \right ] + i \delta u \right ) &amp; \alpha = 1 \end{cases} \]
035 *
036 * <p>The function is continuous with respect to all the parameters; the parameters \( \alpha \)
037 * and \( \beta \) determine the shape and the parameters \( \gamma \) and \( \delta \) determine
038 * the scale and location. The support of the distribution is:
039 *
040 * <p>\[ \text{support} f(x|\alpha,\beta,\gamma,\delta; 0) = \begin{cases} [\delta - \gamma \tan \frac{\pi \alpha}{2}, \infty) &amp; \alpha \lt 1\ and\ \beta = 1 \\
041 * (-\infty, \delta + \gamma \tan \frac{\pi \alpha}{2}] &amp; \alpha \lt 1\ and\ \beta = -1 \\
042 * (-\infty, \infty) &amp; otherwise \end{cases} \]
043 *
044 * <p>The implementation uses the Chambers-Mallows-Stuck (CMS) method as described in:
045 * <ul>
046 *  <li>Chambers, Mallows &amp; Stuck (1976) "A Method for Simulating Stable Random Variables".
047 *      Journal of the American Statistical Association. 71 (354): 340–344.
048 *  <li>Weron (1996) "On the Chambers-Mallows-Stuck method for simulating skewed stable
049 *      random variables". Statistics &amp; Probability Letters. 28 (2): 165–171.
050 * </ul>
051 *
052 * @see <a href="https://en.wikipedia.org/wiki/Stable_distribution">Stable distribution (Wikipedia)</a>
053 * @see <a href="https://link.springer.com/book/10.1007/978-3-030-52915-4">Nolan (2020) Univariate Stable Distributions</a>
054 * @see <a href="https://doi.org/10.1080%2F01621459.1976.10480344">Chambers et al (1976) JOASA 71: 340-344</a>
055 * @see <a href="https://doi.org/10.1016%2F0167-7152%2895%2900113-1">Weron (1996).
056 * Statistics &amp; Probability Letters. 28 (2): 165–171.</a>
057 * @since 1.4
058 */
059public abstract class StableSampler implements SharedStateContinuousSampler {
060    /** pi / 2. */
061    private static final double PI_2 = Math.PI / 2;
062    /** The alpha value for the Gaussian case. */
063    private static final double ALPHA_GAUSSIAN = 2;
064    /** The alpha value for the Cauchy case. */
065    private static final double ALPHA_CAUCHY = 1;
066    /** The alpha value for the Levy case. */
067    private static final double ALPHA_LEVY = 0.5;
068    /** The alpha value for the {@code alpha -> 0} to switch to using the Weron formula.
069     * Note that small alpha requires robust correction of infinite samples. */
070    private static final double ALPHA_SMALL = 0.02;
071    /** The beta value for the Levy case. */
072    private static final double BETA_LEVY = 1.0;
073    /** The gamma value for the normalized case. */
074    private static final double GAMMA_1 = 1.0;
075    /** The delta value for the normalized case. */
076    private static final double DELTA_0 = 0.0;
077    /** The tau value for zero. When tau is zero, this is effectively {@code beta = 0}. */
078    private static final double TAU_ZERO = 0.0;
079    /**
080     * The lower support for the distribution.
081     * This is the lower bound of {@code (-inf, +inf)}
082     * If the sample is not within this bound ({@code lower < x}) then it is either
083     * infinite or NaN and the result should be checked.
084     */
085    private static final double LOWER = Double.NEGATIVE_INFINITY;
086    /**
087     * The upper support for the distribution.
088     * This is the upper bound of {@code (-inf, +inf)}.
089     * If the sample is not within this bound ({@code x < upper}) then it is either
090     * infinite or NaN and the result should be checked.
091     */
092    private static final double UPPER = Double.POSITIVE_INFINITY;
093
094    /** Underlying source of randomness. */
095    private final UniformRandomProvider rng;
096
097    // Implementation notes
098    //
099    // The Chambers-Mallows-Stuck (CMS) method uses a uniform deviate u in (0, 1) and an
100    // exponential deviate w to compute a stable deviate. Chambers et al (1976) published
101    // a formula for alpha = 1 and alpha != 1. The function is discontinuous at alpha = 1
102    // and to address this a trigonmoic rearrangement was provided using half angles that
103    // is continuous with respect to alpha. The original discontinuous formulas were proven
104    // in Weron (1996). The CMS rearrangement creates a deviate in the 0-parameterization
105    // defined by Nolan (2020); the original discontinuous functions create a deviate in the
106    // 1-parameterization defined by Nolan. A shift can be used to convert one parameterisation
107    // to the other. The shift is the magnitude of the zeta term from the 1-parameterisation.
108    // The following table shows how the zeta term -> inf when alpha -> 1 for
109    // different beta (hence the discontinuity in the function):
110    //
111    // Zeta
112    //             Beta
113    // Alpha       1.0         0.5         0.25        0.1         0.0
114    // 0.001       0.001571    0.0007854   0.0003927   0.0001571   0.0
115    // 0.01        0.01571     0.007855    0.003927    0.001571    0.0
116    // 0.05        0.07870     0.03935     0.01968     0.007870    0.0
117    // 0.01        0.01571     0.007855    0.003927    0.001571    0.0
118    // 0.1         0.1584      0.07919     0.03960     0.01584     0.0
119    // 0.5         1.000       0.5000      0.2500      0.1000      0.0
120    // 0.9         6.314       3.157       1.578       0.6314      0.0
121    // 0.95        12.71       6.353       3.177       1.271       0.0
122    // 0.99        63.66       31.83       15.91       6.366       0.0
123    // 0.995       127.3       63.66       31.83       12.73       0.0
124    // 0.999       636.6       318.3       159.2       63.66       0.0
125    // 0.9995      1273        636.6       318.3       127.3       0.0
126    // 0.9999      6366        3183        1592        636.6       0.0
127    // 1.0         1.633E+16   8.166E+15   4.083E+15   1.633E+15   0.0
128    //
129    // For numerical simulation the 0-parameterization is favoured as it is continuous
130    // with respect to all the parameters. When approaching alpha = 1 the large magnitude
131    // of the zeta term used to shift the 1-parameterization results in cancellation and the
132    // number of bits of the output sample is effected. This sampler uses the CMS method with
133    // the continuous function as the base for the implementation. However it is not suitable
134    // for all values of alpha and beta.
135    //
136    // The method computes a value log(z) with z in the interval (0, inf). When z is 0 or infinite
137    // the computation can return invalid results. The open bound for the deviate u avoids
138    // generating an extreme value that results in cancellation, z=0 and an invalid expression.
139    // However due to floating point error this can occur
140    // when u is close to 0 or 1, and beta is -1 or 1. Thus it is not enough to create
141    // u by avoiding 0 or 1 and further checks are required.
142    // The division by the deviate w also results in an invalid expression as the term z becomes
143    // infinite as w -> 0. It should be noted that such events are extremely rare
144    // (frequency in the 1 in 10^15), or will not occur at all depending on the parameters alpha
145    // and beta.
146    //
147    // When alpha -> 0 then the distribution is extremely long tailed and the expression
148    // using log(z) often computes infinity. Certain parameters can create NaN due to
149    // 0 / 0, 0 * inf, or inf - inf. Thus the implementation must check the final result
150    // and perform a correction if required, or generate another sample.
151    // Correcting the original CMS formula has many edge cases depending on parameters. The
152    // alternative formula provided by Weron is easier to correct when infinite values are
153    // created. This correction is made easier by knowing that u is not 0 or 1 as certain
154    // conditions on the intermediate terms can be eliminated. The implementation
155    // thus generates u in the open interval (0,1) but leaves w unchecked and potentially 0.
156    // The sample is generated and the result tested against the expected support. This detects
157    // any NaN and infinite values. Incorrect samples due to the inability to compute log(z)
158    // (extremely rare) and samples where alpha -> 0 has resulted in an infinite expression
159    // for the value d are corrected using the Weron formula and returned within the support.
160    //
161    // The CMS algorithm is continuous for the parameters. However when alpha=1 or beta=0
162    // many terms cancel and these cases are handled with specialised implementations.
163    // The beta=0 case implements the same CMS algorithm with certain terms eliminated.
164    // Correction uses the alternative Weron formula. When alpha=1 the CMS algorithm can
165    // be corrected from infinite cases due to assumptions on the intermediate terms.
166    //
167    // The following table show the failure frequency (result not finite or, when beta=+/-1,
168    // within the support) for the CMS algorithm computed using 2^30 random deviates.
169    //
170    // CMS failure rate
171    //             Beta
172    // Alpha       1.0         0.5         0.25        0.1         0.0
173    // 1.999       0.0         0.0         0.0         0.0         0.0
174    // 1.99        0.0         0.0         0.0         0.0         0.0
175    // 1.9         0.0         0.0         0.0         0.0         0.0
176    // 1.5         0.0         0.0         0.0         0.0         0.0
177    // 1.1         0.0         0.0         0.0         0.0         0.0
178    // 1.0         0.0         0.0         0.0         0.0         0.0
179    // 0.9         0.0         0.0         0.0         0.0         0.0
180    // 0.5         0.0         0.0         0.0         0.0         0.0
181    // 0.25        0.0         0.0         0.0         0.0         0.0
182    // 0.1         0.0         0.0         0.0         0.0         0.0
183    // 0.05        0.0003458   0.0         0.0         0.0         0.0
184    // 0.02        0.009028    6.938E-7    7.180E-7    7.320E-7    6.873E-7
185    // 0.01        0.004878    0.0008555   0.0008553   0.0008554   0.0008570
186    // 0.005       0.1519      0.02896     0.02897     0.02897     0.02897
187    // 0.001       0.6038      0.3903      0.3903      0.3903      0.3903
188    //
189    // The sampler switches to using the error checked Weron implementation when alpha < 0.02.
190    // Unit tests demonstrate the two samplers (CMS or Weron) product the same result within
191    // a tolerance. The switch point is based on a consistent failure rate above 1 in a million.
192    // At this point zeta is small and cancellation leading to loss of bits in the sample is
193    // minimal.
194    //
195    // In common use the sampler will not have a measurable failure rate. The output will
196    // be continuous as alpha -> 1 and beta -> 0. The evaluated function produces symmetric
197    // samples when u and beta are mirrored around 0.5 and 0 respectively. To achieve this
198    // the computation of certain parameters has been changed from the original implementation
199    // to avoid evaluating Math.tan outside the interval (-pi/2, pi/2).
200    //
201    // Note: Chambers et al (1976) use an approximation to tan(x) / x in the RSTAB routine.
202    // A JMH performance test is available in the RNG examples module comparing Math.tan
203    // with various approximations. The functions are faster than Math.tan(x) / x.
204    // This implementation uses a higher accuracy approximation than the original RSTAB
205    // implementation; it has a mean ULP difference to Math.tan of less than 1 and has
206    // a noticeable performance gain.
207
208    /**
209     * Base class for implementations of a stable distribution that requires an exponential
210     * random deviate.
211     */
212    private abstract static class BaseStableSampler extends StableSampler {
213        /** pi/2 scaled by 2^-53. */
214        private static final double PI_2_SCALED = 0x1.0p-54 * Math.PI;
215        /** pi/4 scaled by 2^-53. */
216        private static final double PI_4_SCALED = 0x1.0p-55 * Math.PI;
217        /** -pi / 2. */
218        private static final double NEG_PI_2 = -Math.PI / 2;
219        /** -pi / 4. */
220        private static final double NEG_PI_4 = -Math.PI / 4;
221
222        /** The exponential sampler. */
223        private final ContinuousSampler expSampler;
224
225        /**
226         * @param rng Underlying source of randomness
227         */
228        BaseStableSampler(UniformRandomProvider rng) {
229            super(rng);
230            expSampler = ZigguratSampler.Exponential.of(rng);
231        }
232
233        /**
234         * Gets a random value for the omega parameter ({@code w}).
235         * This is an exponential random variable with mean 1.
236         *
237         * <p>Warning: For simplicity this does not check the variate is not 0.
238         * The calling CMS algorithm should detect and handle incorrect samples as a result
239         * of this unlikely edge case.
240         *
241         * @return omega
242         */
243        double getOmega() {
244            // Note: Ideally this should not have a value of 0 as the CMS algorithm divides
245            // by w and it creates infinity. This can result in NaN output.
246            // Under certain parameterizations non-zero small w also creates NaN output.
247            // Thus output should be checked regardless.
248            return expSampler.sample();
249        }
250
251        /**
252         * Gets a random value for the phi parameter.
253         * This is a uniform random variable in {@code (-pi/2, pi/2)}.
254         *
255         * @return phi
256         */
257        double getPhi() {
258            // See getPhiBy2 for method details.
259            final double x = (nextLong() >> 10) * PI_2_SCALED;
260            // Deliberate floating-point equality check
261            if (x == NEG_PI_2) {
262                return getPhi();
263            }
264            return x;
265        }
266
267        /**
268         * Gets a random value for the phi parameter divided by 2.
269         * This is a uniform random variable in {@code (-pi/4, pi/4)}.
270         *
271         * <p>Note: Ideally this should not have a value of -pi/4 or pi/4 as the CMS algorithm
272         * can generate infinite values when the phi/2 uniform deviate is +/-pi/4. This
273         * can result in NaN output. Under certain parameterizations phi/2 close to the limits
274         * also create NaN output. Thus output should be checked regardless. Avoiding
275         * the extreme values simplifies the number of checks that are required.
276         *
277         * @return phi / 2
278         */
279        double getPhiBy2() {
280            // As per o.a.c.rng.core.utils.NumberFactory.makeDouble(long) but using a
281            // signed shift of 10 in place of an unsigned shift of 11. With a factor of 2^-53
282            // this would produce a double in [-1, 1).
283            // Here the multiplication factor incorporates pi/4 to avoid a separate
284            // multiplication.
285            final double x = (nextLong() >> 10) * PI_4_SCALED;
286            // Deliberate floating-point equality check
287            if (x == NEG_PI_4) {
288                // Sample again using recursion.
289                // A stack overflow due to a broken RNG will eventually occur
290                // rather than the alternative which is an infinite loop
291                // while x == -pi/4.
292                return getPhiBy2();
293            }
294            return x;
295        }
296    }
297
298    /**
299     * Class for implementations of a stable distribution transformed by scale and location.
300     */
301    private static final class TransformedStableSampler extends StableSampler {
302        /** Underlying normalized stable sampler. */
303        private final StableSampler sampler;
304        /** The scale parameter. */
305        private final double gamma;
306        /** The location parameter. */
307        private final double delta;
308
309        /**
310         * @param sampler Normalized stable sampler.
311         * @param gamma Scale parameter. Must be strictly positive.
312         * @param delta Location parameter.
313         */
314        TransformedStableSampler(StableSampler sampler, double gamma, double delta) {
315            // No RNG required
316            super(null);
317            this.sampler = sampler;
318            this.gamma = gamma;
319            this.delta = delta;
320        }
321
322        @Override
323        public double sample() {
324            return gamma * sampler.sample() + delta;
325        }
326
327        @Override
328        public StableSampler withUniformRandomProvider(UniformRandomProvider rng) {
329            return new TransformedStableSampler(sampler.withUniformRandomProvider(rng),
330                                                gamma, delta);
331        }
332
333        @Override
334        public String toString() {
335            // Avoid a null pointer from the unset RNG instance in the parent class
336            return sampler.toString();
337        }
338    }
339
340    /**
341     * Implement the {@code alpha = 2} stable distribution case (Gaussian distribution).
342     */
343    private static final class GaussianStableSampler extends StableSampler {
344        /** sqrt(2). */
345        private static final double ROOT_2 = Math.sqrt(2);
346
347        /** Underlying normalized Gaussian sampler. */
348        private final NormalizedGaussianSampler sampler;
349        /** The standard deviation. */
350        private final double stdDev;
351        /** The mean. */
352        private final double mean;
353
354        /**
355         * @param rng Underlying source of randomness
356         * @param gamma Scale parameter. Must be strictly positive.
357         * @param delta Location parameter.
358         */
359        GaussianStableSampler(UniformRandomProvider rng, double gamma, double delta) {
360            super(rng);
361            this.sampler = ZigguratSampler.NormalizedGaussian.of(rng);
362            // A standardized stable sampler with alpha=2 has variance 2.
363            // Set the standard deviation as sqrt(2) * scale.
364            // Avoid this being infinity to avoid inf * 0 in the sample
365            this.stdDev = Math.min(Double.MAX_VALUE, ROOT_2 * gamma);
366            this.mean = delta;
367        }
368
369        /**
370         * @param rng Underlying source of randomness
371         * @param source Source to copy.
372         */
373        GaussianStableSampler(UniformRandomProvider rng, GaussianStableSampler source) {
374            super(rng);
375            this.sampler = ZigguratSampler.NormalizedGaussian.of(rng);
376            this.stdDev = source.stdDev;
377            this.mean = source.mean;
378        }
379
380        @Override
381        public double sample() {
382            return stdDev * sampler.sample() + mean;
383        }
384
385        @Override
386        public GaussianStableSampler withUniformRandomProvider(UniformRandomProvider rng) {
387            return new GaussianStableSampler(rng, this);
388        }
389    }
390
391    /**
392     * Implement the {@code alpha = 1} and {@code beta = 0} stable distribution case
393     * (Cauchy distribution).
394     */
395    private static final class CauchyStableSampler extends BaseStableSampler {
396        /** The scale parameter. */
397        private final double gamma;
398        /** The location parameter. */
399        private final double delta;
400
401        /**
402         * @param rng Underlying source of randomness
403         * @param gamma Scale parameter. Must be strictly positive.
404         * @param delta Location parameter.
405         */
406        CauchyStableSampler(UniformRandomProvider rng, double gamma, double delta) {
407            super(rng);
408            this.gamma = gamma;
409            this.delta = delta;
410        }
411
412        /**
413         * @param rng Underlying source of randomness
414         * @param source Source to copy.
415         */
416        CauchyStableSampler(UniformRandomProvider rng, CauchyStableSampler source) {
417            super(rng);
418            this.gamma = source.gamma;
419            this.delta = source.delta;
420        }
421
422        @Override
423        public double sample() {
424            // Note:
425            // The CMS beta=0 with alpha=1 sampler reduces to:
426            // S = 2 * a / a2, with a = tan(x), a2 = 1 - a^2, x = phi/2
427            // This is a double angle identity for tan:
428            // 2 * tan(x) / (1 - tan^2(x)) = tan(2x)
429            // Here we use the double angle identity for consistency with the other samplers.
430            final double phiby2 = getPhiBy2();
431            final double a = phiby2 * SpecialMath.tan2(phiby2);
432            final double a2 = 1 - a * a;
433            final double x = 2 * a / a2;
434            return gamma * x + delta;
435        }
436
437        @Override
438        public CauchyStableSampler withUniformRandomProvider(UniformRandomProvider rng) {
439            return new CauchyStableSampler(rng, this);
440        }
441    }
442
443    /**
444     * Implement the {@code alpha = 0.5} and {@code beta = 1} stable distribution case
445     * (Levy distribution).
446     *
447     * Note: This sampler can be used to output the symmetric case when
448     * {@code beta = -1} by negating {@code gamma}.
449     */
450    private static final class LevyStableSampler extends StableSampler {
451        /** Underlying normalized Gaussian sampler. */
452        private final NormalizedGaussianSampler sampler;
453        /** The scale parameter. */
454        private final double gamma;
455        /** The location parameter. */
456        private final double delta;
457
458        /**
459         * @param rng Underlying source of randomness
460         * @param gamma Scale parameter. Must be strictly positive.
461         * @param delta Location parameter.
462         */
463        LevyStableSampler(UniformRandomProvider rng, double gamma, double delta) {
464            super(rng);
465            this.sampler = ZigguratSampler.NormalizedGaussian.of(rng);
466            this.gamma = gamma;
467            this.delta = delta;
468        }
469
470        /**
471         * @param rng Underlying source of randomness
472         * @param source Source to copy.
473         */
474        LevyStableSampler(UniformRandomProvider rng, LevyStableSampler source) {
475            super(rng);
476            this.sampler = ZigguratSampler.NormalizedGaussian.of(rng);
477            this.gamma = source.gamma;
478            this.delta = source.delta;
479        }
480
481        @Override
482        public double sample() {
483            // Levy(Z) = 1 / N(0,1)^2, where N(0,1) is a standard normalized variate
484            final double norm = sampler.sample();
485            // Here we must transform from the 1-parameterization to the 0-parameterization.
486            // This is a shift of -beta * tan(pi * alpha / 2) = -1 when alpha=0.5, beta=1.
487            final double z = (1.0 / (norm * norm)) - 1.0;
488            // In the 0-parameterization the scale and location are a linear transform.
489            return gamma * z + delta;
490        }
491
492        @Override
493        public LevyStableSampler withUniformRandomProvider(UniformRandomProvider rng) {
494            return new LevyStableSampler(rng, this);
495        }
496    }
497
498    /**
499     * Implement the generic stable distribution case: {@code alpha < 2} and
500     * {@code beta != 0}. This routine assumes {@code alpha != 1}.
501     *
502     * <p>Implements the Chambers-Mallows-Stuck (CMS) method using the
503     * formula provided in Weron (1996) "On the Chambers-Mallows-Stuck method for
504     * simulating skewed stable random variables" Statistics &amp; Probability
505     * Letters. 28 (2): 165–171. This method is easier to correct from infinite and
506     * NaN results by boxing intermediate infinite values.
507     *
508     * <p>The formula produces a stable deviate from the 1-parameterization that is
509     * discontinuous at {@code alpha=1}. A shift is used to create the 0-parameterization.
510     * This shift is very large as {@code alpha -> 1} and the output loses bits of precision
511     * in the deviate due to cancellation. It is not recommended to use this sampler when
512     * {@code alpha -> 1} except for edge case correction.
513     *
514     * <p>This produces non-NaN output for all parameters alpha, beta, u and w with
515     * the correct orientation for extremes of the distribution support.
516     * The formulas used are symmetric with regard to beta and u.
517     *
518     * @see <a href="https://doi.org/10.1016%2F0167-7152%2895%2900113-1">Weron, R
519     * (1996). Statistics &amp; Probability Letters. 28 (2): 165–171.</a>
520     */
521    static class WeronStableSampler extends BaseStableSampler {
522        /** Epsilon (1 - alpha). */
523        protected final double eps;
524        /** Alpha (1 - eps). */
525        protected final double meps1;
526        /** Cache of expression value used in generation. */
527        protected final double zeta;
528        /** Cache of expression value used in generation. */
529        protected final double atanZeta;
530        /** Cache of expression value used in generation. */
531        protected final double scale;
532        /** 1 / alpha = 1 / (1 - eps). */
533        protected final double inv1mEps;
534        /** (1 / alpha) - 1 = eps / (1 - eps). */
535        protected final double epsDiv1mEps;
536        /** The inclusive lower support for the distribution. */
537        protected final double lower;
538        /** The inclusive upper support for the distribution. */
539        protected final double upper;
540
541        /**
542         * @param rng Underlying source of randomness
543         * @param alpha Stability parameter. Must be in the interval {@code (0, 2]}.
544         * @param beta Skewness parameter. Must be in the interval {@code [-1, 1]}.
545         */
546        WeronStableSampler(UniformRandomProvider rng, double alpha, double beta) {
547            super(rng);
548            eps = 1 - alpha;
549            // When alpha < 0.5, 1 - eps == alpha is not always true as the reverse is not exact.
550            // Here we store 1 - eps in place of alpha. Thus eps + (1 - eps) = 1.
551            meps1 = 1 - eps;
552
553            // Compute pre-factors for the Weron formula used during error correction.
554            if (meps1 > 1) {
555                // Avoid calling tan outside the domain limit [-pi/2, pi/2].
556                zeta = beta * Math.tan((2 - meps1) * PI_2);
557            } else {
558                zeta = -beta * Math.tan(meps1 * PI_2);
559            }
560
561            // Do not store xi = Math.atan(-zeta) / meps1 due to floating-point division errors.
562            // Directly store Math.atan(-zeta).
563            atanZeta = Math.atan(-zeta);
564            scale = Math.pow(1 + zeta * zeta, 0.5 / meps1);
565            // Note: These terms are used interchangeably in formulas
566            //    1         1
567            // -------  = -----
568            // (1-eps)    alpha
569            inv1mEps = 1.0 / meps1;
570            //    1             eps     (1-alpha)     1
571            // -------  - 1 = ------- = --------- = ----- - 1
572            // (1-eps)        (1-eps)     alpha     alpha
573            epsDiv1mEps = inv1mEps - 1;
574
575
576            // Compute the support. This applies when alpha < 1 and beta = +/-1
577            if (alpha < 1 && Math.abs(beta) == 1) {
578                if (beta == 1) {
579                    // alpha < 0, beta = 1
580                    lower = zeta;
581                    upper = UPPER;
582                } else {
583                    // alpha < 0, beta = -1
584                    lower = LOWER;
585                    upper = zeta;
586                }
587            } else {
588                lower = LOWER;
589                upper = UPPER;
590            }
591        }
592
593        /**
594         * @param rng Underlying source of randomness
595         * @param source Source to copy.
596         */
597        WeronStableSampler(UniformRandomProvider rng, WeronStableSampler source) {
598            super(rng);
599            this.eps = source.eps;
600            this.meps1 = source.meps1;
601            this.zeta = source.zeta;
602            this.atanZeta = source.atanZeta;
603            this.scale = source.scale;
604            this.inv1mEps = source.inv1mEps;
605            this.epsDiv1mEps = source.epsDiv1mEps;
606            this.lower = source.lower;
607            this.upper = source.upper;
608        }
609
610        @Override
611        public double sample() {
612            final double phi = getPhi();
613            final double w = getOmega();
614            return createSample(phi, w);
615        }
616
617        /**
618         * Create the sample. This routine is robust to edge cases and returns a deviate
619         * at the extremes of the support. It correctly handles {@code alpha -> 0} when
620         * the sample is increasingly likely to be +/- infinity.
621         *
622         * @param phi Uniform deviate in {@code (-pi/2, pi/2)}
623         * @param w Exponential deviate
624         * @return x
625         */
626        protected double createSample(double phi, double w) {
627            // Here we use the formula provided by Weron for the 1-parameterization.
628            // Note: Adding back zeta creates the 0-parameterization defined in Nolan (1998):
629            // X ~ S0_alpha(s,beta,u0) with s=1, u0=0 for a standard random variable.
630            // As alpha -> 1 the translation zeta to create the stable deviate
631            // in the 0-parameterization is increasingly large as tan(pi/2) -> infinity.
632            // The max translation is approximately 1e16.
633            // Without this translation the stable deviate is in the 1-parameterization
634            // and the function is not continuous with respect to alpha.
635            // Due to the large zeta when alpha -> 1 the number of bits of the output variable
636            // are very low due to cancellation.
637
638            // As alpha -> 0 or 2 then zeta -> 0 and cancellation is not relevant.
639            // The formula can be modified for infinite terms to compute a result for extreme
640            // deviates u and w when the CMS formula fails.
641
642            // Note the following term is subject to floating point error:
643            // xi = atan(-zeta) / alpha
644            // alphaPhiXi = alpha * (phi + xi)
645            // This is required: cos(phi - alphaPhiXi) > 0 => phi - alphaPhiXi in (-pi/2, pi/2).
646            // Thus we compute atan(-zeta) and use it to compute two terms:
647            // [1] alpha * (phi + xi) = alpha * (phi + atan(-zeta) / alpha) = alpha * phi + atan(-zeta)
648            // [2] phi - alpha * (phi + xi) = phi - alpha * phi - atan(-zeta) = (1-alpha) * phi - atan(-zeta)
649
650            // Compute terms
651            // Either term can be infinite or 0. Certain parameters compute 0 * inf.
652            // t1=inf occurs alpha -> 0.
653            // t1=0 occurs when beta = tan(-alpha * phi) / tan(alpha * pi / 2).
654            // t2=inf occurs when w -> 0 and alpha -> 0.
655            // t2=0 occurs when alpha -> 0 and phi -> pi/2.
656            // Detect zeros and return as zeta.
657
658            // Note sin(alpha * phi + atanZeta) is zero when:
659            // alpha * phi = -atan(-zeta)
660            // tan(-alpha * phi) = -zeta
661            //                   = beta * tan(alpha * pi / 2)
662            // Since |phi| < pi/2 this requires beta to have an opposite sign to phi
663            // and a magnitude < 1. This is possible and in this case avoid a possible
664            // 0 / 0 by setting the result as if term t1=0 and the result is zeta.
665            double t1 = Math.sin(meps1 * phi + atanZeta);
666            if (t1 == 0) {
667                return zeta;
668            }
669            // Since cos(phi) is in (0, 1] this term will not create a
670            // large magnitude to create t1 = 0.
671            t1 /= Math.pow(Math.cos(phi), inv1mEps);
672
673            // Iff Math.cos(eps * phi - atanZeta) is zero then 0 / 0 can occur if w=0.
674            // Iff Math.cos(eps * phi - atanZeta) is below zero then NaN will occur
675            // in the power function. These cases are avoided by phi=(-pi/2, pi/2) and direct
676            // use of arctan(-zeta).
677            final double t2 = Math.pow(Math.cos(eps * phi - atanZeta) / w, epsDiv1mEps);
678            if (t2 == 0) {
679                return zeta;
680            }
681
682            final double x = t1 * t2 * scale + zeta;
683
684            // Check the bounds. Applies when alpha < 1 and beta = +/-1.
685            if (x <= lower) {
686                return lower;
687            }
688            return x < upper ? x : upper;
689        }
690
691        @Override
692        public WeronStableSampler withUniformRandomProvider(UniformRandomProvider rng) {
693            return new WeronStableSampler(rng, this);
694        }
695    }
696
697    /**
698     * Implement the generic stable distribution case: {@code alpha < 2} and
699     * {@code beta != 0}. This routine assumes {@code alpha != 1}.
700     *
701     * <p>Implements the Chambers-Mallows-Stuck (CMS) method from Chambers, et al
702     * (1976) A Method for Simulating Stable Random Variables. Journal of the
703     * American Statistical Association Vol. 71, No. 354, pp. 340-344.
704     *
705     * <p>The formula produces a stable deviate from the 0-parameterization that is
706     * continuous at {@code alpha=1}.
707     *
708     * <p>This is an implementation of the Fortran routine RSTAB. In the event the
709     * computation fails then an alternative computation is performed using the
710     * formula provided in Weron (1996) "On the Chambers-Mallows-Stuck method for
711     * simulating skewed stable random variables" Statistics &amp; Probability
712     * Letters. 28 (2): 165–171. This method is easier to correct from infinite and
713     * NaN results. The error correction path is extremely unlikely to occur during
714     * use unless {@code alpha -> 0}. In general use it requires the random deviates
715     * w or u are extreme. See the unit tests for conditions that create them.
716     *
717     * <p>This produces non-NaN output for all parameters alpha, beta, u and w with
718     * the correct orientation for extremes of the distribution support.
719     * The formulas used are symmetric with regard to beta and u.
720     */
721    static class CMSStableSampler extends WeronStableSampler {
722        /** 1/2. */
723        private static final double HALF = 0.5;
724        /** Cache of expression value used in generation. */
725        private final double tau;
726
727        /**
728         * @param rng Underlying source of randomness
729         * @param alpha Stability parameter. Must be in the interval {@code (0, 2]}.
730         * @param beta Skewness parameter. Must be in the interval {@code [-1, 1]}.
731         */
732        CMSStableSampler(UniformRandomProvider rng, double alpha, double beta) {
733            super(rng, alpha, beta);
734
735            // Compute the RSTAB pre-factor.
736            tau = getTau(alpha, beta);
737        }
738
739        /**
740         * @param rng Underlying source of randomness
741         * @param source Source to copy.
742         */
743        CMSStableSampler(UniformRandomProvider rng, CMSStableSampler source) {
744            super(rng, source);
745            this.tau = source.tau;
746        }
747
748        /**
749         * Gets tau. This is a factor used in the CMS algorithm. If this is zero then
750         * a special case of {@code beta -> 0} has occurred.
751         *
752         * @param alpha Stability parameter. Must be in the interval {@code (0, 2]}.
753         * @param beta Skewness parameter. Must be in the interval {@code [-1, 1]}.
754         * @return tau
755         */
756        static double getTau(double alpha, double beta) {
757            final double eps = 1 - alpha;
758            final double meps1 = 1 - eps;
759            // Compute RSTAB prefactor
760            final double tau;
761
762            // tau is symmetric around alpha=1
763            // tau -> beta / pi/2 as alpha -> 1
764            // tau -> 0 as alpha -> 2 or 0
765            // Avoid calling tan as the value approaches the domain limit [-pi/2, pi/2].
766            if (Math.abs(eps) < HALF) {
767                // 0.5 < alpha < 1.5. Note: This works when eps=0 as tan(0) / 0 == 1.
768                tau = beta / (SpecialMath.tan2(eps * PI_2) * PI_2);
769            } else {
770                // alpha >= 1.5 or alpha <= 0.5.
771                // Do not call tan with alpha > 1 as it wraps in the domain [-pi/2, pi/2].
772                // Since pi is approximate the symmetry is lost by wrapping.
773                // Keep within the domain using (2-alpha).
774                if (meps1 > 1) {
775                    tau = beta * PI_2 * eps * (2 - meps1) * -SpecialMath.tan2((2 - meps1) * PI_2);
776                } else {
777                    tau = beta * PI_2 * eps * meps1 * SpecialMath.tan2(meps1 * PI_2);
778                }
779            }
780
781            return tau;
782        }
783
784        @Override
785        public double sample() {
786            final double phiby2 = getPhiBy2();
787            final double w = getOmega();
788
789            // Compute as per the RSTAB routine.
790
791            // Generic stable distribution that is continuous as alpha -> 1.
792            // This is a trigonomic rearrangement of equation 4.1 from Chambers et al (1976)
793            // as implemented in the Fortran program RSTAB.
794            // Uses the special functions:
795            // tan2 = tan(x) / x
796            // d2 = (exp(x) - 1) / x
797            // The method is implemented as per the RSTAB routine with the exceptions:
798            // 1. The function tan2(x) is implemented with a higher precision approximation.
799            // 2. The sample is tested against the expected distribution support.
800            // Infinite intermediate terms that create infinite or NaN are corrected by
801            // switching the formula and handling infinite terms.
802
803            // Compute some tangents
804            // a in (-1, 1)
805            // bb in [1, 4/pi)
806            // b in (-1, 1)
807            final double a = phiby2 * SpecialMath.tan2(phiby2);
808            final double bb = SpecialMath.tan2(eps * phiby2);
809            final double b = eps * phiby2 * bb;
810
811            // Compute some necessary subexpressions
812            final double da = a * a;
813            final double db = b * b;
814            // a2 in (0, 1]
815            final double a2 = 1 - da;
816            // a2p in [1, 2)
817            final double a2p = 1 + da;
818            // b2 in (0, 1]
819            final double b2 = 1 - db;
820            // b2p in [1, 2)
821            final double b2p = 1 + db;
822
823            // Compute coefficient.
824            // numerator=0 is not possible *in theory* when the uniform deviate generating phi
825            // is in the open interval (0, 1). In practice it is possible to obtain <=0 due
826            // to round-off error, typically when beta -> +/-1 and phiby2 -> -/+pi/4.
827            // This can happen for any alpha.
828            final double z = a2p * (b2 + 2 * phiby2 * bb * tau) / (w * a2 * b2p);
829
830            // Compute the exponential-type expression
831            // Note: z may be infinite, typically when w->0 and a2->0.
832            // This can produce NaN under certain parameterizations due to multiplication by 0.
833            final double alogz = Math.log(z);
834            final double d = SpecialMath.d2(epsDiv1mEps * alogz) * (alogz * inv1mEps);
835
836            // Pre-compute the multiplication factor.
837            // The numerator may be zero. The denominator is not zero as a2 is bounded to
838            // above zero when the uniform deviate that generates phiby2 is not 0 or 1.
839            // The min value of a2 is 2^-52. Assume f cannot be infinite as the numerator
840            // is computed with a in (-1, 1); b in (-1, 1); phiby2 in (-pi/4, pi/4); tau in
841            // [-2/pi, 2/pi]; bb in [1, 4/pi); a2 in (0, 1] limiting the numerator magnitude.
842            final double f = (2 * ((a - b) * (1 + a * b) - phiby2 * tau * bb * (b * a2 - 2 * a))) /
843                    (a2 * b2p);
844
845            // Compute the stable deviate:
846            final double x = (1 + eps * d) * f + tau * d;
847
848            // Test the support
849            if (lower < x && x < upper) {
850                return x;
851            }
852
853            // Error correction path:
854            // x is at the bounds, infinite or NaN (created by 0 / 0,  0 * inf, or inf - inf).
855            // This is caused by extreme parameterizations of alpha or beta, or extreme values
856            // from the random deviates.
857            // Alternatively alpha < 1 and beta = +/-1 and the sample x is at the edge or
858            // outside the support due to floating point error.
859
860            // Here we use the formula provided by Weron which is easier to correct
861            // when deviates are extreme or alpha -> 0. The formula is not continuous
862            // as alpha -> 1 without a shift which reduces the precision of the sample;
863            // for rare edge case correction this has minimal effect on sampler output.
864            return createSample(phiby2 * 2, w);
865        }
866
867        @Override
868        public CMSStableSampler withUniformRandomProvider(UniformRandomProvider rng) {
869            return new CMSStableSampler(rng, this);
870        }
871    }
872
873    /**
874     * Implement the stable distribution case: {@code alpha == 1} and {@code beta != 0}.
875     *
876     * <p>Implements the same algorithm as the {@link CMSStableSampler} with
877     * the {@code alpha} assumed to be 1.
878     *
879     * <p>This sampler specifically requires that {@code beta / (pi/2) != 0}; otherwise
880     * the parameters equal {@code alpha=1, beta=0} as the Cauchy distribution case.
881     */
882    static class Alpha1CMSStableSampler extends BaseStableSampler {
883        /** Cache of expression value used in generation. */
884        private final double tau;
885
886        /**
887         * @param rng Underlying source of randomness
888         * @param beta Skewness parameter. Must be in the interval {@code [-1, 1]}.
889         */
890        Alpha1CMSStableSampler(UniformRandomProvider rng, double beta) {
891            super(rng);
892            tau = beta / PI_2;
893        }
894
895        /**
896         * @param rng Underlying source of randomness
897         * @param source Source to copy.
898         */
899        Alpha1CMSStableSampler(UniformRandomProvider rng, Alpha1CMSStableSampler source) {
900            super(rng);
901            this.tau = source.tau;
902        }
903
904        @Override
905        public double sample() {
906            final double phiby2 = getPhiBy2();
907            final double w = getOmega();
908
909            // Compute some tangents
910            final double a = phiby2 * SpecialMath.tan2(phiby2);
911
912            // Compute some necessary subexpressions
913            final double da = a * a;
914            final double a2 = 1 - da;
915            final double a2p = 1 + da;
916
917            // Compute coefficient.
918
919            // numerator=0 is not possible when the uniform deviate generating phi
920            // is in the open interval (0, 1) and alpha=1.
921            final double z = a2p * (1 + 2 * phiby2 * tau) / (w * a2);
922
923            // Compute the exponential-type expression
924            // Note: z may be infinite, typically when w->0 and a2->0.
925            // This can produce NaN under certain parameterizations due to multiplication by 0.
926            // When alpha=1 the expression
927            // d = d2((eps / (1-eps)) * alogz) * (alogz / (1-eps)) is eliminated to 1 * log(z)
928            final double d = Math.log(z);
929
930            // Pre-compute the multiplication factor.
931            final double f = (2 * (a - phiby2 * tau * (-2 * a))) / a2;
932
933            // Compute the stable deviate:
934            // This does not require correction as f is finite (as per the alpha != 1 case),
935            // tau is non-zero and only d can be infinite due to an extreme w -> 0.
936            return f + tau * d;
937        }
938
939        @Override
940        public Alpha1CMSStableSampler withUniformRandomProvider(UniformRandomProvider rng) {
941            return new Alpha1CMSStableSampler(rng, this);
942        }
943    }
944
945    /**
946     * Implement the generic stable distribution case: {@code alpha < 2} and {@code beta == 0}.
947     *
948     * <p>Implements the same algorithm as the {@link WeronStableSampler} with
949     * the {@code beta} assumed to be 0.
950     *
951     * <p>This routine assumes {@code alpha != 1}; {@code alpha=1, beta=0} is the Cauchy
952     * distribution case.
953     */
954    static class Beta0WeronStableSampler extends BaseStableSampler {
955        /** Epsilon (1 - alpha). */
956        protected final double eps;
957        /** Epsilon (1 - alpha). */
958        protected final double meps1;
959        /** 1 / alpha = 1 / (1 - eps). */
960        protected final double inv1mEps;
961        /** (1 / alpha) - 1 = eps / (1 - eps). */
962        protected final double epsDiv1mEps;
963
964        /**
965         * @param rng Underlying source of randomness
966         * @param alpha Stability parameter. Must be in the interval {@code (0, 2]}.
967         */
968        Beta0WeronStableSampler(UniformRandomProvider rng, double alpha) {
969            super(rng);
970            eps = 1 - alpha;
971            meps1 = 1 - eps;
972            inv1mEps = 1.0 / meps1;
973            epsDiv1mEps = inv1mEps - 1;
974        }
975
976        /**
977         * @param rng Underlying source of randomness
978         * @param source Source to copy.
979         */
980        Beta0WeronStableSampler(UniformRandomProvider rng, Beta0WeronStableSampler source) {
981            super(rng);
982            this.eps = source.eps;
983            this.meps1 = source.meps1;
984            this.inv1mEps = source.inv1mEps;
985            this.epsDiv1mEps = source.epsDiv1mEps;
986        }
987
988        @Override
989        public double sample() {
990            final double phi = getPhi();
991            final double w = getOmega();
992            return createSample(phi, w);
993        }
994
995        /**
996         * Create the sample. This routine is robust to edge cases and returns a deviate
997         * at the extremes of the support. It correctly handles {@code alpha -> 0} when
998         * the sample is increasingly likely to be +/- infinity.
999         *
1000         * @param phi Uniform deviate in {@code (-pi/2, pi/2)}
1001         * @param w Exponential deviate
1002         * @return x
1003         */
1004        protected double createSample(double phi, double w) {
1005            // As per the Weron sampler with beta=0 and terms eliminated.
1006            // Note that if alpha=1 this reduces to sin(phi) / cos(phi) => Cauchy case.
1007
1008            // Compute terms.
1009            // Either term can be infinite or 0. Certain parameters compute 0 * inf.
1010            // Detect zeros and return as 0.
1011
1012            // Note sin(alpha * phi) is only ever zero when phi=0. No value of alpha
1013            // multiplied by small phi can create zero due to the limited
1014            // precision of alpha imposed by alpha = 1 - (1-alpha). At this point cos(phi) = 1.
1015            // Thus 0/0 cannot occur.
1016            final double t1 = Math.sin(meps1 * phi) / Math.pow(Math.cos(phi), inv1mEps);
1017            if (t1 == 0) {
1018                return 0;
1019            }
1020            final double t2 = Math.pow(Math.cos(eps * phi) / w, epsDiv1mEps);
1021            if (t2 == 0) {
1022                return 0;
1023            }
1024            return t1 * t2;
1025        }
1026
1027        @Override
1028        public Beta0WeronStableSampler withUniformRandomProvider(UniformRandomProvider rng) {
1029            return new Beta0WeronStableSampler(rng, this);
1030        }
1031    }
1032
1033    /**
1034     * Implement the generic stable distribution case: {@code alpha < 2} and {@code beta == 0}.
1035     *
1036     * <p>Implements the same algorithm as the {@link CMSStableSampler} with
1037     * the {@code beta} assumed to be 0.
1038     *
1039     * <p>This routine assumes {@code alpha != 1}; {@code alpha=1, beta=0} is the Cauchy
1040     * distribution case.
1041     */
1042    static class Beta0CMSStableSampler extends Beta0WeronStableSampler {
1043        /**
1044         * @param rng Underlying source of randomness
1045         * @param alpha Stability parameter. Must be in the interval {@code (0, 2]}.
1046         */
1047        Beta0CMSStableSampler(UniformRandomProvider rng, double alpha) {
1048            super(rng, alpha);
1049        }
1050
1051        /**
1052         * @param rng Underlying source of randomness
1053         * @param source Source to copy.
1054         */
1055        Beta0CMSStableSampler(UniformRandomProvider rng, Beta0CMSStableSampler source) {
1056            super(rng, source);
1057        }
1058
1059        @Override
1060        public double sample() {
1061            final double phiby2 = getPhiBy2();
1062            final double w = getOmega();
1063
1064            // Compute some tangents
1065            final double a = phiby2 * SpecialMath.tan2(phiby2);
1066            final double b = eps * phiby2 * SpecialMath.tan2(eps * phiby2);
1067            // Compute some necessary subexpressions
1068            final double da = a * a;
1069            final double db = b * b;
1070            final double a2 = 1 - da;
1071            final double a2p = 1 + da;
1072            final double b2 = 1 - db;
1073            final double b2p = 1 + db;
1074            // Compute coefficient.
1075            final double z = a2p * b2 / (w * a2 * b2p);
1076
1077            // Compute the exponential-type expression
1078            final double alogz = Math.log(z);
1079            final double d = SpecialMath.d2(epsDiv1mEps * alogz) * (alogz * inv1mEps);
1080
1081            // Pre-compute the multiplication factor.
1082            // The numerator may be zero. The denominator is not zero as a2 is bounded to
1083            // above zero when the uniform deviate that generates phiby2 is not 0 or 1.
1084            final double f = (2 * ((a - b) * (1 + a * b))) / (a2 * b2p);
1085
1086            // Compute the stable deviate:
1087            final double x = (1 + eps * d) * f;
1088
1089            // Test the support
1090            if (LOWER < x && x < UPPER) {
1091                return x;
1092            }
1093
1094            // Error correction path.
1095            // Here we use the formula provided by Weron which is easier to correct
1096            // when deviates are extreme or alpha -> 0.
1097            return createSample(phiby2 * 2, w);
1098        }
1099
1100        @Override
1101        public Beta0CMSStableSampler withUniformRandomProvider(UniformRandomProvider rng) {
1102            return new Beta0CMSStableSampler(rng, this);
1103        }
1104    }
1105
1106    /**
1107     * Implement special math functions required by the CMS algorithm.
1108     */
1109    static final class SpecialMath {
1110        /** pi/4. */
1111        private static final double PI_4 = Math.PI / 4;
1112        /** 4/pi. */
1113        private static final double FOUR_PI = 4 / Math.PI;
1114        /** tan2 product constant. */
1115        private static final double P0 = -0.5712939549476836914932149599e10;
1116        /** tan2 product constant. */
1117        private static final double P1 = 0.4946855977542506692946040594e9;
1118        /** tan2 product constant. */
1119        private static final double P2 = -0.9429037070546336747758930844e7;
1120        /** tan2 product constant. */
1121        private static final double P3 = 0.5282725819868891894772108334e5;
1122        /** tan2 product constant. */
1123        private static final double P4 = -0.6983913274721550913090621370e2;
1124        /** tan2 quotient constant. */
1125        private static final double Q0 = -0.7273940551075393257142652672e10;
1126        /** tan2 quotient constant. */
1127        private static final double Q1 = 0.2125497341858248436051062591e10;
1128        /** tan2 quotient constant. */
1129        private static final double Q2 = -0.8000791217568674135274814656e8;
1130        /** tan2 quotient constant. */
1131        private static final double Q3 = 0.8232855955751828560307269007e6;
1132        /** tan2 quotient constant. */
1133        private static final double Q4 = -0.2396576810261093558391373322e4;
1134        /**
1135         * The threshold to switch to using {@link Math#expm1(double)}. The following
1136         * table shows the mean (max) ULP difference between using expm1 and exp using
1137         * random +/-x with different exponents (n=2^30):
1138         *
1139         * <pre>
1140         * x        exp  positive x                 negative x
1141         * 64.0      6   0.10004021506756544  (2)   0.0                   (0)
1142         * 32.0      5   0.11177831795066595  (2)   0.0                   (0)
1143         * 16.0      4   0.0986650362610817   (2)   9.313225746154785E-10 (1)
1144         * 8.0       3   0.09863092936575413  (2)   4.9658119678497314E-6 (1)
1145         * 4.0       2   0.10015273280441761  (2)   4.547201097011566E-4  (1)
1146         * 2.0       1   0.14359260816127062  (2)   0.005623611621558666  (2)
1147         * 1.0       0   0.20160607434809208  (2)   0.03312791418284178   (2)
1148         * 0.5      -1   0.3993037799373269   (2)   0.28186883218586445   (2)
1149         * 0.25     -2   0.6307008266448975   (2)   0.5192863345146179    (2)
1150         * 0.125    -3   1.3862918205559254   (4)   1.386285437270999     (4)
1151         * 0.0625   -4   2.772640804760158    (8)   2.772612397558987     (8)
1152         * </pre>
1153         *
1154         * <p>The threshold of 0.5 has a mean ULP below 0.5 and max ULP of 2. The
1155         * transition is monotonic. Neither is true for the next threshold of 0.25.
1156         */
1157        private static final double SWITCH_TO_EXPM1 = 0.5;
1158
1159        /** No instances. */
1160        private SpecialMath() {}
1161
1162        /**
1163         * Evaluate {@code (exp(x) - 1) / x}. For {@code x} in the range {@code [-inf, inf]} returns
1164         * a result in {@code [0, inf]}.
1165         *
1166         * <ul>
1167         * <li>For {@code x=-inf} this returns {@code 0}.
1168         * <li>For {@code x=0} this returns {@code 1}.
1169         * <li>For {@code x=inf} this returns {@code inf}.
1170         * <li>For {@code x=nan} this returns {@code nan}.
1171         * </ul>
1172         *
1173         * <p> This corrects {@code 0 / 0} and {@code inf / inf} division from
1174         * {@code NaN} to either {@code 1} or the upper bound respectively.
1175         *
1176         * @param x value to evaluate
1177         * @return {@code (exp(x) - 1) / x}.
1178         */
1179        static double d2(double x) {
1180            // Here expm1 is only used when use of expm1 and exp consistently
1181            // compute different results by more than 0.5 ULP.
1182            if (Math.abs(x) < SWITCH_TO_EXPM1) {
1183                // Deliberate comparison to floating-point zero
1184                if (x == 0) {
1185                    // Avoid 0 / 0 error
1186                    return 1.0;
1187                }
1188                return Math.expm1(x) / x;
1189            }
1190            // No use of expm1. Accuracy as x moves away from 0 is not required as the result
1191            // is divided by x and the accuracy of the final result is within a few ULP.
1192            if (x < Double.POSITIVE_INFINITY) {
1193                return (Math.exp(x) - 1) / x;
1194            }
1195            // Upper bound (or NaN)
1196            return x;
1197        }
1198
1199        /**
1200         * Evaluate {@code tan(x) / x}.
1201         *
1202         * <p>For {@code x} in the range {@code [0, pi/4]} this returns
1203         * a value in the range {@code [1, 4 / pi]}.
1204         *
1205         * <p>The following properties are desirable for the CMS algorithm:
1206         *
1207         * <ul>
1208         * <li>For {@code x=0} this returns {@code 1}.
1209         * <li>For {@code x=pi/4} this returns {@code 4/pi}.
1210         * <li>For {@code x=pi/4} this multiplied by {@code x} returns {@code 1}.
1211         * </ul>
1212         *
1213         * <p>This method is called by the CMS algorithm when {@code x < pi/4}.
1214         * In this case the method is almost as accurate as {@code Math.tan(x) / x}, does
1215         * not require checking for {@code x=0} and is faster.
1216         *
1217         * @param x the x
1218         * @return {@code tan(x) / x}.
1219         */
1220        static double tan2(double x) {
1221            if (Math.abs(x) > PI_4) {
1222                // Reduction is not supported. Delegate to the JDK.
1223                return Math.tan(x) / x;
1224            }
1225
1226            // Testing with approximation 4283 from Hart et al, as used in the RSTAB
1227            // routine, showed the method was not accurate enough for use with
1228            // double computation. Hart et al state it has max relative error = 1e-10.66.
1229            // For tan(x) / x with x in [0, pi/4] values outside [1, 4/pi] were computed.
1230            // When testing verses Math.tan(x) the mean ULP difference is 93436.3.
1231
1232            // Approximation 4288 from Hart et al (1968, P. 252).
1233            // Max relative error = 1e-26.68 (for tan(x)).
1234            // When testing verses Math.tan(x) the mean ULP difference is 0.590597.
1235
1236            // The approximation is defined as:
1237            // tan(x*pi/4) = x * P(x^2) / Q(x^2)
1238            //   with P and Q polynomials of x squared.
1239            //
1240            // To create tan(x):
1241            // tan(x) = xi * P(xi^2) / Q(xi^2), xi = x * 4/pi
1242            // tan(x) / x = xi * P(xi^2) / Q(xi^2) / x
1243            // tan(x) / x = 4/pi * (P(xi^2) / Q(xi^2))
1244            //            = P(xi^2) / (pi/4 * Q(xi^2))
1245            // The later has a smaller mean ULP difference to Math.tan(x) / x.
1246            final double xi = x * FOUR_PI;
1247
1248            // Use the power form with a reverse summation order to have smaller
1249            // magnitude terms first. Note: x < 1 so greater powers are smaller.
1250            // This has essentially the same accuracy as the nested form of the polynomials
1251            // for a marginal performance increase. See JMH examples for performance tests.
1252            final double x2 = xi * xi;
1253            final double x4 = x2 * x2;
1254            final double x6 = x4 * x2;
1255            final double x8 = x4 * x4;
1256            return          (x8 * P4 + x6 * P3 + x4 * P2 + x2 * P1 + P0) /
1257                    (PI_4 * (x8 * x2 + x8 * Q4 + x6 * Q3 + x4 * Q2 + x2 * Q1 + Q0));
1258        }
1259    }
1260
1261    /**
1262     * @param rng Generator of uniformly distributed random numbers.
1263     */
1264    StableSampler(UniformRandomProvider rng) {
1265        this.rng = rng;
1266    }
1267
1268    /**
1269     * Generate a sample from a stable distribution.
1270     *
1271     * <p>The distribution uses the 0-parameterization: S(alpha, beta, gamma, delta; 0).
1272     */
1273    @Override
1274    public abstract double sample();
1275
1276    /** {@inheritDoc} */
1277    // Redeclare the signature to return a StableSampler not a SharedStateContinuousSampler
1278    @Override
1279    public abstract StableSampler withUniformRandomProvider(UniformRandomProvider rng);
1280
1281    /**
1282     * Generates a {@code long} value.
1283     * Used by algorithm implementations without exposing access to the RNG.
1284     *
1285     * @return the next random value
1286     */
1287    long nextLong() {
1288        return rng.nextLong();
1289    }
1290
1291    /** {@inheritDoc} */
1292    @Override
1293    public String toString() {
1294        // All variations use the same string representation, i.e. no changes
1295        // for the Gaussian, Levy or Cauchy case.
1296        return "Stable deviate [" + rng.toString() + "]";
1297    }
1298
1299    /**
1300     * Creates a standardized sampler of a stable distribution with zero location and unit scale.
1301     *
1302     * <p>Special cases:
1303     *
1304     * <ul>
1305     * <li>{@code alpha=2} returns a Gaussian distribution sampler with
1306     *     {@code mean=0} and {@code variance=2} (Note: {@code beta} has no effect on the distribution).
1307     * <li>{@code alpha=1} and {@code beta=0} returns a Cauchy distribution sampler with
1308     *     {@code location=0} and {@code scale=1}.
1309     * <li>{@code alpha=0.5} and {@code beta=1} returns a Levy distribution sampler with
1310     *     {@code location=-1} and {@code scale=1}. This location shift is due to the
1311     *     0-parameterization of the stable distribution.
1312     * </ul>
1313     *
1314     * <p>Note: To allow the computation of the stable distribution the parameter alpha
1315     * is validated using {@code 1 - alpha} in the interval {@code [-1, 1)}.
1316     *
1317     * @param rng Generator of uniformly distributed random numbers.
1318     * @param alpha Stability parameter. Must be in the interval {@code (0, 2]}.
1319     * @param beta Skewness parameter. Must be in the interval {@code [-1, 1]}.
1320     * @return the sampler
1321     * @throws IllegalArgumentException if {@code 1 - alpha < -1}; or {@code 1 - alpha >= 1};
1322     * or {@code beta < -1}; or {@code beta > 1}.
1323     */
1324    public static StableSampler of(UniformRandomProvider rng,
1325                                   double alpha,
1326                                   double beta) {
1327        validateParameters(alpha, beta);
1328        return create(rng, alpha, beta);
1329    }
1330
1331    /**
1332     * Creates a sampler of a stable distribution. This applies a transformation to the
1333     * standardized sampler.
1334     *
1335     * <p>The random variable \( X \) has
1336     * the stable distribution \( S(\alpha, \beta, \gamma, \delta; 0) \) if:
1337     *
1338     * <p>\[ X = \gamma Z_0 + \delta \]
1339     *
1340     * <p>where \( Z_0 = S(\alpha, \beta; 0) \) is a standardized stable distribution.
1341     *
1342     * <p>Note: To allow the computation of the stable distribution the parameter alpha
1343     * is validated using {@code 1 - alpha} in the interval {@code [-1, 1)}.
1344     *
1345     * @param rng Generator of uniformly distributed random numbers.
1346     * @param alpha Stability parameter. Must be in the interval {@code (0, 2]}.
1347     * @param beta Skewness parameter. Must be in the interval {@code [-1, 1]}.
1348     * @param gamma Scale parameter. Must be strictly positive and finite.
1349     * @param delta Location parameter. Must be finite.
1350     * @return the sampler
1351     * @throws IllegalArgumentException if {@code 1 - alpha < -1}; or {@code 1 - alpha >= 1};
1352     * or {@code beta < -1}; or {@code beta > 1}; or {@code gamma <= 0}; or
1353     * {@code gamma} or {@code delta} are not finite.
1354     * @see #of(UniformRandomProvider, double, double)
1355     */
1356    public static StableSampler of(UniformRandomProvider rng,
1357                                   double alpha,
1358                                   double beta,
1359                                   double gamma,
1360                                   double delta) {
1361        validateParameters(alpha, beta, gamma, delta);
1362
1363        // Choose the algorithm.
1364        // Reuse the special cases as they have transformation support.
1365
1366        if (alpha == ALPHA_GAUSSIAN) {
1367            // Note: beta has no effect and is ignored.
1368            return new GaussianStableSampler(rng, gamma, delta);
1369        }
1370
1371        // Note: As beta -> 0 the result cannot be computed differently to beta = 0.
1372        if (alpha == ALPHA_CAUCHY && CMSStableSampler.getTau(ALPHA_CAUCHY, beta) == TAU_ZERO) {
1373            return new CauchyStableSampler(rng, gamma, delta);
1374        }
1375
1376        if (alpha == ALPHA_LEVY && Math.abs(beta) == BETA_LEVY) {
1377            // Support mirroring for negative beta by inverting the beta=1 Levy sample
1378            // using a negative gamma. Note: The delta is not mirrored as it is a shift
1379            // applied to the scaled and mirrored distribution.
1380            return new LevyStableSampler(rng, beta * gamma, delta);
1381        }
1382
1383        // Standardized sampler
1384        final StableSampler sampler = create(rng, alpha, beta);
1385        // Transform
1386        return new TransformedStableSampler(sampler, gamma, delta);
1387    }
1388
1389    /**
1390     * Creates a standardized sampler of a stable distribution with zero location and unit scale.
1391     *
1392     * @param rng Generator of uniformly distributed random numbers.
1393     * @param alpha Stability parameter. Must be in the interval {@code (0, 2]}.
1394     * @param beta Skewness parameter. Must be in the interval {@code [-1, 1]}.
1395     * @return the sampler
1396     */
1397    private static StableSampler create(UniformRandomProvider rng,
1398                                        double alpha,
1399                                        double beta) {
1400        // Choose the algorithm.
1401        // The special case samplers have transformation support and use gamma=1.0, delta=0.0.
1402        // As alpha -> 0 the computation increasingly requires correction
1403        // of infinity to the distribution support.
1404
1405        if (alpha == ALPHA_GAUSSIAN) {
1406            // Note: beta has no effect and is ignored.
1407            return new GaussianStableSampler(rng, GAMMA_1, DELTA_0);
1408        }
1409
1410        // Note: As beta -> 0 the result cannot be computed differently to beta = 0.
1411        // This is based on the computation factor tau:
1412        final double tau = CMSStableSampler.getTau(alpha, beta);
1413
1414        if (tau == TAU_ZERO) {
1415            // Symmetric case (beta skew parameter is effectively zero)
1416            if (alpha == ALPHA_CAUCHY) {
1417                return new CauchyStableSampler(rng, GAMMA_1, DELTA_0);
1418            }
1419            if (alpha <= ALPHA_SMALL) {
1420                // alpha -> 0 requires robust error correction
1421                return new Beta0WeronStableSampler(rng, alpha);
1422            }
1423            return new Beta0CMSStableSampler(rng, alpha);
1424        }
1425
1426        // Here beta is significant.
1427
1428        if (alpha == 1) {
1429            return new Alpha1CMSStableSampler(rng, beta);
1430        }
1431
1432        if (alpha == ALPHA_LEVY && Math.abs(beta) == BETA_LEVY) {
1433            // Support mirroring for negative beta by inverting the beta=1 Levy sample
1434            // using a negative gamma. Note: The delta is not mirrored as it is a shift
1435            // applied to the scaled and mirrored distribution.
1436            return new LevyStableSampler(rng, beta, DELTA_0);
1437        }
1438
1439        if (alpha <= ALPHA_SMALL) {
1440            // alpha -> 0 requires robust error correction
1441            return new WeronStableSampler(rng, alpha, beta);
1442        }
1443
1444        return new CMSStableSampler(rng, alpha, beta);
1445    }
1446
1447    /**
1448     * Validate the parameters are in the correct range.
1449     *
1450     * @param alpha Stability parameter. Must be in the interval {@code (0, 2]}.
1451     * @param beta Skewness parameter. Must be in the interval {@code [-1, 1]}.
1452     * @throws IllegalArgumentException if {@code 1 - alpha < -1}; or {@code 1 - alpha >= 1};
1453     * or {@code beta < -1}; or {@code beta > 1}.
1454     */
1455    private static void validateParameters(double alpha, double beta) {
1456        // The epsilon (1-alpha) value must be in the interval [-1, 1).
1457        // Logic inversion will identify NaN
1458        final double eps = 1 - alpha;
1459        if (!(-1 <= eps && eps < 1)) {
1460            throw new IllegalArgumentException("alpha is not in the interval (0, 2]: " + alpha);
1461        }
1462        if (!(-1 <= beta && beta <= 1)) {
1463            throw new IllegalArgumentException("beta is not in the interval [-1, 1]: " + beta);
1464        }
1465    }
1466
1467    /**
1468     * Validate the parameters are in the correct range.
1469     *
1470     * @param alpha Stability parameter. Must be in the interval {@code (0, 2]}.
1471     * @param beta Skewness parameter. Must be in the interval {@code [-1, 1]}.
1472     * @param gamma Scale parameter. Must be strictly positive and finite.
1473     * @param delta Location parameter. Must be finite.
1474     * @throws IllegalArgumentException if {@code 1 - alpha < -1}; or {@code 1 - alpha >= 1};
1475     * or {@code beta < -1}; or {@code beta > 1}; or {@code gamma <= 0}; or
1476     * {@code gamma} or {@code delta} are not finite.
1477     */
1478    private static void validateParameters(double alpha, double beta,
1479                                           double gamma, double delta) {
1480        validateParameters(alpha, beta);
1481
1482        InternalUtils.requireStrictlyPositiveFinite(gamma, "gamma");
1483        InternalUtils.requireFinite(delta, "delta");
1484    }
1485}