Monte Carlo Integration

Monte Carlo methods were some of the first methods I ever used for research, and when I learned about them, they seemed like some sort of magic. Their premise is simple: random numbers can be used to integrate arbitrary shapes embedded into other objects. Nowadays, "Monte Carlo" has become a bit of a catch-all term for methods that use random numbers to produce real results, but it all started as a straightforward method to integrate objects. No matter how you slice it, the idea seems a bit crazy at first. After all, random numbers are random. How could they possibly be used to find non-random values?

Well, imagine you have a square. The area of the square is simple, . Since it's a square, the and are the same, so the formula is technically just . If we embed a circle into the square with a radius (shown below), then its area is . For simplicity, we can also say that .

Now, let's say we want to find the area of the circle without an equation. As we said before, it's embedded in the square, so we should be able to find some ratio of the area of the square to the area of the circle:

This means,

So, if we can find the and we know , we should be able to easily find the . The question is, "How do we easily find the ?" Well, one way is with random sampling. We basically just pick a bunch of points randomly in the square, and each point is tested to see whether it's in the circle or not:

function in_circle(x_pos::Float64, y_pos::Float64)

    # Setting radius to 1 for unit circle
    radius = 1
    return x_pos^2 + y_pos^2 < radius^2
end
(defn in-circle? [pv r]
  "take a vector representing point and radius return true if the
  point is inside the circle"
  (< (->>
      pv
      (map #(* % %))
      (reduce +))
     (* r r)))
bool in_circle(double x, double y, double radius) {
    return x * x + y * y < radius * radius;
}
function inCircle(xPos, yPos) {
  // Setting radius to 1 for unit circle
  let radius = 1;
  return xPos * xPos + yPos * yPos < radius * radius;
}
inCircle (x, y) = x^2 + y^2 < 1
fn in_circle(x: f64, y: f64, radius: f64) -> bool {
    x * x + y * y < radius * radius
}
bool inCircle(real x, real y)
{
    return x ^^ 2 + y ^^ 2 < 1.0;
}
func inCircle(x, y float64) bool {
    return x*x+y*y < 1.0 // the radius of an unit circle is 1.0
}
in_circle <- function(x, y, radius = 1){
        # Return True if the point is in the circle and False otherwise.
        return((x*x + y*y) < radius*radius)
}
private static boolean inCircle(double x, double y) {
    return x * x + y * y < 1;
}
func isInCircle(x: Double, y: Double, radius: Double) -> Bool {
    return (x*x) + (y*y) < radius*radius
}
def in_circle(x, y, radius = 1):
    """Return True if the point is in the circle and False otherwise."""
    return (x*x + y*y) < radius*radius
public bool IsInMe(Point point) => Math.Pow(point.X, 2) + Math.Pow(point.Y, 2) < Math.Pow(Radius, 2);

If it's in the circle, we increase an internal count by one, and in the end,

If we use a small number of points, this will only give us a rough approximation, but as we start adding more and more points, the approximation becomes much, much better (as shown below)!

The true power of Monte Carlo comes from the fact that it can be used to integrate literally any object that can be embedded into the square. As long as you can write some function to tell whether the provided point is inside the shape you want (like in_circle() in this case), you can use Monte Carlo integration! This is obviously an incredibly powerful tool and has been used time and time again for many different areas of physics and engineering. I can guarantee that we will see similar methods crop up all over the place in the future!

Example Code

Monte Carlo methods are famous for their simplicity. It doesn't take too many lines to get something simple going. Here, we are just integrating a circle, like we described above; however, there is a small twist and trick. Instead of calculating the area of the circle, we are instead trying to find the value of , and rather than integrating the entire circle, we are only integrating the upper right quadrant of the circle from . This saves a bit of computation time, but also requires us to multiply our output by .

That's all there is to it! Feel free to submit your version via pull request, and thanks for reading!

# function to determine whether an x, y point is in the unit circle
function in_circle(x_pos::Float64, y_pos::Float64)

    # Setting radius to 1 for unit circle
    radius = 1
    return x_pos^2 + y_pos^2 < radius^2
end

# function to integrate a unit circle to find pi via monte_carlo
function monte_carlo(n::Int64)

    pi_count = 0
    for i = 1:n
        point_x = rand()
        point_y = rand()

        if (in_circle(point_x, point_y))
            pi_count += 1
        end
    end

    # This is using a quarter of the unit sphere in a 1x1 box.
    # The formula is pi = (box_length^2 / radius^2) * (pi_count / n), but we
    #     are only using the upper quadrant and the unit circle, so we can use
    #     4*pi_count/n instead
    pi_estimate = 4*pi_count/n
    println("The pi estimate is: ", pi_estimate)
    println("Percent error is: ", signif(100 * abs(pi_estimate - pi) / pi, 3), " %")
end

monte_carlo(10000000)
(ns monte-carlo.core)

(defn in-circle? [pv r]
  "take a vector representing point and radius return true if the
  point is inside the circle"
  (< (->>
      pv
      (map #(* % %))
      (reduce +))
     (* r r)))
(defn rand-point [r]
  "return a random point from (0,0) inclusive to (r,r) exclusive"
  (repeatedly 2 #(rand r)))
(defn monte-carlo [n r]
  "take the number of random points and radius return an estimate to
pi"
  (*' 4 (/ n)
      (loop [i n count 0]
        (if (zero? i)
          count
          (recur (dec i)
                 (if (in-circle? (rand-point r) r)
                   (inc count)
                   count))))))
(defn -main []
  (let [constant-pi Math/PI
        computed-pi (monte-carlo 10000000 2) ;; this may take some time on lower end machines
        difference (Math/abs (- constant-pi computed-pi))
        error (* 100 (/ difference constant-pi))]
    (println "world's PI: " constant-pi
             ",our PI: "    (double computed-pi)
             ",error: " error)))
#include <math.h>
#include <stdio.h>
#include <stdbool.h>
#include <stdlib.h>
#include <time.h>

bool in_circle(double x, double y, double radius) {
    return x * x + y * y < radius * radius;
}

void monte_carlo(int samples) {
    double radius = 1.0;
    int count = 0;

    for (int i = 0; i < samples; ++i) {
        double x = (double)rand() * 2.0 / RAND_MAX - 1.0;
        double y = (double)rand() * 2.0 / RAND_MAX - 1.0;

        if (in_circle(x, y, radius)) {
            count += 1;
        }
    }

    double estimate = 4.0 * count / (samples * radius * radius);

    printf("The estimate of pi is %f\n", estimate);
    printf("Which has an error of %0.2f%\n", 100 * (M_PI - estimate) / M_PI);
}

int main() {
    srand(time(NULL));

    monte_carlo(1000000);

    return 0;
}
// submitted by xam4lor
function inCircle(xPos, yPos) {
  // Setting radius to 1 for unit circle
  let radius = 1;
  return xPos * xPos + yPos * yPos < radius * radius;
}

function monteCarlo(n) {
  let piCount = 0;

  for (let i = 0; i < n; i++) {
    const pointX = Math.random();
    const pointY = Math.random();

    if (inCircle(pointX, pointY)) {
      piCount++;
    }
  }

  // This is using a quarter of the unit sphere in a 1x1 box.
  // The formula is pi = (boxLength^2 / radius^2) * (piCount / n), but we
  // are only using the upper quadrant and the unit circle, so we can use
  // 4*piCount/n instead
  // piEstimate = 4*piCount/n
  const piEstimate = 4 * piCount / n;
  console.log('Percent error is: %s%', 100 * (Math.PI - piEstimate) / Math.PI);
}

monteCarlo(100000000);
import System.Random

monteCarloPi :: RandomGen g => g -> Int -> Float
monteCarloPi g n = count $ filter inCircle $ makePairs
  where makePairs = take n $ toPair (randomRs (0, 1) g :: [Float])
        toPair (x:y:rest) = (x, y) : toPair rest
        inCircle (x, y) = x^2 + y^2 < 1
        count l = 4 * fromIntegral (length l) / fromIntegral n

main = do
  g <- newStdGen
  let p = monteCarloPi g 100000
  putStrLn $ "Estimated pi: " ++ show p
  putStrLn $ "Percent error: " ++ show (100 * abs (pi - p) / pi)
// Submitted by jess 3jane

extern crate rand;

use std::f64::consts::PI;

fn in_circle(x: f64, y: f64, radius: f64) -> bool {
    x * x + y * y < radius * radius
}

fn monte_carlo(n: i64) -> f64 {
    let mut count = 0;

    for _ in 0..n {
        let x = rand::random();
        let y = rand::random();
        if in_circle(x, y, 1.0) {
            count += 1;
        }
    }

    // return our pi estimate
    (4 * count) as f64 / n as f64
}

fn main() {
    let pi_estimate = monte_carlo(10000000);

    println!("Percent error is {:.3}%", (100.0 * (PI - pi_estimate) / PI));
}
///Returns true if a point (x, y) is in the circle with radius r
bool inCircle(real x, real y)
{
    return x ^^ 2 + y ^^ 2 < 1.0;
}

///Calculate pi using monte carlo
real monteCarloPI(ulong n)
{
    import std.algorithm : count;
    import std.random : uniform01;
    import std.range : generate, take;
    import std.typecons : tuple;

    auto piCount =  generate(() => tuple!("x", "y")(uniform01, uniform01))
        .take(n)
        .count!(a => inCircle(a.x, a.y));
    return piCount * 4.0 / n;
}

void main()
{
    import std.math : abs, PI;
    import std.stdio : writeln;

    auto p = monteCarloPI(100_000);
    writeln("Estimated pi: ", p);
    writeln("Percent error: ", abs(p - PI) * 100 / PI);
}
// Submitted by Chinmaya Mahesh (chin123)

package main

import (
    "fmt"
    "math"
    "math/rand"
    "time"
)

func inCircle(x, y float64) bool {
    return x*x+y*y < 1.0 // the radius of an unit circle is 1.0
}

func monteCarlo(samples int) {
    count := 0
    s := rand.NewSource(time.Now().UnixNano())
    r := rand.New(s)

    for i := 0; i < samples; i++ {
        x, y := r.Float64(), r.Float64()

        if inCircle(x, y) {
            count += 1
        }
    }

    estimate := 4.0 * float64(count) / float64(samples)

    fmt.Println("The estimate of pi is", estimate)
    fmt.Printf("Which has an error of %f%%\n", 100*math.Abs(math.Pi-estimate)/math.Pi)
}

func main() {
    monteCarlo(10000000)
}

in_circle <- function(x, y, radius = 1){
        # Return True if the point is in the circle and False otherwise.
        return((x*x + y*y) < radius*radius)
}

monte_carlo <- function(n_samples, radius = 1){
# Return the estimate of pi using the monte carlo algorithm.

        # Sample x, y from the uniform distribution
        x <- runif(n_samples, 0, radius)
        y <- runif(n_samples, 0, radius)

        # Count the number of points inside the circle
        in_circle_count <- sum(in_circle(x, y, radius))

        # Since we've generated points in upper left quadrant ([0,radius], [0,])
        # We need to multiply the number of points by 4    
        pi_estimate <- 4 * in_circle_count / n_samples

        return(pi_estimate)
}

pi_estimate <- monte_carlo(10000000)
percent_error <- abs(pi - pi_estimate)/pi

print(paste("The estimate of pi is: ", formatC(pi_estimate)))
print(paste("The percent error is:: ", formatC(percent_error)))
//submitted by DominikRafacz
import java.util.Random;

public class MonteCarlo {

    public static void main(String[] args) {
        double piEstimation = monteCarlo(1000);
        System.out.println("Estimated pi value: " + piEstimation);
        System.out.printf("Percent error: " + 100 * (Math.PI - piEstimation) / Math.PI);
    }

    //function to check whether point (x,y) is in unit circle
    private static boolean inCircle(double x, double y) {
        return x * x + y * y < 1;
    }

    //function to calculate estimation of pi
    public static double monteCarlo(int samples) {
        int piCount = 0;

        Random random = new Random();

        for (int i = 0; i < samples; i++) {
            double x = random.nextDouble();
            double y = random.nextDouble();
            if (inCircle(x, y)) {
                piCount++;
            }
        }

        double estimation = 4.0 * piCount / samples;
        return estimation;
    }
}
//Double Extension from YannickSteph on StackOverflow: https://stackoverflow.com/questions/25050309/swift-random-float-between-0-and-1

import Foundation

public extension Double {
    public static var random: Double {
        return Double(arc4random()) / 0xFFFFFFFF // Returns a random double between 0.0 and 1.0, inclusive.
    }

    public static func random(min: Double, max: Double) -> Double {
        return Double.random * (max - min) + min
    }
}

func isInCircle(x: Double, y: Double, radius: Double) -> Bool {
    return (x*x) + (y*y) < radius*radius
}

func monteCarlo(n: Int) -> Double {
    let radius: Double = 1
    var piCount = 0
    var randX: Double
    var randY: Double

    for _ in 0...n {
        randX = Double.random(min: 0, max: radius)
        randY = Double.random(min: 0, max: radius)

        if(isInCircle(x: randX, y: randY, radius: radius)) {
            piCount += 1
        }
    }

    let piEstimate = Double(4 * piCount)/(Double(n))
    return piEstimate
}

func main() {
    let piEstimate = monteCarlo(n: 10000)
    print("Pi estimate is: ", piEstimate)
    print("Percent error is: \(100 * abs(piEstimate - Double.pi)/Double.pi)%")
}


main()
import math
import random


def in_circle(x, y, radius = 1):
    """Return True if the point is in the circle and False otherwise."""
    return (x*x + y*y) < radius*radius

def monte_carlo(n_samples, radius = 1):
    """Return the estimate of pi using the monte carlo algorithm."""
    in_circle_count = 0
    for i in range(n_samples):

        # Sample x, y from the uniform distribution
        x = random.uniform(0, radius)
        y = random.uniform(0, radius)

        # Count the number of points inside the circle
        if(in_circle(x, y, radius)):
            in_circle_count += 1

    # Since we've generated points in upper left quadrant ([0,radius], [0, radius])
    # We need to multiply the number of points by 4    
    pi_estimate = 4 * in_circle_count / (n_samples)

    return pi_estimate

if __name__ == '__main__':

    pi_estimate = monte_carlo(100000)
    percent_error = 100*abs(math.pi - pi_estimate)/math.pi

    print("The estimate of pi is: {:.3f}".format(pi_estimate))
    print("The percent error is: {:.3f}".format(percent_error))

MonteCarlo.cs

using System;

namespace MonteCarloIntegration
{
    public class MonteCarlo
    {
        public double Run(int samples)
        {
            var circle = new Circle(1.0);
            var count = 0;
            var random = new Random();

            for (int i = 0; i < samples; i++)
            {
                var point = new Point(random.NextDouble(), random.NextDouble());
                if (circle.IsInMe(point))
                    count++;
            }

            return 4.0 * count / samples;
        }
    }
}

Circle.cs

using System;

namespace MonteCarloIntegration
{
    public struct Point
    {
        public double X { get; set; }
        public double Y { get; set; }

        public Point(double x, double y)
        {
            this.X = x;
            this.Y = y;
        }
    }

    public class Circle
    {
        public double Radius { get; private set; }

        public Circle(double radius) => this.Radius = Math.Abs(radius);

        public bool IsInMe(Point point) => Math.Pow(point.X, 2) + Math.Pow(point.Y, 2) < Math.Pow(Radius, 2);
    }
}

Program.cs

using System;

namespace MonteCarloIntegration
{
    class Program
    {
        static void Main(string[] args)
        {
            var monteCarlo = new MonteCarlo();
            System.Console.WriteLine("Running with 10,000,000 samples.");
            var piEstimate = monteCarlo.Run(10000000);
            System.Console.WriteLine($"The estimate of pi is: {piEstimate}");
            System.Console.WriteLine($"The percent error is: {Math.Abs(piEstimate - Math.PI) / Math.PI * 100}%");
        }
    }
}

results matching ""

    No results matching ""