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 ) & \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 ) & \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) & \alpha \lt 1\ and\ \beta = 1 \\ 041 * (-\infty, \delta + \gamma \tan \frac{\pi \alpha}{2}] & \alpha \lt 1\ and\ \beta = -1 \\ 042 * (-\infty, \infty) & otherwise \end{cases} \] 043 * 044 * <p>The implementation uses the Chambers-Mallows-Stuck (CMS) method as described in: 045 * <ul> 046 * <li>Chambers, Mallows & 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 & 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 & 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 & 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 & 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 & 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}