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.math3.random; 018 019import org.apache.commons.math3.exception.NullArgumentException; 020import org.apache.commons.math3.exception.OutOfRangeException; 021import org.apache.commons.math3.exception.util.LocalizedFormats; 022import 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 * @since 3.0 033 */ 034public class StableRandomGenerator implements NormalizedRandomGenerator { 035 /** Underlying generator. */ 036 private final RandomGenerator generator; 037 038 /** stability parameter */ 039 private final double alpha; 040 041 /** skewness parameter */ 042 private final double beta; 043 044 /** cache of expression value used in generation */ 045 private final double zeta; 046 047 /** 048 * Create a new generator. 049 * 050 * @param generator underlying random generator to use 051 * @param alpha Stability parameter. Must be in range (0, 2] 052 * @param beta Skewness parameter. Must be in range [-1, 1] 053 * @throws NullArgumentException if generator is null 054 * @throws OutOfRangeException if {@code alpha <= 0} or {@code alpha > 2} 055 * or {@code beta < -1} or {@code beta > 1} 056 */ 057 public StableRandomGenerator(final RandomGenerator generator, 058 final double alpha, final double beta) 059 throws NullArgumentException, OutOfRangeException { 060 if (generator == null) { 061 throw new NullArgumentException(); 062 } 063 064 if (!(alpha > 0d && alpha <= 2d)) { 065 throw new OutOfRangeException(LocalizedFormats.OUT_OF_RANGE_LEFT, 066 alpha, 0, 2); 067 } 068 069 if (!(beta >= -1d && beta <= 1d)) { 070 throw new OutOfRangeException(LocalizedFormats.OUT_OF_RANGE_SIMPLE, 071 beta, -1, 1); 072 } 073 074 this.generator = generator; 075 this.alpha = alpha; 076 this.beta = beta; 077 if (alpha < 2d && beta != 0d) { 078 zeta = beta * FastMath.tan(FastMath.PI * alpha / 2); 079 } else { 080 zeta = 0d; 081 } 082 } 083 084 /** 085 * Generate a random scalar with zero location and unit scale. 086 * 087 * @return a random scalar with zero location and unit scale 088 */ 089 public double nextNormalizedDouble() { 090 // we need 2 uniform random numbers to calculate omega and phi 091 double omega = -FastMath.log(generator.nextDouble()); 092 double phi = FastMath.PI * (generator.nextDouble() - 0.5); 093 094 // Normal distribution case (Box-Muller algorithm) 095 if (alpha == 2d) { 096 return FastMath.sqrt(2d * omega) * FastMath.sin(phi); 097 } 098 099 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}