Flip a coin and both outcomes are equally likely. Measure 10,000 heights and most bunch near the middle. Count website visitors per hour and you get bursts around a steady average. Python probability distributions put names, formulas, and SciPy code on those patterns, and this post works through the five you will actually meet: normal, binomial, Poisson, uniform, and exponential.

“Probability is the language of uncertainty.”

Dennis Lindley, statistician

Last Updated: July 2026 | Tested on: Python 3.14.6, SciPy 1.18.0 | Difficulty: Advanced | Reading Time: 19 minutes

Here is the quick mental picture before any math. A distribution is just a recipe that says “outcomes near here are common, outcomes out there are rare.” Think of a dartboard. A beginner who aims for the center will land most darts near the bullseye and a few way out near the rim. Plot where every dart landed and that scatter, dense in the middle and thin at the edges, is a distribution. The rest of this post is about naming the common shapes and reading them with Python.

This post is the bridge between the intuition in the math for data science tutorial and the hypothesis testing tutorial that comes next. Every distribution here shows up again in machine learning: the normal distribution sits under linear regression, the binomial under logistic regression, and the Poisson under event prediction.

Vinay, a 26-year-old data scientist, was building a fraud detection model and assumed transaction amounts were normally distributed. They were not. Money tends to follow a log-normal shape: lots of small payments and a long tail of big ones. Once he log-transformed the amounts so the data looked roughly normal, his model accuracy jumped from 71% to 89%. The distribution mattered more than the algorithm. We will reproduce exactly that fix later in this post with real numbers.

Read the diagram top down. You start with one question, “what kind of data do I have?”, and it splits two ways. Measurements that can be any value (heights, delays, waiting times) point you at the continuous family: normal, uniform, and exponential. Things you count in whole numbers (clicks, defects, tickets) point you at the discrete family: binomial and Poisson. Each box names the shape, its parameters, and a real example, then the bottom row tells you which scipy.stats object to reach for. The rest of this post walks through each of these Python probability distributions with tested code.

Prerequisites

Finish the descriptive statistics tutorial first, since you need to be comfortable with mean and standard deviation. Basic algebra is enough math. Install SciPy with pip install scipy (this pulls in NumPy too). The code here was tested on Python 3.14.6 with SciPy 1.18.0 and NumPy 2.4.6.

Normal Distribution, the Bell Curve

The normal distribution is the most important one in all of statistics, and you already know its shape: the bell curve. Heights, test scores, measurement errors, daily stock returns, they all land roughly normal. The reason is that anything built from many small independent nudges tends to pile up in the middle. It has just two knobs. The mean (the Greek letter μ, mu) sets where the peak sits, and the standard deviation (σ, sigma) sets how wide the bell spreads. The famous 68-95-99.7 rule falls straight out of those two numbers: about 68% of values land within 1σ of the mean, 95% within 2σ, and 99.7% within 3σ.

Think of adult heights: the average sits around 170 cm and the typical wobble around that is 10 cm. Most people cluster near 170 cm. Someone at 200 cm is rare, and 230 cm is almost unheard of. That is a normal distribution doing its thing, and below we model it directly.

📄 normal_dist.py: the bell curve in Python

import numpy as np
from scipy import stats

# Create a normal distribution: mean=170cm, std=10cm (heights)
heights = stats.norm(loc=170, scale=10)

# Probability density at specific points
print(f"PDF at mean (170): {heights.pdf(170):.4f}")
print(f"PDF at 180:        {heights.pdf(180):.4f}")
print(f"PDF at 150:        {heights.pdf(150):.4f}")

# Cumulative probability: P(X <= value)
print(f"\nP(height <= 180): {heights.cdf(180):.4f} ({heights.cdf(180):.1%})")
print(f"P(height <= 170): {heights.cdf(170):.4f} ({heights.cdf(170):.1%})")
print(f"P(height <= 160): {heights.cdf(160):.4f} ({heights.cdf(160):.1%})")

# P(160 < height < 180)
p_range = heights.cdf(180) - heights.cdf(160)
print(f"P(160 < height < 180): {p_range:.4f} ({p_range:.1%})")

# Inverse: what height is at the 95th percentile?
p95 = heights.ppf(0.95)
print(f"\n95th percentile height: {p95:.1f} cm")

# The 68-95-99.7 rule
for n_sigma in [1, 2, 3]:
    p = heights.cdf(170 + n_sigma * 10) - heights.cdf(170 - n_sigma * 10)
    print(f"Within {n_sigma} std: {p:.4f} ({p:.1%})")

# Generate random samples
rng = np.random.default_rng(42)
samples = heights.rvs(size=10000, random_state=42)
print(f"\n10K samples: mean={samples.mean():.1f}, std={samples.std():.1f}")

▶ Output

PDF at mean (170): 0.0399
PDF at 180:        0.0242
PDF at 150:        0.0054

P(height <= 180): 0.8413 (84.1%)
P(height <= 170): 0.5000 (50.0%)
P(height <= 160): 0.1587 (15.9%)
P(160 < height < 180): 0.6827 (68.3%)

95th percentile height: 186.4 cm
Within 1 std: 0.6827 (68.3%)
Within 2 std: 0.9545 (95.4%)
Within 3 std: 0.9973 (99.7%)

10K samples: mean=170.0, std=10.0

What happened here: Four scipy methods do all the work. pdf() gives the height of the curve at a point (it is a density, not a probability, more on that in the FAQ). cdf() answers “what fraction sits at or below this value,” so cdf(170) returns exactly 0.5 because 170 is the mean and the bell is symmetric. Subtract two cdf values and you get the probability of a range, which is why cdf(180) - cdf(160) lands on 0.6827, the 68% from the rule. ppf() is the reverse of cdf(): hand it 0.95 and it tells you the 95th percentile height is 186.4 cm.

Finally, the 68-95-99.7 loop confirms the rule numerically, and the 10,000 random samples come back with a mean of 170.0 and a std of 10.0, matching the distribution we asked for.

The Normal PDF From Scratch

Before you trust a library function, it helps to see what it computes. Peeking at this formula is like lifting the hood of a car you drive every day: you do not need it to get to work, but once you have seen the engine you trust the machine a lot more. The normal probability density has one formula, and it is less scary than it looks. For a value x, with mean μ and standard deviation σ:

📄 normal_scratch.py: the formula by hand, then checked against scipy

import math

# The normal PDF formula, written out by hand (no scipy)
def normal_pdf(x, mu, sigma):
    coefficient = 1 / (sigma * math.sqrt(2 * math.pi))
    exponent = -((x - mu) ** 2) / (2 * sigma ** 2)
    return coefficient * math.exp(exponent)

# Heights: mean 170 cm, std 10 cm
mu, sigma = 170, 10

print(f"PDF at 170 (the peak): {normal_pdf(170, mu, sigma):.4f}")
print(f"PDF at 180 (one std up): {normal_pdf(180, mu, sigma):.4f}")
print(f"PDF at 150 (two std down): {normal_pdf(150, mu, sigma):.4f}")

# Compare with scipy to prove our hand version matches
from scipy import stats
print(f"\nscipy PDF at 170: {stats.norm(170, 10).pdf(170):.4f}")
print(f"scipy PDF at 180: {stats.norm(170, 10).pdf(180):.4f}")

▶ Output

PDF at 170 (the peak): 0.0399
PDF at 180 (one std up): 0.0242
PDF at 150 (two std down): 0.0054

scipy PDF at 170: 0.0399
scipy PDF at 180: 0.0242

What happened here: Our ten-line normal_pdf produces 0.0399 at the peak and 0.0242 one standard deviation up, the exact same numbers SciPy returns. So stats.norm.pdf() is not magic, it is this one formula evaluated for you, with all the edge cases and array support handled. Once you have seen the formula match by hand, you can lean on the library with full confidence and never type the math again.

Binomial Distribution, Success and Failure Counts

Switch from measuring to counting. The binomial distribution answers one question: if you run the same yes-or-no trial n times, and each trial succeeds with probability p, how many successes do you get? Coin flips are the textbook case (n flips, p = 0.5 for heads), but the real action is in click rates, defect rates, and conversion rates. Picture a basketball player who sinks 15% of long shots. Take 100 attempts and the binomial tells you how likely it is to make exactly 20, or 10, or 30. Each shot is independent, and only two things can happen: in or out.

📄 binomial_dist.py: counting successes in fixed trials

from scipy import stats

# Email campaign: 100 emails sent, 15% click rate
# What's the probability of getting exactly 20 clicks?
email = stats.binom(n=100, p=0.15)

print(f"P(exactly 20 clicks): {email.pmf(20):.4f} ({email.pmf(20):.1%})")
print(f"P(15 or fewer clicks): {email.cdf(15):.4f} ({email.cdf(15):.1%})")
print(f"P(more than 20 clicks): {1 - email.cdf(20):.4f}")
print(f"Expected clicks: {email.mean():.0f}")
print(f"Std dev: {email.std():.1f}")

# Quality control: 1000 products, 2% defect rate
# P(more than 30 defects)?
defects = stats.binom(n=1000, p=0.02)
p_over_30 = 1 - defects.cdf(30)
print(f"\nP(> 30 defects in 1000 items): {p_over_30:.4f} ({p_over_30:.1%})")
print(f"Expected defects: {defects.mean():.0f}")

▶ Output

P(exactly 20 clicks): 0.0402 (4.0%)
P(15 or fewer clicks): 0.5683 (56.8%)
P(more than 20 clicks): 0.0663
Expected clicks: 15
Std dev: 3.6

P(> 30 defects in 1000 items): 0.0126 (1.3%)
Expected defects: 20

What happened here: Notice the new method. For counts you use pmf() (probability mass function), not pdf(). Since you can only get a whole number of clicks, each count has a real probability you can read off directly: exactly 20 clicks happens about 4.0% of the time. cdf() still means “this many or fewer,” so 15 or fewer clicks comes in at 56.8%, and “more than 20” is just 1 - cdf(20).

The expected value lands at 15 (which is simply n times p, 100 × 0.15), and the spread is 3.6 clicks. The quality-control example shows the same idea at scale: with a 2% defect rate on 1000 items you expect 20 defects, and the chance of seeing more than 30 is a small but real 1.3%.

Poisson Distribution, Events Per Time Period

The Poisson distribution is the binomial's cousin for rare events spread over time or space. You do not have a fixed number of trials, you just have an average rate: 5 support tickets per hour, 2 server crashes per month, 3 typos per page. Given that average (the Greek letter λ, lambda), Poisson tells you the probability of seeing any exact count in one period. The everyday version: a busy coffee shop gets about 5 walk-ins every 10 minutes. Some 10-minute stretches see 2, some see 8, but they bunch around 5. Poisson is the math behind that bunching.

📄 poisson_dist.py: counting rare events

from scipy import stats

# Website gets average 5 support tickets per hour
tickets = stats.poisson(mu=5)

print(f"P(exactly 3 tickets): {tickets.pmf(3):.4f}")
print(f"P(0 tickets, quiet hour): {tickets.pmf(0):.4f}")
print(f"P(10+ tickets, busy hour): {1 - tickets.cdf(9):.4f}")
print(f"Expected: {tickets.mean()}, Std dev: {tickets.std():.2f}")

# Server crashes: average 2 per month
crashes = stats.poisson(mu=2)
p_zero = crashes.pmf(0)
print(f"\nP(0 crashes this month): {p_zero:.4f} ({p_zero:.1%})")
print(f"P(5+ crashes): {1 - crashes.cdf(4):.4f}")

▶ Output

P(exactly 3 tickets): 0.1404
P(0 tickets, quiet hour): 0.0067
P(10+ tickets, busy hour): 0.0318
Expected: 5.0, Std dev: 2.24

P(0 crashes this month): 0.1353 (13.5%)
P(5+ crashes): 0.0527

What happened here: Poisson takes a single parameter, mu, which is the average rate λ. With 5 tickets per hour as the average, the chance of an exact slow hour with 0 tickets is tiny (0.67%), while a flood of 10 or more happens about 3.2% of the time. One neat fact pops out of the numbers: for a Poisson distribution the variance equals the mean, so the std is the square root of the mean.

Here mean is 5.0 and std is 2.24, which is √5. The crash example reads the same way: with an average of 2 crashes per month, a clean zero-crash month has a 13.5% chance, and a bad month of 5 or more sits around 5.3%.

Uniform and Exponential Distributions

Two more Python probability distributions round out the continuous toolkit, and they are close partners. The uniform distribution is the flat one: every value in a range is equally likely, like a fair spinner. The exponential distribution models the waiting time between rare random events, and it is the natural companion to Poisson. If crashes arrive at a Poisson rate, the gaps between them follow an exponential curve: short gaps are common, long quiet stretches are rare but possible.

A homely way to feel the difference: the uniform distribution is a bus that is equally likely to be 0 to 10 minutes late, no minute favored. The exponential distribution is waiting for the next text message, where most gaps are short but every so often nobody texts you for hours.

📄 uniform_exponential.py: flat odds and waiting times

from scipy import stats

# Uniform: a random delay between 2 and 5 seconds before a retry
# scipy.stats.uniform takes loc (start) and scale (width), so here loc=2, scale=3
retry_delay = stats.uniform(loc=2, scale=3)

print("Uniform distribution (retry delay, 2 to 5 seconds)")
print(f"P(delay <= 3 seconds): {retry_delay.cdf(3):.4f} ({retry_delay.cdf(3):.1%})")
print(f"P(3 < delay < 4):      {retry_delay.cdf(4) - retry_delay.cdf(3):.4f}")
print(f"Mean delay: {retry_delay.mean():.1f} s, Std dev: {retry_delay.std():.2f} s")

# Exponential: time between server crashes, average one crash every 30 days
# The rate is 1 crash per 30 days, so the average gap (scale) is 30 days
gap = stats.expon(scale=30)

print("\nExponential distribution (days between crashes, avg gap 30 days)")
print(f"P(next crash within 10 days): {gap.cdf(10):.4f} ({gap.cdf(10):.1%})")
print(f"P(no crash for 60 days):      {1 - gap.cdf(60):.4f} ({1 - gap.cdf(60):.1%})")
print(f"Mean gap: {gap.mean():.0f} days, Median gap: {gap.median():.1f} days")

▶ Output

Uniform distribution (retry delay, 2 to 5 seconds)
P(delay <= 3 seconds): 0.3333 (33.3%)
P(3 < delay < 4):      0.3333
Mean delay: 3.5 s, Std dev: 0.87 s

Exponential distribution (days between crashes, avg gap 30 days)
P(next crash within 10 days): 0.2835 (28.3%)
P(no crash for 60 days):      0.1353 (13.5%)
Mean gap: 30 days, Median gap: 20.8 days

What happened here: The uniform setup is the one that trips people up. SciPy's uniform takes a start (loc) and a width (scale), not a start and end, so a 2-to-5 second range is loc=2, scale=3. Because every value is equally likely, the odds of landing in the first third (2 to 3 seconds) are exactly 33.3%, and any 1-second slice has the same chance. For the exponential gap, the average wait is 30 days, yet the median is only 20.8 days.

That gap between mean and median is the signature of a right-skewed curve: a few very long quiet spells drag the average up above the typical wait. Notice too that the chance of a 60-day quiet stretch (13.5%) is exactly e^-2, since 60 days is two average gaps.

The Catch: When Your Data Is Not Normal

Here is the mistake that bit Vinay, and it bites almost everyone once. Picture the average net worth of everyone in a small cafe: it stays modest until one billionaire walks in and yanks the average way up, even though almost every person in the room is still perfectly ordinary. That long thin tail is exactly what skewed data looks like. It is tempting to assume your data is normal and move on.

Real-world data often is not, and money is the classic offender. Transaction amounts, salaries, and house prices are usually right-skewed: a big pile of small values and a long thin tail of large ones. Run a tool that expects a bell curve on skewed data and your results quietly go wrong. The fix is often a single log transform, which squashes that long tail back into a bell shape. Let us reproduce Vinay's exact fix and prove it with a normality test.

📄 not_normal.py: detect skew, then fix it with a log transform

import numpy as np
from scipy import stats

# Vinay's fraud-detection story: transaction amounts are NOT normal.
# Real money data is usually skewed right (lots of small amounts, a few huge ones).
rng = np.random.default_rng(7)
amounts = rng.lognormal(mean=3.0, sigma=1.0, size=5000)  # log-normal, skewed

# Normality test: Shapiro-Wilk. If p > 0.05, the data looks normal.
stat_raw, p_raw = stats.shapiro(amounts[:500])  # shapiro caps the sample size
print(f"Raw amounts:  Shapiro p-value = {p_raw:.6f}  -> normal? {p_raw > 0.05}")
print(f"Raw skewness: {stats.skew(amounts):.2f} (0 means symmetric)")

# Log-transform, then test again
log_amounts = np.log(amounts)
stat_log, p_log = stats.shapiro(log_amounts[:500])
print(f"\nLogged amounts: Shapiro p-value = {p_log:.6f}  -> normal? {p_log > 0.05}")
print(f"Logged skewness: {stats.skew(log_amounts):.2f}")

▶ Output

Raw amounts:  Shapiro p-value = 0.000000  -> normal? False
Raw skewness: 6.79 (0 means symmetric)

Logged amounts: Shapiro p-value = 0.767848  -> normal? True
Logged skewness: 0.03

What happened here: The raw transaction amounts fail the Shapiro-Wilk normality test hard. The p-value rounds to 0.000000 (it is actually about 1e-28), which is far below 0.05, so we reject the idea that the data is normal. The skewness of 6.79 confirms a heavy right tail (a perfectly symmetric bell would read 0). After taking the natural log of every amount, the picture flips: the p-value jumps to 0.77, comfortably above 0.05, and the skewness drops to 0.03, basically symmetric. Same data, one transform, and now a model that assumes normality has a fair chance. That single line, np.log(amounts), is the whole reason Vinay's accuracy climbed from 71% to 89%.

Choosing Between Python Probability Distributions

When you are staring at a new dataset, this table is the 30-second decision guide. Ask whether you are measuring something continuous or counting whole numbers, then match the rest of the row.

One closing tip from experience: do not guess. Plot a histogram of your data first, then run a quick Shapiro-Wilk test like we did above. The shape of the data tells you the distribution, not the other way around. Picking the wrong one is the kind of silent bug that makes a model look fine in testing and fall apart in production.

Conclusion

You now have the five workhorse Python probability distributions and the SciPy calls to drive them. The normal bell curve models anything built from many small nudges, the binomial counts successes in a fixed number of yes-or-no trials, and Poisson counts rare events over time. Uniform gives flat odds across a range and exponential models the waiting time between those rare events.

Along the way you saw the four methods every one of them shares, pdf or pmf, cdf, ppf, and rvs, plus the single most valuable habit: check whether your data is actually normal before you assume it, and reach for a log transform when it is not. That one check is what turned Vinay's 71% model into an 89% one.

Next up is the hypothesis testing tutorial, where these distributions become the engine behind t-tests and chi-square tests that tell you whether a result is real or just noise. For the full learning path from basics to machine learning, head back to the Python + AI/ML tutorial series home.

Frequently Asked Questions

Interview Questions on Probability Distributions

If you can walk through these without peeking, you are ready for this topic in an interview.

Q: Why do you use pmf() for a binomial or Poisson distribution but pdf() for a normal one?

Binomial and Poisson are discrete, so each whole-number outcome has a real probability you can read off directly with pmf() (probability mass function). The normal is continuous, so the chance of landing on any exact value is effectively zero; pdf() returns a density, not a probability. To get an actual probability from a continuous distribution you integrate over a range, which in SciPy means subtracting two cdf() values.

Q: What is the difference between cdf() and ppf()?

They are inverses of each other. cdf(x) takes a value and returns the cumulative probability of landing at or below it, so cdf(180) might return 0.84. ppf(q) takes a probability and returns the value at that percentile, so ppf(0.84) returns roughly 180. You use cdf() to answer "what fraction is below this?" and ppf() to answer "what value marks this percentile?", which is how you compute things like the 95th percentile cutoff.

Q: For a Poisson distribution, what is the relationship between the mean and the variance?

For Poisson they are equal: the variance equals the mean, both equal to lambda. That means the standard deviation is the square root of lambda. So with an average of 5 events per period, the mean is 5, the variance is 5, and the standard deviation is the square root of 5, about 2.24. If real count data shows variance far larger than its mean, that is overdispersion and a sign Poisson is the wrong fit; a negative binomial is often better.

Q: A colleague sets up a uniform delay between 2 and 5 seconds with stats.uniform(2, 5) and the numbers look wrong. What is the bug?

SciPy's uniform takes loc (the start) and scale (the width), not the start and end. stats.uniform(2, 5) actually produces a range from 2 to 2 plus 5, which is 2 to 7, not 2 to 5. The correct call for a 2-to-5 second range is stats.uniform(loc=2, scale=3). This loc-and-scale convention trips up almost everyone the first time and is worth double-checking whenever a continuous distribution gives surprising results.

Q: Your fraud model trains fine but performs poorly in production, and you notice the transaction-amount feature is heavily right-skewed. What do you check and try first?

First confirm the skew with stats.skew() and a Shapiro-Wilk test; a large positive skewness and a tiny p-value tell you the feature is far from normal. Many models and distance metrics implicitly assume roughly symmetric, bell-shaped features, so a long tail of large amounts distorts them. The first fix to try is a log transform with np.log() (or np.log1p() if zeros are present), which compresses the tail and often pulls the feature close to normal. Re-run the normality check after transforming to confirm it worked before retraining.

Q: When is it safe to approximate a binomial distribution with a Poisson one, and why would you bother?

When the number of trials n is large and the success probability p is small, the binomial is well approximated by a Poisson with lambda equal to n times p. A common rule of thumb is n above 20 and p below 0.05. It is worth doing because Poisson needs only one parameter instead of two and is cheaper to compute for rare events like defects, crashes, or fraud hits, where n can be huge and p tiny.

Q: Why does the normal distribution show up so often in real-world data and machine learning?

Because of the Central Limit Theorem: the average of many independent samples from almost any distribution tends toward normal as the sample grows, regardless of the original shape. Anything built from many small independent effects, like measurement errors or combined noise sources, piles up in a bell curve. This is why so many statistical tests assume normality of the mean and why the normal sits under models like linear regression.

Further reading: the official Python documentation is the authoritative source on this.

Related Posts

Previous: Statistics: Descriptive (Mean, Median, Mode, Std Dev)

Next: Statistics: Hypothesis Testing with t-test and chi-square

Series Home: Python + AI/ML Tutorial Series