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 0x11a7c2fd0>, <matplotlib.axis.XTick object at 0x11a9187d0>, <matplotlib.axis.XTick object at 0x11a918b90>, <matplotlib.axis.XTick object at 0x11a918f50>, <matplotlib.axis.XTick object at 0x11a919310>, <matplotlib.axis.XTick object at 0x11a9196d0>, <matplotlib.axis.XTick object at 0x11a919a90>, <matplotlib.axis.XTick object at 0x11a919e50>, <matplotlib.axis.XTick object at 0x11a91a210>, <matplotlib.axis.XTick object at 0x11a91a5d0>, <matplotlib.axis.XTick object at 0x11a91a990>, <matplotlib.axis.XTick object at 0x11a91ad50>, <matplotlib.axis.XTick object at 0x11a91b110>, <matplotlib.axis.XTick object at 0x11a91b4d0>, <matplotlib.axis.XTick object at 0x11a91b890>, <matplotlib.axis.XTick object at 0x11a91bc50>, <matplotlib.axis.XTick object at 0x11a954050>, <matplotlib.axis.XTick object at 0x11a954410>, <matplotlib.axis.XTick object at 0x11a9547d0>, <matplotlib.axis.XTick object at 0x11a954b90>], [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 .