Monte-Carlo integration as a usage example.
Project: http://git-wip-us.apache.org/repos/asf/commons-rng/repo Commit: http://git-wip-us.apache.org/repos/asf/commons-rng/commit/a9092d38 Tree: http://git-wip-us.apache.org/repos/asf/commons-rng/tree/a9092d38 Diff: http://git-wip-us.apache.org/repos/asf/commons-rng/diff/a9092d38 Branch: refs/heads/master Commit: a9092d38d9d00e44cf178a8f164fa693705e87b7 Parents: 6ee7951 Author: Gilles <er...@apache.org> Authored: Fri Aug 12 18:03:18 2016 +0200 Committer: Gilles <er...@apache.org> Committed: Fri Aug 12 18:03:18 2016 +0200 ---------------------------------------------------------------------- .../apache/commons/rng/userguide/ComputePi.java | 91 ++++++++++++++++++++ .../rng/userguide/MonteCarloIntegration.java | 90 +++++++++++++++++++ 2 files changed, 181 insertions(+) ---------------------------------------------------------------------- http://git-wip-us.apache.org/repos/asf/commons-rng/blob/a9092d38/src/userguide/java/org/apache/commons/rng/userguide/ComputePi.java ---------------------------------------------------------------------- diff --git a/src/userguide/java/org/apache/commons/rng/userguide/ComputePi.java b/src/userguide/java/org/apache/commons/rng/userguide/ComputePi.java new file mode 100644 index 0000000..5ab4c8a --- /dev/null +++ b/src/userguide/java/org/apache/commons/rng/userguide/ComputePi.java @@ -0,0 +1,91 @@ +/* + * 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.userguide; + +import org.apache.commons.rng.UniformRandomProvider; +import org.apache.commons.rng.RandomSource; + +/** + * Computation of \( \pi \) using Monte-Carlo integration. + * + * The computation estimates the value by computing the probability that + * a point \( p = (x, y) \) will lie in the circle of radius \( r = 1 \) + * inscribed in the square of side \( r = 1 \). + * The probability could be computed by \[ area_{circle} / area_{square} \], + * where \( area_{circle} = \pi * r^2 \) and \( area_{square} = 4 r^2 \). + * Hence, the probability is \( \frac{\pi}{4} \). + * + * The Monte Carlo simulation will produce \( N \) points. + * Defining \( N_c \) as the number of point that satisfy \( x^2 + y^2 <= 1 \), + * we will have \( \frac{N_c}{N} \approx \frac{\pi}{4} \). + */ +public class ComputePi extends MonteCarloIntegration { + /** Domain dimension. */ + private static final int DIMENSION = 2; + + /** + * Simulation constructor. + * + * @param source RNG algorithm. + */ + public ComputePi(RandomSource source) { + super(source, DIMENSION); + } + + /** + * Arguments: + * <ul> + * <li> + * Error on the approximation: computation will stop when the + * absolute difference between successive approximations is lower + * than the given value. + * </li> + * <li> + * {@link RandomSource Random source identifier}. + * </li> + * </ul> + */ + public static void main(String[] args) { + if (args.length != 2) { + throw new IllegalStateException("Missing arguments"); + } + + final long numPoints = Long.valueOf(args[0]); + final RandomSource randomSource = RandomSource.valueOf(args[1]); + + final ComputePi piApp = new ComputePi(randomSource); + final double piMC = piApp.compute(numPoints); + + System.out.println("After generating " + (DIMENSION * numPoints) + + " random numbers, the error on ð is " + Math.abs(piMC - Math.PI)); + } + + /** + * @param n Number of random points to generate. + * @return the approximate value of pi. + */ + public double compute(long numPoints) { + return 4 * integrate(numPoints); + } + + /** {@inheritDoc} */ + @Override + protected boolean isInside(double ... rand) { + final double r2 = rand[0] * rand[0] + rand[1] * rand[1]; + return r2 <= 1; + } +} http://git-wip-us.apache.org/repos/asf/commons-rng/blob/a9092d38/src/userguide/java/org/apache/commons/rng/userguide/MonteCarloIntegration.java ---------------------------------------------------------------------- diff --git a/src/userguide/java/org/apache/commons/rng/userguide/MonteCarloIntegration.java b/src/userguide/java/org/apache/commons/rng/userguide/MonteCarloIntegration.java new file mode 100644 index 0000000..e7e17a8 --- /dev/null +++ b/src/userguide/java/org/apache/commons/rng/userguide/MonteCarloIntegration.java @@ -0,0 +1,90 @@ +/* + * 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.userguide; + +import org.apache.commons.rng.UniformRandomProvider; +import org.apache.commons.rng.RandomSource; + +/** + * <a href="https://en.wikipedia.org/wiki/Monte_Carlo_integration";>Monte-Carlo method</a> + * for approximating an integral on a n-dimensional unit cube. + */ +public abstract class MonteCarloIntegration { + /** RNG. */ + private final UniformRandomProvider rng; + /** Integration domain dimension. */ + private final int dimension; + + /** + * Simulation constructor. + * + * @param source RNG algorithm. + * @param dimension Integration domain dimension. + */ + public MonteCarloIntegration(RandomSource source, + int dimension) { + this.rng = RandomSource.create(source); + this.dimension = dimension; + } + + /** + * Run the Monte-Carlo integration. + * + * @param n Number of random points to generate. + * @return the integral. + */ + public double integrate(long n) { + double result = 0; + long inside = 0; + long total = 0; + while (total < n) { + if (isInside(generateU01())) { + ++inside; + } + + ++total; + result = inside / (double) total; + } + + return result; + } + + /** + * Indicates whether the given points is inside the region whose + * integral is computed. + * + * @param point Point whose coordinates are random numbers uniformly + * distributed in the unit interval. + * @return {@code true} if the {@code point} is inside. + */ + protected abstract boolean isInside(double ... point); + + /** + * @return a value from a random sequence uniformly distributed + * in the {@code [0, 1)} interval. + */ + private double[] generateU01() { + final double[] rand = new double[dimension]; + + for (int i = 0; i < dimension; i++) { + rand[i] = rng.nextDouble(); + } + + return rand; + } +}
