Code #3402

# monte_carlo.py
# Monte Carlo Simulation for Pi Estimation
# A complete implementation with extensive documentation

import random
import math
from typing import Tuple, List, Optional


def monte_carlo_pi_estimation(
 num_samples: int = 1000000,
 seed: Optional[int] = None,
 return_points: bool = False,
 verbose: bool = False
) -> Tuple[float, float, Optional[List[Tuple[float, float, bool]]]]:
 """
 Estimate the value of pi using the Monte Carlo method.
 
 This function uses a probabilistic approach to estimate pi by randomly
 generating points within a unit square [0, 1] x [0, 1] and determining
 what fraction of those points fall within a quarter circle of radius 1
 centered at the origin.
 
 The mathematical basis:
 - Area of quarter circle = (pi * r^2) / 4 = pi / 4 (when r = 1)
 - Area of unit square = 1
 - Ratio of points inside circle to total points ≈ pi / 4
 - Therefore, pi ≈ 4 * (points inside circle / total points)
 
 Parameters
 ----------
 num_samples : int, optional
 The number of random points to generate for the simulation.
 Higher values produce more accurate estimates but require more
 computation time. Default is 1,000,000.
 
 Must be a positive integer greater than 0.
 
 Recommended values:
 - Quick estimate: 10,000 - 100,000
 - Good accuracy: 1,000,000
 - High precision: 10,000,000+
 
 seed : int or None, optional
 Random seed for reproducibility. If None, the random number
 generator is not seeded, producing different results each run.
 Default is None.
 
 Useful for:
 - Testing and debugging
 - Reproducible experiments
 - Comparing different sample sizes
 
 return_points : bool, optional
 If True, returns a list of all generated points along with
 whether each point fell inside the quarter circle.
 Default is False.
 
 Warning: Setting this to True with large num_samples values
 will consume significant memory.
 
 verbose : bool, optional
 If True, prints progress updates during the simulation.
 Updates are printed at 10%, 25%, 50%, 75%, and 100% completion.
 Default is False.
 
 Returns
 -------
 Tuple[float, float, Optional[List[Tuple[float, float, bool]]]]
 A tuple containing:
 
 pi_estimate : float
 The estimated value of pi based on the Monte Carlo simulation.
 
 error : float
 The absolute error between the estimate and the actual value
 of pi (math.pi). Calculated as |pi_estimate - math.pi|.
 
 points : List[Tuple[float, float, bool]] or None
 If return_points is True, a list of tuples where each tuple
 contains (x, y, is_inside) for each generated point.
 If return_points is False, this value is None.
 
 Raises
 ------
 ValueError
 If num_samples is not a positive integer.
 TypeError
 If num_samples is not an integer or seed is not an integer/None.
 
 Examples
 --------
 Basic usage with default parameters:
 
 >>> pi_est, error, _ = monte_carlo_pi_estimation()
 >>> print(f'Estimated pi: {pi_est:.6f}')
 Estimated pi: 3.141592 # Approximate, varies each run
 
 Using a seed for reproducible results:
 
 >>> pi_est, error, _ = monte_carlo_pi_estimation(num_samples=100000, seed=42)
 >>> print(f'Estimated pi: {pi_est:.6f}, Error: {error:.6f}')
 Estimated pi: 3.139880, Error: 0.001713
 
 Getting the generated points for visualization:
 
 >>> pi_est, error, points = monte_carlo_pi_estimation(
 ... num_samples=1000,
 ... seed=123,
 ... return_points=True
 ... )
 >>> inside_count = sum(1 for _, _, inside in points if inside)
 >>> print(f'Points inside circle: {inside_count} / 1000')
 Points inside circle: 780 / 1000
 
 Using verbose mode to track progress:
 
 >>> pi_est, error, _ = monte_carlo_pi_estimation(
 ... num_samples=1000000,
 ... verbose=True
 ... )
 Progress: 10% complete (100000 / 1000000 samples)
 Progress: 25% complete (250000 / 1000000 samples)
 Progress: 50% complete (500000 / 1000000 samples)
 Progress: 75% complete (750000 / 1000000 samples)
 Progress: 100% complete (1000000 / 1000000 samples)
 
 Comparing accuracy with different sample sizes:
 
 >>> for n in [1000, 10000, 100000, 1000000]:
 ... pi_est, error, _ = monte_carlo_pi_estimation(n, seed=42)
 ... print(f'n={n:>7}: pi≈{pi_est:.6f}, error={error:.6f}')
 n= 1000: pi≈3.164000, error=0.022407
 n= 10000: pi≈3.153600, error=0.012007
 n= 100000: pi≈3.139880, error=0.001713
 n=1000000: pi≈3.142628, error=0.001035
 
 Notes
 -----
 The Monte Carlo method for estimating pi works as follows:
 
 1. Generate random points (x, y) where both x and y are uniformly
 distributed in the range [0, 1].
 
 2. For each point, check if it falls inside the quarter circle
 by testing if x^2 + y^2 <= 1.
 
 3. The ratio of points inside the circle to total points
 approximates the ratio of the quarter circle area to the
 square area, which equals pi/4.
 
 4. Multiply this ratio by 4 to get the estimate of pi.
 
 The accuracy of the estimate improves with the square root of the
 number of samples (Law of Large Numbers). To halve the error,
 you need approximately 4 times as many samples.
 
 Time Complexity: O(n) where n is num_samples
 Space Complexity: O(1) if return_points is False, O(n) otherwise
 
 References
 ----------
 .. [1] Metropolis, N., & Ulam, S. (1949). 'The Monte Carlo Method'.
 Journal of the American Statistical Association, 44(247), 335-341.
 
 .. [2] Robert, C., & Casella, G. (2004). 'Monte Carlo Statistical Methods'.
 Springer.
 
 See Also
 --------
 random.random : Generates random floats in [0.0, 1.0)
 math.pi : The mathematical constant pi
 
 """
 # =========================================================================
 # Input Validation
 # =========================================================================
 
 # Validate num_samples type
 if not isinstance(num_samples, int):
 raise TypeError(
 f'num_samples must be an integer, got {type(num_samples).__name__}'
 )
 
 # Validate num_samples value
 if num_samples <= 0:
 raise ValueError(
 f'num_samples must be a positive integer, got {num_samples}'
 )
 
 # Validate seed type
 if seed is not None and not isinstance(seed, int):
 raise TypeError(
 f'seed must be an integer or None, got {type(seed).__name__}'
 )
 
 # =========================================================================
 # Initialize Random Number Generator
 # =========================================================================
 
 # Set the random seed if provided for reproducibility
 if seed is not None:
 random.seed(seed)
 
 # =========================================================================
 # Monte Carlo Simulation
 # =========================================================================
 
 # Counter for points that fall inside the quarter circle
 points_inside_circle: int = 0
 
 # Optional list to store all generated points
 points_list: Optional[List[Tuple[float, float, bool]]] = []
 if not return_points:
 points_list = None
 
 # Define progress milestones for verbose output
 progress_milestones: List[float] = [0.10, 0.25, 0.50, 0.75, 1.00]
 next_milestone_index: int = 0
 
 # Main simulation loop
 for i in range(num_samples):
 # Generate random x and y coordinates in the unit square [0, 1]
 x: float = random.random()
 y: float = random.random()
 
 # Calculate the squared distance from the origin
 # We use squared distance to avoid the costly sqrt operation
 # Point is inside quarter circle if x^2 + y^2 <= r^2 (r=1)
 distance_squared: float = x * x + y * y
 
 # Check if the point is inside the quarter circle
 is_inside: bool = distance_squared <= 1.0
 
 # Increment counter if point is inside
 if is_inside:
 points_inside_circle += 1
 
 # Store point if requested
 if return_points and points_list is not None:
 points_list.append((x, y, is_inside))
 
 # Print progress if verbose mode is enabled
 if verbose and next_milestone_index < len(progress_milestones):
 current_progress: float = (i + 1) / num_samples
 if current_progress >= progress_milestones[next_milestone_index]:
 percentage: int = int(progress_milestones[next_milestone_index] * 100)
 print(
 f'Progress: {percentage}% complete '
 f'({i + 1} / {num_samples} samples)'
 )
 next_milestone_index += 1
 
 # =========================================================================
 # Calculate Pi Estimate
 # =========================================================================
 
 # The ratio of points inside the quarter circle to total points
 # approximates pi/4, so we multiply by 4 to get pi
 ratio: float = points_inside_circle / num_samples
 pi_estimate: float = 4.0 * ratio
 
 # Calculate the absolute error from the true value of pi
 absolute_error: float = abs(pi_estimate - math.pi)
 
 # =========================================================================
 # Return Results
 # =========================================================================
 
 return pi_estimate, absolute_error, points_list


# =============================================================================
# Main Execution Block - Demonstration
# =============================================================================

if __name__ == '__main__':
 """
 Demonstration of the Monte Carlo pi estimation function.
 
 This block runs when the script is executed directly and showcases
 various features of the monte_carlo_pi_estimation function.
 """
 
 print('=' * 70)
 print('MONTE CARLO PI ESTIMATION - DEMONSTRATION')
 print('=' * 70)
 print()
 
 # -------------------------------------------------------------------------
 # Demo 1: Basic Usage
 # -------------------------------------------------------------------------
 print('Demo 1: Basic Usage (1 million samples)')
 print('-' * 40)
 
 pi_est, error, _ = monte_carlo_pi_estimation(num_samples=1000000, seed=42)
 
 print(f'Estimated value of pi: {pi_est:.10f}')
 print(f'Actual value of pi: {math.pi:.10f}')
 print(f'Absolute error: {error:.10f}')
 print(f'Relative error: {(error / math.pi) * 100:.6f}%')
 print()
 
 # -------------------------------------------------------------------------
 # Demo 2: Accuracy vs Sample Size
 # -------------------------------------------------------------------------
 print('Demo 2: Accuracy vs Sample Size')
 print('-' * 40)
 print(f'{"Samples":>12} | {"Pi Estimate":>14} | {"Error":>12} | {"Rel. Error":>10}')
 print('-' * 55)
 
 sample_sizes = [100, 1000, 10000, 100000, 1000000, 10000000]
 
 for n in sample_sizes:
 pi_est, error, _ = monte_carlo_pi_estimation(num_samples=n, seed=42)
 rel_error = (error / math.pi) * 100
 print(f'{n:>12,} | {pi_est:>14.10f} | {error:>12.10f} | {rel_error:>9.6f}%')
 
 print()
 
 # -------------------------------------------------------------------------
 # Demo 3: Reproducibility with Seeds
 # -------------------------------------------------------------------------
 print('Demo 3: Reproducibility with Seeds')
 print('-' * 40)
 
 print('Running 3 times with seed=123:')
 for run in range(1, 4):
 pi_est, _, _ = monte_carlo_pi_estimation(num_samples=100000, seed=123)
 print(f' Run {run}: pi ≈ {pi_est:.10f}')
 
 print()
 print('Running 3 times without seed (different each time):')
 for run in range(1, 4):
 pi_est, _, _ = monte_carlo_pi_estimation(num_samples=100000, seed=None)
 print(f' Run {run}: pi ≈ {pi_est:.10f}')
 
 print()
 
 # -------------------------------------------------------------------------
 # Demo 4: Verbose Mode
 # -------------------------------------------------------------------------
 print('Demo 4: Verbose Mode')
 print('-' * 40)
 
 pi_est, error, _ = monte_carlo_pi_estimation(
 num_samples=500000,
 seed=42,
 verbose=True
 )
 print(f'Final estimate: pi ≈ {pi_est:.10f}')
 print()
 
 # -------------------------------------------------------------------------
 # Demo 5: Getting Points for Analysis
 # -------------------------------------------------------------------------
 print('Demo 5: Getting Points for Analysis')
 print('-' * 40)
 
 pi_est, error, points = monte_carlo_pi_estimation(
 num_samples=10000,
 seed=42,
 return_points=True
 )
 
 if points is not None:
 inside_points = [(x, y) for x, y, inside in points if inside]
 outside_points = [(x, y) for x, y, inside in points if not inside]
 
 print(f'Total points generated: {len(points)}')
 print(f'Points inside quarter circle: {len(inside_points)}')
 print(f'Points outside quarter circle: {len(outside_points)}')
 print(f'Ratio inside: {len(inside_points) / len(points):.6f}')
 print(f'Expected ratio (pi/4): {math.pi / 4:.6f}')
 
 print()
 print('=' * 70)
 print('DEMONSTRATION COMPLETE')
 print('=' * 70)