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 */
017package org.apache.commons.math4.random;
018
019import org.apache.commons.rng.UniformRandomProvider;
020import org.apache.commons.math4.exception.NullArgumentException;
021import org.apache.commons.math4.exception.OutOfRangeException;
022import org.apache.commons.math4.exception.util.LocalizedFormats;
023import org.apache.commons.math4.util.FastMath;
024
025/**
026 * <p>This class provides a stable normalized random generator. It samples from a stable
027 * distribution with location parameter 0 and scale 1.</p>
028 *
029 * <p>The implementation uses the Chambers-Mallows-Stuck method as described in
030 * <i>Handbook of computational statistics: concepts and methods</i> by
031 * James E. Gentle, Wolfgang H&auml;rdle, Yuichi Mori.</p>
032 *
033 * @since 3.0
034 */
035public class StableRandomGenerator implements NormalizedRandomGenerator {
036    /** Underlying generator. */
037    private final UniformRandomProvider generator;
038    /** stability parameter */
039    private final double alpha;
040    /** skewness parameter */
041    private final double beta;
042    /** cache of expression value used in generation */
043    private final double zeta;
044
045    /**
046     * Create a new generator.
047     *
048     * @param generator Underlying random generator
049     * @param alpha Stability parameter. Must be in range (0, 2]
050     * @param beta Skewness parameter. Must be in range [-1, 1]
051     * @throws NullArgumentException if generator is null
052     * @throws OutOfRangeException if {@code alpha <= 0} or {@code alpha > 2}
053     * or {@code beta < -1} or {@code beta > 1}
054     */
055    public StableRandomGenerator(final UniformRandomProvider generator,
056                                 final double alpha, final double beta)
057        throws NullArgumentException, OutOfRangeException {
058        if (generator == null) {
059            throw new NullArgumentException();
060        }
061
062        if (!(alpha > 0d && alpha <= 2d)) {
063            throw new OutOfRangeException(LocalizedFormats.OUT_OF_RANGE_LEFT,
064                    alpha, 0, 2);
065        }
066
067        if (!(beta >= -1d && beta <= 1d)) {
068            throw new OutOfRangeException(LocalizedFormats.OUT_OF_RANGE_SIMPLE,
069                    beta, -1, 1);
070        }
071
072        this.generator = generator;
073        this.alpha = alpha;
074        this.beta = beta;
075        if (alpha < 2d && beta != 0d) {
076            zeta = beta * FastMath.tan(FastMath.PI * alpha / 2);
077        } else {
078            zeta = 0d;
079        }
080    }
081
082    /**
083     * Generate a random scalar with zero location and unit scale.
084     *
085     * @return a random scalar with zero location and unit scale
086     */
087    @Override
088    public double nextNormalizedDouble() {
089        // we need 2 uniform random numbers to calculate omega and phi
090        double omega = -FastMath.log(generator.nextDouble());
091        double phi = FastMath.PI * (generator.nextDouble() - 0.5);
092
093        // Normal distribution case (Box-Muller algorithm)
094        if (alpha == 2d) {
095            return FastMath.sqrt(2d * omega) * FastMath.sin(phi);
096        }
097
098        double x;
099        // when beta = 0, zeta is zero as well
100        // Thus we can exclude it from the formula
101        if (beta == 0d) {
102            // Cauchy distribution case
103            if (alpha == 1d) {
104                x = FastMath.tan(phi);
105            } else {
106                x = FastMath.pow(omega * FastMath.cos((1 - alpha) * phi),
107                    1d / alpha - 1d) *
108                    FastMath.sin(alpha * phi) /
109                    FastMath.pow(FastMath.cos(phi), 1d / alpha);
110            }
111        } else {
112            // Generic stable distribution
113            double cosPhi = FastMath.cos(phi);
114            // to avoid rounding errors around alpha = 1
115            if (FastMath.abs(alpha - 1d) > 1e-8) {
116                double alphaPhi = alpha * phi;
117                double invAlphaPhi = phi - alphaPhi;
118                x = (FastMath.sin(alphaPhi) + zeta * FastMath.cos(alphaPhi)) / cosPhi *
119                    (FastMath.cos(invAlphaPhi) + zeta * FastMath.sin(invAlphaPhi)) /
120                     FastMath.pow(omega * cosPhi, (1 - alpha) / alpha);
121            } else {
122                double betaPhi = FastMath.PI / 2 + beta * phi;
123                x = 2d / FastMath.PI * (betaPhi * FastMath.tan(phi) - beta *
124                    FastMath.log(FastMath.PI / 2d * omega * cosPhi / betaPhi));
125
126                if (alpha != 1d) {
127                    x += beta * FastMath.tan(FastMath.PI * alpha / 2);
128                }
129            }
130        }
131        return x;
132    }
133}