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   * @version $Id: StableRandomGenerator.java 1394763 2012-10-05 19:54:00Z psteitz $
33   * @since 3.0
34   */
35  public class StableRandomGenerator implements NormalizedRandomGenerator {
36      /** Underlying generator. */
37      private final RandomGenerator generator;
38  
39      /** stability parameter */
40      private final double alpha;
41  
42      /** skewness parameter */
43      private final double beta;
44  
45      /** cache of expression value used in generation */
46      private final double zeta;
47  
48      /**
49       * Create a new generator.
50       *
51       * @param generator underlying random generator to use
52       * @param alpha Stability parameter. Must be in range (0, 2]
53       * @param beta Skewness parameter. Must be in range [-1, 1]
54       * @throws NullArgumentException if generator is null
55       * @throws OutOfRangeException if {@code alpha <= 0} or {@code alpha > 2}
56       * or {@code beta < -1} or {@code beta > 1}
57       */
58      public StableRandomGenerator(final RandomGenerator generator,
59                                   final double alpha, final double beta)
60          throws NullArgumentException, OutOfRangeException {
61          if (generator == null) {
62              throw new NullArgumentException();
63          }
64  
65          if (!(alpha > 0d && alpha <= 2d)) {
66              throw new OutOfRangeException(LocalizedFormats.OUT_OF_RANGE_LEFT,
67                      alpha, 0, 2);
68          }
69  
70          if (!(beta >= -1d && beta <= 1d)) {
71              throw new OutOfRangeException(LocalizedFormats.OUT_OF_RANGE_SIMPLE,
72                      beta, -1, 1);
73          }
74  
75          this.generator = generator;
76          this.alpha = alpha;
77          this.beta = beta;
78          if (alpha < 2d && beta != 0d) {
79              zeta = beta * FastMath.tan(FastMath.PI * alpha / 2);
80          } else {
81              zeta = 0d;
82          }
83      }
84  
85      /**
86       * Generate a random scalar with zero location and unit scale.
87       *
88       * @return a random scalar with zero location and unit scale
89       */
90      public double nextNormalizedDouble() {
91          // we need 2 uniform random numbers to calculate omega and phi
92          double omega = -FastMath.log(generator.nextDouble());
93          double phi = FastMath.PI * (generator.nextDouble() - 0.5);
94  
95          // Normal distribution case (Box-Muller algorithm)
96          if (alpha == 2d) {
97              return FastMath.sqrt(2d * omega) * FastMath.sin(phi);
98          }
99  
100         double x;
101         // when beta = 0, zeta is zero as well
102         // Thus we can exclude it from the formula
103         if (beta == 0d) {
104             // Cauchy distribution case
105             if (alpha == 1d) {
106                 x = FastMath.tan(phi);
107             } else {
108                 x = FastMath.pow(omega * FastMath.cos((1 - alpha) * phi),
109                     1d / alpha - 1d) *
110                     FastMath.sin(alpha * phi) /
111                     FastMath.pow(FastMath.cos(phi), 1d / alpha);
112             }
113         } else {
114             // Generic stable distribution
115             double cosPhi = FastMath.cos(phi);
116             // to avoid rounding errors around alpha = 1
117             if (FastMath.abs(alpha - 1d) > 1e-8) {
118                 double alphaPhi = alpha * phi;
119                 double invAlphaPhi = phi - alphaPhi;
120                 x = (FastMath.sin(alphaPhi) + zeta * FastMath.cos(alphaPhi)) / cosPhi *
121                     (FastMath.cos(invAlphaPhi) + zeta * FastMath.sin(invAlphaPhi)) /
122                      FastMath.pow(omega * cosPhi, (1 - alpha) / alpha);
123             } else {
124                 double betaPhi = FastMath.PI / 2 + beta * phi;
125                 x = 2d / FastMath.PI * (betaPhi * FastMath.tan(phi) - beta *
126                     FastMath.log(FastMath.PI / 2d * omega * cosPhi / betaPhi));
127 
128                 if (alpha != 1d) {
129                     x = x + beta * FastMath.tan(FastMath.PI * alpha / 2);
130                 }
131             }
132         }
133         return x;
134     }
135 }