Python for loops are widely considered interpretable, but they are also remarkably slow when running statistical simulations on data collections containing thousands of elements.

This article shows three examples of using NumPy array operations — known as vectorized operations — to run simulations like Monte Carlo or Bootstrap models on large datasets nearly instantly.

We will use real datasets from the gakudo-ai/open-datasets repository to illustrate the three strategies and tricks.

1. Bootstrapping Means

Suppose we take the penguins dataset, known for its rather small size, pick its body_mass_g attribute, and want to calculate a 95% confidence interval using 10,000 samples obtained through bootstrapping.

Instead of using a for loop, the vectorized solution is as shown below. We first generate a large 2D array of random indices associated with penguin instances in the base dataset, and compute means across the axis: all in no time, virtually!

import pandas as pd
import numpy as np

# Loading Penguins dataset directly
url = "https://raw.githubusercontent.com/gakudo-ai/open-datasets/main/penguins.csv"
masses = pd.read_csv(url)['body_mass_g'].dropna().values

# Generating a 2D array of 10,000 samples in one go
idx = np.random.choice(len(masses), size=(10000, len(masses)), replace=True)
boot_means = masses[idx].mean(axis=1)

What matters most in the previous example is not printing an output, but realizing how nearly instantaneous the execution time is.

2. Monte Carlo Shocks

Here’s a larger, dataset-driven simulation. Consider the California housing dataset, where the target attribute, median_house_value, contains average house prices for each of the over 20K districts in the state of California. Our goal is to simulate 5,000 different market scenarios, all of which have housing prices randomly fluctuating by up to ±10%.

NumPy broadcasting can also be easily exploited in this scenario. The trick boils down to broadcasting the 1D array of original prices against a 2D array of randomly generated price modifiers ranging between 0.9 and 1.1. This approach to Monte Carlo simulation uses zero loops and may take up to 2 seconds to run, including the vectorized product applied to the resulting massive 2D matrix:

url = "https://raw.githubusercontent.com/gakudo-ai/open-datasets/main/housing.csv"
prices = pd.read_csv(url)['median_house_value'].dropna().values

# Broadcasting 1D prices against a massive, 5000xN matrix of random modifiers
mods = np.random.uniform(0.9, 1.1, size=(5000, len(prices)))
sim_markets = prices * mods

3. Probabilistic Boolean Evaluation

The last example uses the 50 Startups dataset, fetching one of its predictor attributes, R&D Spend, to run 1,000 simulations aimed at predicting every startup’s survival rate if their R&D investments dropped by 5–15%. For simplicity, we use a $50K threshold to define a condition that will be evaluated through boolean array operations on an entire, massive matrix of simulated values. Neither row-by-row nor if-else evaluations are needed to perform this experiment in virtually no time:

url = "https://raw.githubusercontent.com/gakudo-ai/open-datasets/main/50_Startups.csv"
spend = pd.read_csv(url)['R&D Spend'].values
threshold = 50000 

# Creating modifiers and evaluating logical conditions globally on 50 startups
reductions = np.random.uniform(0.85, 0.95, size=(1000, len(spend)))
new_spend = spend * reductions
survival_rates = (new_spend > threshold).mean(axis=1)

Note that the resulting survival_rates variable will contain a 1D array where each element represents the proportion of startups in a given simulation that achieved an R&D spend above the threshold after reductions — in other words, a “survival rate”.