Distributions you’ll actually meet: Gaussian, Bernoulli, and the long tail

Most people meet the Gaussian distribution in school, draw a bell curve, and quietly assume that real-world data behaves something like it.

Then they start working with real data.

Website response times occasionally spike to several seconds. A handful of words appear thousands of times while most appear only once. A yes-or-no event cannot possibly follow a continuous bell curve. Even some measurements from nature produce two distinct peaks instead of one.

The problem isn’t that the Gaussian distribution is wrong. It’s that different data-generating processes produce different distribution shapes, and choosing the wrong model can lead to seriously misleading predictions.

In this tutorial, we’ll build a practical gallery of probability distributions using Python, NumPy, SciPy, and Matplotlib. Rather than stopping at textbook formulas, we’ll fit distributions to real observations, visualize their behavior, and investigate what happens in the tails.

By the end, you’ll understand:

• Why some continuous measurements are approximately Gaussian.

• How the Bernoulli distribution models binary outcomes.

• Why word frequencies follow a radically different pattern from heights.

• How to model positive, right-skewed measurements such as HTTP latency.

• Why averages can hide the events that matter most.

• How to recognize when a dataset cannot reasonably be described by one bell curve.

We’ll finish with a particularly revealing real-world dataset: a famous natural phenomenon whose measurements form two distinct groups.

1. What is a probability distribution?

A probability distribution describes how likely different values of a random variable are.

The random variable could represent:

• The height of a randomly selected adult.

• Whether a user clicks a button.

• The number of orders received in one minute.

• The response time of an HTTP request.

• The number of times a particular word appears in a document.

The important distinction is between the values a variable can take and the probabilities assigned to those values.

Discrete versus continuous distributions

A discrete distribution assigns probabilities to individual outcomes.

For example, consider a button click:

• Click: 1

• No click: 0

We can meaningfully ask for the probability of exactly 1.

A continuous distribution instead describes probability density across a range of values.

Suppose an adult’s height is modeled by a continuous distribution. The probability of one mathematically exact height is zero, but the probability of falling between two heights can be positive.

This distinction produces two important functions.

Probability mass function (PMF)

Used for discrete random variables. It assigns a probability to each possible value.

Probability density function (PDF)

Used for continuous random variables. The area under the density curve over an interval gives the probability of falling within that interval.

A density value itself is not necessarily a probability and can exceed 1. The total area under a valid continuous density curve must equal 1.

Two additional functions are especially useful in practice:

• Cumulative distribution function (CDF): the probability of observing a value less than or equal to a threshold.

• Survival function (SF): the probability of observing a value strictly greater than a threshold.

For a threshold t:

P(X > t) = 1 − CDF(t)

In code, SciPy’s survival function is often preferable to subtracting the CDF from 1 because it can be more numerically accurate for extreme tail probabilities.

We’ll use these functions repeatedly.

2. Set up the Python environment

We’ll use four tools:

Library	Purpose
NumPy	Numerical arrays, sampling, and percentiles
SciPy	Probability distributions, parameter fitting, and statistical diagnostics
Matplotlib	Histograms, density curves, and diagnostic plots
Python standard library	Downloading datasets, parsing text, and measuring request durations

Use Python 3.12 or newer. As of October 2026, SciPy 1.18.1 is the current stable bug-fix release, published on August 21, 2026. The examples use documented APIs supported by this release family.

News+1

Create a project directory and virtual environment:

mkdir distribution-gallery
cd distribution-gallery
python -m venv .venv

Activate it on macOS or Linux:

source .venv/bin/activate

Install the dependencies:

python -m pip install --upgrade pip
python -m pip install numpy scipy matplotlib

We’ll create four Python scripts. Each is self-contained and can be run independently. The dataset examples need an internet connection for their initial downloads; the latency experiment also needs network access during measurement.

3. The Gaussian distribution: a useful starting point

The Gaussian distribution, also called the normal distribution, is one of the most widely used continuous probability models.

It has the familiar symmetric bell shape.

Two parameters define it:

• Mean (μ): the center of the distribution.

• Standard deviation (σ): the spread around that center.

A standard deviation must be positive. Increasing it makes the distribution wider and flatter. Changing the mean shifts the bell curve without changing its shape.

A Gaussian distribution is symmetric around its mean. Moving away from the center reduces probability density.

For an exactly Gaussian population:

• Approximately 68.3% of observations fall within one standard deviation of the mean.

• Approximately 95.4% fall within two standard deviations.

• Approximately 99.7% fall within three standard deviations.

These familiar percentages are properties of the mathematical model, not guarantees about every real dataset.

Why does this distribution appear so often?

One important reason is the central limit theorem.

Under suitable conditions, an average or suitably standardized sum of many random contributions approaches a normal distribution as the number of contributions increases.

For example, a measurement affected by many individually small influences may have an approximately Gaussian shape.

However, the theorem does not say every raw measurement is normally distributed. Nor does it guarantee that a sum of strongly dependent, extremely heavy-tailed, or heterogeneous variables will be well approximated by a Gaussian at practical sample sizes.

That distinction matters when fitting distributions to real observations.

4. Gallery #1: fitting a Gaussian to real height measurements

Let’s use a historical dataset collected by Francis Galton.

The dataset contains 928 records of adult children’s heights and corresponding mid-parent heights, compiled from Galton’s work. We’ll use the adult children’s heights, measured in inches.

The values were tabulated into approximately one-inch intervals, so the data also demonstrates how measurement precision can affect the apparent shape of a distribution.

Galton • HistData

Create a file named heights.py:

import csv
import io
from urllib.request import urlopen

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

URL = (
    "https://raw.githubusercontent.com/vincentarelbundock/"
    "Rdatasets/master/csv/HistData/Galton.csv"
)

with urlopen(URL, timeout=30) as response:
    content = response.read().decode("utf-8-sig")

reader = csv.DictReader(io.StringIO(content))

heights = np.array(
    [float(row["child"]) for row in reader],
    dtype=float,
)

# Estimate the Gaussian parameters by maximum likelihood.
mu, sigma = stats.norm.fit(heights)

print(f"Observations: {len(heights)}")
print(f"Mean height: {mu:.2f} inches")
print(f"Standard deviation: {sigma:.2f} inches")

# Compare empirical and Gaussian model percentiles.
empirical_95 = np.percentile(heights, 95)
model_95 = stats.norm.ppf(0.95, loc=mu, scale=sigma)

print(f"Empirical 95th percentile: {empirical_95:.2f}")
print(f"Gaussian 95th percentile: {model_95:.2f}")

fig, axes = plt.subplots(1, 2, figsize=(12, 4))

# Histogram and fitted density.
axes[0].hist(
    heights,
    bins=20,
    density=True,
    alpha=0.65,
    label="Observed heights",
)

x = np.linspace(heights.min(), heights.max(), 400)

axes[0].plot(
    x,
    stats.norm.pdf(x, loc=mu, scale=sigma),
    linewidth=2,
    label="Fitted Gaussian",
)

axes[0].set_xlabel("Height (inches)")
axes[0].set_ylabel("Probability density")
axes[0].legend()

# Normal Q-Q plot.
stats.probplot(heights, dist="norm", plot=axes[1])
axes[1].set_title("Normal Q-Q plot")

fig.tight_layout()
fig.savefig("heights.png", dpi=160)
plt.close(fig)

Run it:

python heights.py

The script creates heights.png, containing a histogram with a fitted Gaussian curve and a normal Q-Q plot.

What exactly did we fit?

The call to stats.norm.fit() estimates the Gaussian parameters using maximum likelihood.

For normally modeled observations, the fitted mean is the arithmetic mean of the data. The fitted variance uses the average squared deviation divided by the number of observations, not the sample variance correction used for an unbiased variance estimator.

The fitted curve is a model. It does not force the observations to follow that curve.

Compare two quantities:

• The empirical 95th percentile, calculated directly from observed heights.

• The Gaussian 95th percentile, predicted from the estimated mean and standard deviation.

If the numbers differ, that difference tells you something about the model’s approximation of the observed distribution.

Reading a Q-Q plot

A quantile-quantile plot compares the ordered observations with theoretical quantiles from a normal distribution.

If the Gaussian model is a reasonable approximation, the points should lie roughly along a straight line.

Look for:

• Curvature: possible skewness or a different tail shape.

• Strong departures at the ends: extreme observations differ from normal-model predictions.

• Steps or clusters: discretization, rounding, or repeated measurements.

• An unusual overall shape: possibly mixed populations rather than one homogeneous distribution.

Galton’s measurements are historical, rounded, and drawn from a particular collection of families. They should not be interpreted as a representative height distribution for today’s global population.

That’s an important statistical lesson: even a reasonable-looking fit cannot repair an unrepresentative sample.

5. Bernoulli: the distribution for yes-or-no outcomes

Many datasets do not contain continuous measurements at all.

Consider:

• Did a customer purchase an item?

• Did a scheduled job complete successfully?

• Did a user click a link?

• Does a text block contain a particular word?

Each observation has two possible outcomes.

A Bernoulli random variable is conventionally coded as:

• 1 for success.

• 0 for failure.

The model has one parameter, p, representing the probability of success.

Its probability mass function is:

• P(X = 1) = p

• P(X = 0) = 1 − p

The expected value is p, and its variance is p(1 − p).

For independent, identically distributed Bernoulli trials, the sample proportion of successes estimates p.

Notice that the meaning of “success” is entirely contextual. A success might mean a click, an order, or simply a matching word.

Bernoulli versus binomial

This distinction is easy to remember:

• Bernoulli: one binary trial.

• Binomial: the number of successes across a fixed number of independent Bernoulli trials with the same success probability.

If we record whether each of 100 independently sampled users clicks a button, each individual outcome can be Bernoulli.

The total click count can be modeled as binomial if the independence and common-probability assumptions hold.

Without those assumptions, the binomial model may underestimate or misrepresent variability.

6. Gallery #2: Bernoulli outcomes and the long tail of language

We’ll use a real text corpus: Pride and Prejudice by Jane Austen, available through Project Gutenberg.

The downloadable text is publicly available, and its electronic edition was updated in September 2026.

Project Gutenberg

This single dataset lets us examine two different statistical phenomena:

1. Whether a block of text contains a particular word.

2. How frequently different words occur throughout an entire book.

These are different random variables constructed from the same source data.

Part A: construct real Bernoulli observations

We’ll divide the book into non-overlapping blocks of 100 tokens.

For each block, define an event:

• 1 if the block contains the word “darcy”.

• 0 otherwise.

We can estimate the Bernoulli parameter from one portion of the book and compare it with observed outcomes in another portion.

Part B: examine word-frequency ranks

We’ll then count every word and rank words by descending frequency.

The most frequent words should appear many times. Farther down the ranking, frequencies decrease sharply.

Create words.py:

import re
from collections import Counter
from urllib.request import urlopen

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

URL = "https://www.gutenberg.org/cache/epub/1342/pg1342.txt"

with urlopen(URL, timeout=30) as response:
    raw_text = response.read().decode("utf-8-sig")

# Remove Project Gutenberg's surrounding metadata and license.
start_marker = "*** START OF THE PROJECT GUTENBERG EBOOK"
end_marker = "*** END OF THE PROJECT GUTENBERG EBOOK"

start = raw_text.index(start_marker)
body_start = raw_text.index("\n", start)
body_end = raw_text.index(end_marker, body_start)

novel = raw_text[body_start:body_end]

# Simple English-oriented tokenization.
tokens = re.findall(
    r"[a-z]+(?:'[a-z]+)?",
    novel.lower(),
)

print(f"Total tokens: {len(tokens)}")

# --------------------------------------------------
# Bernoulli: does a 100-token block contain "darcy"?
# --------------------------------------------------

block_size = 100

blocks = [
    tokens[i:i + block_size]
    for i in range(0, len(tokens) - block_size + 1, block_size)
]

events = np.array(
    [int("darcy" in block) for block in blocks],
    dtype=int,
)

split = len(events) // 2

train = events[:split]
test = events[split:]

p_hat = float(train.mean())

outcomes = np.array([0, 1])
predicted = stats.bernoulli.pmf(outcomes, p_hat)

observed = np.bincount(
    test,
    minlength=2,
) / len(test)

print(f"Fitted Bernoulli p: {p_hat:.3f}")
print(f"Held-out positive rate: {test.mean():.3f}")

fig, ax = plt.subplots(figsize=(7, 4))

ax.bar(
    outcomes - 0.18,
    predicted,
    width=0.36,
    label="Model fitted on first half",
)

ax.bar(
    outcomes + 0.18,
    observed,
    width=0.36,
    label="Observed in second half",
)

ax.set_xticks([0, 1])
ax.set_xticklabels(["Absent", "Present"])
ax.set_ylabel("Probability")
ax.set_title("Bernoulli word-presence model")
ax.legend()

fig.tight_layout()
fig.savefig("bernoulli.png", dpi=160)
plt.close(fig)

# --------------------------------------------------
# Word-frequency rank distribution.
# --------------------------------------------------

counts = Counter(tokens)

frequencies = np.array(
    sorted(counts.values(), reverse=True)
)

ranks = np.arange(1, len(frequencies) + 1)

fig, ax = plt.subplots(figsize=(8, 5))

ax.loglog(
    ranks,
    frequencies,
    ".",
    markersize=3,
)

ax.set_xlabel("Word-frequency rank (log scale)")
ax.set_ylabel("Occurrences (log scale)")
ax.set_title("Word-frequency ranks")
ax.grid(True, which="both", alpha=0.2)

fig.tight_layout()
fig.savefig("words_zipf.png", dpi=160)
plt.close(fig)

print(f"Distinct word types: {len(counts)}")
print("Ten most common words:")

for word, count in counts.most_common(10):
    print(f"{word}: {count}")

Run the script:

python words.py

You should now have two plots: bernoulli.png and words_zipf.png.

Understanding the Bernoulli result

The fitted parameter is simply the fraction of training blocks containing the target word.

For instance, if the estimated probability were 0.30, our model would assign:

• 30% probability to a block containing the word.

• 70% probability to a block not containing it.

Those percentages are illustrative; the actual values come from the downloaded text.

The held-out comparison is more interesting than plotting the training proportions alone.

Different portions of a novel feature different characters, conversations, and scenes. Consequently, the probability of mentioning a particular character may change across the book.

A Bernoulli distribution can describe the marginal yes-or-no outcome at one point, while the assumption of a constant probability across all observations may still be unrealistic.

This is a useful introduction to distribution shift: the frequency of an event can change with context.

Our blocks also come from one continuous narrative, so we should not casually assume independent observations.

Understanding the word-frequency result

Word-frequency rankings display a very different pattern.

A small number of words occur repeatedly, while many distinct words are rare.

This behavior is often described by Zipf’s law, which approximately relates a word’s frequency to its rank:

frequency(rank) ≈ C / rank^s

Here, C is a positive scaling constant and s is a positive exponent. In many language datasets, s is roughly near 1 over parts of the ranking, although the precise relationship depends on the corpus and analysis.

Illustrative Zipf-like rank-frequency relationship. Logarithmic axes make multiplicative changes visible.

There is an important technical subtlety here.

Zipf’s law describes a rank-frequency relationship. It is not automatically the same thing as fitting a continuous Pareto distribution to individual observations.

The ranked frequencies are also derived from a finite text. Common words, rare words, punctuation handling, capitalization, and tokenization choices all affect the curve.

A roughly straight line on logarithmic axes is suggestive, not proof of an exact power law.

Nevertheless, the pattern is unmistakable: language contains enormous differences between common and uncommon words.

Why this matters beyond linguistics

Long-tailed frequency patterns appear in many technology applications:

• Product popularity in e-commerce.

• Content views on social platforms.

• Requests for individual URLs.

• Search-query frequencies.

• File-access patterns.

• Customer activity levels.

The important consequence is that an average item can be a very poor description of a typical or popular item.

For caching, recommendation systems, and search infrastructure, this distinction can fundamentally affect system design.

7. The long tail: why rare values deserve attention

So far, we’ve seen two different ideas of a tail.

In the Gaussian case, the tails describe observations far from the mean.

In the word-frequency case, the long tail consists of many low-frequency word types when ranked by popularity.

In operational data, engineers often talk about another phenomenon: a small fraction of observations with unusually large positive values.

Consider response times:

Request group	Illustrative response time
Typical	80–150 ms
Somewhat slow	200–400 ms
Very slow	1,000–2,000 ms
Exceptional	Several seconds or more

Illustrative numbers, not measurements from a production service.

The distribution is constrained at zero: a valid duration cannot be negative.

Meanwhile, high values can extend considerably farther from the center than low values.

That asymmetry rules out an exact Gaussian distribution.

Some common positive-valued models include:

• Exponential: useful for certain waiting-time processes.

• Gamma: a flexible family for positive, skewed values.

• Lognormal: appropriate when the logarithm of a positive variable is approximately normally distributed.

• Weibull: useful in reliability and duration modeling.

The exact choice depends on the process and the observed data.

A long right tail does not automatically mean the data follows a power law. The term “long tail” is often used informally, while rigorous definitions of heavy-tailed distributions are more specific.

8. Gallery #3: measuring and fitting real HTTP response times

Now let’s collect our own real measurements.

We’ll make a small, rate-limited series of HTTP requests to https://example.com/ and record elapsed wall-clock time.

The measurement includes more than server processing: connection establishment, TLS negotiation, network transfer, and client-side effects can all influence the result.

We’ll compare:

1. A fitted Gaussian.

2. A fitted lognormal distribution.

3. The empirical response-time percentiles.

SciPy currently exposes the required norm.fit, lognorm.fit, pdf, and ppf methods.

SciPy v1.18.0 Manual+1

Create latency.py:

from time import perf_counter, sleep
from urllib.error import URLError
from urllib.request import urlopen

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

URL = "https://example.com/"
N_REQUESTS = 50

measurements = []
failures = 0

# Small, sequential sample with a pause between requests.
for _ in range(N_REQUESTS):
    start = perf_counter()

    try:
        with urlopen(URL, timeout=10) as response:
            response.read()
    except (URLError, TimeoutError):
        failures += 1
    else:
        elapsed_ms = (perf_counter() - start) * 1000
        measurements.append(elapsed_ms)

    sleep(0.15)

latencies = np.array(measurements, dtype=float)

if len(latencies) < 20:
    raise RuntimeError(
        "Too few successful requests to fit distributions."
    )

np.savetxt(
    "latency_ms.csv",
    latencies,
    delimiter=",",
    header="latency_ms",
    comments="",
    fmt="%.3f",
)

print(f"Successful requests: {len(latencies)}")
print(f"Failed requests: {failures}")

# Gaussian fit.
mu, sigma = stats.norm.fit(latencies)

# Lognormal fit, constrained to positive support.
shape, loc, scale = stats.lognorm.fit(
    latencies,
    floc=0,
)

# Empirical percentiles.
empirical_p50 = np.percentile(latencies, 50)
empirical_p95 = np.percentile(latencies, 95)

# Model-based 95th percentiles.
normal_p95 = stats.norm.ppf(
    0.95,
    loc=mu,
    scale=sigma,
)

lognormal_p95 = stats.lognorm.ppf(
    0.95,
    shape,
    loc=0,
    scale=scale,
)

print(f"Empirical p50: {empirical_p50:.1f} ms")
print(f"Empirical p95: {empirical_p95:.1f} ms")
print(f"Gaussian p95: {normal_p95:.1f} ms")
print(f"Lognormal p95: {lognormal_p95:.1f} ms")

# Compare fitted probability densities.
x = np.linspace(0, latencies.max() * 1.1, 500)

fig, ax = plt.subplots(figsize=(9, 5))

ax.hist(
    latencies,
    bins="auto",
    density=True,
    alpha=0.5,
    label="Measured latencies",
)

ax.plot(
    x,
    stats.norm.pdf(x, loc=mu, scale=sigma),
    linewidth=2,
    label="Gaussian fit",
)

ax.plot(
    x,
    stats.lognorm.pdf(
        x,
        shape,
        loc=0,
        scale=scale,
    ),
    linewidth=2,
    label="Lognormal fit",
)

ax.axvline(
    empirical_p95,
    linestyle="--",
    label="Empirical p95",
)

ax.set_xlabel("Request duration (ms)")
ax.set_ylabel("Probability density")
ax.set_title("Measured HTTP response times")
ax.legend()

fig.tight_layout()
fig.savefig("latencies.png", dpi=160)
plt.close(fig)

Run the experiment:

python latency.py

You’ll get a plot and a CSV containing the successfully measured durations.

The exact results will differ by location, network conditions, proxies, and time of measurement. If outbound requests are blocked, point the script at an HTTP endpoint you control and are authorized to benchmark.

Why might a lognormal fit work?

A positive random variable is lognormally distributed when its logarithm follows a normal distribution.

In plain terms:

• The original values are strictly positive.

• Their distribution can be strongly right-skewed.

• Extremely large values are possible, though increasingly unlikely.

The lognormal model is not guaranteed to fit HTTP latency. It is simply a useful candidate to compare with the Gaussian model.

Notice that we constrained its location parameter to zero.

This gives the fitted lognormal a natural support of positive response times and avoids the additional flexibility of fitting an arbitrary horizontal offset.

The critical lesson: p95 is different from the average

Suppose a service has the following illustrative latency measurements:

• Mean: 120 ms

• Median: 90 ms

• 95th percentile: 280 ms

• 99th percentile: 950 ms

The mean might sound excellent.

Yet approximately 1% of requests take at least 950 ms, depending on the percentile definition and empirical measurement method.

For a popular service receiving one million requests, 1% represents 10,000 requests.

That illustrates why production service-level objectives often focus on latency percentiles instead of averages alone.

There is also an estimation problem.

With only 50 successful observations, the empirical 99th percentile is determined by the very largest measurements. It is not a reliable estimate of a rare production tail event.

A serious latency study requires many more observations, representative traffic, consistent measurement conditions, and attention to dependency and time variation.

Compare models quantitatively

Looking at a histogram is useful, but visual inspection alone can be misleading.

One common comparison tool is the Akaike information criterion, or AIC.

AIC = 2k − 2 × log-likelihood

Here, k is the number of fitted parameters.

Lower AIC values indicate a better relative balance between in-sample likelihood and model complexity among the candidate models being compared.

Both models in our latency example have two free parameters: Gaussian mean and standard deviation, or lognormal shape and scale with location fixed at zero.

Add this code at the end of latency.py to calculate their AIC values:

normal_loglik = np.sum(
    stats.norm.logpdf(
        latencies,
        loc=mu,
        scale=sigma,
    )
)

lognormal_loglik = np.sum(
    stats.lognorm.logpdf(
        latencies,
        shape,
        loc=0,
        scale=scale,
    )
)

aic_normal = 2 * 2 - 2 * normal_loglik
aic_lognormal = 2 * 2 - 2 * lognormal_loglik

print(f"Gaussian AIC: {aic_normal:.2f}")
print(f"Lognormal AIC: {aic_lognormal:.2f}")

Whichever model produces the lower AIC is preferred under that specific comparison.

However, AIC does not prove that either model is correct. It also does not guarantee good predictions of the extreme tail.

Check held-out observations and compare empirical versus predicted exceedance probabilities when the intended application depends on rare events.

For performance engineering, you should also investigate whether your data represents multiple underlying situations, such as different network paths, cache states, or server workloads.

9. Cherry on the cake: Old Faithful isn’t one bell curve

One of the most memorable counterexamples to the assumption that continuous natural measurements should look Gaussian comes from Yellowstone National Park.

Old Faithful is a famous geyser that erupts repeatedly. You might reasonably expect eruption duration to vary around some typical value, forming one bell-shaped distribution.

But the classic Old Faithful dataset tells a different story.

The standard R dataset contains 272 observations of eruption durations and waiting times until the next eruption. The measurements are real, although the historical eruption-duration values include rounding.

Seminar for Statistics

When you plot the eruption durations, you find two prominent regions of concentration rather than one symmetric central peak.

Broadly speaking, the observations separate into shorter eruptions and longer eruptions.

Fit the data yourself

Create faithful.py:

import csv
import io
from urllib.request import urlopen

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

URL = (
    "https://raw.githubusercontent.com/vincentarelbundock/"
    "Rdatasets/master/csv/datasets/faithful.csv"
)

with urlopen(URL, timeout=30) as response:
    content = response.read().decode("utf-8-sig")

rows = list(csv.DictReader(io.StringIO(content)))

eruptions = np.array(
    [float(row["eruptions"]) for row in rows]
)

waiting = np.array(
    [float(row["waiting"]) for row in rows]
)

mu, sigma = stats.norm.fit(eruptions)
kde = stats.gaussian_kde(eruptions)

x = np.linspace(eruptions.min(), eruptions.max(), 500)

fig, axes = plt.subplots(1, 2, figsize=(12, 4))

axes[0].hist(
    eruptions,
    bins=22,
    density=True,
    alpha=0.55,
    label="Observed durations",
)

axes[0].plot(
    x,
    stats.norm.pdf(x, loc=mu, scale=sigma),
    linewidth=2,
    label="Single Gaussian",
)

axes[0].plot(
    x,
    kde(x),
    linewidth=2,
    label="Kernel density estimate",
)

axes[0].set_xlabel("Eruption duration (minutes)")
axes[0].set_ylabel("Probability density")
axes[0].legend()

axes[1].scatter(
    eruptions,
    waiting,
    alpha=0.6,
    s=18,
)

axes[1].set_xlabel("Eruption duration (minutes)")
axes[1].set_ylabel("Waiting until next eruption (minutes)")
axes[1].set_title("Eruption duration vs waiting time")

fig.tight_layout()
fig.savefig("old_faithful.png", dpi=160)
plt.close(fig)

short = eruptions < 3

print(f"Total observations: {len(eruptions)}")
print(f"Durations under 3 minutes: {short.sum()}")
print(f"Durations at least 3 minutes: {(~short).sum()}")

Run it:

python faithful.py

What should you notice?

In the left-hand plot, the single Gaussian distributes density around one central location.

But the empirical data has two distinct areas of high concentration, with relatively fewer observations between them.

The kernel density estimate can reveal those multiple peaks without assuming one specific parametric distribution.

A KDE is itself a smoothed estimate, however. Its apparent number and sharpness of peaks depend on bandwidth selection, so it should be interpreted alongside the histogram.

In the right-hand plot, the scatterplot shows how eruption duration relates to the subsequent waiting time.

The observations reveal structure that a single one-dimensional Gaussian fit cannot explain well.

The deeper explanation: mixture distributions

One way to describe multimodal data is with a mixture model.

A two-component Gaussian mixture has the form:

f(x) = w × f₁(x) + (1 − w) × f₂(x)

Here:

• f₁ and f₂ are two component probability densities.

• w is the mixing weight, between 0 and 1.

• Each component can have its own mean and standard deviation.

This does not mean the geyser literally has two independent Gaussian-generating mechanisms.

A mixture is a statistical model that can represent heterogeneity in the observations; deciding whether it corresponds to distinct physical processes requires additional evidence.

A similar mistake happens in technology datasets.

Suppose you combine response times for cache hits and cache misses. One group might cluster around tens of milliseconds and the other around hundreds.

A single Gaussian average would hide the operational distinction between those populations.

Before fitting a distribution, ask whether you are measuring one process or accidentally combining several different processes.

That question is frequently more valuable than choosing a more complicated statistical formula.

10. How to choose a distribution in practice

You now have several examples of why the support and shape of the data matter.

Use this reference as a starting point, not a rigid classification.

Data pattern	Candidate model	Example
Binary 0/1 outcome	Bernoulli	Purchase or no purchase
Number of successes in a fixed number of independent trials	Binomial	Clicks across a fixed number of exposures
Counts in a fixed interval	Poisson	Arrivals per minute under suitable assumptions
Symmetric continuous values	Gaussian	Approximately normal measurement noise
Positive right-skewed values	Lognormal or gamma	HTTP response times
Waiting times from a memoryless process	Exponential	Time until the next event under a constant-rate Poisson process
Strongly unequal ranked frequencies	Zipf-like rank law	Word or product popularity
Multiple distinct peaks	Mixture model	Measurements combining different operating conditions

For example, a Poisson model has the property that its mean equals its variance.

Real event counts sometimes have much larger variance, especially when arrival rates change over time or events occur in bursts.

In that situation, alternatives such as a negative binomial model or a time-varying rate model may be more appropriate.

A distribution is useful because its assumptions are plausible and its predictions are adequate, not merely because its curve looks familiar.

11. Five practical checks before trusting a fitted distribution

1. Inspect the support

Can the variable take negative values?

If not, a Gaussian model may assign probability to impossible outcomes.

Can the variable only be 0 or 1?

Then a continuous Gaussian density is not an appropriate generative model for individual observations.

2. Look for skewness and multiple peaks

A histogram is a good starting point.

But the appearance of a histogram depends on bin width, so inspect multiple reasonable choices.

KDEs can supplement histograms for continuous data, but excessive smoothing can conceal important modes.

3. Compare empirical and predicted quantiles

A fitted model is useful when it predicts quantities you care about.

For performance data, that might mean p90, p95, or p99 response times.

For product-demand data, it might mean the threshold exceeded by the most popular 5% of products.

A reasonable mean estimate does not guarantee reasonable tail predictions.

4. Investigate whether observations are independent

Examples of dependence include:

• Repeated measurements of the same server.

• Transactions from the same customer.

• Measurements recorded close together in time.

• Text blocks from the same narrative.

• Measurements from related family members.

Dependence does not necessarily prevent describing the marginal distribution.

But it changes how much independent information the data contains and can invalidate naïve uncertainty estimates.

5. Separate model fit from prediction quality

A model can fit a training dataset beautifully yet perform poorly on new observations.

Use representative validation data when possible.

For time-dependent data, evaluate on later periods rather than relying exclusively on randomly shuffled observations.

For rare-event modeling, gather enough observations to estimate the tail meaningfully.

In some applications, the right answer may be an empirical distribution or quantile estimate rather than a simple named probability family.

12. Try these exercises

The fastest way to build intuition about distributions is to make them fail in controlled experiments.

Exercise A: change the height model

Take the real height dataset and compare its histogram with Gaussian curves using different standard deviations.

Try the fitted standard deviation multiplied by 0.5 and by 2.

Observe how a curve can have the correct center while describing the spread badly.

Then compare the empirical 5th and 95th percentiles with the model’s predictions.

Exercise B: change the word-presence event

Replace the word “darcy” with another character name or a common word.

Try block sizes of 25, 100, and 500 tokens.

Ask yourself:

• How does the Bernoulli success probability change?

• Does the fitted parameter remain stable between the two halves of the book?

• Does the proportion of positive blocks approach 1 for very common words?

• Why does block size change the definition of the random variable?

You can also tokenize another public-domain text to see how much the word-frequency ranks depend on the corpus.

Exercise C: examine latency tails more seriously

Increase the sample size of the latency experiment gradually, using an endpoint you control and an appropriate request rate.

Compare the empirical median, p90, and p95 across successive batches.

Plot the fraction of requests exceeding different thresholds.

Then compare those empirical exceedance rates with predictions from the Gaussian and lognormal models.

Pay particular attention to the threshold range where the models disagree most.

Do not assume that one afternoon’s measurements describe the same system at peak traffic.

13. Final lessons to carry into real projects

The most valuable thing you can learn about probability distributions is not how to memorize their formulas.

It’s how to recognize the structure of a dataset before choosing a mathematical model.

Our gallery demonstrated four distinct patterns:

• Measured heights: a setting where a Gaussian can be a useful approximate model, but rounding, sampling, and population heterogeneity matter.

• Binary word-presence events: a natural use of Bernoulli probabilities, with context-dependent success rates.

• Word frequencies: a sharply unequal rank-frequency pattern that a bell curve cannot describe.

• HTTP response times: positive measurements for which right-skewed models and tail percentiles can be more informative than averages.

• Old Faithful eruption durations: a real-world demonstration that even natural continuous measurements can be strongly multimodal.

The common theme is that distributions describe data-generating behavior, not just plot shapes.

A good statistical workflow starts with the variable’s possible values, inspects the observations, proposes plausible models, and evaluates the predictions that matter for the actual task.

Once you develop that habit, distribution fitting stops being an exercise in drawing curves and becomes a tool for making better decisions.

Further reading

For deeper exploration, these resources provide authoritative documentation and reusable examples:

• SciPy statistical distributions documentation — distribution functions, fitting, probability calculations, and diagnostics.

• NumPy random Generator documentation — reproducible simulation with modern random-number generators.

• HistData: Galton height measurements — provenance and limitations of the historical height dataset.

• Project Gutenberg: Pride and Prejudice — the public-domain text used for word-frequency experiments.

• R documentation: Old Faithful dataset — dataset definitions, measurements, and original research references.

Your next step

Run all four scripts and build your own distribution gallery. Keep the original observations, generated plots, fitted parameters, and empirical percentiles together so you can reproduce your comparisons.

Then take one dataset from a project of your own—retail orders, server telemetry, household energy usage, or another everyday measurement—and ask three questions:

1. What values are possible?

2. What shape does the observed data actually have?

3. Which distribution makes useful predictions without hiding important behavior?

If you can answer those questions and test the answers with code, you’ve learned something far more practical than simply recognizing a bell curve.

You’ve learned how to start thinking statistically.