可靠性工程常用分布_二参数威布尔分布计算实例

威布尔分布为一种常用的连续分布,通常用于机械零件及应力强度疲劳寿命分析。说到连续分布就自然会想到离散型分布，这里顺便总结分类下常用的离散型分布和连续型分布。
1）离散型分布：二项分布，泊松分布；
2）连续型分布：正态分布，对数正态分布，指数分布，威布尔分布，伽马分布。

关于威布尔分布，标准的威布尔分布有二参数和三参数两种形式。这里主要分析下两种计算二参数威布尔分布形状参数β,以及尺度参数η的方法(概率纸法,以及Excel计算法）。
如下面第一张截图可知二参数威布尔累积分布（故障概率）函数以及可靠性函数分别如下：
F(t)=1-exp[-(t/η)^β]————公式1；
R(t)=1-F(t)=exp[-(t/η)^β]—–公式2；
由公式也可得知威布尔分布形状参数β与浴盆曲线的关系如下（参考如下第二张截图）：
1）β<1时，处于浴盆曲线早夭期，故障率逐渐降低；
2）β=1时，威布尔分布为指数分布，处于浴盆曲线底部，此时故障率恒定；
3）β>1时，处于浴盆曲线磨损期，故障率逐渐增高；

举一实例如下截图1，20个产品测试1000小时，其中有14个产品在测试过程中分别在如下时间失效：70, 128, 204, 291, 312, 377, 473, 549, 591, 663, 748, 827, 903, 955小时，假设符合二参数威布尔分布，计算威布尔分布的形状参数β,以及尺度参数η。

解题方法1——威布尔概率纸法；
1）将失效时间依次从小到大排序作为第一行（用作X轴）—如下截图2；
2）查询中位秩表（如下截图4），样品总量为20，则中位秩值依次为：0.034,0.083,0.132……0.672作为第二行（用作Y轴）; 将各个中位秩值依次对应第一步排序的失效时间（如下截图2）
这里中位秩值也可以通过公式：Median Rank rj=(j-0.3)/(N+0.4)（如下截图5）计算获得，这里j为失效次序（1,2….),N为样本总量20，例如第一个失效的中位秩值为（1-0.3）/(20+0.4)=0.0343;第二个失效的中位秩值为（2-0.3）/(20+0.4)=0.0833,等等。
3）在威布尔概率纸上依次对应画点如下截图3；
4）对这些点用线性拟合的方法画出一条直线（图上右下方实线）；
5）在威布尔概率纸右上方，经过纵坐标故障率0.623（这里0.623取值如何得出？解答如下：F(t)=1-exp[-(t/η)^β]， 假设t=η,这F(η)=1-exp(-1)=1-0.367=0.623，参考如下截图6)，画出一条平行直线（图上左上方虚线），与威布尔概率纸上方的交叉点即为β值，得出该例β值约为1.3.
6）由步骤5的解释，经过纵坐标0.623画一条平行于X轴的直线，与拟合直线的交叉点即为η的取值，得出该例η值约为900.

解题方法2——Excel拟合计算法；
1）这里第一步还是先写出故障概率函数如下：
F(t)=1-exp[-(t/η)^β]————公式1；
2）将公式1转为包含（X,Y关系）的线性方程，依次推导过程如下，可参考如下第一张截图；
F(t)=1-exp[-(t/η)^β]，
1-F(t)=exp[-(t/η)^β]；
1/(1-F(t))=exp[(t/η)^β]；
Ln[1/(1-F(t))]=(t/η)^β;
Ln{Ln[1/(1-F(t))]}=Ln[(t/η)^β];
Ln{Ln[1/(1-F(t))]}=β[Ln(t)-Ln(η)];—-公式3
假设Ln{Ln[1/(1-F(t))]}=Y, Ln(t)=X,则公式3可转为：Y=βX-βLn(η)—-公式4；
3）将失效时间70, 128, …955依次从小到大排序,见如下第二张截图B列；
4）通过第一种方法中的中位秩值计算法：Median Rank rj=(j-0.3)/(N+0.4)，依次得出F(t)的值为0.0343,0.0833，….0.6716.见如下第二张截图C列;
5）用Excel依次计算出Y以及X的值，见如下第二张截图F列和G列；
6）在Excel中对如下第二张截图F列和G列的数值进行线性拟合，线性拟合步骤入下:
1: 如下截图3,选中X,Y对应区域数据,插入图表X,Y(Scatter)，接下一步；
2: 如下截图4,对X，Y轴坐标对应坐标点进行赋值，接下一步；
3: 如下截图5,得出X，Y坐标轴上的散点图，接下一步；
4: 如下截图6，选中散点进行线性拟合，接下一步；
5: 如下截图7，得出拟合直线，显示线性方程Y=1.2818X-8.7644。
7）由上线性方程Y=1.2818X-8.7644，得出β=1.2818；βLn(η)=8.7644；
8）算出η=exp{8.7644/β)=exp{8.7644/1.2818)=932.22

顺便拓展下：1)经常听说风扇的L10寿命,为何MTTF=7*L10？：
L10为测试有10%的的不良率的产品寿命，通常风扇失效模型采用威布尔分布，则F(t)=0.1，F(t)=1-exp[-(t/η)^β]=0.1， 
则exp[-(t/η)^β]=0.9，t=L10,η为尺度参数，也可以表述成特征寿命，
exp[(L10/η)^β]=1/0.9=1.111;—-公式5；
MTTF为产品寿命期望值为63.2%不良时对应的预期时间，由以上分析，可设η=MTTF（这里未找到出处….）
则(L10/MTTF)^β=Ln1.111=0.10536;
这里MTTF与L10的关系为与β相关的值，而不是简单写成7，当然如果风扇厂商有多年数据统计，可能得出β值，这里发现只有当β=1.155时，才约为7倍关系。如果有不同推导方式的朋友可以分享下~~

2）当产品寿命服从威布尔分布时MTTF=ηΓ(1/β+1)—-公式6；这里Γ为伽马函数（参考如下截图）；但这里在Excel中用什么函数计算呢？
答案：Excel中未提供直接计算伽马函数的公式，但提供了GAMMALN函数，它给出了伽马函数的自然对数，因此可以通过e的幂函数计算，为GAMMA(X)=EXP(GAMMLN(X))；
这里可以举个特例，当β取值为1时，Γ(1/1+1)=1，可自行到Excel中计算。

写到这里发现第一个问题与第二个问题应该是关联的，应该可以通过第二个问题反推出MTTF与L10的关系：
例如由公式6可得—η=MTTF/Γ(1/β+1)—-公式7;
将公式7代入公式5，得exp{L10/[MTTF/Γ(1/β+1)]}^β=1/0.9=1.111—-公式8；
由公式8得{L10/[MTTF/Γ(1/β+1)]}^β=Ln(1/0.9);—-公式9；
公式9才是真正能反应MTTF与L10的关系，与β相关。

个人笔记，欢迎交流~~