基于贝叶斯理论的定时截尾加速验收抽样试验优化设计方法

本发明是一种基于贝叶斯理论的加速寿命验收抽样试验优化设计方法，属于加速寿命试 验和可靠性验证试验技术领域，用于解决可靠性与系统工程领域的技术问题。

寿命验证试验是为了确定产品的寿命特征是否达到了研制要求而进行的试验，它包括寿 命鉴定试验和寿命验收试验，不同的阶段展开不同的试验。

寿命验收抽样试验在实施前必须要有一个详细的试验计划，如何对待试产品进行抽样， 如何根据抽样结果进行决策推断，如何对于试验方案中的两类风险计进行算等，都要有明确 的说明，一个好的试验设计方案可以根据待试产品的抽样样本，运用严谨的数学理论，对产 品的寿命参数进行分析，完成对寿命参数的统计推断，进而对产品的寿命进行验证，辅助决 策者对产品的寿命做出正确的决策，不仅如此，一个好的试验设计方法，还能在先进试验理 论的指导下，运用最小的试验样本量或最短的试验时间，投入最少的人力，物力和财力，进 行试验方案设计，并且不会影响寿命统计验证的准确度。随着科技的发展，产品成本越来越 高，用于验证产品寿命的试验费用不断增长，如何在节省费用的条件下不降低验收试验的准 确度，这些都体现了寿命验收抽样试验设计的重要性。

随着科技和工艺的发展，对于拥有高可靠长寿命的产品而言，利用传统的寿命证试验对 寿命指标进行验证时经常需要花费较长的试验时间，成本非常高，一般的企业无法承受。尤 其是对于产品型号种类多而批量小的生产厂家，用传统的寿命试验方法进行寿命验证是一件 极其让人头疼的事。

基于上述原因，通过高量级应力的试验方法来快速地评估产品在实际的使用环境中的寿 命水平越来越成为人们关注的焦点，因此，近年来一些学者提出一些加速环境的寿命验证试 验来快速评估寿命指标的设想，希望利用这种加速试验能达到缩短寿命验证试验时间、降低 试验成本的目的。一直以来，加速试验在快速激发产品缺陷，有效改善产品设计和制造方面 都有着巨大的优势，目前加速寿命试验、可靠性强化试验等都得到广泛的应用，并不断走向 成熟。那么，如何利用加速试验来进行寿命指标的验证是当前亟待解决的一个问题。

目前在加速验证试验方案方面，针对加速因子未知的情况，基于指数分布，Yum和Kim (1990)在两个不同的应力水平下设计了定数截尾验收抽样试验，但方案计算十分复杂，并且 误差较大，Hsieh(1994)丰富了Yum和Kim(1990)的方法，给出了更好的简化计算公式，且最 小化了总截尾数量。基于Weibull和对数正态分布，Bai(1993)等在两个高于使用条件的应力 水平下设计了定数截尾的加速寿命抽样验收方案，为了进一步研究定时截尾试验方案， Bai(1995)等后来又在Bai(1993)研究的基础上加入了预期的时间约束；Seo(2009)等在这些方法 的基础上，对Weibull分布的形状参数是非常数的情况进行了定数加速寿命抽样验收方案设 计。对于加速因子已知的情况，Kim和Yum(2009)假定Weibull分布的形状参数未知，设计了 定时截尾下的加速寿命抽样方案，随后他们通过假定Weibull分布的形状参数已知，又开展了 混合截尾加速寿命抽样验收方案设计的研究。

传统的产品寿命统计验收抽样试验设计方法是以数理统计中的大样本统计分析理论为基 础的统计决策方法，主要是根据现场试验样本提供的信息对所考虑的寿命指标进行假设检验， 做出接受或拒绝统计假设的决策。Bayes方法是寿命验收抽样试验设计的最佳选择之一，Bayes 产品寿命统计验收抽样试验设计方法是以Bayes理论为基础，在充分利用可以利用的信息（如 产品历史信息、研制各阶段的试验信息、类似型号产品的信息、专家信息等）和只进行少量 现场试验的情况下，对产品的寿命指标进行统计决策的方法。

对于Bayes方法应用在寿命验证试验技术方面，N.Balakrishnan(2007)等在寿命试验截断 于先前固定时间的情况下运用Bayes方法设计了验收抽样方案，Chien-Tai Lin(2008)等基于指 数分布下的I型和II型混合结尾抽样设计了Bayes可变抽样方案，TaChen Liang和 Ming-ChungYang(2011)对基于混合截尾样本的指数寿命分布进行了最优Bayes抽样方案设计， Zeinab Amin(2012)等基于Pareto寿命模型设计了验收抽样方案，Mohammad Saber Fallah  Nezhad(2012)等在验收试验中考虑产品的检验误差，运用Bayes方法设计了抽样验收试验。 贝叶斯试验设计一般很难获得后验分布的数学表达式，因此衍生了解决该问题的两种方法： 一种是基于仿真，一种是基于大样本理论。Gladys D.C.Barriga等人提出指数-威布尔寿命分 布和Arrhenius模型下，基于马尔可夫链蒙特卡罗（MCMC）的ALT的贝叶斯方法。在《贝 叶斯可靠性》一书中，作者Michael S.Hamada,Alyson G.Wilson等基于平均风险准则和验后 风险准则，利用MCMC方法分别设计了二项分布，泊松分布和Weibull分布下的Bayes验收 抽样方案。

但是，随着加速应力的施加，加速因子的不确定性会改变原有的使用方和生产方的双方 风险和所需要的试验样本数量，目前，考虑加速因子的不确定性的验收抽样试验设计还很少 有人在研究，为了使得寿命验收抽样试验的设计在加速因子不确定性的影响下足够强健，并 考虑定时截尾方案在工程应用中的广泛和实用性，我们可以结合加速寿命试验技术和Bayes 方法设计基于贝叶斯理论的定时截尾加速验收抽样试验方案。

本发明的目的是为了解决在加速试验中，由于加速因子的估计和计算不准确而影响寿命 验证试验的结果。因此，将加速因子的波动以概率分布的形式进行描述，然后在此基础上设 计寿命验收抽样试验可以削弱波动带来的影响，结合贝叶斯方法对于先验信息的利用，可以 提高验证试验的准确程度。由于引入了加速因子以及产品寿命的先验分布，故使用MCMC方 法进行贝叶斯定时截尾的加速验收抽样试验方案优化设计。

本发明是基于贝叶斯理论的定时截尾加速验收抽样试验优化设计方法，包括以下几个步 骤：

步骤一、确定产品的寿命分布和验证指标参数，然后建立统计假设；

步骤二、对选定的产品设计试验决策法则，确定其试验和接收流程；

步骤三、基于验后风险准则推导满足双方风险关于验证指标参数的约束条件；

步骤四、基于历史数据给出验证指标参数的先验分布，基于马尔科夫蒙特卡洛(Markov  Chain Monte Carlo,MCMC)方法，利用WinBUGS14计算满足条件的方案集合；

步骤五、确定试验的费用约束，计算得到最优方案。

本发明的优点在于：

（1）本发明方法首次将贝叶斯理论引入到加速验收抽样试验的优化设计中，基于贝叶斯 理论，可充分利用试验前的历史数据，相似产品信息等，在确定验证指标参数先验分布的情 况下对试验进行优化设计，可以在保证试验精度的前提下降低试验费用；

（2）本发明将验后风险准则应用在试验方案的求解过程中，验后风险准则借助于现场的 试验数据，反向推导指标参数的先验分布中与现场试验得出的寿命水平不同的概率，其计算 方法侧重对于参数的先验分布主观认可。

（3）本发明方法首次考虑了试验中加速因子不确定性的影响，并将影响量化表示，在此 基础上设计寿命验收抽样试验可以削弱波动带来的影响，提高试验的准确程度。

图1是本发明基于贝叶斯理论的定时截尾加速验收抽样试验优化设计方法的流程图；

图2是本发明步骤一中的抽样特性曲线图；

图3是本发明步骤二中的定时截尾抽验规则方框图；

图4是本发明步骤三中的定时截尾方案求解流程图；

下面将结合附图和实施例对本发明作进一步的详细说明。

本发明将验后风险准则应用在试验方案的求解过程中，建立基于贝叶斯理论的定时截尾 加速验收抽样试验优化设计框架，给出了基于贝叶斯理论的加速验收抽样试验优化设计的具 体步骤。实施例中选用指数分布为产品的寿命分布来阐述本发明提出的贝叶斯优化设计方法。

本发明是一种基于贝叶斯理论的定时截尾加速验收抽样试验优化设计方法，如图1所示， 包括以下几个步骤：

步骤一、确定产品的寿命分布，寿命验证指标参数，然后建立统计假设；

产品的寿命决定于设计与制造中对其功能、结构、原材料等的选择及质量控制过程中各 种随机因素的影响。它是一个服从一定统计规律的随机变量，一般用寿命的分布函数(也称累 积分布函数)来描述。

从寿命试验中得到的数据，是从某批产品(总体)中得到的一个样本，用统计推断的理论， 可以判断出产品的寿命分布，得到累积分布函数。由此可以计算产品的可靠性参数，如可靠 度、失效率、概率密度函数，以及各种寿命特征量，如平均寿命、可靠寿命、特征寿命等。

在选择统计试验方案时，首先要对产品失效前寿命分布进行假设，因为统计试验方案是 用在某一寿命分布情况下。常用的分布形式主要分为连续型和离散型两种。

多数产品寿命服从连续型随机变量的概率分布，常用的有指数分布、正态分布、威布尔 分布等。某些产品以工作次数、循环周期数作为其寿命度量单位，如开关的开关次数等，这 时可以用离散型随机变量的概率分布来描述其寿命分布的规律，如二项分布、泊松分布等。

例如，在可靠性理论中，指数分布是最基本、最常用的分布，电子产品的寿命和复杂系 统的故障时间均可用指数分布来描述。指数分布寿命的失效率为常数。很多电子产品在早期 失效之后及损耗故障期之前，产品的失效率基本上是稳定的。

指数分布的密度函数有两种表达形式：

f(t)=λexp(-λt)与 $f\left(t\right)=\frac{1}{\theta }\exp (-t/\theta )---\left(1\right)$

式中，λ为指数分布的失效率；θ为指数分布的平均寿命。

两个表达式实质相同，在使用条件下，参数λU和θU(产品使用条件下的失效率和平均寿 命)的关系为λU=1/θU，相应的指数分布函数的形式为：

F(t)=1-exp(-t/θ)与F(t)=1-exp(-λt)                     (2)

假设试验的加速因子为AF，则根据指数分布场合，其加速应力条件下对使用条件下的加 速因子定义为：

$\mathrm{AF}=\frac{{\theta }_U}{{\theta }_A}=\frac{{\lambda }_A}{{\lambda }_U}---\left(3\right)$

为了方便先验分布的选取，本发明选取加速条件下的失效率λA作为其寿命验证指标参数， 根据协定的双方风险，在产品的抽样特性曲线(OC Curve)上选择对应的检验上下限λ0(=1/θ0) 和(λ1=1/θ1)，如图2所示，然后建立统计假设如下：

H0：λA≤λ0·AF H1：λA>λ1·AF                     (4)

双方商定，当产品批的加速条件下的失效率λA≤λ0·AF时，以大概率接收这批产品，即规 定生产方风险为α，则接收概率L(λ)≥1-α；当产品批的加速应力下失效率λA≥λ1·AF时，以 小概率接收(高概率拒收)这批产品，即规定使用方风险为β，则接收概率L(λ)≤β。

步骤二、对选定的产品设计试验决策法则，确定其试验和接收流程；

确定了产品的分布及统计假设后，需要确定定时截尾方案的抽验规则，定时截尾试验是 指对n个样品进行试验，事先规定试验截尾t0，到了时刻t0所有试验样品停止试验，利用试 验数据评估产品的可靠性特征量。按试验过程中对发生故障的产品所采取的措施，又可分为 有替换和无替换两种方案。有替换指的是试验中某产品发生了故障，立即用新产品代替，保 持整个实验过程中样本数不变，而无替换是指当产品发生故障就立刻撤去，在试验过程中， 随着故障产品的增加而使样本减少，在本发明中使用的是无替换定时截尾试验。

定时截尾试验方案的优点在于最大积累试验时间是事先确定的，因此在试验以前就可以 确定试验设备，人力物力的最大需要量，便于计划管理，因此得到广泛应用。

以指数分布型产品的定时截尾方案为例，方案通常记为(c，T)，其中T为截尾试验时间， c为接收拒绝数。在本发明中，我们使用恒定应力加速试验(CSALT)，如图3所示，寿命试验 的决策法则即为：

（1）在产品批中选择n个样品进行CSALT，加速因子为AF，试验为无替换试验；

（2）试验进行到试验累计时间到达预定值T时停止试验，记录试验过程中的失效数；

（3）设在试验过程中出现了r故障，如果r≤c，则认为产品批合格，接收；如果r>c， 则认为产品不合格，拒收批产品。

因此，定时截尾试验设计的主要任务就是选择合适的c和T。

步骤三、基于验后风险准则推导满足双方风险关于验证指标参数的约束条件；

本发明将验后风险准则应用在试验方案的求解过程中，验后风险准则借助于现场的试验 数据，反向推导指标参数的先验分布中与现场试验得出的寿命水平不同的概率，其计算方法 侧重对于参数的先验分布主观认可。

下面介绍验后风险准则(Posterior Risk Criteria)下的双方风险计算原理。

基于概率论与数理统计的原理，基于验后风险准则原理，弃真风险(生产方风险)α的计算 公式如式(5)所示：

$\alpha =P(R\in {\theta }_0|Z\in D_1)={\int }_{{\theta }_0}p\left(R\right|Z\in D_1)\mathrm{dR}=\frac{{\int }_{{\theta }_0}P(Z\in D_1|R)\pi \left(R\right)\mathrm{dR}}{{\int }_{\theta }P(Z\in D_1|R)\pi \left(R\right)\mathrm{dR}}---\left(5\right)$

式(5)中,θ表示产品寿命参数R的全集，即取值范围p(R|Z∈D1)表示在给定Z∈D1的条件下 R的概率密度函数，π(R)是寿命参数R的先验分布。

由(5)式所见，弃真风险(生产方风险)α的物理含义为：根据决策法则，依据抽样结果拒绝 原假设的情况下，但是，产品总体的寿命却是满足要求的。其数学含义为：在拒绝原假设的 前提下，产品寿命参数的验后分布满足要求的概率。

在验后风险准则中，采伪风险(使用方风险)β的定义如式(6)所示：

$\beta =P(R\in {\theta }_1|Z\in D_0)={\int }_{{\theta }_1}p\left(R\right|Z\in D_0)\mathrm{dR}=\frac{{\int }_{{\theta }_1}P(Z\in D_0|R)\pi \left(R\right)\mathrm{dR}}{{\int }_{\theta }P(Z\in D_0|R)\pi \left(R\right)\mathrm{dR}}---\left(6\right)$

式(6)中，p(R|Z∈D0)表示在现场试验数据得出Z∈D0基础上的R的概率密度函数。

由式(6)可见，采伪风险β的物理含义为：根据决策法则，依据抽样结果接受原假设的情 况下，但是，产品总体的寿命却是不满足要求的。其数学含义为：在接受原假设的前提下， 产品寿命参数的验后分布不满足要求的概率。

以指数分布型产品为例，首先推导指数分布定时截尾下产品的接收概率L(λ)。

根据指数分布的累积分布函数F(t)=1-exp(-t/θ)可知,产品的可靠度R(t)=exp(-λt)，到时间 t时，n个产品中出现r个故障的概率为

$C_n^{r}F{\left(t\right)}^{r}R{\left(t\right)}^{n-r}---\left(7\right)$

到时间t时，产品的故障率r≤c，从而产品被接受的概率为

$L\left(\lambda \right)=\underset{r=0}{\overset{c}{\Sigma }}{C}_n^{r}F{\left(t\right)}^{r}R{\left(t\right)}^{n-r}---\left(8\right)$

由于λ的值一般都很小，故将R(t)=exp(-λt)泰勒展开可得

$R\left(t\right)=e^{-\mathrm{\lambda t}}=1-\mathrm{\lambda t}+\frac{1}{2!}{\lambda }^2{t}^2-···\approx 1-\mathrm{\lambda t}$

F(t)=1-R(t)=λt               (9)

即可得接受概率

$L\left(\lambda \right)=\underset{r=0}{\overset{c}{\Sigma }}{C}_n^r{\left(\mathrm{\lambda t}\right)}^r{(1-\mathrm{\lambda t})}^{n-r}---\left(10\right)$

在nλt≤5，F(t)≤10%的条件下，二项概率可用泊松概率近似，于是得到：

$L\left(\lambda \right)=\underset{r=0}{\overset{c}{\Sigma }}{e}^{-\mathrm{n\lambda t}}\frac{{\left(\mathrm{n\lambda t}\right)}^r}{r!}---\left(11\right)$

一般情况下n都较小，故T≈nt，从而

$L\left(\lambda \right)=\underset{r=0}{\overset{c}{\Sigma }}{e}^{-\mathrm{\lambda T}}\frac{{\left(\mathrm{\lambda T}\right)}^r}{r!}---\left(12\right)$

对于验证指标失效率λ，根据Bayes理论，取其共轭先验分布为Gamma分布，记为 Gamma(a,b),即：

$\pi \left(\lambda \right|a,b)=\frac{b^a}{\Gamma \left(a\right)}{\lambda }^{a-1}{e}^{-\mathrm{b\lambda }}---\left(13\right)$

其中，Γ(a)为Gamma函数，其定义为：

$\Gamma \left(a\right)={\int }_0^{\infty }{e}^{-x}{x}^{a-1}\mathrm{dx}---\left(14\right)$

根据验后风险准则和对指数型产品先验分布和接受概率的表达式的推导以及加速条件 λA=λU·AF，对于生产方后验风险α(c,T)的计算公式为：

$\begin{array}{c}\alpha (c,T)=P(\lambda \le \mathrm{AF}·{\lambda }_0|t\ge T)={\int }_0^{\mathrm{AF}·{\lambda }_0}p\left(\lambda \right|t\ge T)\mathrm{d\lambda }\\ =\frac{{\int }_0^{\mathrm{AF}·{\lambda }_0}P(t<T|\lambda )\pi \left(\lambda \right)\mathrm{d\lambda }}{{\int }_0^{\infty }P(t<T|\lambda )\pi \left(\lambda \right)\mathrm{d\lambda }}=\frac{{\int }_0^{\mathrm{AF}·{\lambda }_0}(1-\underset{r=0}{\overset{c}{\Sigma }}{e}^{-\mathrm{\lambda T}}\frac{{\left(\mathrm{\lambda T}\right)}^r}{r!})\pi \left(\lambda \right)\mathrm{d\lambda }}{{\int }_0^{\infty }(1-\underset{r=0}{\overset{c}{\Sigma }}{e}^{-\mathrm{\lambda T}}\frac{{\left(\mathrm{\lambda T}\right)}^r}{r!})\pi \left(\lambda \right)\mathrm{d\lambda }}\end{array}---\left(15\right)$

使用方后验风险β(c,T)的计算公式为：

$\begin{array}{c}\beta (r,T)=P(\lambda >\mathrm{AF}·{\lambda }_1|t\ge T)={\int }_{\mathrm{AF}·{\lambda }_1}^{\infty }p\left(\lambda \right|t\ge T)\mathrm{d\lambda }\\ =\frac{{\int }_{\mathrm{AF}·{\lambda }_1}^{\infty }P(t\ge T|\lambda )\pi \left(\lambda \right)\mathrm{d\lambda }}{{\int }_0^{\infty }P(t\ge T|\lambda )\pi \left(\lambda \right)\mathrm{d\lambda }}=\frac{{\int }_{\mathrm{AF}·{\lambda }_1}^{\infty }(1-\underset{r=0}{\overset{c}{\Sigma }}{e}^{-\mathrm{\lambda T}}\frac{{\left(\mathrm{\lambda T}\right)}^r}{r!})\pi \left(\lambda \right)\mathrm{d\lambda }}{{\int }_0^{\infty }(1-\underset{r=0}{\overset{c}{\Sigma }}{e}^{-\mathrm{\lambda T}}\frac{{\left(\mathrm{\lambda T}\right)}^r}{r!})\pi \left(\lambda \right)\mathrm{d\lambda }}\end{array}---\left(16\right)$

为了获得验证试验方案，需要解下列满足双方风险的约束条件，流程图如图4：

α(c,T)=P(λ≤AF·λ0|t≥T)≥1-α              (17)

β(c,T)=P(λ>AF·λ1|t≥T)≤β                (18)

步骤四、基于历史数据给出验证指标参数的先验分布，基于马尔科夫蒙特卡洛(Markov  Chain Monte Carlo,MCMC)方法，利用WinBUGS14计算满足条件的方案集合；

贝叶斯试验设计一般很难获得后验分布的数学表达式，为了解决这个问题，本发明使用 基于马尔科夫蒙特卡洛(Markov Chain Monte Carlo,MCMC)方法，并利用WinBUGS14计算满 足条件的方案。

具体的，以指数分布型产品为例，由步骤三可知，指数分布型产品对于验证指标失效率 λ，根据Bayes理论，一般取其共轭先验分布即Gamma分布，记为Gamma(a,b)，分别定义参 数a和b均服从分层先验分布Gamma(η,κ).

为了描述加速因子的不确定性，将加速因子的波动以概率分布的形式进行描述，为了简 单起见，定义加速因子服从的概率分布为均匀分布U(A,B).

则由MCMC方法并利用WinBUGS14可得到参数的后验均值，假设蒙特卡洛模拟的次数 为N次，λ的N次后验抽样计为λ(j)，先验分布两参数a，b的N次后验抽样分别计为(a(j),b(j))， 加速因子AF的N次抽样计为AF(j),令

$g^{\left(j\right)}\left(r\right)=\frac{{\left(b^{\left(j\right)}\right)}^{a^{\left(j\right)}}{T}^r}{r!\Gamma \left(a^{\left(j\right)}\right){(T+b^{\left(j\right)})}^{r+a^{\left(j\right)}}}---\left(19\right)$

则生产方的后验风险(15)可写为：

$\alpha (c,T)=\frac{\underset{j=1}{\overset{N}{\Sigma }}[1-\underset{r=0}{\overset{c}{\Sigma }}{e}^{-{\lambda }^{\left(j\right)}T}\frac{{\left({\lambda }^{\left(j\right)}T\right)}^r}{r!}]I({\lambda }^{\left(j\right)}\le {\lambda }_0·{\mathrm{AF}}^{\left(j\right)})}{\underset{j=1}{\overset{N}{\Sigma }}[1-\underset{r=0}{\overset{c}{\Sigma }}{e}^{-{\lambda }^{\left(j\right)}T}\frac{{\left({\lambda }^{\left(j\right)}T\right)}^r}{r!}]}$

$=\frac{\underset{j=1}{\overset{N}{\Sigma }}\left[\underset{r=0}{\overset{c}{\Sigma }}{g}^{\left(j\right)}\left(r\right)\gamma (r+a^{\left(j\right)},(T+b^{\left(j\right)}){\lambda }_0{\mathrm{AF}}^{\left(j\right)})\right]}{\underset{j=1}{\overset{N}{\Sigma }}\left[\underset{r=0}{\overset{c}{\Sigma }}{g}^{\left(j\right)}\left(r\right)\Gamma (r+a^{\left(j\right)})\right]}---\left(20\right)$

使用方的后验风险为：

$\begin{array}{c}\beta (r,T)=\frac{\underset{j=1}{\overset{N}{\Sigma }}[1-\underset{r=0}{\overset{c}{\Sigma }}{e}^{-{\lambda }^{\left(j\right)}T}\frac{{\left({\lambda }^{\left(j\right)}T\right)}^r}{r!}]I({\lambda }^{\left(j\right)}\ge {\lambda }_1·{\mathrm{AF}}^{\left(j\right)})}{\underset{j=1}{\overset{N}{\Sigma }}[1-\underset{r=0}{\overset{c}{\Sigma }}{e}^{-{\lambda }^{\left(j\right)}}\frac{{\left({\lambda }^{\left(j\right)}T\right)}^r}{r!}]}\\ =\frac{\underset{j=1}{\overset{N}{\Sigma }}\left[\underset{r=0}{\overset{c}{\Sigma }}{g}^{\left(j\right)}\left(r\right)(\Gamma (r+a^{\left(j\right)})-\gamma (r+a^{\left(j\right)},(T+b^{\left(j\right)}){\lambda }_1{\mathrm{AF}}^{\left(j\right)}))\right]}{\underset{j=1}{\overset{N}{\Sigma }}\left[\underset{r=0}{\overset{c}{\Sigma }}{g}^{\left(j\right)}\left(r\right)\Gamma (r+a^{\left(j\right)})\right]}\end{array}---\left(21\right)$

式中的是不完全gamma函数。这时候求解关于双方风险的约束条件
(17)和(18)，可以得到满足约束条件的方案集合(c,T)。

其中，具体的抽样模拟过程采用下面的步骤：

子步骤1.在参数λ对应的先验分布Gamma分布中抽取i次参数λi，在抽样分布f(x|λi) 中抽取j次仿真失效数据xij。

子步骤2.λi对应的Gamma分布的两参数分别为ai和bi，结合分层先验分布Gamma(η,κ) 和仿真失效数据xij，利用Winbugs软件和MCMC方法可得到参数ai和bi的后验平均值 E(π(ai|xij))和E(π(bi|xij))。

子步骤3.设定蒙特卡洛模拟次数，然后重复子步骤1和子步骤2仿真得到足够的数据。

子步骤4.利用仿真数据求解(20)和(21)方程组，得到满足条件的方案集合(c,T)。

步骤五、确定试验的费用约束，计算得到最优方案。

由(17)和(18)可获得满足双方风险条件的试验方案(c,T)集合，以及实际的双方风险，在考 虑试验费用的情况下，我们可以对方案进行进一步优化，考虑优化参数为c,T，α，β，假设 下列参数的定义如下：

a1：与试验时间相关的试验费用，包括试验过程中的电力，物力及人力损耗，试验时间 越长损耗越高。

a21：与试验样品有关的试验费用，失效的样品越多则损耗越高。

a22：经过试验但未失效的样本造成的试验损失，这些样本虽未失效，但是已经过了一定 时间的使用损耗，不具有最初的性能。

a3：与生产方风险有关的试验费用，包括产品的重新设计与生产延误、销售市场的萎缩 等。

a4：与使用方风险有关的试验费用,包括产品使用过程当中的维护保障、任务延迟损失等。

由以上参数可得考虑产品试验费用为：

f(c,T,α,β)=a1·T+a12·c+a22·(n-c)+a3·α+a4·β     (22)

优化目标即使上述试验费用最小。

实施例：

以某寿命分布服从指数分布的电子产品为例，采用本发明提出的基于贝叶斯理论的定时 截尾加速验收抽样试验优化设计方案进行验收试验，应用步骤和方法如下：

步骤一、确定产品的寿命分布，寿命验证指标参数，然后建立统计假设；

此电子产品寿命服从指数分布，假设根据试验资料以及历史数据等先验信息可确定试验 的加速因子大概为AF=3，选取失效率λA作为其寿命验证指标参数，根据协定的双方风险 α=0.2，β=0.2，在产品的抽样特性曲线(OC Curve)上选择对应的检验上下限λ0=0.004和 λ1=0.005，然后建立统计假设如下：

H0：λA≤λ0·AF H1：λA>λ1·AF

步骤二、对选定的产品设计试验决策法则，确定其试验和接收流程；

使用恒定应力加速试验(CSALT)，寿命试验的决策法则即为：

（1）在产品批中选择n=30个样品进行CSALT，加速因子为AF=3，试验为无替换试验；

（2）试验进行到试验累计时间到达预定值T时停止试验，记录试验过程中的失效数；

（3）设在试验过程中出现了r故障，如果r≤c，则认为产品批合格，接收；如果r>c， 则认为产品不合格，拒收批产品。

步骤三、基于验后风险准则推导满足双方风险关于验证指标参数的约束条件；

基于验后风险准则，指数分布型产品满足双方风险关于验证指标参数的约束条件即为式 (17)和(18)。

步骤四、基于历史数据给出验证指标参数的先验分布，基于马尔科夫蒙特卡洛(Markov  Chain Monte Carlo,MCMC)方法，利用WinBUGS14计算满足条件的方案集合；

根据历史数据等信息，确定验证指标参数λ服从的先验分布为λ～Gamma(1,1000)，确定 λ的先验分布伽马分布Gamma(a,b)中两参数a和b的分层先验分布分别为a～Gamma(1,1), b～Gamma(1,1000)，确定加速因子的概率分布为均值为3的平均分布U(2.5,3.5)。

在λ的先验分布中抽取i=100次λ的值，然后分别在抽样分布为指数分布f(x|λi)中获取 j=100次仿真失效数据xij，由WinBUGS14软件，输入选定的产品寿命的先验分布，可得到 100个参数a和b的后验分布平均值E(π(ai|xij))和E(π(bi|xij))。选定N=100，分别利用a和b 的后验平均值和加速因子的概率分布U(2.5,3.5)抽样100次带入式(20)和(21)，运算得到满足 条件的方案结果部分列表如下：

表1满足双方风险条件的试验方案集合

试验时间T(h) 试验接收数c 生产方实际风险α 使用方实际风险β

35 0 0.1983 0.1631

69 1 0.1997 0.1592

105 2 0.1996 0.1551

141 3 0.1990 0.1534

177 4 0.1996 0.1527

214 5 0.1994 0.1514

步骤五、确定试验的费用约束，计算得到最优方案。

由表格可发现生产方实际风险基本无变化，而使用方实际风险随着实际的试验时间变大 而变小。所以在比较试验费用时可以忽略生产方实际风险的影响，假设根据实际情况，取 a1=100元，a21=3000元，a22=1000元，a3=3000000元，则

试验费用计算如下：

表2试验方案费用列表

试验时间T(h) 试验接收数c 生产方实际风险α 使用方实际风险β 试验费用(元)

35 0 0.1983 0.1631 522800

69 1 0.1997 0.1592 516500

105 2 0.1996 0.1551 509800

141 3 0.1990 0.1534 510300

177 4 0.1996 0.1527 513800

214 5 0.1994 0.1514 515600

由此可得到最优方案应该为(c,T)=(2,105)。