How to determine Weibull Law Parameters with Python? (Part 2)

If we start from the Weibull Probability that we determined previously:

After a few simple mathematical operations (take the log of both sides), we can convert this expression into a linear expression, such as the following one:

This means that we can pose:

and

And the previous equation can be written:

Let’s import first the python modules we will need for the study:

• os is a classic module always useful to handle the link with files and the system
• numpy is here for the numerical calculations
• matplotlib will be useful to draw the graphs
• scipy will provide us with an useful function to do regression of the curve and fit the parameters
• pandas will facilitate the interaction with the data

import os
import numpy as np
import matplotlib.pyplot as pl
from scipy import stats
import pandas as pd

Let’s define now a function to calculate the Y=ln[-ln(1-P_F)] value mentionned previously

def weibull_inv(p):
    return np.log(-np.log(1.0-p))

Now let’s suppose that we have the following data into a data.csv file

Let’s use pandas to import those datas into a DataFrame

#importing the CSV file
Samples_DF = pd.read_csv('data.csv')

#Sorting the data according to stress values and re-indexing
Samples_DF = Samples_DF.sort_values(by='Sigma Failure(MPa)')
Samples_DF = Samples_DF.reset_index(drop=True)

Let’s plot that now

pl.plot(Samples_DF['Sigma Failure(MPa)'],'o')
pl.grid()
pl.xlabel('# Sample')
pl.ylabel('Stress at fracture (Principal)')
pl.title("Failure Stress results of experiment")

Let’s add to the dataframe table 2 new columns corresponding respectively to:

• The law of probability (an empirical distribution function)
• The inversed Weibull distribution law

Samples_DF['Proba'] = (Samples_DF.index-Samples_DF.index[0]+1) / (len((Samples_DF.index))+1)
Samples_DF['Weibull'] = weibull_inv(Samples_DF['Proba'])

We obtain the following dataframe:

Sample N  failure load(N)  Sigma Failure(MPa)     Proba   Weibull
0          1             1168                17.7  0.047619 -3.020227
1          2             1175                17.8  0.095238 -2.301751
2          3             1195                18.1  0.142857 -1.869825
3          4             1287                19.5  0.190476 -1.554433
4          5             1360                20.6  0.238095 -1.302197
5          6             1366                20.7  0.285714 -1.089240
6          7             1426                21.6  0.333333 -0.902720
7          8             1445                21.9  0.380952 -0.734859
8          9             1452                22.0  0.428571 -0.580505
9         10             1491                22.6  0.476190 -0.435985
10        11             1524                23.1  0.523810 -0.298490
11        12             1537                23.3  0.571429 -0.165703
12        13             1551                23.5  0.619048 -0.035543
13        14             1564                23.7  0.666667  0.094048
14        15             1564                23.7  0.714286  0.225351
15        16             1577                23.9  0.761905  0.361224
16        17             1637                24.8  0.809524  0.505750
17        18             1656                25.1  0.857143  0.665730
18        19             1676                25.4  0.904762  0.855000
19        20             1682                25.5  0.952381  1.113344

In order to perform a linear regression, we have to extract w and lnsw as below:

w = Samples_DF['Weibull']
lnsw = np.log(Samples_DF['Sigma Failure(MPa)'])

Let’s calculate now the parameters of the linear regression:

with:

and

m, lnsm0, *t = stats.linregress(lnsw,w)
sigma0 = np.exp(- lnsm0 / m)

print('m=', m)
print('sigma0=',sigma0)

Here’s what I get:

m= 9.23254658432
sigma0= 23.3758678631

m expresses the dispersion around the Weibull Law

sigma0 is the scale factor of the Weibull Law.

It has the same order than the stress at failure.

Let’s draw all of that now with mathplotlib:

pl.figure()
pl.plot(lnsw,w)
pl.plot(lnsw,w,'*')

x = np.arange(lnsw.iloc[0], lnsw.iloc[-1]+0.1, 0.1)
y = (lambda x : m * x + lnsm0)(x)
pl.plot(x, y)
pl.plot()
pl.grid()
pl.ylabel('log(-log(1 - Probability of fracture))')
pl.ylabel('-log(Stress at fracture)/m')
pl.title("Weibull Analysis of experiment data")