NumPy Broadcasting, Statistics and Simulation

The previous lesson showed that NumPy applies operations to whole arrays. This lesson goes one step deeper into what makes NumPy so expressive. Broadcasting is the rule that lets you combine arrays of different shapes: subtracting a column mean from every row of a table, scaling each store's numbers by its own factor, or adding a delivery fee to every cell of a grid. Understanding it removes the need for most nested loops.

We then build a toolkit for descriptive statistics: summaries along axes, percentiles, cumulative sums, rankings, correlations and histograms. We finish with a Monte Carlo simulation, a technique where you answer a hard question by running thousands of random experiments on the computer. This is the bridge from "describing what happened" to "estimating what might happen", and it uses nothing more than the random generator and array maths you already know.

Worked examples

Example 1: Broadcasting a scalar and a row

Question: Apply a store-specific multiplier to a table.

import numpy as np

orders = np.array([[310, 295, 330, 342],
                   [220, 245, 238, 260],
                   [405, 390, 421, 433]])

print(orders + 10)

week_factor = np.array([1.0, 1.0, 1.1, 1.1])
print(orders * week_factor)

store_factor = np.array([[0.9], [1.0], [1.2]])
print((orders * store_factor).round(0))
print(week_factor.shape, store_factor.shape, (orders * store_factor).shape)

[[320 305 340 352]
 [230 255 248 270]
 [415 400 431 443]]
[[310.  295.  363.  376.2]
 [220.  245.  261.8 286. ]
 [405.  390.  463.1 476.3]]
[[279. 266. 297. 308.]
 [220. 245. 238. 260.]
 [486. 468. 505. 520.]]
(4,) (3, 1) (3, 4)

How it works: When two arrays have different shapes, NumPy compares their shapes from the right. Two dimensions are compatible if they are equal or if one of them is 1; a size-1 dimension is "stretched" to match. A scalar stretches to everything. A row of shape (4,) is stretched down all three rows. A column of shape (3, 1) is stretched across all four columns. No data is copied; the stretching is virtual.

Reading the result: The row factor multiplied each column (week) by its own value, while the (3, 1) column multiplied each row (store). The final result still has shape (3, 4).

Example 2: Standardise columns

Question: Convert each column to z-scores.

rng = np.random.default_rng(0)
data = rng.normal(loc=[25, 500, 4], scale=[5, 120, 1], size=(200, 3))

col_mean = data.mean(axis=0)
col_std = data.std(axis=0)
z = (data - col_mean) / col_std

print(col_mean.round(1), col_std.round(1))
print(z.mean(axis=0).round(6), z.std(axis=0).round(6))
print(z.shape)

[ 24.7 504.8   4. ] [  5.1 119.8   1. ]
[-0. -0. -0.] [1. 1. 1.]
(200, 3)

How it works: data has shape (200, 3) and col_mean has shape (3,). Broadcasting subtracts each column's mean from that whole column, then divides by that column's standard deviation. After standardising, every column has mean 0 and standard deviation 1, which puts different units (minutes, rupees, ratings) on a comparable scale. This is a standard step before many machine-learning algorithms.

Reading the result: The original column means were [ 24.7 504.8 4. ], and after standardising they are all ≈ 0.0. The one-line formula replaced a double loop.

Example 3: Pairwise differences with np.newaxis

Question: Compute the distance between every pair of dark stores.

coords = np.array([[12.97, 77.59],   # Bengaluru
                   [19.07, 72.88],   # Mumbai
                   [28.61, 77.21],   # Delhi
                   [13.08, 80.27]])  # Chennai
names = ["Bengaluru", "Mumbai", "Delhi", "Chennai"]

diff = coords[:, np.newaxis, :] - coords[np.newaxis, :, :]
dist = np.sqrt((diff ** 2).sum(axis=2))

print(diff.shape, dist.shape)
print(dist.round(1))
i, j = np.unravel_index(np.argmax(dist), dist.shape)
print(names[i], names[j], dist[i, j].round(1))

(4, 4, 2) (4, 4)
[[ 0.   7.7 15.6  2.7]
 [ 7.7  0.  10.5  9.5]
 [15.6 10.5  0.  15.8]
 [ 2.7  9.5 15.8  0. ]]
Delhi Chennai 15.8

How it works: Inserting a new axis with np.newaxis (or None) turns shape (4, 2) into (4, 1, 2) and (1, 4, 2). Subtracting them broadcasts to (4, 4, 2): every pair of rows. Squaring, summing over the last axis and taking the root gives a 4×4 matrix of straight-line distances (in degrees here, which is only for illustration). This trick replaces a double loop over pairs.

Reading the result: The two stores farthest apart (in coordinate space) are Delhi and Chennai, at 15.8 degrees. The diagonal is all zeros because every store is zero distance from itself.

Example 4: Statistics along axes

Question: Summarise a 2-D table in different directions.

weekly = np.array([[310, 295, 330, 342],
                   [220, 245, 238, 260],
                   [405, 390, 421, 433]])

print(weekly.mean(axis=0), weekly.mean(axis=1))
print(weekly.std(axis=1).round(1))
print(np.median(weekly, axis=0))
print(weekly.argmax(axis=1), weekly.argmax())
print(weekly.max(axis=1) - weekly.min(axis=1))
print(np.percentile(weekly, 75, axis=1))
print(weekly.cumsum(axis=1))

[311.66666667 310.         329.66666667 345.        ] [319.25 240.75 412.25]
[18.1 14.4 16.2]
[310. 295. 330. 342.]
[3 3 3] 11
[47 40 43]
[333.   248.75 424.  ]
[[ 310  605  935 1277]
 [ 220  465  703  963]
 [ 405  795 1216 1649]]

How it works: Almost every reduction takes an axis. argmax returns the position of the maximum, per row when you give axis=1, or in the flattened array when you do not. cumsum runs a cumulative total along the axis. Percentiles accept the same argument.

Reading the result: The best week for each store is week index [3 3 3] (counting from 0), and the biggest swing between best and worst week is 47 orders.

Example 5: Sorting and ranking

Question: Rank products by revenue and extract the top three.

products = np.array(["Milk", "Eggs", "Curd", "Paneer", "Bread", "Butter"])
revenue = np.array([5400, 8200, 3100, 7600, 2900, 6100])

order = np.argsort(revenue)[::-1]
print(products[order])
print(revenue[order])

top3 = products[np.argsort(revenue)[-3:][::-1]]
print(top3)

ranks = np.empty_like(order)
ranks[order] = np.arange(1, len(order) + 1)
print({str(p): int(r) for p, r in zip(products, ranks)})

share = revenue / revenue.sum()
print(share.round(3), np.cumsum(share[order]).round(2))

['Eggs' 'Paneer' 'Butter' 'Milk' 'Curd' 'Bread']
[8200 7600 6100 5400 3100 2900]
['Eggs' 'Paneer' 'Butter']
{'Milk': 4, 'Eggs': 1, 'Curd': 5, 'Paneer': 2, 'Bread': 6, 'Butter': 3}
[0.162 0.246 0.093 0.228 0.087 0.183] [0.25 0.47 0.66 0.82 0.91 1.  ]

How it works: argsort returns the indices that would sort the array; indexing another array with them reorders it consistently (here, the product names). Reversing gives descending order. The rank trick assigns rank numbers back to each original position. Cumulative share of the sorted values is how a Pareto (80/20) analysis is built.

Reading the result: The top three products are ['Eggs', 'Paneer', 'Butter'], and together they hold 66% of revenue.

Example 6: Differences, growth and rolling windows

Question: Compute day-to-day change and a 3-day moving average.

daily = np.array([120, 135, 128, 150, 162, 158, 171, 190, 185, 200])

change = np.diff(daily)
pct = change / daily[:-1] * 100
print(change)
print(pct.round(1))

window = 3
kernel = np.ones(window) / window
moving = np.convolve(daily, kernel, mode="valid")
print(moving.round(1))

cum = daily.cumsum()
print(cum[-1], np.cumprod(1 + pct[:3] / 100).round(3))

[15 -7 22 12 -4 13 19 -5 15]
[12.5 -5.2 17.2  8.  -2.5  8.2 11.1 -2.6  8.1]
[127.7 137.7 146.7 156.7 163.7 173.  182.  191.7]
1599 [1.125 1.067 1.25 ]

How it works: np.diff gives differences between neighbours (one value shorter), and dividing by the previous value gives percentage growth. A moving average is a convolution with a kernel of equal weights; mode="valid" keeps only positions where the window fully fits, so the output is window - 1 values shorter. Pandas hides this behind .rolling(3).mean().

Reading the result: Sales rose from 120 to 200 over 10 days, an average daily change of 8.9. The smoothed series has 8 points, from 127.7 to 191.7.

Example 7: Correlation

Question: Do orders with more items take longer to deliver?

import csv

orders_rows = list(csv.DictReader(open("orders.csv", newline="", encoding="utf-8")))
items_rows = list(csv.DictReader(open("order_items.csv", newline="", encoding="utf-8")))

n_items = {}
for r in items_rows:
    n_items[r["order_id"]] = n_items.get(r["order_id"], 0) + int(r["quantity"])

x, y = [], []
for r in orders_rows:
    if r["delivery_minutes"] and r["order_id"] in n_items:
        x.append(n_items[r["order_id"]])
        y.append(float(r["delivery_minutes"]))

x, y = np.array(x), np.array(y)
r = np.corrcoef(x, y)
print(x.size)
print(r.round(3))
print(np.cov(x, y).round(2))

5630
[[1.    0.041]
 [0.041 1.   ]]
[[ 8.99  0.86]
 [ 0.86 48.02]]

How it works: np.corrcoef returns a correlation matrix; the off-diagonal value is the Pearson correlation between the two arrays, always between -1 and 1. np.cov gives the covariance matrix. A value near 0 means no linear relationship.

Reading the result: Using 5,630 orders, the correlation between basket size and delivery time is 0.041. Remember that correlation measures straight-line association, not causation, and a tiny value does not rule out non-linear patterns.

Example 8: Histograms and bins

Question: Count deliveries in time buckets.

counts, edges = np.histogram(y, bins=[0, 15, 20, 25, 30, 40, 100])
for lo, hi, c in zip(edges[:-1], edges[1:], counts):
    print(f"{lo:>3.0f}-{hi:<3.0f}{c:>6,}  {'#' * int(c / 60)}")

labels = np.array(["<15", "15-20", "20-25", "25-30", "30-40", "40+"])
bucket = np.digitize(y, [15, 20, 25, 30, 40])
print(labels[bucket][:6])
print(np.bincount(bucket))

  0-15    565  #########
 15-20    975  ################
 20-25  1,550  #########################
 25-30  1,417  #######################
 30-40  1,059  #################
 40-100    64  #
['30-40' '20-25' '30-40' '25-30' '20-25' '30-40']
[ 565  975 1550 1417 1059   64]

How it works: np.histogram counts values in each bin and returns the counts and the edges. np.digitize assigns each value the index of its bin, which we can use to label values or to count with bincount. The text bar chart is a quick way to see a distribution without plotting.

Reading the result: The most common delivery bucket is 20-25 minutes with 1,550 orders, and 1,123 deliveries (20%) took 30 minutes or more.

Example 9: Random sampling

Question: Draw random samples and shuffle reproducibly.

rng = np.random.default_rng(2026)

customers = np.array([f"C{i:04d}" for i in range(1, 601)])
sample = rng.choice(customers, size=5, replace=False)
print(sample)

winners = rng.choice(customers, size=3, replace=False, p=None)
print(len(set(winners)))

weights = np.array([0.5, 0.3, 0.15, 0.05])
draw = rng.choice(["UPI", "Card", "COD", "Wallet"], size=10000, p=weights)
values, counts = np.unique(draw, return_counts=True)
print({str(v): int(c) for v, c in zip(values, counts)})

perm = rng.permutation(10)
print(perm, np.sort(perm))

['C0384' 'C0220' 'C0107' 'C0508' 'C0016']
3
{'COD': 1514, 'Card': 3029, 'UPI': 4932, 'Wallet': 525}
[7 2 5 3 1 8 4 0 6 9] [0 1 2 3 4 5 6 7 8 9]

How it works: default_rng(seed) creates a generator whose results are repeatable. choice samples with or without replacement, and p supplies probabilities. permutation returns a shuffled copy. Drawing 10,000 values with known weights shows the law of large numbers: the observed proportions land near the intended ones.

Reading the result: The observed UPI share was 49.3% against the intended 50%. With more draws the match would be even closer.

Example 10: A Monte Carlo simulation

Question: If each delivery is on time with probability 0.9, how likely is it that a rider with 12 deliveries has at least 2 late ones?

rng = np.random.default_rng(7)
trials = 100_000

late = rng.random((trials, 12)) > 0.9
late_per_rider = late.sum(axis=1)

prob = (late_per_rider >= 2).mean()
print(late_per_rider[:10])
print(round(prob, 4))
print(np.bincount(late_per_rider)[:6] / trials)

from math import comb
exact = 1 - sum(comb(12, k) * 0.1**k * 0.9**(12 - k) for k in range(2))
print(round(exact, 4))

[0 2 1 0 1 1 3 0 1 2]
0.3399
[0.28523 0.37485 0.2299  0.08442 0.02145 0.00355]
0.341

How it works: We simulate 100,000 riders, each making 12 deliveries that are late with probability 0.1. A single call creates a 100,000 × 12 matrix of random numbers; comparing with 0.9 gives booleans, and summing along axis=1 counts the late deliveries per rider. The mean of "two or more late" estimates the probability. Because this particular problem has an exact binomial formula, we can check the answer.

Reading the result: The simulation estimates 0.3399, close to the exact 0.3410, so about 34% of such riders would have two or more late deliveries. Monte Carlo is valuable precisely when no neat formula exists.

Common mistakes

1. Mismatched shapes in broadcasting.

np.ones((3, 4)) + np.ones((3,))

Traceback (most recent call last):
  File "example.py", line 1, in <module>
    np.ones((3, 4)) + np.ones((3,))
    ~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~
ValueError: operands could not be broadcast together with shapes (3,4) (3,) 

2. Confusing axis=0 and axis=1. Remember: the named axis is the one that disappears.

3. Forgetting keepdims or newaxis when you need a column result, such as (x - x.mean(axis=1)) failing; use x.mean(axis=1, keepdims=True).

4. Using a different random seed per run and being unable to reproduce results.

5. Reading a correlation as causation.

Practice questions

Practice 1 (Easy): Centre the data

Subtract each column's mean from weekly and verify that the column means become zero.

centred = weekly - weekly.mean(axis=0)
print(centred.round(1))
print(centred.mean(axis=0).round(10))

[[ -1.7 -15.    0.3  -3. ]
 [-91.7 -65.  -91.7 -85. ]
 [ 93.3  80.   91.3  88. ]]
[-0.  0. -0.  0.]

Practice 2 (Medium): Row percentages

Convert weekly into the percentage each week contributes to its own store's total, so each row sums to 100.

pct = weekly / weekly.sum(axis=1, keepdims=True) * 100
print(pct.round(1))
print(pct.sum(axis=1))

[[24.3 23.1 25.8 26.8]
 [22.8 25.4 24.7 27. ]
 [24.6 23.7 25.5 26.3]]
[100. 100. 100.]

Practice 3 (Hard): Simulate the chance of a stock-out

A store holds 60 units of milk. Daily demand is normally distributed with mean 50 and standard deviation 12 (rounded to whole units, never negative). Estimate the probability that demand on a single day exceeds stock, using 200,000 simulated days.

rng = np.random.default_rng(11)
demand = np.maximum(0, rng.normal(50, 12, 200_000).round())

p_out = (demand > 60).mean()
print(round(p_out, 4))
print(np.percentile(demand, 95))

0.1904
70.0

Key takeaways

• Broadcasting aligns shapes from the right: dimensions must match or be 1.
• Use axis for reductions and keepdims=True or np.newaxis to keep shapes broadcastable.
• Standardising columns is one line: (x - x.mean(axis=0)) / x.std(axis=0).
• argsort, cumsum, diff, histogram and digitize cover ranking, growth and binning.
• np.corrcoef measures linear association, not causation.
• A seeded default_rng gives reproducible random draws and simulations.
• Monte Carlo simulation estimates probabilities by repeating random experiments.