Monte Carlo Simulation

Author

Robert W. Walker

Published

February 16, 2026

Monte Carlo Business Case

Gemini is a powerful tool for Monte Carlo simulation.

Prompt:

Can you generate a Monte Carlo simulation of Potential.Customers, distributed Poisson with arrival rate 7, Customers making a purchase with probability 0.6 from the set of Potential.Customers, and each Customer spending normal with mean 500 dollars and standard deviation of 100 dollars. Try 1000 iterations.

How’s that done?

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# Set seed for reproducibility
np.random.seed(42)

# Parameters
iterations = 1000
lambda_poisson = 7
p_purchase = 0.6
mean_spending = 500
std_spending = 100

results = []

for i in range(iterations):
    # Step 1: Sample Potential Customers
    potential_customers = np.random.poisson(lambda_poisson)
    
    # Step 2: Sample Actual Customers (Binomial)
    # Alternatively: check each potential customer with a 0.6 prob
    actual_customers = np.random.binomial(n=potential_customers, p=p_purchase)
    
    # Step 3: Sample Spending for each customer
    if actual_customers > 0:
        individual_spendings = np.random.normal(mean_spending, std_spending, actual_customers)
        total_spending = np.sum(individual_spendings)
    else:
        total_spending = 0.0
        
    results.append({
        'Iteration': i + 1,
        'Potential_Customers': potential_customers,
        'Actual_Customers': actual_customers,
        'Total_Spending': total_spending
    })

df = pd.DataFrame(results)

# Save to CSV
df.to_csv('monte_carlo_simulation_results.csv', index=False)

# Summary statistics
summary_stats = df.describe()
print(summary_stats)

         Iteration  Potential_Customers  Actual_Customers  Total_Spending
count  1000.000000          1000.000000       1000.000000     1000.000000
mean    500.500000             6.940000          4.192000     2091.394737
std     288.819436             2.492913          2.063345     1035.160073
min       1.000000             0.000000          0.000000        0.000000
25%     250.750000             5.000000          3.000000     1346.837689
50%     500.500000             7.000000          4.000000     2013.612058
75%     750.250000             9.000000          5.000000     2773.369027
max    1000.000000            16.000000         13.000000     6465.587864

How’s that done?

# Plotting the distribution of Total Spending
plt.hist(df['Total_Spending'], bins=30, color='skyblue', edgecolor='black')
plt.title('Distribution of Total Spending (1000 Iterations)')
plt.xlabel('Total Spending ($)')
plt.ylabel('Frequency')
plt.grid(axis='y', alpha=0.75)
plt.savefig('total_spending_distribution.png')

# Plotting Potential vs Actual Customers
plt.clf()
plt.figure(figsize=(10, 6))
plt.hist([df['Potential_Customers'], df['Actual_Customers']], label=['Potential Customers', 'Actual Customers'], bins=np.arange(0, 20) - 0.5, alpha=0.7)

(array([[  1.,   3.,  11.,  60.,  90., 132., 165., 151., 131., 105.,  61.,
         47.,  23.,  13.,   4.,   1.,   2.,   0.,   0.],
       [ 16.,  60., 142., 180., 189., 168., 113.,  70.,  35.,  15.,   8.,
          0.,   3.,   1.,   0.,   0.,   0.,   0.,   0.]]), array([-0.5,  0.5,  1.5,  2.5,  3.5,  4.5,  5.5,  6.5,  7.5,  8.5,  9.5,
       10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5]), <a list of 2 BarContainer objects>)

How’s that done?

plt.title('Distribution of Potential vs Actual Customers')
plt.xlabel('Number of Customers')
plt.ylabel('Frequency')
plt.legend()
plt.xticks(range(0, 20))

([<matplotlib.axis.XTick object at 0x11a0a2fd0>, <matplotlib.axis.XTick object at 0x11a1f87d0>, <matplotlib.axis.XTick object at 0x11a1f8b90>, <matplotlib.axis.XTick object at 0x11a1f8f50>, <matplotlib.axis.XTick object at 0x11a1f9310>, <matplotlib.axis.XTick object at 0x11a1f96d0>, <matplotlib.axis.XTick object at 0x11a1f9a90>, <matplotlib.axis.XTick object at 0x11a1f9e50>, <matplotlib.axis.XTick object at 0x11a1fa210>, <matplotlib.axis.XTick object at 0x11a1fa5d0>, <matplotlib.axis.XTick object at 0x11a1fa990>, <matplotlib.axis.XTick object at 0x11a1fad50>, <matplotlib.axis.XTick object at 0x11a1fb110>, <matplotlib.axis.XTick object at 0x11a1fb4d0>, <matplotlib.axis.XTick object at 0x11a1fb890>, <matplotlib.axis.XTick object at 0x11a1fbc50>, <matplotlib.axis.XTick object at 0x11a234050>, <matplotlib.axis.XTick object at 0x11a234410>, <matplotlib.axis.XTick object at 0x11a2347d0>, <matplotlib.axis.XTick object at 0x11a234b90>], [Text(0, 0, '0'), Text(1, 0, '1'), Text(2, 0, '2'), Text(3, 0, '3'), Text(4, 0, '4'), Text(5, 0, '5'), Text(6, 0, '6'), Text(7, 0, '7'), Text(8, 0, '8'), Text(9, 0, '9'), Text(10, 0, '10'), Text(11, 0, '11'), Text(12, 0, '12'), Text(13, 0, '13'), Text(14, 0, '14'), Text(15, 0, '15'), Text(16, 0, '16'), Text(17, 0, '17'), Text(18, 0, '18'), Text(19, 0, '19')])

How’s that done?

plt.savefig('customer_counts.png')

How’s that done?

my.df <- read.csv("monte_carlo_simulation_results.csv")
library(DT)
DT::datatable(my.df)

How’s that done?

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

# 1. Visualization of Linkage: Potential vs Actual
plt.figure(figsize=(8, 6))
# Adding jitter to show density better
sns.regplot(data=df, x='Potential_Customers', y='Actual_Customers', x_jitter=0.2, y_jitter=0.2, 
            scatter_kws={'alpha':0.3}, line_kws={'color':'red'})
plt.title('Linkage: Potential Customers vs Actual Customers')
plt.xlabel('Potential Customers (Poisson $\\lambda=7$)')
plt.ylabel('Actual Customers (Binomial $p=0.6$)')
plt.savefig('linkage_potential_vs_actual.png')

# 2. Visualization of Linkage: Actual Customers vs Total Spending
plt.clf()
plt.figure(figsize=(8, 6))
sns.regplot(data=df, x='Actual_Customers', y='Total_Spending', x_jitter=0.2, 
            scatter_kws={'alpha':0.3}, line_kws={'color':'green'})
plt.title('Linkage: Actual Customers vs Total Spending')
plt.xlabel('Actual Customers')
plt.ylabel('Total Spending ($)')
plt.savefig('linkage_actual_vs_spending.png')


# Correlations
correlations = df[['Potential_Customers', 'Actual_Customers', 'Total_Spending']].corr()
print("Correlation Matrix:")

Correlation Matrix:

How’s that done?

print(correlations)

                     Potential_Customers  Actual_Customers  Total_Spending
Potential_Customers             1.000000          0.758090        0.742648
Actual_Customers                0.758090          1.000000        0.979867
Total_Spending                  0.742648          0.979867        1.000000

Correlation Matrix:
                     Potential_Customers  Actual_Customers  Total_Spending
Potential_Customers             1.000000          0.758090        0.742648
Actual_Customers                0.758090          1.000000        0.979867
Total_Spending                  0.742648          0.979867        1.000000

I have refined the Monte Carlo simulation to explicitly demonstrate and visualize how the variables are interlinked.

Implementation Logic

Step 1 (Arrival): A random number of Potential Customers is drawn from a Poisson distribution ().
Step 2 (Conversion): The number of Actual Customers is then drawn from a Binomial distribution where the number of trials () is specifically the result from Step 1, and the success probability is .
Step 3 (Spending): For every actual customer identified in Step 2, a unique spending amount is drawn from a Normal distribution (). The Total Spending is the sum of these individual draws.

Visualizing the Linkages

Potential vs. Actual Customers: The first plot shows a direct linear relationship. Because the Binomial “n” is the Potential Customer count, you can see that the number of actual customers is constrained by the number of arrivals.
Actual Customers vs. Total Spending: This shows an extremely strong linear link. Since each additional customer adds an average of to the total, the volume of customers is the strongest predictor of total revenue.
Individual Customer Spending: To confirm that each customer is an individual draw, the third plot shows the distribution of every single purchase across all iterations. It forms a perfect normal bell curve centered at .