StableSampler.java

  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;

  19. import org.apache.commons.rng.UniformRandomProvider;

  21. /**
  22.  * Samples from a stable distribution.
  23.  *
  24.  * <p>Several different parameterizations exist for the stable distribution.
  25.  * This sampler uses the 0-parameterization distribution described in Nolan (2020) "Univariate Stable
  26.  * Distributions: Models for Heavy Tailed Data". Springer Series in Operations Research and
  27.  * Financial Engineering. Springer. Sections 1.7 and 3.3.3.
  28.  *
  29.  * <p>The random variable \( X \) has
  30.  * the stable distribution \( S(\alpha, \beta, \gamma, \delta; 0) \) if its characteristic
  31.  * function is given by:
  32.  *
  33.  * <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 \\
  34.  * \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} \]
  35.  *
  36.  * <p>The function is continuous with respect to all the parameters; the parameters \( \alpha \)
  37.  * and \( \beta \) determine the shape and the parameters \( \gamma \) and \( \delta \) determine
  38.  * the scale and location. The support of the distribution is:
  39.  *
  40.  * <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 \\
  41.  * (-\infty, \delta + \gamma \tan \frac{\pi \alpha}{2}] &amp; \alpha \lt 1\ and\ \beta = -1 \\
  42.  * (-\infty, \infty) &amp; otherwise \end{cases} \]
  43.  *
  44.  * <p>The implementation uses the Chambers-Mallows-Stuck (CMS) method as described in:
  45.  * <ul>
  46.  *  <li>Chambers, Mallows &amp; Stuck (1976) "A Method for Simulating Stable Random Variables".
  47.  *      Journal of the American Statistical Association. 71 (354): 340–344.
  48.  *  <li>Weron (1996) "On the Chambers-Mallows-Stuck method for simulating skewed stable
  49.  *      random variables". Statistics &amp; Probability Letters. 28 (2): 165–171.
  50.  * </ul>
  51.  *
  52.  * @see <a href="https://en.wikipedia.org/wiki/Stable_distribution">Stable distribution (Wikipedia)</a>
  53.  * @see <a href="https://link.springer.com/book/10.1007/978-3-030-52915-4">Nolan (2020) Univariate Stable Distributions</a>
  54.  * @see <a href="https://doi.org/10.1080%2F01621459.1976.10480344">Chambers et al (1976) JOASA 71: 340-344</a>
  55.  * @see <a href="https://doi.org/10.1016%2F0167-7152%2895%2900113-1">Weron (1996).
  56.  * Statistics &amp; Probability Letters. 28 (2): 165–171.</a>
  57.  * @since 1.4
  58.  */
  59. public abstract class StableSampler implements SharedStateContinuousSampler {
  60.     /** pi / 2. */
  61.     private static final double PI_2 = Math.PI / 2;
  62.     /** The alpha value for the Gaussian case. */
  63.     private static final double ALPHA_GAUSSIAN = 2;
  64.     /** The alpha value for the Cauchy case. */
  65.     private static final double ALPHA_CAUCHY = 1;
  66.     /** The alpha value for the Levy case. */
  67.     private static final double ALPHA_LEVY = 0.5;
  68.     /** The alpha value for the {@code alpha -> 0} to switch to using the Weron formula.
  69.      * Note that small alpha requires robust correction of infinite samples. */
  70.     private static final double ALPHA_SMALL = 0.02;
  71.     /** The beta value for the Levy case. */
  72.     private static final double BETA_LEVY = 1.0;
  73.     /** The gamma value for the normalized case. */
  74.     private static final double GAMMA_1 = 1.0;
  75.     /** The delta value for the normalized case. */
  76.     private static final double DELTA_0 = 0.0;
  77.     /** The tau value for zero. When tau is zero, this is effectively {@code beta = 0}. */
  78.     private static final double TAU_ZERO = 0.0;
  79.     /**
  80.      * The lower support for the distribution.
  81.      * This is the lower bound of {@code (-inf, +inf)}
  82.      * If the sample is not within this bound ({@code lower < x}) then it is either
  83.      * infinite or NaN and the result should be checked.
  84.      */
  85.     private static final double LOWER = Double.NEGATIVE_INFINITY;
  86.     /**
  87.      * The upper support for the distribution.
  88.      * This is the upper bound of {@code (-inf, +inf)}.
  89.      * If the sample is not within this bound ({@code x < upper}) then it is either
  90.      * infinite or NaN and the result should be checked.
  91.      */
  92.     private static final double UPPER = Double.POSITIVE_INFINITY;

  94.     /** Underlying source of randomness. */
  95.     private final UniformRandomProvider rng;

  97.     // Implementation notes
  98.     //
  99.     // 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.

  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;

  222.         /** The exponential sampler. */
  223.         private final ContinuousSampler expSampler;

  225.         /**
  226.          * @param rng Underlying source of randomness
  227.          */
  228.         BaseStableSampler(UniformRandomProvider rng) {
  229.             super(rng);
  230.             expSampler = ZigguratSampler.Exponential.of(rng);
  231.         }

  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.         }

  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.         }

  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.     }

  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;

  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.         }

  322.         @Override
  323.         public double sample() {
  324.             return gamma * sampler.sample() + delta;
  325.         }

  327.         @Override
  328.         public StableSampler withUniformRandomProvider(UniformRandomProvider rng) {
  329.             return new TransformedStableSampler(sampler.withUniformRandomProvider(rng),
  330.                                                 gamma, delta);
  331.         }

  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.     }

  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);

  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;

  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.         }

  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.         }

  380.         @Override
  381.         public double sample() {
  382.             return stdDev * sampler.sample() + mean;
  383.         }

  385.         @Override
  386.         public GaussianStableSampler withUniformRandomProvider(UniformRandomProvider rng) {
  387.             return new GaussianStableSampler(rng, this);
  388.         }
  389.     }

  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;

  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.         }

  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.         }

  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.         }

  437.         @Override
  438.         public CauchyStableSampler withUniformRandomProvider(UniformRandomProvider rng) {
  439.             return new CauchyStableSampler(rng, this);
  440.         }
  441.     }

  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;

  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.         }

  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.         }

  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.         }

  492.         @Override
  493.         public LevyStableSampler withUniformRandomProvider(UniformRandomProvider rng) {
  494.             return new LevyStableSampler(rng, this);
  495.         }
  496.     }

  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;

  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;

  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.             }

  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;


  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.         }

  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.         }

  610.         @Override
  611.         public double sample() {
  612.             final double phi = getPhi();
  613.             final double w = getOmega();
  614.             return createSample(phi, w);
  615.         }

  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.

  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.

  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)

  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.

  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);

  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.             }

  682.             final double x = t1 * t2 * scale + zeta;

  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.         }

  691.         @Override
  692.         public WeronStableSampler withUniformRandomProvider(UniformRandomProvider rng) {
  693.             return new WeronStableSampler(rng, this);
  694.         }
  695.     }

  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;

  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);

  735.             // Compute the RSTAB pre-factor.
  736.             tau = getTau(alpha, beta);
  737.         }

  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.         }

  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;

  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.             }

  781.             return tau;
  782.         }

  784.         @Override
  785.         public double sample() {
  786.             final double phiby2 = getPhiBy2();
  787.             final double w = getOmega();

  789.             // Compute as per the RSTAB routine.

  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.

  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;

  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;

  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);

  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);

  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);

  845.             // Compute the stable deviate:
  846.             final double x = (1 + eps * d) * f + tau * d;

  848.             // Test the support
  849.             if (lower < x && x < upper) {
  850.                 return x;
  851.             }

  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.

  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.         }

  867.         @Override
  868.         public CMSStableSampler withUniformRandomProvider(UniformRandomProvider rng) {
  869.             return new CMSStableSampler(rng, this);
  870.         }
  871.     }

  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;

  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.         }

  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.         }

  904.         @Override
  905.         public double sample() {
  906.             final double phiby2 = getPhiBy2();
  907.             final double w = getOmega();

  909.             // Compute some tangents
  910.             final double a = phiby2 * SpecialMath.tan2(phiby2);

  912.             // Compute some necessary subexpressions
  913.             final double da = a * a;
  914.             final double a2 = 1 - da;
  915.             final double a2p = 1 + da;

  917.             // Compute coefficient.

  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);

  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);

  930.             // Pre-compute the multiplication factor.
  931.             final double f = (2 * (a - phiby2 * tau * (-2 * a))) / a2;

  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.         }

  939.         @Override
  940.         public Alpha1CMSStableSampler withUniformRandomProvider(UniformRandomProvider rng) {
  941.             return new Alpha1CMSStableSampler(rng, this);
  942.         }
  943.     }

  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;

  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.         }

  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.         }

  988.         @Override
  989.         public double sample() {
  990.             final double phi = getPhi();
  991.             final double w = getOmega();
  992.             return createSample(phi, w);
  993.         }

  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.

  1008.             // Compute terms.
  1009.             // Either term can be infinite or 0. Certain parameters compute 0 * inf.
  1010.             // Detect zeros and return as 0.

  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.         }

  1027.         @Override
  1028.         public Beta0WeronStableSampler withUniformRandomProvider(UniformRandomProvider rng) {
  1029.             return new Beta0WeronStableSampler(rng, this);
  1030.         }
  1031.     }

  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.         }

  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.         }

  1059.         @Override
  1060.         public double sample() {
  1061.             final double phiby2 = getPhiBy2();
  1062.             final double w = getOmega();

  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);

  1077.             // Compute the exponential-type expression
  1078.             final double alogz = Math.log(z);
  1079.             final double d = SpecialMath.d2(epsDiv1mEps * alogz) * (alogz * inv1mEps);

  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);

  1086.             // Compute the stable deviate:
  1087.             final double x = (1 + eps * d) * f;

  1089.             // Test the support
  1090.             if (LOWER < x && x < UPPER) {
  1091.                 return x;
  1092.             }

  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.         }

  1100.         @Override
  1101.         public Beta0CMSStableSampler withUniformRandomProvider(UniformRandomProvider rng) {
  1102.             return new Beta0CMSStableSampler(rng, this);
  1103.         }
  1104.     }

  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;

  1159.         /** No instances. */
  1160.         private SpecialMath() {}

  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.         }

  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.             }

  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.

  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.

  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;

  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.     }

  1261.     /**
  1262.      * @param rng Generator of uniformly distributed random numbers.
  1263.      */
  1264.     StableSampler(UniformRandomProvider rng) {
  1265.         this.rng = rng;
  1266.     }

  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();

  1276.     /** {@inheritDoc} */
  1277.     // Redeclare the signature to return a StableSampler not a SharedStateContinuousSampler
  1278.     @Override
  1279.     public abstract StableSampler withUniformRandomProvider(UniformRandomProvider rng);

  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.     }

  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.     }

  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.     }

  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);

  1363.         // Choose the algorithm.
  1364.         // Reuse the special cases as they have transformation support.

  1366.         if (alpha == ALPHA_GAUSSIAN) {
  1367.             // Note: beta has no effect and is ignored.
  1368.             return new GaussianStableSampler(rng, gamma, delta);
  1369.         }

  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.         }

  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.         }

  1383.         // Standardized sampler
  1384.         final StableSampler sampler = create(rng, alpha, beta);
  1385.         // Transform
  1386.         return new TransformedStableSampler(sampler, gamma, delta);
  1387.     }

  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.

  1405.         if (alpha == ALPHA_GAUSSIAN) {
  1406.             // Note: beta has no effect and is ignored.
  1407.             return new GaussianStableSampler(rng, GAMMA_1, DELTA_0);
  1408.         }

  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);

  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.         }

  1426.         // Here beta is significant.

  1428.         if (alpha == 1) {
  1429.             return new Alpha1CMSStableSampler(rng, beta);
  1430.         }

  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.         }

  1439.         if (alpha <= ALPHA_SMALL) {
  1440.             // alpha -> 0 requires robust error correction
  1441.             return new WeronStableSampler(rng, alpha, beta);
  1442.         }

  1444.         return new CMSStableSampler(rng, alpha, beta);
  1445.     }

  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.     }

  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);

  1482.         InternalUtils.requireStrictlyPositiveFinite(gamma, "gamma");
  1483.         InternalUtils.requireFinite(delta, "delta");
  1484.     }
  1485. }
Created with JaCoCo 0.8.12.202403310830