Monte Carlo Pi Estimation

Description

Monte Carlo Pi Estimation is a randomized algorithm that approximates the value of π using random sampling. The method uses the relationship between the area of a circle and its circumscribing square to estimate π through random point generation and probability.

How It Works

Consider a quarter circle of radius r inside a square of side r
Generate random points within the square
Count points that fall inside the quarter circle
Ratio of points inside circle to total points ≈ π/4
Multiply ratio by 4 to estimate π

The accuracy of the estimation improves with the number of points generated, following a rate of approximately 1/√n where n is the number of points.

Visualization

Speed:

Points Inside: 0

Total Points: 0

π Estimate: 0

Error: 0%

Enter number of points and click Start

Complexity Analysis

Time Complexity

Point Generation: O(n)
- Each point requires constant time to generate
- n is the number of points requested
Point Checking: O(n)
- Each point requires constant time to check
- Uses Pythagorean theorem for distance

Space Complexity

Variables: O(1)
- Only counters and current point stored
- No additional data structures needed

Error Analysis

Expected Error: O(1/√n)
- Error decreases with square root of points
- Doubling accuracy requires 4x points

Advantages and Disadvantages

Advantages

Simple to understand and implement
Naturally parallelizable
Can be stopped at any time
Demonstrates probability concepts

Disadvantages

Slow convergence rate
Results vary between runs
Not suitable for high precision
Requires good random number generator

Implementation


import random
import math
from typing import Tuple, List
import matplotlib.pyplot as plt

class MonteCarloPi:
    def __init__(self, seed: int = None):
        """Initialize the Monte Carlo Pi estimator."""
        if seed is not None:
            random.seed(seed)
        
    def estimate_pi(self, num_points: int) -> Tuple[float, List[float]]:
        """
        Estimate π using Monte Carlo method.
        
        Args:
            num_points: Number of random points to generate
            
        Returns:
            Tuple containing:
            - Final π estimate
            - List of intermediate π estimates
        """
        points_inside = 0
        estimates = []
        
        for i in range(1, num_points + 1):
            # Generate random point in [0, 1] x [0, 1]
            x = random.random()
            y = random.random()
            
            # Check if point lies inside quarter circle
            if math.sqrt(x*x + y*y) <= 1:
                points_inside += 1
            
            # Calculate current estimate
            pi_estimate = 4 * points_inside / i
            estimates.append(pi_estimate)
        
        return estimates[-1], estimates
    
    def visualize_convergence(self, estimates: List[float]) -> None:
        """Plot convergence of π estimates."""
        plt.figure(figsize=(10, 6))
        plt.plot(estimates, label='π estimate')
        plt.axhline(y=math.pi, color='r', linestyle='--', label='Actual π')
        plt.xlabel('Number of Points')
        plt.ylabel('Estimated π')
        plt.title('Convergence of Monte Carlo π Estimation')
        plt.legend()
        plt.grid(True)
        plt.show()

# Example usage
if __name__ == "__main__":
    # Create estimator with fixed seed for reproducibility
    mc = MonteCarloPi(seed=42)
    
    # Run estimation with different numbers of points
    for points in [1000, 10000, 100000]:
        pi_estimate, estimates = mc.estimate_pi(points)
        error = abs(pi_estimate - math.pi) / math.pi * 100
        
        print(f"\nNumber of points: {points}")
        print(f"π estimate: {pi_estimate:.6f}")
        print(f"Actual π: {math.pi:.6f}")
        print(f"Error: {error:.2f}%")
        
        # Visualize convergence
        mc.visualize_convergence(estimates)


#include <iostream>
#include <random>
#include <chrono>
#include <vector>
#include <cmath>
#include <iomanip>

class MonteCarloPi {
private:
    std::mt19937 gen;
    std::uniform_real_distribution<double> dis;

public:
    MonteCarloPi(int seed = -1) : dis(0.0, 1.0) {
        // Use provided seed or generate from time
        if (seed == -1) {
            seed = std::chrono::system_clock::now()
                .time_since_epoch()
                .count();
        }
        gen.seed(seed);
    }
    
    std::pair<double, std::vector<double>> estimate_pi(int num_points) {
        int points_inside = 0;
        std::vector<double> estimates;
        estimates.reserve(num_points);
        
        for (int i = 1; i <= num_points; ++i) {
            // Generate random point
            double x = dis(gen);
            double y = dis(gen);
            
            // Check if point is inside quarter circle
            if (std::sqrt(x*x + y*y) <= 1.0) {
                points_inside++;
            }
            
            // Calculate current estimate
            double pi_estimate = 4.0 * points_inside / i;
            estimates.push_back(pi_estimate);
        }
        
        return {estimates.back(), estimates};
    }
    
    void print_results(int num_points, double pi_estimate) {
        const double PI = 3.14159265358979323846;
        double error = std::abs(pi_estimate - PI) / PI * 100;
        
        std::cout << std::fixed << std::setprecision(6);
        std::cout << "\nNumber of points: " << num_points << std::endl;
        std::cout << "π estimate: " << pi_estimate << std::endl;
        std::cout << "Actual π: " << PI << std::endl;
        std::cout << "Error: " << error << "%" << std::endl;
    }
};

int main() {
    // Create estimator with fixed seed for reproducibility
    MonteCarloPi mc(42);
    
    // Run estimation with different numbers of points
    for (int points : {1000, 10000, 100000}) {
        auto [pi_estimate, estimates] = mc.estimate_pi(points);
        mc.print_results(points, pi_estimate);
    }
    
    return 0;
}


using System;
using System.Collections.Generic;
using System.Linq;

public class MonteCarloPi
{
    private readonly Random rng;
    
    public MonteCarloPi(int? seed = null)
    {
        rng = seed.HasValue ? new Random(seed.Value) : new Random();
    }
    
    public (double estimate, List<double> estimates) EstimatePi(int numPoints)
    {
        int pointsInside = 0;
        var estimates = new List<double>(numPoints);
        
        for (int i = 1; i <= numPoints; i++)
        {
            // Generate random point
            double x = rng.NextDouble();
            double y = rng.NextDouble();
            
            // Check if point is inside quarter circle
            if (Math.Sqrt(x*x + y*y) <= 1.0)
            {
                pointsInside++;
            }
            
            // Calculate current estimate
            double piEstimate = 4.0 * pointsInside / i;
            estimates.Add(piEstimate);
        }
        
        return (estimates.Last(), estimates);
    }
    
    public void PrintResults(int numPoints, double piEstimate)
    {
        const double PI = Math.PI;
        double error = Math.Abs(piEstimate - PI) / PI * 100;
        
        Console.WriteLine($"\nNumber of points: {numPoints}");
        Console.WriteLine($"π estimate: {piEstimate:F6}");
        Console.WriteLine($"Actual π: {PI:F6}");
        Console.WriteLine($"Error: {error:F2}%");
    }
    
    public static void Main()
    {
        // Create estimator with fixed seed for reproducibility
        var mc = new MonteCarloPi(seed: 42);
        
        // Run estimation with different numbers of points
        foreach (int points in new[] {1000, 10000, 100000})
        {
            var (piEstimate, estimates) = mc.EstimatePi(points);
            mc.PrintResults(points, piEstimate);
            
            // Calculate statistics
            double mean = estimates.Average();
            double stdDev = Math.Sqrt(
                estimates.Sum(x => Math.Pow(x - mean, 2)) / estimates.Count
            );
            
            Console.WriteLine($"Standard Deviation: {stdDev:F6}");
        }
    }
}

Quick Info