Density Function Estimation
Histogram to estimate density function
Whether the goal is prediction or statistical inference, density function helps us in either case. To understand how the mass of a random variable is distributed across the datapoints — estimating an accurate density function of the random variable would be a real help. Before seeing how we can estimate a density function, let’s take look at how density would help us with simple example.
Say you have a data of general and malnourished children heights in normal distribution:
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import norm, expon, binom
import seaborn as sns
# General population heights
pop_heights = np.sort(np.random.normal(loc=170, scale=10, size=500))
# Malnourished children heights
mal_heights = np.sort(np.random.normal(loc=150, scale=12, size=500))To get density function(PDF) of these random variables:
pdf_pop = norm.pdf(pop_heights, loc=np.mean(pop_heights), scale=np.std(pop_heights))
pdf_mal = norm.pdf(mal_heights, loc=np.mean(mal_heights), scale=np.std(mal_heights))norm.pdf() takes a datapoint and return it’s density in Standard Normal Distribution. We can change the SND to our data’s scale by giving the mean loc and standard deviation scale
plt.figure(figsize=(10,5))
plt.plot(pop_heights, pdf_pop, color='#4C78A8', lw=2.5, label='General Children PDF')
plt.plot(mal_heights, pdf_mal, color='#E15759', lw=2.5, label='Malnourished Children PDF')
plt.fill_between(pop_heights, pdf_pop, color='#4C78A8', alpha=0.1)
plt.fill_between(mal_heights, pdf_mal, color='#E15759', alpha=0.1)
plt.title("Estimating Children Heights")
plt.xlabel("Height (cm)")
plt.ylabel("Density")
plt.legend()
plt.show()
Observations from the Graph:
Shift in Mean Heights:
- The peak of the general children is centered around ~170 cm.
- The peak of the malnourished children is centered around ~150 cm.
- This shows that malnourished children are, on average, shorter than the general children.
Overlap Between Distributions:
- There is some overlap between the two distributions in the 150–165 cm range.
- This means that within this height range, it is harder to distinguish whether a child belongs to the malnourished group or the general group based on height alone.
Spread (Variance) Difference:
- The red curve is wider, showing greater variability in heights.
- The blue curve is narrower, showing less variability.
Tail Behavior:
- The general children distribution extends more towards taller heights ~200 cm.
- The malnourished distribution is more spread but with very few children beyond ~185 cm, showing a ceiling effect in growth.
Practical Implication:
- Height can be used as a proxy indicator for nutritional status, but due to overlap, statistical inference is needed to make reliable group distinctions rather than relying on raw comparison.
From these density functions with predictions we can infer probabilities like :
What fraction of children are taller than 180 cm?
What’s the most common height range?
also performing hypothesis tests, finding confidence intervals, doing interpretation, etc. This helps in applications like public health or clothing industry sizing.
Some other examples like:
Time until next customer arrives in queue management, call centers, and traffic modeling.
Probability of heads in coin flips, the foundation of hypothesis testing and fairness checks in experiments.
For all these above examples we created the data on specific distribution in mind. We can plot histogram for the data to observe how the mass is distributed just like I did for the above two plots. The real-world data may won’t follow these text book distributions. We can see some weird shapes like:
Plots help intuition, but density functions give numbers we can compute, compare, and infer from.
These density functions allow us to move beyond pictures and actually quantify aspects of the distribution, like where it tends to center, how spread out it is, or whether it has heavy tails. This is where density estimation comes in — instead of just plotting data, we aim to mathematically estimate the underlying density function that generated it. That way, we can use the data not only to visualize but also to reason, compare, and make inferences. A simple method to estimate the density of a random variable is Histogram density estimation.
Histogram Density Estimation :
Let X = { x1, x2, x3, …, xn} be a random variable with n observations sampled from some population density function f .
Now divide the interval [a, b] into m equal bins
each bin has an equal width of:
Our goal is to estimate the density at each data point xi. Instead of calculating the density for every single point, we can group nearby points into bins and use the density of each bin as a good approximation for the density of the points it contains. Now assign each data point xi to a bin according to the value it contains. Let vj be the number of observations in bin j, then the probability of an observation getting into bin j would be :
This is only an estimate of the true probability of an observation
falling into bin j. If the underlying population has density f, then the actual probability of landing inside bin j is
By the law of large numbers:
Now we know each bin(rectangle) has a width of h and we know the probability(area) of getting into that bin. Our goal is to find the density(height) of the bin, we know
that would make our bin’s density
The histogram estimator is defined by :
which we can write more clearly as
Note that this density function is an estimation, that means
The bin-width h is smoothing parameter,
Get Mohith’s stories in your inbox
Join Medium for free to get updates from this writer.
As the bin-width h increases, each bin averages over more data points, so individual points have less influence on the histogram. (high-bias, over-smoothing)
As the bin-width h decreases, each bin contains fewer data points, so individual points have more influence on the histogram. (high-variance, under-smoothing)
This is where the trade-off between bias and variance comes in, high-bias results in overly smooth, inaccurate density estimates, while high-variance makes the histogram overly sensitive to individual points, leading to a noisy estimate.
Optimal Smoothing Parameter(h):
As a loss function we will use the Integrated Squared Error(ISE):
The Risk or Mean Integrated Squared Error(MISE) for a bin-width h:
The Risk can be written as,
Minimizing risk corresponds to balancing bias and variance. In histogram estimation we do that with smoothing parameter h.
To find out optimal h value, we can use cross-validation to compute the risk of the density function with bin-width of h and compare the risk with other h values risks.
Normally, cross-validation for histograms would mean doing leave-one-out CV. That is — for every single data point, we throw it out, rebuild the histogram with the remaining n−1 samples, then check how well the left-out point is predicted. Doing this n times for each candidate bin-width h quickly becomes painful.
The nice thing about the formula
is that it’s a shortcut. Instead of refitting the histogram n times, you only need to compute it once per h. So the shortcut doesn’t eliminate the need to try different bin widths — it just saves you from the heavy lifting of recomputing the histogram n times for each candidate h.
Expectation and Variance of Histogram Density Function :
Consider fixed x and fixed m, and let Bj be the bin containing x. Then,
Enough of theory :
let’s see the HDE working on Old Faithful Dataset, particularly the column named waiting , which represent the waiting times for the next volcanic eruption.
df = pd.read_csv('data.tsv', sep='\t')
waiting_times = df['waiting'].to_numpy()
m = 50 # Initial no of bins
min_value = waiting_times.min()
max_value = waiting_times.max()
bins = np.linspace(min_value, max_value, m+1) # m+1 for edges
h = bins[1] - bins[0] # Bin width
# Counting the values with in each bin
lengths = np.zeros(m)
for i in range(m):
if i < m-1:
lengths[i] = len(waiting_times[(waiting_times >= bins[i]) & (waiting_times < bins[i+1])])
else:
lengths[i] = len(waiting_times[(waiting_times >= bins[i]) & (waiting_times <= bins[i+1])])
areas = lengths / len(waiting_times)
densities = areas / h
plt.figure(figsize=(10, 6))
plt.plot(bins[:-1], densities, color='#E15759', linewidth=2, label='Density Estimate')
plt.xlabel("Waiting Time (minutes)", fontsize=14)
plt.ylabel("Density", fontsize=14)
plt.title("Density Estimate of Waiting Times", fontsize=16, pad=14)
plt.legend()
plt.tight_layout()
plt.savefig("density_estimate.png", dpi=250)
plt.show()
Let’s make a bar plot and see how each bins density looks like:
plt.figure(figsize=(10, 6))
plt.bar(bins[:-1], densities, width=h, align='edge', color='#4C78A8', alpha=0.7, edgecolor='black')
plt.xlabel("Waiting Time (minutes)", fontsize=14)
plt.ylabel("Estimated PDF", fontsize=14)
plt.title("Histogram-Based PDF Estimation", fontsize=16, pad=14)
plt.tight_layout()
plt.savefig("histogram_pdf.png", dpi=250)
plt.show()
We got this histogram plot from the idea of Histogram Density Function, then it should be similar to matplotlib’s histogram plot given the waiting_times ,bins = mand density = True
plt.figure(figsize=(10, 6))
plt.hist(waiting_times, bins=m, density=True, color='#4C78A8', alpha=0.7, edgecolor='black')
plt.xlabel("Waiting Time (minutes)", fontsize=14)
plt.ylabel("Density", fontsize=14)
plt.title("Normalized Histogram of Waiting Times", fontsize=16, pad=14)
plt.tight_layout()
plt.savefig("normalized_histogram.png", dpi=250)
plt.show()
Let’s take a look at how over-smoothed and under-smoothed density functions looks like
Over-Smoothing :
m = 80
min_value = waiting_times.min()
max_value = waiting_times.max()
h = (max_value - min_value) / m
print("h :",h) # 0.6625
bins = np.linspace(min_value, max_value, m+1)
lengths, _ = np.histogram(waiting_times, bins=bins) #np.histogram does the same thing as our counting loop is doing
areas = lengths / len(waiting_times)
density = areas / h
plt.figure(figsize=(10, 6))
plt.plot(bins[:-1], density, color='#E15759', linewidth=2, label='Density Estimate (Overfitting)')
plt.xlabel("Waiting Time (minutes)", fontsize=14)
plt.ylabel("Density", fontsize=14)
plt.title("Overfitting Density Estimate (m=80)", fontsize=16, pad=14)
plt.legend()
plt.tight_layout()
plt.savefig("overfitting_density.png", dpi=250)
plt.show()
See the oscillations and how the data is too much focused on the sample data.
Under-Smoothing :
m = 10
min_value = waiting_times.min()
max_value = waiting_times.max()
h = (max_value - min_value) / m
print('h :', h) #5.3
bins = np.linspace(min_value, max_value, m+1)
lengths, _ = np.histogram(waiting_times, bins=bins)
areas = lengths / len(waiting_times)
density = areas / h
plt.figure(figsize=(10, 6))
plt.plot(bins[:-1], density, color='#E15759', linewidth=2, label='Density Estimate (Underfitting)')
plt.xlabel("Waiting Time (minutes)", fontsize=14)
plt.ylabel("Density", fontsize=14)
plt.title("Underfitting Density Estimate (m=10)", fontsize=16, pad=14)
plt.legend()
plt.tight_layout()
plt.savefig("underfitting_density.png", dpi=250)
plt.show()
Each bin of width 5.3 is averaging more data, reducing the influence of the individual datapoints.
Confidence Band for density function :
Remember that our density function is a random variable here ? Since we estimated the density function for the waiting_times (m = 50 PDF)and we also agree with the law of large numbers, we can estimate where the population density function f(x) is using Confidence Bands. To know more about Confidence Intervals you can see Confidence Is The New Truth
The 1-α confidence band for a histogram density function looks like:
alpha = 1 - 0.95 # 95% confidence level
z = norm.ppf(1 - alpha/2)
n = len(waiting_times)
SE = np.sqrt(densities * (1 - densities * h) / n)
lower = np.maximum(densities - z * SE, 0)
upper = densities + z * SE
plt.figure(figsize=(10, 6))
plt.plot(bins[:-1], densities, color='#E15759', linewidth=2, label='Density Estimate')
plt.plot(bins[:-1], upper, color='#4C78A8', linestyle='--', label='Upper Bound')
plt.plot(bins[:-1], lower, color='#4C78A8', linestyle='--', label='Lower Bound')
plt.fill_between(bins[:-1], lower, upper, color='gray', alpha=0.3, label='95% Confidence Interval')
plt.xlabel("Waiting Time (minutes)", fontsize=14)
plt.ylabel("Density", fontsize=14)
plt.title("Density Estimate with Confidence Intervals", fontsize=16, pad=14)
plt.legend()
plt.tight_layout()
plt.savefig("confidence_intervals.png", dpi=250)
plt.show()
Conclusion :
We started with the idea that density functions aren’t just pretty plots — they give us a mathematical handle on how data is spread out. Histograms gave us the first step into density estimation, where the choice of bin-width h controls everything: too small and the estimate becomes noisy, too large and it gets oversmoothed. We also saw how to attach confidence bands to remind ourselves that the estimated density is still a random variable.
Histograms are simple and intuitive, but they’re also discontinuous, with abrupt jumps at bin edges. In the next blog, we’ll look at a really cool Kernel Density Estimation (KDE), which smooths out those jumps and gives a more flexible way of estimating densities.
