1 /* 2 * Licensed to the Apache Software Foundation (ASF) under one or more 3 * contributor license agreements. See the NOTICE file distributed with 4 * this work for additional information regarding copyright ownership. 5 * The ASF licenses this file to You under the Apache License, Version 2.0 6 * (the "License"); you may not use this file except in compliance with 7 * the License. You may obtain a copy of the License at 8 * 9 * http://www.apache.org/licenses/LICENSE-2.0 10 * 11 * Unless required by applicable law or agreed to in writing, software 12 * distributed under the License is distributed on an "AS IS" BASIS, 13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14 * See the License for the specific language governing permissions and 15 * limitations under the License. 16 */ 17 18 package org.apache.commons.rng.examples.quadrature; 19 20 import org.apache.commons.rng.UniformRandomProvider; 21 import org.apache.commons.rng.simple.RandomSource; 22 23 /** 24 * <a href="https://en.wikipedia.org/wiki/Monte_Carlo_integration">Monte-Carlo method</a> 25 * for approximating an integral on a n-dimensional unit cube. 26 */ 27 public abstract class MonteCarloIntegration { 28 /** RNG. */ 29 private final UniformRandomProvider rng; 30 /** Integration domain dimension. */ 31 private final int dimension; 32 33 /** 34 * Simulation constructor. 35 * 36 * @param source RNG algorithm. 37 * @param dimension Integration domain dimension. 38 */ 39 public MonteCarloIntegration(RandomSource source, 40 int dimension) { 41 this.rng = RandomSource.create(source); 42 this.dimension = dimension; 43 } 44 45 /** 46 * Run the Monte-Carlo integration. 47 * 48 * @param n Number of random points to generate. 49 * @return the integral. 50 */ 51 public double integrate(long n) { 52 double result = 0; 53 long inside = 0; 54 long total = 0; 55 while (total < n) { 56 if (isInside(generateU01())) { 57 ++inside; 58 } 59 60 ++total; 61 result = inside / (double) total; 62 } 63 64 return result; 65 } 66 67 /** 68 * Indicates whether the given points is inside the region whose 69 * integral is computed. 70 * 71 * @param point Point whose coordinates are random numbers uniformly 72 * distributed in the unit interval. 73 * @return {@code true} if the {@code point} is inside. 74 */ 75 protected abstract boolean isInside(double... point); 76 77 /** 78 * @return a value from a random sequence uniformly distributed 79 * in the {@code [0, 1)} interval. 80 */ 81 private double[] generateU01() { 82 final double[] rand = new double[dimension]; 83 84 for (int i = 0; i < dimension; i++) { 85 rand[i] = rng.nextDouble(); 86 } 87 88 return rand; 89 } 90 }