๐Ÿ“Monte Carlo Math โ‰…

๐Ÿ“Monte Carlo Math โ‰…#

Just math with no plotting (much faster!)#

Rendering graphics is a very slow process. If we just want to estimate ฯ€, we can do so much faster with Monte Carlo integration. Here is an example using a single core:

import random
import math
import time
import numpy as np

def estimate_pi(num_points):
    points_inside_circle = 0
    for _ in range(num_points):
        x = random.uniform(-1, 1)
        y = random.uniform(-1, 1)
        if x**2 + y**2 <= 1:
            points_inside_circle += 1
    return 4 * points_inside_circle / num_points

sample_sizes_sequential = np.logspace(1, 8, num=8, base=10.0, dtype=int)
time_results = {}
accuracy_results = {}
for num_samples in sample_sizes_sequential:
    start_time = time.time()
    estimated_val = estimate_pi(num_samples)
    end_time = time.time()
    elapsed_time = end_time - start_time
    time_results[num_samples] = elapsed_time
    accuracy_results[num_samples] = abs(math.pi - estimated_val)
    print(f"Estimated value with {num_samples} samples: {estimated_val} (Time taken: {elapsed_time:.4f} seconds)")

Estimated value with 10 samples: 2.0 (Time taken: 0.0000 seconds)
Estimated value with 100 samples: 3.24 (Time taken: 0.0001 seconds)
Estimated value with 1000 samples: 3.184 (Time taken: 0.0004 seconds)
Estimated value with 10000 samples: 3.1308 (Time taken: 0.0040 seconds)
Estimated value with 100000 samples: 3.13476 (Time taken: 0.0402 seconds)

Estimated value with 1000000 samples: 3.142468 (Time taken: 0.4037 seconds)

Estimated value with 10000000 samples: 3.1410452 (Time taken: 3.7168 seconds)

Estimated value with 100000000 samples: 3.14159168 (Time taken: 39.7457 seconds)

Explanation#

  • Area of square:

    • \(2 * 2 = 4\)

  • Area of circle:

    • \(ฯ€ * 1 * 1 = ฯ€\)

  • Every point generated in the random.uniform(-1, 1) by random.uniform(-1, 1) grid has equal probability of being anywhere inside of 4 unit area square

  • Probability of falling inside the circle, then, must be:

    • \(P = \frac{\text{area of circle}}{\text{area of square}} = \frac{\pi}{4}\)

  • Since the likelihood of points falling inside or outside the circle is effectively a direct function of their respective areas, the following is also true:

    • \(P = \frac{\text{number points inside circle}}{\text{number of total points}}\) = \(\frac{\pi}{4}\)

    • Therefore, if we want to estimate ฯ€, we multiply P by 4

Parallelizing the Monte Carlo ฯ€ estimation#

We can parallelize the Monte Carlo ฯ€ estimation by splitting the total number of points into multiple chunks and estimating ฯ€ for each chunk in parallel. This significantly reduces the time it takes to estimate ฯ€.

Here is an example using 128 processes:

import random
import math
import time
from multiprocessing import Pool
import numpy as np

def estimate_pi(num_points):
    points_inside_circle = 0
    for _ in range(num_points):
        x = random.uniform(-1, 1)
        y = random.uniform(-1, 1)
        if x**2 + y**2 <= 1:
            points_inside_circle += 1
    return points_inside_circle

def parallel_estimate_pi(num_points, num_processes):
    pool = Pool(processes=num_processes)
    points_per_process = [num_points // num_processes] * num_processes
    results = pool.map(estimate_pi, points_per_process)
    pool.close()
    pool.join()
    total_points_inside_circle = sum(results)
    return 4 * total_points_inside_circle / num_points

sample_sizes_parallel = np.logspace(1, 9, num=9, base=10.0, dtype=int)
num_processes = 128
parallel_time_results = {}
parallel_accuracy_results = {}

for num_samples in sample_sizes_parallel:
    start_time = time.time()
    estimated_val = parallel_estimate_pi(num_samples, num_processes)
    end_time = time.time()
    elapsed_time = end_time - start_time
    parallel_time_results[num_samples] = elapsed_time
    parallel_accuracy_results[num_samples] = abs(math.pi - estimated_val)
    print(f"Estimated value with {num_samples} samples: {estimated_val} (Time taken: {elapsed_time:.4f} seconds)")

Estimated value with 10 samples: 0.0 (Time taken: 0.3635 seconds)

Estimated value with 100 samples: 0.0 (Time taken: 0.3284 seconds)

Estimated value with 1000 samples: 2.856 (Time taken: 0.3214 seconds)

Estimated value with 10000 samples: 3.1444 (Time taken: 0.3743 seconds)

Estimated value with 100000 samples: 3.14068 (Time taken: 0.3624 seconds)

Estimated value with 1000000 samples: 3.1448 (Time taken: 0.4140 seconds)

Estimated value with 10000000 samples: 3.1426908 (Time taken: 0.5384 seconds)

Estimated value with 100000000 samples: 3.1412786 (Time taken: 1.8293 seconds)

Estimated value with 1000000000 samples: 3.141500728 (Time taken: 14.1222 seconds)

Numba#

Numba is a just-in-time compiler for Python that can significantly speed up numerical computations. Here is an example using Numba:

from numba import njit
import numpy as np
import time

@njit(fastmath=True)
def numba_estimate_pi(num_points):
    points_inside_circle = 0
    for _ in range(num_points):
        x = np.random.uniform(-1, 1)
        y = np.random.uniform(-1, 1)
        if x**2 + y**2 <= 1:
            points_inside_circle += 1
    return 4 * points_inside_circle / num_points

# Sample sizes for testing
sample_sizes_numba = np.logspace(1, 9, num=9, base=10.0, dtype=int)
numba_time_results = {}
numba_accuracy_results = {}

for num_samples in sample_sizes_numba:
    start_time = time.time()
    estimated_val = numba_estimate_pi(num_samples)
    end_time = time.time()
    elapsed_time = end_time - start_time
    
    # Store results using Python int as keys
    numba_time_results[int(num_samples)] = elapsed_time
    numba_accuracy_results[int(num_samples)] = abs(math.pi - estimated_val)
    print(f"Estimated value with {num_samples} samples: {estimated_val} (Time taken: {elapsed_time:.4f} seconds)")

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[3], line 1
----> 1 from numba import njit
      2 import numpy as np
      3 import time

ModuleNotFoundError: No module named 'numba'

Analysis of Results#

We can see that Numba is significantly faster than the sequential and parallel approaches.

import matplotlib.pyplot as plt

# Convert all sample sizes to Python int
sample_sizes_parallel = [int(s) for s in sample_sizes_parallel]
sample_sizes_sequential = [int(s) for s in sample_sizes_sequential]  # Assuming same sizes for sequential
sample_sizes_numba = [int(s) for s in sample_sizes_numba]

# Plotting the time results
plt.figure(figsize=(14, 6))

# Time taken plot (linear scale)
plt.subplot(1, 2, 1)
plt.plot(sample_sizes_parallel, [parallel_time_results[int(s)] for s in sample_sizes_parallel], marker='o', label='Parallel Time')
plt.plot(sample_sizes_sequential, [time_results[int(s)] for s in sample_sizes_sequential], marker='o', label='Sequential Time')
plt.plot(sample_sizes_numba, [numba_time_results[int(s)] for s in sample_sizes_numba], marker='o', label='Numba Time')
plt.xlabel('Number of Samples')
plt.ylabel('Time (seconds)')
plt.title('Time taken to estimate ฯ€')
plt.legend()
plt.grid(True, linestyle='--')

# Accuracy plot (linear scale)
plt.subplot(1, 2, 2)
plt.yscale('log')
plt.plot(sample_sizes_parallel, [parallel_accuracy_results[int(s)] for s in sample_sizes_parallel], marker='o', label='Parallel Accuracy')
plt.plot(sample_sizes_sequential, [accuracy_results[int(s)] for s in sample_sizes_sequential], marker='o', label='Sequential Accuracy')
plt.plot(sample_sizes_numba, [numba_accuracy_results[int(s)] for s in sample_sizes_numba], marker='o', label='Numba Accuracy')
plt.xlabel('Number of Samples')
plt.ylabel('Error')
plt.title('Accuracy of ฯ€ estimation')
plt.legend()
plt.grid(True, linestyle='--')

plt.tight_layout()
plt.show()

Examples of Monte Carlo Simulations#

Monte Carlo simulations are used in a variety of fields to solve problems that might be deterministic in principle but are too complex to solve directly.

  • Finance:

    • Model the probability things like the pricing of options and other derivatives, risk management, and portfolio optimization.

  • Physics:

    • Simulate the behavior of particles, such as in the study of statistical mechanics, quantum mechanics, and particle physics.

  • Engineering:

    • Assess the reliability and performance of systems and components, such as in the analysis of structural reliability, network performance, and manufacturing processes.

  • Medicine:

    • Model the spread of diseases, the effectiveness of treatments, and the planning of radiation therapy for cancer patients.

  • Environmental Science:

    • Model complex environmental systems, such as climate change predictions, the dispersion of pollutants, and the assessment of ecological risks.

  • Computer Graphics:

    • Simulate the interaction of light with surfaces, which helps in rendering realistic images and animations.

  • Game Development:

    • Create realistic scenarios and behaviors in games, such as simulating the roll of dice, the shuffle of cards, or the movement of characters.

Too many more to list.