韦伯尔分布/威布尔分布/Weibull distribution 置信区间确定

 韦布尔分布理论介绍见：

Equations of supported distributions — reliability 0.8.3 documentation

韦伯分布（威布尔分布，Weibull distribution)_心态与做事习惯决定人生高度的博客-CSDN博客_威布尔分布

其pdf, cdf为：

 那么可以根据散点的方法确定置信区间。其思路为采用多个散点，基于分布函数（Distribution Function，CDF），找到相应的分位数就可以。比如95%的置信区间，可以确定0.25 和0.975所对应的x值，进而确定相应的置信区间。 

这里需要注意的是，威布尔分布是非对称分布，不可以根据解析式确定精确的置信区间。在部分文章里有关于非对称分布的近似求解方法，但求得的结果与置信区间有误差，故本文没采用。感兴趣的读者可参考 https://en.wikipedia.org/wiki/Confidence_interval。 
相关代码见：

import numpy as np
import matplotlib.pyplot as plt


# define the pdf of weibull distribution
def weibpdf(x, scale, shape):
    return (shape / scale) * (x / scale)**(shape - 1) * np.exp(-(x / scale) ** shape)
# define the cdf of weibull distribution
def weibcdf(x, scale, shape):
    return 1 - np.exp( -1 * (x/scale) ** shape)
# define the confidence interval of weibull distribution
def weib_confidence_interval(scale, shape, confidence_level = 0.95, sampleNum = 50000):
    assert(sampleNum >= 1000)
    x = np.linspace(1, scale * 3, sampleNum)
    ycdf = np.array( list( map(weibcdf, x, np.ones(len(x)) * scale, np.ones(len(x)) * shape) ))
    lower_level, upper_level = 0.5 * (1 - confidence_level), confidence_level + 0.5 * (1 - confidence_level)
    x_lower_level = x[np.min( np.where(ycdf >= lower_level))]
    x_upper_level = x[np.min( np.where(ycdf >= upper_level))]
    print('scale, shape, confidence level', scale, shape, confidence_level)
    print('confidence interval',x_lower_level, x_upper_level, 'corresponding cdf', weibcdf(x_lower_level, scale, shape), weibcdf(x_upper_level, scale, shape))

    return x_lower_level, x_upper_level

scales, shapes = [50, 50, 50, 30, 70], [1.5, 2.5, 3.5, 1.5, 1.5]
colors = ['b', 'g', 'r', 'c', 'm', 'y', 'k', 'w'][:len(shapes)]
for scale, shape, color in zip(scales, shapes, colors):
    #x = np.arange(1, scale*3)
    x = np.linspace(1, scale * 3, 1000)
    ypdf = np.array( list( map(weibpdf, x, np.ones(len(x)) * scale, np.ones(len(x)) * shape) ))
    ycdf = np.array( list( map(weibcdf, x, np.ones(len(x)) * scale, np.ones(len(x)) * shape) ))

    x_lower_level, x_upper_level = weib_confidence_interval(scale, shape)

    plt.subplot(2,1,1)
    plt.plot(x, ypdf, color = color , label='scale=%d, shape=%0.2f'%(scale, shape))
    plt.legend()
    plt.subplot(2,1,2)
    plt.plot(x, ycdf, color = color ,  label='scale=%d, shape=%0.2f'%(scale, shape))
    for xi in [x_lower_level, x_upper_level]:
        plt.vlines(xi, ymin= 0, ymax=1, color = color , linestyles='--', linewidth = 0.5)
    plt.legend()
plt.show()

结果：