Probability Distributions

A probability distribution describes how likely each possible value of a random quantity is. "Marks are roughly bell-shaped around 60", "each delivery has a 5% chance of being late", "a shop gets about 3 orders a day": each sentence is a distribution in disguise. Once you know which distribution fits, you can answer questions such as "How likely is a day with no orders?" without waiting a year to find out.

In this lesson you will use scipy.stats to work with four distributions every data scientist meets, and you will check whether a distribution fits real data from the practice files.

Random variables, discrete and continuous

A random variable is a quantity whose value depends on chance, such as tomorrow's number of orders.

• A discrete variable takes separate values you can count: 0, 1, 2 orders. Its distribution gives the probability of each value, the probability mass function (PMF).
• A continuous variable can take any value in a range: a temperature of 31.27 °C. The chance of any exact value is zero, so its distribution is described by a probability density function (PDF), and probabilities are areas under the curve.

Every distribution in scipy.stats has the same set of methods:

The normal distribution

The normal distribution is the bell curve. It is defined by its mean (loc) and standard deviation (scale). About 68% of values fall within 1 SD of the mean, 95% within 2 SD and 99.7% within 3 SD. Many measurements, such as heights and exam marks for a large group, are roughly normal.

from scipy import stats

marks = stats.norm(loc=62, scale=14)       # a "frozen" distribution: mean 62, SD 14

print("P(mark < 33):  ", round(marks.cdf(33), 4))
print("P(mark > 90):  ", round(marks.sf(90), 4))
print("P(50 to 75):   ", round(marks.cdf(75) - marks.cdf(50), 4))
print("Top 10% cut-off:", round(marks.ppf(0.90), 1))
print("Within 2 SD:   ", round(marks.cdf(62 + 28) - marks.cdf(62 - 28), 4))

P(mark < 33):   0.0192
P(mark > 90):   0.0228
P(50 to 75):    0.6278
Top 10% cut-off: 79.9
Within 2 SD:    0.9545

stats.norm(loc=62, scale=14) creates a frozen distribution with fixed parameters, so you do not repeat them in every call. If marks followed this curve, about 2% of students would fail with less than 33, and a student needs about 80 to be in the top 10%. ppf() answers "what value?" questions; cdf() and sf() answer "what probability?" questions.

Does real data look normal?

A model is only useful if it fits. Compare the English marks in students.csv with a normal curve using the same mean and SD.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats

students = pd.read_csv("students.csv").drop_duplicates()
english = students["english"]
model = stats.norm(loc=english.mean(), scale=english.std())

fig, ax = plt.subplots(figsize=(7, 4))
ax.hist(english, bins=range(15, 101, 5), density=True, color="#cbd5e1", edgecolor="white", label="Students")
x = np.linspace(15, 100, 200)
ax.plot(x, model.pdf(x), color="#c2410c", linewidth=2, label="Normal model")
ax.set_title("English marks vs a normal curve")
ax.set_xlabel("Marks")
ax.legend()
plt.show()

print("Share within 1 SD, data: ", round(english.between(model.mean() - model.std(), model.mean() + model.std()).mean(), 3))
print("Share within 1 SD, model:", round(model.cdf(model.mean() + model.std()) - model.cdf(model.mean() - model.std()), 3))

Share within 1 SD, data:  0.715
Share within 1 SD, model: 0.683

The grey histogram follows the orange bell curve closely, rising to a peak near 60 marks and falling away symmetrically. density=True scales the histogram so its total area is 1, like the curve's. The numbers agree well too: about 72% of students fall within one standard deviation of the mean, close to the 68% the normal model predicts (the real marks are slightly more bunched in the middle). So for English marks, the normal distribution is a reasonable model.

The binomial distribution

The binomial distribution counts successes in a fixed number of independent tries, each with the same probability. Its parameters are n (number of tries) and p (chance of success each time). Typical questions: how many of 20 deliveries will be late, how many of 10 penalty kicks will score?

from scipy import stats

late = stats.binom(n=20, p=0.05)       # 20 deliveries, each 5% likely to be late

print("Expected late deliveries:", late.mean())
print("P(none late):     ", round(late.pmf(0), 4))
print("P(exactly 2 late):", round(late.pmf(2), 4))
print("P(3 or more late):", round(late.sf(2), 4))
for k in range(5):
    print(k, "#" * round(late.pmf(k) * 100))

Expected late deliveries: 1.0
P(none late):      0.3585
P(exactly 2 late): 0.1887
P(3 or more late): 0.0755
0 ####################################
1 ######################################
2 ###################
3 ######
4 #

On an average day, 1 of 20 deliveries is late, and there is about a 36% chance that none is. "3 or more" is written as sf(2), the probability of a value greater than 2, because the binomial is discrete. The little text bar chart shows the shape of the PMF: most likely 0 or 1, rarely more than 3.

The Poisson distribution

The Poisson distribution counts events in a fixed period when they happen independently at a steady average rate mu: customers per hour, calls per minute, orders per day. Let us test it on the practice shop. How many orders arrive per day, and does a Poisson model match?

import pandas as pd
from scipy import stats

sales = pd.read_csv("sales.csv", parse_dates=["order_date"]).drop_duplicates()
per_day = sales.groupby("order_date").size().reindex(pd.date_range("2025-01-01", "2025-12-31"), fill_value=0)
quiet = per_day[~per_day.index.month.isin([10, 11])]      # leave out the festive months

mu = quiet.mean()
model = stats.poisson(mu)
print("Average orders per day (Jan-Sep, Dec):", round(mu, 2))

table = pd.DataFrame({
    "observed": quiet.value_counts(normalize=True).sort_index(),
    "poisson": [model.pmf(k) for k in quiet.value_counts().sort_index().index],
}).round(3)
print(table.head(8))
print("Days with 7+ orders, observed:", round((quiet >= 7).mean(), 3), " model:", round(model.sf(6), 3))

Average orders per day (Jan-Sep, Dec): 3.0
   observed  poisson
0     0.059    0.050
1     0.168    0.149
2     0.168    0.224
3     0.253    0.224
4     0.161    0.168
5     0.105    0.101
6     0.049    0.051
7     0.023    0.022
Days with 7+ orders, observed: 0.036  model: 0.034

reindex() adds the days with no orders as 0, which groupby().size() alone would leave out, a common mistake when counting events. The observed share of days with 0, 1, 2 orders and so on is close to the Poisson prediction, so a Poisson model with about 3 orders a day describes normal days well. The festive months were left out because their rate is much higher; one Poisson model cannot describe two different rates.

The uniform distribution and sampling

In a uniform distribution every value in a range is equally likely, like the arrival minute of a bus that comes "sometime in the next 15 minutes". Every distribution can also generate random samples with rvs(), which is how simulations are built.

from scipy import stats

wait = stats.uniform(loc=0, scale=15)      # any time from 0 to 15 minutes
print("Mean wait:", wait.mean(), " P(wait > 10):", round(wait.sf(10), 3))

samples = stats.norm(loc=62, scale=14).rvs(size=5, random_state=42)
print(samples.round(1))

orders = stats.poisson(3.1).rvs(size=7, random_state=1)
print("A simulated week of orders:", orders)

Mean wait: 7.5  P(wait > 10): 0.333
[69.  60.1 71.1 83.3 58.7]
A simulated week of orders: [2 1 1 3 3 3 4]

Note SciPy's naming: for uniform, loc is the start and scale is the width of the range, not its end. random_state makes samples repeatable, like a NumPy seed.

Common mistakes

• Assuming everything is normal. Order values, incomes and rainfall are right-skewed. Check with a histogram before using a normal model.
• Off-by-one with discrete cdf and sf. For "3 or more", use sf(2) or 1 - cdf(2), not sf(3).
• Reading a PDF value as a probability. For continuous distributions, probabilities come from cdf() differences, not pdf() values.
• Ignoring changing rates. A Poisson model needs a steady rate; split periods such as festive and normal months.

Key points

• A distribution describes how likely each value of a random variable is; discrete ones have a PMF, continuous ones a PDF.
• scipy.stats distributions share pmf/pdf, cdf, sf, ppf, mean, std and rvs.
• Normal (norm): bell curve set by mean and SD; about 68%, 95% and 99.7% within 1, 2 and 3 SD.
• Binomial (binom): successes in n tries; Poisson (poisson): counts of events at a steady rate; uniform: equal chances in a range.
• Always compare a model with real data before trusting it.

from scipy import stats

calls = stats.poisson(4.5)
print("P(0 calls):", round(calls.pmf(0), 4))
print("P(more than 8):", round(calls.sf(8), 4))
print("95% of minutes have at most", int(calls.ppf(0.95)), "calls")
print(calls.rvs(size=10, random_state=7))

P(0 calls): 0.0111
P(more than 8): 0.0403
95% of minutes have at most 8 calls
[5 2 4 5 5 6 6 6 6 1]