You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by ah...@apache.org on 2021/06/29 13:47:52 UTC
[commons-rng] branch master updated: RNG-147: Levy distribution
sampler
This is an automated email from the ASF dual-hosted git repository.
aherbert pushed a commit to branch master
in repository https://gitbox.apache.org/repos/asf/commons-rng.git
The following commit(s) were added to refs/heads/master by this push:
new 5f4346c RNG-147: Levy distribution sampler
5f4346c is described below
commit 5f4346ca54c57440d0aee345263ed9678e14be2e
Author: aherbert <ah...@apache.org>
AuthorDate: Tue Jun 29 14:47:48 2021 +0100
RNG-147: Levy distribution sampler
JMH test shows it is 5-10x faster than using an inverse transform.
---
.../ContinuousSamplersPerformance.java | 8 +-
.../distribution/LevySamplersPerformance.java | 135 +++++++++++++++++++++
.../rng/sampling/distribution/LevySampler.java | 99 +++++++++++++++
.../distribution/ContinuousSamplersList.java | 5 +
.../rng/sampling/distribution/LevySamplerTest.java | 98 +++++++++++++++
src/changes/changes.xml | 3 +
6 files changed, 346 insertions(+), 2 deletions(-)
diff --git a/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/sampling/distribution/ContinuousSamplersPerformance.java b/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/sampling/distribution/ContinuousSamplersPerformance.java
index ece69a2..9eb3754 100644
--- a/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/sampling/distribution/ContinuousSamplersPerformance.java
+++ b/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/sampling/distribution/ContinuousSamplersPerformance.java
@@ -26,6 +26,7 @@ import org.apache.commons.rng.sampling.distribution.ChengBetaSampler;
import org.apache.commons.rng.sampling.distribution.ContinuousSampler;
import org.apache.commons.rng.sampling.distribution.ContinuousUniformSampler;
import org.apache.commons.rng.sampling.distribution.InverseTransformParetoSampler;
+import org.apache.commons.rng.sampling.distribution.LevySampler;
import org.apache.commons.rng.sampling.distribution.LogNormalSampler;
import org.apache.commons.rng.sampling.distribution.MarsagliaNormalizedGaussianSampler;
import org.apache.commons.rng.sampling.distribution.ZigguratNormalizedGaussianSampler;
@@ -78,6 +79,7 @@ public class ContinuousSamplersPerformance {
"AhrensDieterExponentialSampler",
"AhrensDieterGammaSampler",
"MarsagliaTsangGammaSampler",
+ "LevySampler",
"LogNormalBoxMullerNormalizedGaussianSampler",
"LogNormalMarsagliaNormalizedGaussianSampler",
"LogNormalZigguratNormalizedGaussianSampler",
@@ -115,8 +117,10 @@ public class ContinuousSamplersPerformance {
// This tests the Ahrens-Dieter algorithm since alpha < 1
sampler = AhrensDieterMarsagliaTsangGammaSampler.of(rng, 0.76, 9.8);
} else if ("MarsagliaTsangGammaSampler".equals(samplerType)) {
- // This tests the Marsaglia-Tsang algorithm since alpha > 1
- sampler = AhrensDieterMarsagliaTsangGammaSampler.of(rng, 12.34, 9.8);
+ // This tests the Marsaglia-Tsang algorithm since alpha > 1
+ sampler = AhrensDieterMarsagliaTsangGammaSampler.of(rng, 12.34, 9.8);
+ } else if ("LevySampler".equals(samplerType)) {
+ sampler = LevySampler.of(rng, 1.23, 4.56);
} else if ("LogNormalBoxMullerNormalizedGaussianSampler".equals(samplerType)) {
sampler = LogNormalSampler.of(BoxMullerNormalizedGaussianSampler.of(rng), 12.3, 4.6);
} else if ("LogNormalMarsagliaNormalizedGaussianSampler".equals(samplerType)) {
diff --git a/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/sampling/distribution/LevySamplersPerformance.java b/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/sampling/distribution/LevySamplersPerformance.java
new file mode 100644
index 0000000..ee3ba4b
--- /dev/null
+++ b/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/sampling/distribution/LevySamplersPerformance.java
@@ -0,0 +1,135 @@
+/*
+ * 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.rng.examples.jmh.sampling.distribution;
+
+import org.apache.commons.math3.distribution.LevyDistribution;
+import org.apache.commons.rng.UniformRandomProvider;
+import org.apache.commons.rng.sampling.distribution.ContinuousInverseCumulativeProbabilityFunction;
+import org.apache.commons.rng.sampling.distribution.ContinuousSampler;
+import org.apache.commons.rng.sampling.distribution.InverseTransformContinuousSampler;
+import org.apache.commons.rng.sampling.distribution.LevySampler;
+import org.apache.commons.rng.simple.RandomSource;
+import org.openjdk.jmh.annotations.Benchmark;
+import org.openjdk.jmh.annotations.BenchmarkMode;
+import org.openjdk.jmh.annotations.Fork;
+import org.openjdk.jmh.annotations.Measurement;
+import org.openjdk.jmh.annotations.Mode;
+import org.openjdk.jmh.annotations.OutputTimeUnit;
+import org.openjdk.jmh.annotations.Param;
+import org.openjdk.jmh.annotations.Scope;
+import org.openjdk.jmh.annotations.Setup;
+import org.openjdk.jmh.annotations.State;
+import org.openjdk.jmh.annotations.Warmup;
+
+import java.util.concurrent.TimeUnit;
+
+/**
+ * Executes a benchmark to compare the speed of generation of Levy distributed random numbers
+ * using different methods.
+ */
+@BenchmarkMode(Mode.AverageTime)
+@OutputTimeUnit(TimeUnit.NANOSECONDS)
+@Warmup(iterations = 5, time = 1, timeUnit = TimeUnit.SECONDS)
+@Measurement(iterations = 5, time = 1, timeUnit = TimeUnit.SECONDS)
+@State(Scope.Benchmark)
+@Fork(value = 1, jvmArgs = {"-server", "-Xms128M", "-Xmx128M"})
+public class LevySamplersPerformance {
+ /**
+ * The value.
+ *
+ * <p>This must NOT be final!</p>
+ */
+ private double value;
+
+ /**
+ * The samplers's to use for testing. Defines the RandomSource and the type of Levy sampler.
+ */
+ @State(Scope.Benchmark)
+ public static class Sources {
+ /**
+ * RNG providers.
+ *
+ * <p>Use different speeds.</p>
+ *
+ * @see <a href="https://commons.apache.org/proper/commons-rng/userguide/rng.html">
+ * Commons RNG user guide</a>
+ */
+ @Param({"SPLIT_MIX_64",
+ "MWC_256",
+ "JDK"})
+ private String randomSourceName;
+
+ /**
+ * The sampler type.
+ */
+ @Param({"LevySampler", "InverseTransformDiscreteSampler"})
+ private String samplerType;
+
+ /** The sampler. */
+ private ContinuousSampler sampler;
+
+ /**
+ * @return the sampler.
+ */
+ public ContinuousSampler getSampler() {
+ return sampler;
+ }
+
+ /** Instantiates sampler. */
+ @Setup
+ public void setup() {
+ final RandomSource randomSource = RandomSource.valueOf(randomSourceName);
+ final UniformRandomProvider rng = randomSource.create();
+ if ("LevySampler".equals(samplerType)) {
+ sampler = LevySampler.of(rng);
+ } else {
+ final ContinuousInverseCumulativeProbabilityFunction levyFunction =
+ new ContinuousInverseCumulativeProbabilityFunction() {
+ /** Use CM for the inverse CDF. null is for the unused RNG. */
+ private final LevyDistribution dist = new LevyDistribution(null, 0.0, 1.0);
+ @Override
+ public double inverseCumulativeProbability(double p) {
+ return dist.inverseCumulativeProbability(p);
+ }
+ };
+ sampler = InverseTransformContinuousSampler.of(rng, levyFunction);
+ }
+ }
+ }
+
+ /**
+ * Baseline for the JMH timing overhead for production of an {@code double} value.
+ *
+ * @return the {@code double} value
+ */
+ @Benchmark
+ public double baseline() {
+ return value;
+ }
+
+ /**
+ * Run the sampler.
+ *
+ * @param sources Source of randomness.
+ * @return the sample value
+ */
+ @Benchmark
+ public double sample(Sources sources) {
+ return sources.getSampler().sample();
+ }
+}
diff --git a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/LevySampler.java b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/LevySampler.java
new file mode 100644
index 0000000..8a5b02c
--- /dev/null
+++ b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/LevySampler.java
@@ -0,0 +1,99 @@
+/*
+ * 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.rng.sampling.distribution;
+
+import org.apache.commons.rng.UniformRandomProvider;
+
+/**
+ * Sampling from a Lévy distribution.
+ *
+ * @since 1.4
+ * @see <a href="https://en.wikipedia.org/wiki/L%C3%A9vy_distribution">Lévy distribution</a>
+ */
+public final class LevySampler implements SharedStateContinuousSampler {
+ /** Gaussian sampler. */
+ private final NormalizedGaussianSampler gaussian;
+ /** Location. */
+ private final double location;
+ /** Scale. */
+ private final double scale;
+ /** RNG (used for the toString() method). */
+ private final UniformRandomProvider rng;
+
+ /**
+ * @param rng Generator of uniformly distributed random numbers.
+ * @param location Location of the Lévy distribution.
+ * @param scale Scale of the Lévy distribution.
+ */
+ private LevySampler(UniformRandomProvider rng,
+ double location,
+ double scale) {
+ this.gaussian = ZigguratNormalizedGaussianSampler.of(rng);
+ this.location = location;
+ this.scale = scale;
+ this.rng = rng;
+ }
+
+ /**
+ * @param rng Generator of uniformly distributed random numbers.
+ * @param source Source to copy.
+ */
+ private LevySampler(UniformRandomProvider rng,
+ LevySampler source) {
+ this.gaussian = ZigguratNormalizedGaussianSampler.of(rng);
+ this.location = source.location;
+ this.scale = source.scale;
+ this.rng = rng;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double sample() {
+ final double n = gaussian.sample();
+ return scale / (n * n) + location;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public String toString() {
+ return "Lévy deviate [" + rng.toString() + "]";
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public LevySampler withUniformRandomProvider(UniformRandomProvider rng) {
+ return new LevySampler(rng, this);
+ }
+
+ /**
+ * Create a new Lévy distribution sampler.
+ *
+ * @param rng Generator of uniformly distributed random numbers.
+ * @param location Location of the Lévy distribution.
+ * @param scale Scale of the Lévy distribution.
+ * @return the sampler
+ * @throws IllegalArgumentException if {@code scale <= 0}
+ */
+ public static LevySampler of(UniformRandomProvider rng,
+ double location,
+ double scale) {
+ if (scale <= 0) {
+ throw new IllegalArgumentException("scale is not strictly positive: " + scale);
+ }
+ return new LevySampler(rng, location, scale);
+ }
+}
diff --git a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/ContinuousSamplersList.java b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/ContinuousSamplersList.java
index 499a071..6802b49 100644
--- a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/ContinuousSamplersList.java
+++ b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/ContinuousSamplersList.java
@@ -138,6 +138,11 @@ public final class ContinuousSamplersList {
final double cLevy = 0.76;
add(LIST, new org.apache.commons.math3.distribution.LevyDistribution(unusedRng, muLevy, cLevy),
RandomSource.TWO_CMRES.create());
+ // Levy sampler
+ add(LIST, new org.apache.commons.math3.distribution.LevyDistribution(unusedRng, muLevy, cLevy),
+ LevySampler.of(RandomSource.JSF_64.create(), muLevy, cLevy));
+ add(LIST, new org.apache.commons.math3.distribution.LevyDistribution(unusedRng, 0.0, 1.0),
+ LevySampler.of(RandomSource.JSF_64.create(), 0.0, 1.0));
// Log normal ("inverse method").
final double scaleLogNormal = 2.345;
diff --git a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/LevySamplerTest.java b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/LevySamplerTest.java
new file mode 100644
index 0000000..4ae4802
--- /dev/null
+++ b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/LevySamplerTest.java
@@ -0,0 +1,98 @@
+/*
+ * 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.rng.sampling.distribution;
+
+import org.apache.commons.rng.UniformRandomProvider;
+import org.apache.commons.rng.core.source64.SplitMix64;
+import org.apache.commons.rng.sampling.RandomAssert;
+import org.apache.commons.rng.simple.RandomSource;
+
+import org.junit.Assert;
+import org.junit.Test;
+
+/**
+ * Test for the {@link LevySampler}.
+ */
+public class LevySamplerTest {
+ /**
+ * Test the constructor with a negative scale.
+ */
+ @Test(expected = IllegalArgumentException.class)
+ public void testConstructorThrowsWithNegativeScale() {
+ final UniformRandomProvider rng = RandomSource.SPLIT_MIX_64.create(0L);
+ final double location = 1;
+ final double scale = -1e-6;
+ LevySampler.of(rng, location, scale);
+ }
+
+ /**
+ * Test the constructor with a zero scale.
+ */
+ @Test(expected = IllegalArgumentException.class)
+ public void testConstructorThrowsWithZeroScale() {
+ final UniformRandomProvider rng = RandomSource.SPLIT_MIX_64.create(0L);
+ final double location = 1;
+ final double scale = 0;
+ LevySampler.of(rng, location, scale);
+ }
+
+ /**
+ * Test the SharedStateSampler implementation.
+ */
+ @Test
+ public void testSharedStateSampler() {
+ final UniformRandomProvider rng1 = RandomSource.SPLIT_MIX_64.create(0L);
+ final UniformRandomProvider rng2 = RandomSource.SPLIT_MIX_64.create(0L);
+ final double location = 4.56;
+ final double scale = 1.23;
+ final LevySampler sampler1 = LevySampler.of(rng1, location, scale);
+ final LevySampler sampler2 = sampler1.withUniformRandomProvider(rng2);
+ RandomAssert.assertProduceSameSequence(sampler1, sampler2);
+ }
+
+ /**
+ * Test the support of the standard distribution is {@code [0, inf)}.
+ */
+ @Test
+ public void testSupport() {
+ final double location = 0.0;
+ final double scale = 1.0;
+ // Force the underlying ZigguratNormalizedGaussianSampler to create 0
+ final LevySampler s1 = LevySampler.of(
+ new SplitMix64(0L) {
+ @Override
+ public long next() {
+ return 0L;
+ }
+ }, location, scale);
+ Assert.assertEquals(Double.POSITIVE_INFINITY, s1.sample(), 0.0);
+
+ // Force the underlying ZigguratNormalizedGaussianSampler to create +inf
+ final LevySampler s2 = LevySampler.of(
+ new SplitMix64(0L) {
+ @Override
+ public long next() {
+ return (Long.MAX_VALUE << 7) & Long.MAX_VALUE;
+ }
+ @Override
+ public double nextDouble() {
+ return 0.0;
+ }
+ }, location, scale);
+ Assert.assertEquals(0.0, s2.sample(), 0.0);
+ }
+}
diff --git a/src/changes/changes.xml b/src/changes/changes.xml
index b848b7f..f74e736 100644
--- a/src/changes/changes.xml
+++ b/src/changes/changes.xml
@@ -77,6 +77,9 @@ re-run tests that fail, and pass the build if they succeed
within the allotted number of reruns (the test will be marked
as 'flaky' in the report).
">
+ <action dev="aherbert" type="add" issue="147">
+ New "LevySampler" to sample from a Levy distribution.
+ </action>
<action dev="aherbert" type="add" issue="145">
"ContinuousUniformSampler": Add optional support for an open interval: (lower, upper).
</action>