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