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