title: How to get the confidence interval of a Weibull distribution using Python?
body: <p>I want to perform a probability Weibull fit with 0.95% confidence bounds by means of Python. As test data, I use fail cycles of a measurement which are plotted against the reliability R(t).</p>
  
  <p>So far, I found a way to perform the Weibull fit, however, I still do not manage to get the confidence bounds. The Weibull plot with the same test data set was already performed with origin, therfore I know which shape I would "expect" for the confidence interval. But I do not understand how to get there.</p>
  
  <p><a href="https://i.stack.imgur.com/Wbn5F.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/Wbn5F.png" alt="example for confidence interval"></a></p>
  
  <p>I found information about Weibull confidence intervals on <a href="http://reliawiki.org/index.php/The_Weibull_Distribution" rel="nofollow noreferrer">reliawiki</a>(cf. Bounds on Reliability based on Fisher Matrix confidence bounds) and used the description there to calculate the variance and the upper and lower confidence bound (R_U and R_L).</p>
  
  <p>Here is a working code example for my Weibull fit and my confidence bounds with the test data set based on the discription of <a href="http://reliawiki.org/index.php/The_Weibull_Distribution" rel="nofollow noreferrer">reliawiki</a> (cf. Bounds on Reliability). For the fit, I used a OLS model fit.</p>
  
  <pre class="lang-py prettyprint-override"><code>import os, sys
  import pandas as pd
  import numpy as np
  import matplotlib.pyplot as plt
  import matplotlib as mpl
  from scipy.optimize import curve_fit
  import math
  import statsmodels.api as sm
  
  def weibull_ticks(y, pos):
      return "{:.0f}%".format(100 * (1 - np.exp(-np.exp(y))))
  
  def loglog(x):
      return np.log(-np.log(1 - np.asarray(x)))
  
  class weibull_example(object):
  
      def __init__(self, dat):
          self.fits = {}
          dat.index = np.arange(1, len(dat) + 1)
          dat.sort_values('data', inplace=True)
          #define yaxis-values
          dat['percentile'] = dat.index*1/len(dat)
          self.data = dat
  
          self.fit()
          self.plot_data()
  
      def fit(self):
          #fit the data points with a the OLS model
          self.data=self.data[:-1]
          x0 = np.log(self.data.dropna()['data'].values)
          Y = loglog(self.data.dropna()['percentile'])
          Yx = sm.add_constant(Y)
          model = sm.OLS(x0, Yx)
          results = model.fit()
          yy = loglog(np.linspace(.001, .999, 100))
          YY = sm.add_constant(yy)
          XX = np.exp(results.predict(YY))
          self.eta = np.exp(results.params[0])
          self.beta = 1 / results.params[1]
          self.fits['syx'] = {'results': results, 'model': model,
                              'line': np.row_stack([XX, yy]),
                              'beta': self.beta,
                              'eta': self.eta}
  
          cov = results.cov_params()
          #get variance and covariance
          self.beta_var = cov[1, 1]
          self.eta_var = cov[0, 0]
          self.cov = cov[1, 0]
  
      def plot_data(self, fit='yx'):
          dat = self.data
          #plot data points
          plt.semilogx(dat['data'], loglog(dat['percentile']), 'o')
          fit = 's' + fit
          self.plot_fit(fit)
  
          ax = plt.gca()
          formatter = mpl.ticker.FuncFormatter(weibull_ticks)
          ax.yaxis.set_major_formatter(formatter)
          yt_F = np.array([0.001, 0.01, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5,
                           0.6, 0.7, 0.8, 0.9, 0.95, 0.99])
          yt_lnF = loglog(yt_F)
          plt.yticks(yt_lnF)
  
          plt.ylim(loglog([.01, .99]))
  
      def plot_fit(self, fit='syx'):
          dat = self.fits[fit]['line']
          plt.plot(dat[0], dat[1])
  
          #calculate variance to get confidence bound
          def variance(x):
              return (math.log(x) - math.log(self.eta)) ** 2 * self.beta_var + \
                     (self.beta/self.eta) ** 2 * self.eta_var - \
                     2 * (math.log(x) - math.log(self.eta)) * (-self.beta/self.eta) * self.cov
  
          #calculate confidence bounds
          def confidence_upper(x):
              return 1-np.exp(-np.exp(self.beta*(math.log(x)-math.log(self.eta)) - 0.95*np.sqrt(variance(x))))
          def confidence_lower(x):
              return 1-np.exp(-np.exp(self.beta*(math.log(x)-math.log(self.eta)) + 0.95*np.sqrt(variance(x))))
  
          yvals_1 = list(map(confidence_upper, dat[0]))
          yvals_2 = list(map(confidence_lower, dat[0]))
  
          #plot confidence bounds
          plt.semilogx(dat[0], loglog(yvals_1), linestyle="solid", color="black", linewidth=2,
                   label="fit_u_1", alpha=0.8)
          plt.semilogx(dat[0], loglog(yvals_2), linestyle="solid", color="green", linewidth=2,
                   label="fit_u_1", alpha=0.8)
  
  def main():
      fig, ax1 = plt.subplots()
      ax1.set_xlabel("$Cycles\ til\ Failure$")
      ax1.set_ylabel("$Weibull\ Percentile$")
  
      #my data points
      data = pd.DataFrame({'data': [1556, 2595, 11531, 38079, 46046, 57357]})
      weibull_example(data)
  
      plt.savefig("Weibull.png")
      plt.close(fig)
  
  if __name__ == "__main__":
      main()
  </code></pre>
  
  <p>The confidence bounds in my plot look not like I expected. I tried a lot of different 'variances', just to understand the function and to check, if the problem is just a typing error. Meanwhile, I am convinced that the problem is more general and that I understood something false from the description on <a href="http://reliawiki.org/index.php/The_Weibull_Distribution" rel="nofollow noreferrer">reliawiki</a>. Unfortunately, I really do not get what's the problem and I do not know anyone else I can ask. In the internet and on different forums, I did not find an appropriate answer.</p>
  
  <p>That's why I decided to ask this question here. It's the first time I ask a question in a forum. Therefore, I hope that I explained everything sufficiently and that the code example is useful. 
  Thank you very much :)</p>