Fit a curve to a binned distribution

fit_with_integral.py

# %%
from typing import Iterable

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

if __name__ == "__main__":
    rng = np.random.default_rng(123456)
    # %% generate a skewed beta distribution with a=2, b=8
    X = rng.beta(2, 8, size=1000)
    plt.hist(X)

    # %% bin random values into NON-equidistant bins
    y, bins = np.histogram(
        X, bins=[0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.4, 0.5, 0.6, 1]
    )
    x = (bins[1:] + bins[:-1]) / 2

    plt.plot(x, y)

    # %% default fitting w/o accounting for the integral
    def f_distr(x: Iterable[float], a: float, b: float, c: float) -> float:
        """Beta function with area normalisation `c`"""
        _y = c * beta.pdf(x, a, b)
        return _y

    popt, pcov = curve_fit(f_distr, x, y)
    print(f"Fit WITHOUT bin width: {popt[0]:.3f}, {popt[1]:.3f}")

    _ = plt.plot(x, y, marker="o")
    _ = plt.plot(x, f_distr(x, popt[0], popt[1], popt[2]), ls="--", c="red")

    # %% fitting using integrals in bins
    def f_integral(x: Iterable[float], *args) -> float:
        """Integral of the beta distribution"""
        result = [
            integrate.quad(f_distr, v_min, v_max, args=args)
            for v_min, v_max in zip(bins[:-1], bins[1:])
        ]

        return [v[0] for v in result]

    popt, pcov = curve_fit(f_integral, x, y, p0=(1, 2, 3))
    print(f"Fit WITH bin width: {popt[0]:.3f}, {popt[1]:.3f}")

    _ = plt.plot(x, y, marker="o")
    y_fitted = f_integral(x, *popt)
    _ = plt.plot(x, y_fitted, ls="--", c="red")