Should histograms be normalized first before fitting them?

Question

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()

Answer 1

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.