001 /*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements. See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License. You may obtain a copy of the License at
008 *
009 * http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017 package org.apache.commons.math3.random;
018
019 import org.apache.commons.math3.exception.NullArgumentException;
020 import org.apache.commons.math3.exception.OutOfRangeException;
021 import org.apache.commons.math3.exception.util.LocalizedFormats;
022 import org.apache.commons.math3.util.FastMath;
023
024 /**
025 * <p>This class provides a stable normalized random generator. It samples from a stable
026 * distribution with location parameter 0 and scale 1.</p>
027 *
028 * <p>The implementation uses the Chambers-Mallows-Stuck method as described in
029 * <i>Handbook of computational statistics: concepts and methods</i> by
030 * James E. Gentle, Wolfgang Härdle, Yuichi Mori.</p>
031 *
032 * @version $Id: StableRandomGenerator.java 1394763 2012-10-05 19:54:00Z psteitz $
033 * @since 3.0
034 */
035 public class StableRandomGenerator implements NormalizedRandomGenerator {
036 /** Underlying generator. */
037 private final RandomGenerator generator;
038
039 /** stability parameter */
040 private final double alpha;
041
042 /** skewness parameter */
043 private final double beta;
044
045 /** cache of expression value used in generation */
046 private final double zeta;
047
048 /**
049 * Create a new generator.
050 *
051 * @param generator underlying random generator to use
052 * @param alpha Stability parameter. Must be in range (0, 2]
053 * @param beta Skewness parameter. Must be in range [-1, 1]
054 * @throws NullArgumentException if generator is null
055 * @throws OutOfRangeException if {@code alpha <= 0} or {@code alpha > 2}
056 * or {@code beta < -1} or {@code beta > 1}
057 */
058 public StableRandomGenerator(final RandomGenerator generator,
059 final double alpha, final double beta)
060 throws NullArgumentException, OutOfRangeException {
061 if (generator == null) {
062 throw new NullArgumentException();
063 }
064
065 if (!(alpha > 0d && alpha <= 2d)) {
066 throw new OutOfRangeException(LocalizedFormats.OUT_OF_RANGE_LEFT,
067 alpha, 0, 2);
068 }
069
070 if (!(beta >= -1d && beta <= 1d)) {
071 throw new OutOfRangeException(LocalizedFormats.OUT_OF_RANGE_SIMPLE,
072 beta, -1, 1);
073 }
074
075 this.generator = generator;
076 this.alpha = alpha;
077 this.beta = beta;
078 if (alpha < 2d && beta != 0d) {
079 zeta = beta * FastMath.tan(FastMath.PI * alpha / 2);
080 } else {
081 zeta = 0d;
082 }
083 }
084
085 /**
086 * Generate a random scalar with zero location and unit scale.
087 *
088 * @return a random scalar with zero location and unit scale
089 */
090 public double nextNormalizedDouble() {
091 // we need 2 uniform random numbers to calculate omega and phi
092 double omega = -FastMath.log(generator.nextDouble());
093 double phi = FastMath.PI * (generator.nextDouble() - 0.5);
094
095 // Normal distribution case (Box-Muller algorithm)
096 if (alpha == 2d) {
097 return FastMath.sqrt(2d * omega) * FastMath.sin(phi);
098 }
099
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 }