# Should histograms be normalized first before fitting them?

Assume I have some data that follow a power law and I would like to estimate the exponent $$a$$. An obvious way to do so is to bin the data and then fit the power law to the histogram. However I found that the result can be unreliable as it not only depends on the number of free parameters, but also on the binning.

For example, if I want to fit the histogram using the actual number of data in each bin, I also need to fit a scaling factor $$c$$ as well as the exponent:

$$f(x) = c \cdot x^a$$

If I'm not mistaken however, to improve the results I should generally reduce the number of free parameters as much as possible. I this case I can, simply by normalizing the histogram (i.e. by dividing the height of each bin by the maximum height). This way I can drop the scaling factor again (as $$c$$ should be unity):

$$f(x) = x^a$$

While this actually allows me to get closer to the truth, this also makes my fit much more susceptible to the number of bins:

When I allow for the scaling factor to be fit as well, the result is pretty independent of the number of bins, but the exponent is pretty far off; when I normalize my data, I can achive better results for the exponent, but the number of bins matters quite a lot. What is the best way to approach this?

Code for reproduction (Python):

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

power = 4
data = np.random.power(power, 1000)
bins = [3, 5, 20, 50, 100, 900]

fig, axs = plt.subplots(2, 3, figsize=(9, 4), constrained_layout=True)
for i, ax in enumerate(axs.flatten()):
bin_heights, bin_borders = np.histogram(data, bins=bins[i])
bin_centers = bin_borders[:-1] + np.diff(bin_borders) / 2
x_interval_for_fit = np.linspace(bin_borders[0], bin_borders[-1], 1000)

# fit using the 'raw' heights of the bins
popt, _ = curve_fit(
lambda x, c, a: c*x**a, bin_centers, bin_heights, p0=[1, power]
)
ax.plot(
x_interval_for_fit, x_interval_for_fit ** popt[1],
label=f'numbers: c = {popt[0]:.1f}, a = {popt[1]:.1f}', c='green'
)

# fit using the normalized heights of the bins
bin_heights = bin_heights/max(bin_heights)
popt, _ = curve_fit(
lambda x, a: x**a, bin_centers, bin_heights, p0=[4]
)
ax.plot(
x_interval_for_fit, x_interval_for_fit**popt[0],
label=f'normalized: a = {popt[0]:.1f}', c='orange'
)

# plotting the histogram etc.
ax.bar(
bin_borders[:-1], height=bin_heights, width=np.diff(bin_borders),
align='edge', facecolor='c', alpha=0.5, zorder=0,
)
ax.text(0.2, 0.3, f'# bins: {bins[i]}', transform=ax.transAxes)
ax.plot(
x_interval_for_fit, x_interval_for_fit**power, '--', c='r',
label=f'actual: a = {power}'
)
ax.legend()

fig.show()

• curve_fit is not correct for this purpose. One of the most salient reasons is that the variance of any histogram bin is directly proportional to its area, so when the bar areas vary substantially across the histogram--as they do in all your cases--it often produces terrible estimates. You can fit curves to histograms, by using maximum likelihood estimates based on the multinomial distribution of the bin counts.
– whuber
Commented Feb 9, 2022 at 17:04

Generally speaking, you shouldn't be using histograms to estimate parameters (unless you really have to; see bottom part of this answer for details). It would be much better to use a more formal method of inference, such as Maximum likelihood, to estimate $$\alpha$$.

### Defining the Problem

From the numpy manual

np.random.power refers to the "Power distribution" with density function

$$f(x|\theta) = \theta x^{\theta - 1}, \ 0 \leq x \leq 1, \ \theta > 0$$

Note that the power you have chosen (power = 4) corresponds to this $$\theta$$, which is related to your $$\alpha$$ as $$\alpha = \theta - 1$$. Note also that we don't have to worry about estimating $$c$$ in this formulation.

Finally, we can note that this distribution is the inverse of the Pareto distribution and can also be viewed as a special case of the Beta distribution (e.g., $$\text{Beta}(\theta, 1)$$).

### The MLE

Without going into the derivation, it is easy to show that the Maximum Likelihood Estimator of $$\theta$$ here is simply

$$\hat\theta = \frac{-n}{\sum_{i=1}^n \log x_i}$$

where $$n$$ is the sample size and $$x_i$$ refers to the $$i^{th}$$ data point. Equivalently (by invariance), we can put this in terms of $$\alpha$$ as

$$\hat\alpha = \frac{-n}{\sum_{i=1}^n \log x_i} + 1.$$

This approach has two major advantages

1. There is no binning parameter to worry about (and thus you don't have to worry about normalizing the histogram)
2. This approach will generally lead to better estimation of $$\alpha$$ than your current approach.

### Censored Likelihood

In some rare cases, you may actually have to use the histogram to estimate $$\alpha$$. For instance, you may be using a dataset for which you don't have access to the specific numbers, but only previously binned values. In this setting, you still don't have to worry about choosing bin size.

Here, I would recommend looking into interval censored likelihood methods, which is going to be much better than simply fitting a curve to the histogram.

• Hey, thanks a lot for the excellent answer! Unfortunately it shows that I've made the problem unnecessarily more complicated by using numpy.random.power, making the naive assumption that it simply draws from a function in the form of $f(x|\theta) = x^\theta$, which is why that's what I'm trying to fit. That's on me for not checking first. Sorry. The real data that I'm trying to fit is actually much closer to $x^{-4}$, but I wanted to make my example a bit simpler, and numpy.random.power doesn't allow you to draw from functions with $\theta < 1$.
– mapf
Commented Feb 9, 2022 at 19:16
• That being said, would $\hat\theta = \frac{-n}{\sum_{i=1}^n \log x_i}$ still give me the correct result in the $f(x|\theta) = x^\theta$, $\theta \approx -4$ case?
– mapf
Commented Feb 9, 2022 at 19:19
• @mapf, I'm not sure, but I think you are looking for the Pareto distribution (numpy.random.pareto). I would spend some time researching the Pareto distribution. The maximum likelihood estimator may be different, but you can find details online or post another question. Commented Feb 9, 2022 at 19:51
• I see, I will look into it!
– mapf
Commented Feb 9, 2022 at 21:07