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    *      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.math3.random;
18  
19  import org.apache.commons.math3.exception.NullArgumentException;
20  import org.apache.commons.math3.exception.OutOfRangeException;
21  import org.apache.commons.math3.exception.util.LocalizedFormats;
22  import org.apache.commons.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&auml;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 }