View Javadoc
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&auml;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 }