Apache
/
commons-math
mirror of https://github.com/apache/commons-math.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

MATH-1615: Functionality is in "Commons RNG" (cf. RNG-137).

This commit is contained in:
Gilles Sadowski 2021-07-11 02:44:10 +02:00
parent 71bfa2daeb
commit 5fee542f82
2 changed files with 0 additions and 263 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

133
commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/random/StableRandomGenerator.java
View File

@ -1,133 +0,0 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math4.legacy.random;
import org.apache.commons.rng.UniformRandomProvider;
import org.apache.commons.math4.legacy.exception.NullArgumentException;
import org.apache.commons.math4.legacy.exception.OutOfRangeException;
import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
import org.apache.commons.math4.legacy.core.jdkmath.AccurateMath;
/**
* <p>This class provides a stable normalized random generator. It samples from a stable
* distribution with location parameter 0 and scale 1.</p>
*
* <p>The implementation uses the Chambers-Mallows-Stuck method as described in
* <i>Handbook of computational statistics: concepts and methods</i> by
* James E. Gentle, Wolfgang H&auml;rdle, Yuichi Mori.</p>
*
* @since 3.0
*/
public class StableRandomGenerator implements NormalizedRandomGenerator {
/** Underlying generator. */
private final UniformRandomProvider generator;
/** stability parameter. */
private final double alpha;
/** skewness parameter. */
private final double beta;
/** cache of expression value used in generation. */
private final double zeta;
/**
* Create a new generator.
*
* @param generator Underlying random generator
* @param alpha Stability parameter. Must be in range (0, 2]
* @param beta Skewness parameter. Must be in range [-1, 1]
* @throws NullArgumentException if generator is null
* @throws OutOfRangeException if {@code alpha <= 0} or {@code alpha > 2}
* or {@code beta < -1} or {@code beta > 1}
*/
public StableRandomGenerator(final UniformRandomProvider generator,
final double alpha, final double beta)
throws NullArgumentException, OutOfRangeException {
if (generator == null) {
throw new NullArgumentException();
}
if (!(alpha > 0d && alpha <= 2d)) {
throw new OutOfRangeException(LocalizedFormats.OUT_OF_RANGE_LEFT,
alpha, 0, 2);
}
if (!(beta >= -1d && beta <= 1d)) {
throw new OutOfRangeException(LocalizedFormats.OUT_OF_RANGE_SIMPLE,
beta, -1, 1);
}
this.generator = generator;
this.alpha = alpha;
this.beta = beta;
if (alpha < 2d && beta != 0d) {
zeta = beta * AccurateMath.tan(AccurateMath.PI * alpha / 2);
} else {
zeta = 0d;
}
}
/**
* Generate a random scalar with zero location and unit scale.
*
* @return a random scalar with zero location and unit scale
*/
@Override
public double nextNormalizedDouble() {
// we need 2 uniform random numbers to calculate omega and phi
double omega = -AccurateMath.log(generator.nextDouble());
double phi = AccurateMath.PI * (generator.nextDouble() - 0.5);
// Normal distribution case (Box-Muller algorithm)
if (alpha == 2d) {
return AccurateMath.sqrt(2d * omega) * AccurateMath.sin(phi);
}
double x;
// when beta = 0, zeta is zero as well
// Thus we can exclude it from the formula
if (beta == 0d) {
// Cauchy distribution case
if (alpha == 1d) {
x = AccurateMath.tan(phi);
} else {
x = AccurateMath.pow(omega * AccurateMath.cos((1 - alpha) * phi),
1d / alpha - 1d) *
AccurateMath.sin(alpha * phi) /
AccurateMath.pow(AccurateMath.cos(phi), 1d / alpha);
}
} else {
// Generic stable distribution
double cosPhi = AccurateMath.cos(phi);
// to avoid rounding errors around alpha = 1
if (AccurateMath.abs(alpha - 1d) > 1e-8) {
double alphaPhi = alpha * phi;
double invAlphaPhi = phi - alphaPhi;
x = (AccurateMath.sin(alphaPhi) + zeta * AccurateMath.cos(alphaPhi)) / cosPhi *
(AccurateMath.cos(invAlphaPhi) + zeta * AccurateMath.sin(invAlphaPhi)) /
AccurateMath.pow(omega * cosPhi, (1 - alpha) / alpha);
} else {
double betaPhi = AccurateMath.PI / 2 + beta * phi;
x = 2d / AccurateMath.PI * (betaPhi * AccurateMath.tan(phi) - beta *
AccurateMath.log(AccurateMath.PI / 2d * omega * cosPhi / betaPhi));
if (alpha != 1d) {
x += beta * AccurateMath.tan(AccurateMath.PI * alpha / 2);
}
}
}
return x;
}
}

130
commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/random/StableRandomGeneratorTest.java
View File

@ -1,130 +0,0 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math4.legacy.random;
import org.apache.commons.math4.legacy.exception.OutOfRangeException;
import org.apache.commons.rng.UniformRandomProvider;
import org.apache.commons.rng.simple.RandomSource;
import org.apache.commons.math4.legacy.stat.StatUtils;
import org.apache.commons.math4.legacy.stat.descriptive.DescriptiveStatistics;
import org.junit.Assert;
import org.junit.Test;
/**
* Tests for the class {@link StableRandomGenerator}.
*/
public class StableRandomGeneratorTest {
private final UniformRandomProvider rg = RandomSource.create(RandomSource.WELL_19937_C, 100);
private static final int sampleSize = 10000;
/**
* Run the double nextDouble() method test Due to leptokurtic property the
* acceptance range is widened.
*
* TODO: verify that tolerance this wide is really OK
*/
@Test
public void testNextDouble() {
StableRandomGenerator generator = new StableRandomGenerator(rg, 1.3,
0.1);
double[] sample = new double[2 * sampleSize];
for (int i = 0; i < sample.length; ++i) {
sample[i] = generator.nextNormalizedDouble();
}
Assert.assertEquals(0.0, StatUtils.mean(sample), 0.3);
}
/**
* If alpha = 2, than it must be Gaussian distribution
*/
@Test
public void testGaussianCase() {
StableRandomGenerator generator = new StableRandomGenerator(rg, 2d, 0.0);
double[] sample = new double[sampleSize];
for (int i = 0; i < sample.length; ++i) {
sample[i] = generator.nextNormalizedDouble();
}
Assert.assertEquals(0.0, StatUtils.mean(sample), 0.02);
Assert.assertEquals(1.0, StatUtils.variance(sample), 0.02);
}
/**
* If alpha = 1, than it must be Cauchy distribution
*/
@Test
public void testCauchyCase() {
StableRandomGenerator generator = new StableRandomGenerator(rg, 1d, 0.0);
DescriptiveStatistics summary = new DescriptiveStatistics();
for (int i = 0; i < sampleSize; ++i) {
double sample = generator.nextNormalizedDouble();
summary.addValue(sample);
}
// Standard Cauchy distribution should have zero median and mode
double median = summary.getPercentile(50);
Assert.assertEquals(0.0, median, 0.2);
}
/**
* Input parameter range tests
*/
@Test
public void testAlphaRangeBelowZero() {
try {
new StableRandomGenerator(rg,
-1.0, 0.0);
Assert.fail("Expected OutOfRangeException");
} catch (OutOfRangeException e) {
Assert.assertEquals(-1.0, e.getArgument());
}
}
@Test
public void testAlphaRangeAboveTwo() {
try {
new StableRandomGenerator(rg,
3.0, 0.0);
Assert.fail("Expected OutOfRangeException");
} catch (OutOfRangeException e) {
Assert.assertEquals(3.0, e.getArgument());
}
}
@Test
public void testBetaRangeBelowMinusOne() {
try {
new StableRandomGenerator(rg,
1.0, -2.0);
Assert.fail("Expected OutOfRangeException");
} catch (OutOfRangeException e) {
Assert.assertEquals(-2.0, e.getArgument());
}
}
@Test
public void testBetaRangeAboveOne() {
try {
new StableRandomGenerator(rg,
1.0, 2.0);
Assert.fail("Expected OutOfRangeException");
} catch (OutOfRangeException e) {
Assert.assertEquals(2.0, e.getArgument());
}
}
}