Monte Carlo Simulation for estimating π

The Circle in the Square

Consider a square with side length 2, centered at the origin, with a circle of radius 1 inscribed within it.

Key points:

• Square area: 4 square units
• Circle area: π square units

The ratio of these areas is π/4.

The Monte Carlo Approach

Imagine randomly throwing darts at this square. The proportion of darts landing inside the circle should approximate the ratio of the circle’s area to the square’s area: (Darts inside circle) / (Total darts) ≈ π/4 Multiply both sides by 4, and we have our π estimate.

We use random number generators to create (x, y) coordinates within the square, then check if each point falls within the circle using x^2 + y^2 ≤ 1.

def generate_data(NUM_SAMPLES: int, min_: int, max_: int):
    min_, max_ = -1, 1
    circle = 0

    Xs = []
    Ys = []
    for i in range(NUM_SAMPLES):
        X = np.random.uniform(min_, max_)
        Y = np.random.uniform(min_, max_)

        Xs.append(X)
        Ys.append(Y)

        if X ** 2 + Y **2 < 1**2:
            circle += 1

    df = pd.DataFrame({"x": Xs, "y": Ys})
    df['in_circle'] = np.where(df['x']**2 + df['y']**2 < 1**2, 1, 0)
    rolling_est = 4 * np.cumsum(df['in_circle'][1:]) / np.arange(1, len(df))
    rolling_est = np.insert(rolling_est, 0, 0)
    df['rolling_est'] = rolling_est
    return df

df = generate_data(NUM_SAMPLES=2500, min_=-1, max_=1)

We add a column in_circle to indicate whether the point falls inside the circle or outside. Additionally, we create another column rolling_est which is a rolling estimate of π. This is what the dataset looks like.

df.head()

	x	y	in_circle	rolling_est
0	0.455318	0.499045	1	0.0
1	-0.326996	0.797670	1	4.0
2	-0.672441	-0.206805	1	4.0
3	0.630730	-0.050815	1	4.0
4	-0.890021	0.124752	1	4.0

Let’s run the simulation.

import plotly.express as px
import pandas as pd
import numpy as np

# Generate sample data
n_frames = int(len(df) / 10)
N = 10

# Create a DataFrame with cumulative data
df_list = []
for i in range(1, n_frames + 1):
    df_temp = pd.DataFrame({
        'x': df['x'][:i*N], 
        'y': df['y'][:i*N],
        'in_circle': df['in_circle'][:i*N],
        'rolling_est': df['rolling_est'][:i*N],
        'Iteration': i
    })
    df_list.append(df_temp)

df_with_iterations_for_viz = pd.concat(df_list, ignore_index=True)

# Create the animated scatter plot
fig = px.scatter(
    df_with_iterations_for_viz,
    x='x',
    y='y',
    animation_frame='Iteration',
    range_x=[-2, 2], range_y=[-1, 1]
    ) 

# Add a circle
fig = fig.add_shape(type="circle",
    xref="x", yref="y",
    x0=-1, y0=-1, x1=1, y1=1,
    line_color="red",
)

# Add a square
fig = fig.add_shape(type="rect",
    xref="x", yref="y",
    x0=-1, y0=-1,
    x1=1, y1=1,
    line_color="Black"
)

# Fix axes
fig = fig.update_yaxes(
    scaleanchor = "x",
    scaleratio = 1,
)

# Update layout for better appearance
fig = fig.update_layout(
    title='Cumulative Scatter Plot Animation',
    xaxis_title='X Axis',
    yaxis_title='Y Axis',
    showlegend=False
)

fig.show()

Rolling estimate of π converges over many iterations.

import pandas as pd
import plotly.express as px

# Assuming df is your DataFrame
# Calculate the rolling average
window_size = 20
rolling_avg = df['rolling_est'].rolling(window=window_size).mean()

# Create a new DataFrame to store the rolling averages
rolling_avg_df = df.copy()
rolling_avg_df['rolling_avg'] = rolling_avg

# Drop NaN values that were created by the rolling mean
rolling_avg_df = rolling_avg_df.dropna()

# Prepare the DataFrame for visualization
n_frames = int(len(rolling_avg_df) / 10)
N = 10

df_list = []
for i in range(1, n_frames + 1):
    df_temp = pd.DataFrame({
        'index': rolling_avg_df.index[:i*N],  # Use the index for x-axis
        'rolling_avg': rolling_avg_df['rolling_avg'][:i*N],
        'Iteration': i
    })
    df_list.append(df_temp)

# Combine all the frames into one DataFrame
rolling_avg_with_iterations = pd.concat(df_list, ignore_index=True)

# Create the animated line plot with fixed x-axis
fig = px.scatter(
    rolling_avg_with_iterations,
    x='index',
    y='rolling_avg',
    animation_frame='Iteration',
    range_x=[0, len(rolling_avg_df)],  # Fix the x-axis range to cover the entire index range
    range_y=[1, 3.5],
    title='Animated Rolling Average Estimation Time Series',
    labels={'index': 'DataFrame Index', 'rolling_avg': 'Rolling Average Estimation'},
)

fig = fig.add_shape(
        type='line',
        x0=0,
        y0=3.1415926535897932,
        x1=3000,
        y1=3.1415926535897932,
        line=dict(
            color='Red',
        )
)

# Update layout to improve appearance
fig = fig.update_traces(line=dict(width=2))
fig = fig.update_layout(
    updatemenus=[{
        'buttons': [
            {
                'args': [None, {'frame': {'duration': 50, 'redraw': True}, 'fromcurrent': True}],
                'label': 'Play',
                'method': 'animate'
            },
            {
                'args': [[None], {'frame': {'duration': 0, 'redraw': True}, 'mode': 'immediate', 'transition': {'duration': 0}}],
                'label': 'Pause',
                'method': 'animate'
            }
        ],
        'direction': 'left',
        'pad': {'r': 10, 't': 87},
        'showactive': False,
        'type': 'buttons',
        'x': 0.1,
        'xanchor': 'right',
        'yanchor': 'top'
    }]
)
fig.show()