Fitting a Weibull distribution in python with stats.exponweib.fit

Question

I've been trying to fit a Weibull distribution with stats.exponweib.fit - there is not a fit in Scipy for just Weibull, so, one needs to utilize the fit for exponential Weibull and set the first shape parameter to 1. However, when the stats.exponweib.fit function is fed with data from a Weibull distribution with known shape parameters - the fit returns a different set of shape parameters. Some example code to display this behavior would be:

from numpy import random, exp, log
import matplotlib.pyplot as plt
from scipy import stats
import csv

# Expoential Weibull PDF 
def expweibPDF(x, k, lam, alpha):
    return (alpha * (k/lam) *
            ((x/lam)**(k-1))  *
            ((1 - exp(-(x/lam)**k))**(alpha-1)) *
            exp(-(x/lam)**k))

# Expoential Weibull CDF
def exp_cdf(x, k, lam, alpha):
    return (1 - exp(-(x / lam)**k))**alpha

# Expoential Weibull Inverse CDF
def exp_inv_cdf(p, k, lam, alpha):
    return lam * ( - log( (1 - p)**(1/alpha) ))**(1/k)

# parameters for the fit - alpha = 1.0 reduces to normal Webull
# the shape parameters k = 5.0 and lam = 1.0 are demonstrated on Wikipedia:
# https://en.wikipedia.org/wiki/Weibull_distribution

alpha = 1.0
k0 = 5.0
lam0 = 1.0
x = []
y = []

# create a Weibull distribution
random.seed(123)
n = 1000  
for i in range(1,n) :
    p = random.random()
    x0 = exp_inv_cdf(p,k0,lam0,alpha)
    x += [ x0 ]
    y += [ expweibPDF(x0,k0,lam0,alpha) ]


# now fit the Weibull using python library
# setting f0=1 should set alpha = 1.0
# so, shape parameters should be the k0 = 5.0 and lam = 1.0

(exp1, k1, loc1, lam1)  = stats.exponweib.fit(y,floc=0, f0=1)

print (exp1, k1, loc1, lam1)

The output here is:

(1, 2.8146777019890856, 0, 1.4974049126907345)

I would have expected:

(1, 5.0, 0, 1.0)

When we plot the curves:

# plotting the two curves
fig, ax = plt.subplots(2, 1)
ax[0].plot(x,y, 'ro', lw=2)
ax[1].plot(x,stats.exponweib.pdf(x,exp1,k1,loc1,lam1), 'ro', lw=2)
plt.show()

We get the following curves showing the input data from a known Weibull distribution with shape factors k=5 and lambda=1 and output from the exponweib.fit with different shape factors.

Input Weibull data and output from exponweib.fit

First post on stackoverflow - so, hopefully the above is the right way to frame a question. Welcome any ideas on the above and any pointers on posting :)

Check the answer to this (duplicate?) question: [How to fit a weibull distribution to data using python](http://stackoverflow.com/questions/34908492/how-to-fit-a-weibull-distribution-to-data-using-python) . In your case, your `x` variable contains a random sample of values from the original distribution, so this is what you should pass to `stats.exponweib.fit` — Pablo Arnalte-Mur, Feb 23 '16 at 14:39
Thanks Pablo for the quick response. Exactly right - i was fitting the PDF instead of the samples. And yes, similar issue is addressed in the previous question: [How to fit a weibull distribution to data using python](http://stackoverflow.com/questions/34908492/how-to-fit-a-weibull-distribution-to-data-using-python) — TurboToad, Feb 23 '16 at 18:47
Possible duplicate of [How to fit a weibull distribution to data using python?](https://stackoverflow.com/questions/34908492/how-to-fit-a-weibull-distribution-to-data-using-python) — Him, Aug 08 '18 at 14:30

score 1 · Answer 1 · answered Nov 04 '20 at 06:18

Scipy does offer the "standard" Weibull distribution that you will find on Wikipedia. The function you should use for this is scipy.stats.weibull_min

Scipy's implementation of Weibull can be a little confusing, and its ability to fit 3 parameter Weibull distributions sometimes gives wild results. You're also unable to fit censored data using Scipy. I suggest that you might want to check out the Python reliability library which makes the process of creating, fitting, and using probability distributions fairly simple compared to Scipy.

score 0 · Answer 2 · answered Oct 28 '20 at 01:06

In my notebook I tried the WeibullMaxFactory of OpenTURNS to fit a Weibull distribution on your x

import openturns as ot
from openturns.viewer import View
sample = ot.Sample(x, 1) # formats your x into a 'Sample' of dimension = 1
distribution = ot.WeibullMaxFactory().build(sample) # fits a Weibull to your data 
graph = distribution.drawPDF() # build the PDF
graph.setLegends(['Weibull'])
View(graph)

To obtain the Weibull parameters:

print(distribution)
>>> WeibullMax(beta = 0.618739, alpha = 2.85518, gamma = 1.48269)

Is there a `scipy` or `numpy` way of doing this instead of using a not-so-well-known package? — Rajesh Swarnkar, Apr 05 '22 at 16:51

Fitting a Weibull distribution in python with stats.exponweib.fit

2 Answers2