NumPy’s random module is often recognized for generating arrays of random numbers, but its capabilities go far beyond basic random number generation. With NumPy, we can create simulations that allow us to model different scenarios, test hypotheses, and gain insights into various systems. In this article, we will use NumPy to create random numbers and build simulations, covering examples from estimating π with Monte Carlo methods to simulating ecosystem dynamics.
Example 1: Monte Carlo Estimation of π
One of the classic applications of random number generation is estimating π using Monte Carlo simulation. The idea is simple: if we randomly throw darts at a square with an inscribed quarter-circle, the ratio of darts landing inside the quarter-circle to total throws will approximate π/4.
import numpy as np
# Set seed for reproducibility
np.random.seed(42)
def estimate_pi(n_points):
# Generate random points
x = np.random.uniform(0, 1, n_points)
y = np.random.uniform(0, 1, n_points)
# Calculate distances from origin
distances = np.sqrt(x**2 + y**2)
# Count points inside quarter circle
points_inside = np.sum(distances <= 1)
# Estimate pi
pi_estimate = 4 * points_inside / n_points
return pi_estimate, x, y, distances
# Run simulation with increasing number of points
point_counts = [1000, 10000, 100000, 1000000]
results = []
for n in point_counts:
pi_est, _, _, _ = estimate_pi(n)
results.append(pi_est)
print(f"π estimate with {n:,} points: {pi_est:.6f}")
print(f"Error: {abs(pi_est - np.pi):.6f}")
print()
Here are the results of our π estimation:
π estimate with 1,000 points: 3.112000 Error: 0.029593 π estimate with 10,000 points: 3.154400 Error: 0.012807 π estimate with 100,000 points: 3.136400 Error: 0.005193 π estimate with 1,000,000 points: 3.141132 Error: 0.000461
Let’s break down what’s happening in this simulation:
1. Random Number Type: We use np.random.uniform() to generate uniformly distributed random numbers between 0 and 1. This gives us random (x,y) coordinates within a unit square.
2. Simulation Process:
– Each point represents a “dart throw” onto a 1×1 square
– We check if each point falls within a quarter circle of radius 1
– The ratio of points inside the quarter circle to total points, multiplied by 4, approximates π
3. Results Analysis:
– Notice how the error decreases roughly with the square root of the number of points
– With 1,000 points, we’re off by about 0.03 (1% error)
– With 1,000,000 points, we’re off by only 0.0005 (0.015% error)
– This demonstrates the “law of large numbers” in action
For a visual simulation of this example, please take a look at this video tutorial.
Example 2: Customer Service Queue Simulation
def simulate_service_desk(n_hours, arrival_rate, service_rate):
# Convert rates to same time unit (minutes)
arrival_rate = arrival_rate / 60 # customers per minute
service_rate = service_rate / 60 # customers served per minute
# Simulate arrival times (using exponential distribution)
n_minutes = n_hours * 60
arrival_times = np.random.exponential(1/arrival_rate, size=int(n_minutes*arrival_rate*1.5))
arrival_times = np.cumsum(arrival_times)
arrival_times = arrival_times[arrival_times <= n_minutes]
# Simulate service times
service_times = np.random.exponential(1/service_rate, size=len(arrival_times))
# Calculate waiting times and queue lengths
wait_times = []
queue_lengths = np.zeros(n_minutes)
current_time = 0
for i, (arrival, service) in enumerate(zip(arrival_times, service_times)):
# Update current time and calculate wait
current_time = max(current_time, arrival)
wait = current_time - arrival
wait_times.append(wait)
# Update queue lengths array
start_minute = int(arrival)
end_minute = int(min(n_minutes, current_time + service))
queue_lengths[start_minute:end_minute] += 1
# Move time forward
current_time += service
return wait_times, queue_lengths
# Run simulation
hours = 8
arrival_rate = 30 # customers per hour
service_rate = 40 # customers served per hour
wait_times, queue_lengths = simulate_service_desk(hours, arrival_rate, service_rate)
Our simulation results show:
Average wait time: 2.48 minutes Maximum wait time: 11.22 minutes Maximum queue length: 8 customers
Let’s analyze this simulation in detail:
1. Random Number Types Used:
– np.random.exponential() for both arrival and service times
– We use exponential distribution because it models the time between independent events well
– The arrival rate parameter (λ) determines the average time between customer arrivals
– The service rate parameter determines how quickly customers are served
2. Simulation Process:
– Models an 8-hour day with an average of 30 customers arriving per hour
– Service rate is 40 customers per hour (when someone is being served)
– Each customer’s arrival and service time is randomly generated
3. Results Analysis:
– The average wait time of 2.48 minutes is quite reasonable
– Maximum wait of 11.22 minutes indicates occasional busy periods
– Queue never exceeded 8 people, suggesting adequate staffing
– These metrics could help determine if another service desk is needed
Example 3: Stock Price Movement Simulation
def simulate_stock_prices(initial_price, days, mu, sigma, n_simulations):
# Parameters
dt = 1/252 # Trading days in a year
n_steps = days
# Generate random walks
random_walks = np.random.normal(
(mu - sigma**2/2)*dt,
sigma*np.sqrt(dt),
size=(n_simulations, n_steps)
)
# Calculate price paths
price_paths = initial_price * np.exp(np.cumsum(random_walks, axis=1))
# Add initial price at the beginning
price_paths = np.hstack([
initial_price * np.ones((n_simulations, 1)),
price_paths
])
return price_paths
# Simulation parameters
initial_price = 100 # Starting stock price
days = 252 # One trading year
mu = 0.10 # Expected annual return
sigma = 0.20 # Annual volatility
n_sims = 1000 # Number of simulations
# Run simulation
prices = simulate_stock_prices(initial_price, days, mu, sigma, n_sims)
Our simulation produced these results:
Expected final price: $111.55 95% confidence interval: $72.95 to $160.55
Let’s examine this simulation more closely:
1. Random Number Type:
– np.random.normal() generates normally distributed random numbers
– This models the random “shocks” that affect stock prices
– The normal distribution is used because stock returns are often assumed to be normally distributed (in log space)
2. Simulation Process:
– Uses Geometric Brownian Motion, the standard model for stock price movements
– Incorporates both drift (expected return) and volatility
– Runs 1,000 different possible price paths over one trading year
3. Results Analysis:
– Starting from $100, the average ending price is $111.55 (11.55% return)
– The wide confidence interval ($72.95 to $160.55) shows the significant uncertainty
– This matches our input of 10% expected return (μ) and 20% volatility (σ)
– Demonstrates why diversification is important in investing
Example 4: Ecosystem Population Dynamics
def simulate_ecosystem(n_days, initial_rabbits, initial_foxes, params):
rabbit_pop = np.zeros(n_days)
fox_pop = np.zeros(n_days)
# Set initial populations
rabbit_pop[0] = initial_rabbits
fox_pop[0] = initial_foxes
# Simulation parameters
rabbit_growth = params['rabbit_growth']
fox_death = params['fox_death']
predation = params['predation']
fox_efficiency = params['fox_efficiency']
# Run simulation
for day in range(1, n_days):
# Calculate population changes
rabbits = rabbit_pop[day-1]
foxes = fox_pop[day-1]
# Add randomness to parameters
random_factor = np.random.normal(1, 0.1)
# Update populations
new_rabbits = rabbits + (
rabbit_growth * rabbits * random_factor -
predation * rabbits * foxes
)
new_foxes = foxes + (
fox_efficiency * predation * rabbits * foxes -
fox_death * foxes
) * random_factor
# Ensure populations don't go negative
rabbit_pop[day] = max(0, new_rabbits)
fox_pop[day] = max(0, new_foxes)
return rabbit_pop, fox_pop
The simulation ended with:
Final rabbit population: 1129 Final fox population: 1
Let’s analyze this ecological simulation:
1. Random Number Type:
– np.random.normal() with mean 1 and standard deviation 0.1
– Models random environmental factors affecting population growth
– The normal distribution is appropriate for modeling small random variations around an expected value
2. Simulation Process:
– Based on the Lotka-Volterra predator-prey equations
– Adds random environmental factors to make it more realistic
– Runs for 365 days starting with 100 rabbits and 20 foxes
3. Results Analysis:
– The final populations (1129 rabbits, 1 fox) show an interesting outcome
– The fox population has nearly died out
– The rabbit population has grown significantly
– This suggests our parameters might need adjustment for a more stable ecosystem
– Demonstrates how random variations can lead to unexpected outcomes in complex systems
Conclusion
These examples demonstrate the ability of NumPy to create practical simulations. From financial modeling to ecosystem dynamics, the key is combining NumPy’s random number generation with domain-specific logic to create meaningful simulations.
Some key takeaways:
1. Always set a random seed for reproducibility
2. Consider using multiple simulations to understand variation
3. Validate your simulations against known results when possible
4. Use NumPy’s vectorized operations for better performance
Remember that while these examples are simplified, the same principles apply to more complex simulations. The combination of NumPy’s efficient array operations and random number generation makes it a suitable choice for simulation work.
For more advanced simulations, you might want to explore specialized libraries like SimPy for discrete event simulation or Mesa for agent-based modeling. However, NumPy provides an excellent foundation for understanding and implementing basic simulation concepts.
