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 * https://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
18 /*
19 * This is not the original file distributed by the Apache Software Foundation
20 * It has been modified by the Hipparchus project
21 */
22 package org.hipparchus.random;
23
24 import org.hipparchus.exception.LocalizedCoreFormats;
25 import org.hipparchus.exception.MathIllegalArgumentException;
26 import org.hipparchus.exception.NullArgumentException;
27 import org.hipparchus.util.FastMath;
28 import org.hipparchus.util.MathUtils;
29 import org.hipparchus.util.SinCos;
30
31 /**
32 * <p>This class provides a stable normalized random generator. It samples from a stable
33 * distribution with location parameter 0 and scale 1.</p>
34 *
35 * <p>The implementation uses the Chambers-Mallows-Stuck method as described in
36 * <i>Handbook of computational statistics: concepts and methods</i> by
37 * James E. Gentle, Wolfgang Härdle, Yuichi Mori.</p>
38 *
39 */
40 public class StableRandomGenerator implements NormalizedRandomGenerator {
41 /** Underlying generator. */
42 private final RandomGenerator generator;
43
44 /** stability parameter */
45 private final double alpha;
46
47 /** skewness parameter */
48 private final double beta;
49
50 /** cache of expression value used in generation */
51 private final double zeta;
52
53 /**
54 * Create a new generator.
55 *
56 * @param generator underlying random generator to use
57 * @param alpha Stability parameter. Must be in range (0, 2]
58 * @param beta Skewness parameter. Must be in range [-1, 1]
59 * @throws NullArgumentException if generator is null
60 * @throws MathIllegalArgumentException if {@code alpha <= 0} or {@code alpha > 2}
61 * or {@code beta < -1} or {@code beta > 1}
62 */
63 public StableRandomGenerator(final RandomGenerator generator,
64 final double alpha, final double beta)
65 throws MathIllegalArgumentException, NullArgumentException {
66 if (generator == null) {
67 throw new NullArgumentException();
68 }
69
70 if (!(alpha > 0d && alpha <= 2d)) {
71 throw new MathIllegalArgumentException(LocalizedCoreFormats.OUT_OF_RANGE_LEFT,
72 alpha, 0, 2);
73 }
74
75 MathUtils.checkRangeInclusive(beta, -1, 1);
76
77 this.generator = generator;
78 this.alpha = alpha;
79 this.beta = beta;
80 if (alpha < 2d && beta != 0d) {
81 zeta = beta * FastMath.tan(FastMath.PI * alpha / 2);
82 } else {
83 zeta = 0d;
84 }
85 }
86
87 /**
88 * Generate a random scalar with zero location and unit scale.
89 *
90 * @return a random scalar with zero location and unit scale
91 */
92 @Override
93 public double nextNormalizedDouble() {
94 // we need 2 uniform random numbers to calculate omega and phi
95 double omega = -FastMath.log(generator.nextDouble());
96 double phi = FastMath.PI * (generator.nextDouble() - 0.5);
97
98 // Normal distribution case (Box-Muller algorithm)
99 if (alpha == 2d) {
100 return FastMath.sqrt(2d * omega) * FastMath.sin(phi);
101 }
102
103 double x;
104 // when beta = 0, zeta is zero as well
105 // Thus we can exclude it from the formula
106 if (beta == 0d) {
107 // Cauchy distribution case
108 if (alpha == 1d) {
109 x = FastMath.tan(phi);
110 } else {
111 x = FastMath.pow(omega * FastMath.cos((1 - alpha) * phi),
112 1d / alpha - 1d) *
113 FastMath.sin(alpha * phi) /
114 FastMath.pow(FastMath.cos(phi), 1d / alpha);
115 }
116 } else {
117 // Generic stable distribution
118 double cosPhi = FastMath.cos(phi);
119 // to avoid rounding errors around alpha = 1
120 if (FastMath.abs(alpha - 1d) > 1e-8) {
121 final SinCos scAlphaPhi = FastMath.sinCos(alpha * phi);
122 final SinCos scInvAlphaPhi = FastMath.sinCos(phi * (1.0 - alpha));
123 x = (scAlphaPhi.sin() + zeta * scAlphaPhi.cos()) / cosPhi *
124 (scInvAlphaPhi.cos() + zeta * scInvAlphaPhi.sin()) /
125 FastMath.pow(omega * cosPhi, (1 - alpha) / alpha);
126 } else {
127 double betaPhi = FastMath.PI / 2 + beta * phi;
128 x = 2d / FastMath.PI * (betaPhi * FastMath.tan(phi) - beta *
129 FastMath.log(FastMath.PI / 2d * omega * cosPhi / betaPhi));
130
131 if (alpha != 1d) {
132 x += beta * FastMath.tan(FastMath.PI * alpha / 2);
133 }
134 }
135 }
136 return x;
137 }
138 }