How to make a Weibull fit on an existing histogram?

Question

I created a histogram plot based on my datasets. I would like to create a Weibull fit for this histogram. I used scipy and the stats.weibull function, but unfortunately, it does not work.

Do you have an idea of how to use the stats.weibull in this case?

Here is the code:

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import scipy.stats as stats

data = 'Figures/Histogram/Histogram.xlsx'
hist= pd.read_excel('Histogram/Histogram.xlsx')

# x= hist['DeltaT_value'] 
x= hist['DeltaT_-250_2017']
X=x[(x>0)]

plt.figure(figsize=(15,4))
plt.hist(X, bins= np.arange (0,1500,25), color='#0504aa', edgecolor ='red', rwidth= 0.8)
plt.ylabel('Number of EL')
plt.xlabel('Delta T (years CE) between EL')
plt.xlim(0, 401)
plt.xticks(np.arange(0,401,25))
plt.yticks(np.arange(0,2.2,1))`


# Weibull 
####

shape, loc, scale = stats.weibull_min.fit(X)
x = np.linspace(stats.weibull_min.ppf(0.01, shape, loc=loc, scale=scale), stats.weibull_min.ppf(0.99, shape, loc=loc, scale=scale), 100)
plt.plot(x, stats.weibull_min.pdf(x, shape, loc=loc, scale=scale), 'r-', lw=5, alpha=0.6, label='weibull')

I tried this:

shape, loc, scale = stats.weibull_min.fit(X)
x = np.linspace(stats.weibull_min.ppf(0.01, shape, loc=loc, scale=scale), stats.weibull_min.ppf(0.99, shape, loc=loc, scale=scale), 100)
plt.plot(x, stats.weibull_min.pdf(x, shape, loc=loc, scale=scale), 'r-', lw=5, alpha=0.6, label='weibull')

Unfortunately, it seems another graph is created on top of the histogram instead of a fit.

Answer 1

After working on it for a little while, I found a solution:

shape, loc, scale = stats.weibull_min.fit(X, floc = 0, f0 = 1)
W = np.linspace(stats.weibull_min.ppf(0.001, shape, loc=loc, scale=scale), stats.weibull_min.ppf(0.99, shape, loc=loc, scale=scale), 1000)

p =  stats.weibull_min.pdf(W, shape, scale = scale)

Robert Dodier's comments on my original post are very interesting and I will try to see how I can use a "log likelihood function" on the original dataset rather than a fit on the histogram. This is the first time I will use a log likelihood function. If anyone has some advice, I would gladly take it.

Thanks.

Answer 2

Mat, here's what I get with the approach I was arguing for, which is to work with the likelihood function as derived from quantile data (i.e., the cumulative summation of histogram bars).

I've devised some code for the computer algebra system Maxima (https://maxima.sourceforge.io) to do the calculations. I'll link to the code, but it's probably difficult to use at this point; mostly what I want to say is that this business about working with the likelihood function is a good path forward. See also my comments on this other question about deriving the likelihood function for binned data: https://stats.stackexchange.com/questions/11176/can-anyone-explain-quantile-maximum-probability-estimation-qmpe/442966#442966

The code I'm using is the package robert-dodier/qmpe under https://github.com/maxima-project-on-github/maxima-packages . That in turn makes use of as-yet-unreleased functions (to be released in the next version of Maxima, namely 5.47) from the Maxima package distrib; that unreleased version may be found at https://sourceforge.net/p/maxima/code/ci/master/tree/share/distrib/ .

Looks like the data shown in the histograms is just:

q: [0, 25, 50, 75, 100, 125, 150, 175, 200, 225, 250, 275, 300, 325];
n: [0,  1,  3,  2,   2,   3,   0,   2,   1,   0,   1,   1,   0,   1];

with q being the quantiles (0, 25, 50, ... 325 years) and n being the number of events in each bin (0, 1, 2, or 3 being the only observed values).

I tried a lognormal distribution and a Weibull distribution. Other distributions are possible by the same method, I just haven't written the code for them yet; I might still try it.

For the lognormal, I get mu = 4.546, sigma = 0.7794, with negative log likelihood nll = 2.528. For the Weibull, I get shape = 1.488, scale = 136.2, with nll = 2.491. Here are figures comparing lognormal to data and Weibull to data. The green line is the final fit; the red line, barely visible in the Weibull plot, is from the initial parameters for the optimization.

As you can see, the Weibull fit is a little better as determined from nll, but neither one is really very good, and in fact it's likely there's no distribution that gives a good fit -- there are only a few data per bin and they're pretty noisy.

As a postscript, here is the code I am working with to generate the figures shown. As I was saying, I don't expect anyone to run this code, it's just to present the general method about working with the likelihood function as derive from the histogram data.

/* inspired by: https://stackoverflow.com/questions/75090377/how-to-make-a-weibull-fit-on-an-existing-histogram */

/* data shown in histogram https://i.stack.imgur.com/uYt70.png
 * values appear to be just 0, 1, 2, or 3
 */

q: [0, 25, 50, 75, 100, 125, 150, 175, 200, 225, 250, 275, 300, 325];
n: [0,  1,  3,  2,   2,   3,   0,   2,   1,   0,   1,   1,   0,   1];
p_unnormalized: makelist (lsum (n1, n1, firstn (n, k)), k, 1, length (n));
p: p_unnormalized / last (p_unnormalized);

load (distrib);

load ("qmpe.mac");

quux: construct_qmpe (cdf_weibull (u, shape, scale), 'u, '[shape, scale], mle_weibull, 1e-4, [1, 0]);

mumble: quux (q, p);

set_plot_option ([svg_file, "./SO-75090377-weibull-fit-compare-cdf.svg"]);
plot_qmpe_comparison_cdf (q, p, cdf_weibull (u, shape, scale), u, '[shape, scale], quantile_weibull (0.99, shape, scale), assoc ('initial, mumble), assoc ('final, mumble));

set_plot_option ([svg_file, "./SO-75090377-weibull-fit-compare-pdf.svg"]);
plot_qmpe_comparison_pdf (q, p, pdf_weibull (u, shape, scale), u, '[shape, scale], quantile_weibull (0.99, shape, scale), assoc ('initial, mumble), assoc ('final, mumble));

foo: construct_qmpe (cdf_lognormal (u, location, scale), 'u, '[location, scale], mle_lognormal, 1e-4, [1, 0]);

baz: foo (q, p);

set_plot_option ([svg_file, "./SO-75090377-lognormal-fit-compare-cdf.svg"]);
plot_qmpe_comparison_cdf (q, p, cdf_lognormal (u, location, scale), u, '[location, scale], quantile_lognormal (0.99, location, scale), assoc ('initial, baz), assoc ('final, baz));

set_plot_option ([svg_file, "./SO-75090377-lognormal-fit-compare-pdf.svg"]);
plot_qmpe_comparison_pdf (q, p, pdf_lognormal (u, location, scale), u, '[location, scale], quantile_lognormal (0.99, location, scale), assoc ('initial, baz), assoc ('final, baz));

ev (assoc ('nll, mumble), assoc ('final, mumble), nouns, numer);
ev (assoc ('nll, baz), assoc ('final, baz), nouns, numer);