CN102411537A - 一种基于混合贝叶斯先验分布的可靠性验证测试方法 - Google Patents

Abstract

本发明公开了一种基于混合贝叶斯先验分布的可靠性验证测试方法，采用共轭先验分布法确定未知参数的先验分布，通过先验矩方法和最大熵方法分别求出两组不同参数，即得到不同的先验分布，再根据第二类极大似然方法确定以上两种先验分布的权重，按权重将先验矩方法和最大熵方法得出的参数融合，进而得到最终的先验分布比单纯使用其中一种方法得到的先验分布更加准确，与真实的分布拟合度更好。通过最终的贝叶斯先验分布及试验信息，计算出可靠性验证测试中所需的最小无失效用例数，这种方法相对于无先验知识的情况，可以有效的减少测试用例量。

Description

一种基于混合贝叶斯先验分布的可靠性验证测试方法

技术领域

本发明涉及软件测试领域，特别涉及一种基于混合贝叶斯先验分布的可靠性验证测试方法。

背景技术

随着计算机在民用和军用产品上应用的日益增多，软件缺陷引发的产品故障以及灾难性事故也越来越严重。软件可靠性是整个系统可靠性的重要保证，甚至是决定性的因素。软件可靠性测试是软件可靠性保证过程中非常重要的一步。经过可靠性测试的软件并不能确保该软件中剩余的错误数量最小，但是可以确保该软件达到较高的可靠性要求。从工程上来看，软件的可靠性高即意味着该软件的失效率低，又意味着一旦该软件发生失效，由此造成的危害也较小。软件可靠性测试的目的就是要保证软件中没有对可靠性影响较大的错误。

学术界和工业界一直努力在试图解决，如何准确、客观、高可信地验证软件是否达到了可靠性要求的问题，经过几十年的发展，形成了多种软件可靠性验证测试方法。利用基于经典统计理论的可靠性验证测试方法进行可靠性测评，从理论上讲，其评价结果的可信度是很高的，但是，随着软件规模、软件复杂度的日益增大，以及对可靠性要求的不断提高，由于经典统计理论的可靠性验证测试方法没有考虑已经存在的先验信息，因此所需要的测试用例数量非常大，即造成测试用例开销大，又导致测试的持续期长，从而使得基于传统方法的软件验证测试工作效率极低。然而在实际的工程实践中，为了保证高可信软件的可靠性，开发方既要进行可靠性设计，又要在开发过程中对软件进行严格的可靠性增长测试。因此，当对软件进行可靠性验证测试时，软件已经具备了较高的可靠性，可靠性增长测试过程是可以当成先验知识加以利用的。因此利用先验贝叶斯统计原理，从统计学的角度推断如何减少验证测试用例数，缩短测试时间，降低测试开销，是解决高可靠性软件验证测试的有效途径。

利用贝叶斯方法进行统计推断的最基本问题便是如何确定统计量的先验分布，经过多年的研究，已出现多种先验分布求解方法，包括共轭先验分布法、最大熵原则、林德莱原则、Jeffreys原则、最大数据信息原则、不变测度等方法，这些方法中，各有理论优点，需要结合具体应用选取相应的方法。其中，共轭先验分布具有良好的数学表达和很广泛的应用基础。

贝叶斯统计学中先验分布中所含的未知参数称为超参数。关于超参数的计算经对现有文献检索发现，文献：覃志东等《安全关键软件可靠性验证测试方法研究》，介绍了先验矩方法求解软件失效概率的概率密度函数的先验分布。文献：杜小翔，钱红兵《基于先验贝叶斯推断验证法的高可靠性评估》也采用了先验矩方法求解超参数。文献：Savchuk Vladimir P.Bayes Reliability Estimation Using Multiple Sources of Prior Information：BinmialSampling，介绍了求解超参数的新方法，提出了如何运用最大熵的方法来确定二项式分布的情况下的先验分布。文献：Yoon Won Hyo.Systematic Bayes Prior-Assignment byCoupling the Mini-Max Entropy and Moment-Matching Methods，研究了用正态分布最大熵方法来确定先验分布的超参数值。文献：詹吴可，姜礼平《共轭最大熵先验下的贝叶斯估计》，进一步扩充了最大熵方法，研究了指数的情况下，最大熵先验下的贝叶斯估计方法。但是这些方法只是提供了求先验分布的手段缺少对先验分布的合理性和准确性的研究，即使在相同的共轭先验分布簇中，不同的计算方法得到的参数值也可能不同。在文献《软件可靠性验证测试中降低测试用例量方法研究》中，对比了两种不同方法求解先验分布参数，在特定的实验条件下得出基于共轭最大熵方法的先验分布比基于先验矩方法的先验分布的可信度更高的结论。但对于其他的先验数据，基于共轭最大熵方法的超参数计算并不一定优于基于先验矩的计算方法。

发明内容

本发明的目的是为了解决上述问题，提供一种基于混合贝叶斯先验分布的可靠性验证测试方法。本发明采用共轭先验分布法确定未知参数的先验分布，通过先验矩方法和最大熵方法分别求出两组不同参数，即得到不同的先验分布，再根据第二类极大似然方法确定以上两种先验分布的置信因子，将置信因子看成对应方法的权重，按权重将先验矩方法和最大熵方法得出的参数融合，进而得到最终的先验分布比单纯使用其中一种方法得到的先验分布更加准确，与真实的分布拟合度更好。通过最终的先验分布和试验信息得到一个后验分布，利用此后验分布，得到可靠性验证测试中的最小用例数。

本发明的一种基于混合贝叶斯先验分布的可靠性验证测试方法，包含以下几个步骤：

步骤一、采用共轭分布确定未知参数的先验分布。

步骤二、根据可靠性先验信息，利用先验矩方法确定先验分布中的超参数。

步骤三、根据可靠性先验信息，利用最大熵方法确定先验分布中的超参数。

步骤四、采用第二类极大似然方法确定步骤二和步骤三中得出的两种的先验分布的权重。

步骤五、按不同的权重将两种方法估计出的超参数值进行融合，确定最终先验分布的超参数。

步骤六、根据贝叶斯公式考虑先验分布信息和试验信息，得到一个后验分布，根据要求的置信度和失效概率，利用后验分布，计算出可靠性验证测试中的最小用例数。

本发明的优点在于：

本发明充分利用了贝叶斯先验特性，采用共轭先验分布法确定未知参数的先验分布，通过先验矩方法和最大熵方法分别求出两组不同参数，根据第二类极大似然方法确定以上两种先验分布的置信因子，将置信因子看成对应方法的权重，按权重将先验矩方法和最大熵方法得出的参数融合，进而得到最终的先验分布更合理，比单纯使用其中一种方法得到的先验分布更加准确，与真实的分布拟合度更好。通过最终的先验分布和试验信息得到一个后验分布，利用后验分布，得到可靠性验证测试中的最小用例数，最终达到不影响可靠度的前提下降低测试量的作用。

附图说明

图1为本发明的方法流程图；

图2为本发明实施例中先验分布对比图。

具体实施方式

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

本发明的一种基于混合贝叶斯先验分布的可靠性验证测试方法，流程如图1所示，具体包括以下几个步骤：

步骤一、采用共轭先验分布的方法，确定未知失效概率p的先验分布；

共轭先验分布具有良好的数学表达和很广泛的应用基础，因此，本发明采用共轭分布方法确定先验分布的形式。

假定某软件在任一运行时刻、任意选择输入的失效概率为p，且每次输入操作执行是均满足贝努利实验的独立统计性，则在n次执行中失效r次的概率为：

$p\left(r\right)=C_n^r{p}^r{(1-p)}^{n-r}$

因为二项式分布的共轭先验分布服从beta分布，所以失效概率p的概率密度函数的先验分布为：

f(p)＝p^a-1(1-p)^b-1/B(a，b)    (1.0)

其中：a＞0，b＞0，为beta函数的超参数，且 $B(a,b)={\∫}_0^1{p}^{a-1}{(1-p)}^{b-1}\mathrm{dp}.$

步骤二、利用先验矩方法确定先验分布中的超参数；

首先收集测试增长阶段的测试信息，假设选取最后的m组测试信息作为先验信息，每组测试中含有n个测试用例，其中m组测试信息中导致失效的测试用例数的个数分别为k₁，k₂，...，k_m，k₁，k₂，...，k_m组成样本x，则样本x的边缘分布为：

$h\left(x\right)={\∫}_0^{1}f\left(p\right)\mathrm{\π}\left(x\right|p)\mathrm{dp}---\left(2.0\right)$

其中：π(x|p)为样本x对p条件分布，

将公式(1.0)带入公式(2.0)得出公式(2.1)：

$h\left(x\right)={\∫}_0^1\frac{1}{B(a,b)}{p}^{a-1}{(1-p)}^{b-1}{C}_n^x{p}^x{(1-p)}^{n-x}\mathrm{dp}={\∫}_0^1\frac{1}{B(a,b)}{C}_n^x{p}^{a+x-1}{(1-p)}^{b+n-x-1}\mathrm{dp}$

$=\frac{C_n^x}{B(a,b)}{\∫}_0^1{p}^{a+x-1}{(1-p)}^{b+n-x-1}\mathrm{dp}$

因为 $B(a,b)={\∫}_0^1{p}^{a-1}{(1-p)}^{b-1}\mathrm{dp},$ 所以 ${\∫}_0^1{p}^{a+x-1}{(1-p)}^{b+n-x-1}\mathrm{dp}=B(a+x,n+b-x),$ 所以得到：

$h\left(x\right)=C_n^{x}B(a+x,n+b-x)/B(a,b)---\left(2.1\right)$

则h(x)的一阶矩，二阶矩具体为：

$E\left(x\right)=\∫\mathrm{xh}\left(x\right)\mathrm{dx}=\underset{x=0}{\overset{n}{\mathrm{\Σ}}}\mathrm{xh}\left(x\right)={\∫}_0^{1}f\left(p\right)\underset{x=0}{\overset{n}{\mathrm{\Σ}}}{C}_n^x{\mathrm{xp}}^x{(1-p)}^{n-x}\mathrm{dp}=\mathrm{na}/(a+b)---\left(2.2\right)$

$E\left[(x-1)x\right]=\∫x(x-1)h\left(x\right)\mathrm{dx}=\underset{x=0}{\overset{n}{\mathrm{\Σ}}}x(x-1)h\left(x\right)$

$={\∫}_0^{1}f\left(p\right)\underset{x=0}{\overset{n}{\mathrm{\Σ}}}{C}_n^{x}x(x-1)p^x{(1-p)}^{n-x}\mathrm{dp}=\frac{n(n-1)a(a+1)}{(a+b)(a+b+1)}---\left(2.3\right)$

通过公式(2.2)(2.3)获取a，b的值如下：

a＝Ex(Ex²-nEx)/[(n-1)(Ex)²+n(Ex-Ex²)]        (2.4)

b＝(n-Ex)(Ex²-nEx)/[(n-1)(Ex)²+n(Ex-Ex²)]    (2.5)

其中：Ex为h(x)的一阶矩，用样本均值

来估计超参数a和b，Ex²为h(x)的二阶矩，用来估计超参数a和b，最后，通过先验矩方法得到的超参数a和b的值分别记为a₁和b₁。

步骤三、利用最大熵方法确定先验分布中的超参数；

一个分布的熵越大，就表示该分布的不确定性越大。因此要在先验信息的约束之下，尽可能少的引入主观成分，可以将先验信息看成约束条件，通过最大化熵，以确定约束条件下的先验分布。

将失效概率p的先验信息用f(p)如式(3.0)所述的约束方式表示：

$E\left[g_k\left(p\right)\right]={\∫}_0^{1}f\left(p\right)g_k\left(p\right)\mathrm{dp}={\mathrm{\μ}}_k---\left(3.0\right)$

其中：E[g_k(p)]为f(p)的k阶原点矩，即g_k(p)＝p^k，将E[g_k(p)]记为μ_k。本发明中E[p^k]采用f(p)的二阶原点矩，因此得到共轭最大熵求解公式如公式(3.1)所示。

$\left\{\begin{array}{c}\max H\left(p\right)=-{\∫}_0^{1}f\left(p\right)\ln \left(\frac{f\left(p\right)}{f_0\left(p\right)}\mathrm{dp}\right)\\ s.t.E\left[g_k\left(p\right)\right]={\∫}_0^{1}f\left(p\right)g_k\left(p\right)\mathrm{dp}={\mathrm{\μ}}_k\end{array}\right\}---\left(3.1\right)$

其中：H(p)表示p的熵，f₀(p)为问题自然的不变的无信息先验分布，取值为1。则将公式(1.0)代入公式(3.1)，得到公式(3.2)：

$H\left(p\right)=-\frac{1}{B(a,b)}{\∫}_0^1{p}^{a-1}{(1-p)}^{b-1}[\ln p^{a-1}{(1-p)}^{b-1}-\ln B(a,b)]\mathrm{dp}$

$=\frac{\ln B(a,b)}{B(a,b)}-\frac{1}{B(a,b)}{\∫}_0^1{p}^{a-1}{(1-p)}^{b-1}\ln p^{a-1}{(1-p)}^{b-1}\mathrm{dp}---\left(3.2\right)$

取k＝2，则：

$E\left[g_k\left(p\right)\right]=\frac{1}{B(a,b)}{\∫}_0^1{p}^2{p}^{a-1}{(1-p)}^{b-1}\mathrm{dp}$

$=\frac{B(a+2,b)}{B(a,b)}=\frac{a(a+1)}{(a+b)(a+b+1)}={\mathrm{\μ}}_2---\left(3.3\right)$

公式(3.1)为条件极值问题，利用拉格朗日乘数法将条件极值转换为无条件极值。具体为：

令：

获取

的无条件极值，则。

式中：λ是拉格朗日求极值过程的一个参数，F_a(a，b，λ)表示F(a，b，λ)对a求偏导，F_b(a，b，λ)表示F(a，b，λ)对b求偏导，F_λ(a，b，λ)表示F(a，b，λ)对λ求偏导，h_a表示对a求偏导、h_b表示对b求偏导。根据方程组(3.5)获取得到a和b的值，分别记为a₂和b₂。

步骤四、采用第二类极大似然方法确定步骤二和步骤三中得出的两种的先验分布的权重；

第二类似然估计方法是把样本信息看成是由边缘分布产生的，依据在不同的先验分布中样本的似然值的大小，判定不同先验分布的可信度。似然值越大，可信度越高，置信因子也就越大，对应的先验分布就越真实。因此可以依据置信因子的大小确定参数的权重。具体步骤如下：

步骤4.1：将先验矩方法得出的分布函数记为π₁(p)， ${\mathrm{\π}}_1\left(p\right)=p^{a_1-1}{(1-p)}^{b_1-1}/B(a_1,b_1),$ 最大熵方法得出的先验分布函数记为π₂(p)， ${\mathrm{\π}}_2\left(p\right)=p^{a_2-1}{(1-p)}^{b_2-1}/B(a_2,b_2),$ 根据公式(4.0)和(4.1)分别取出确定这两种不同方法得出的先验分布的似然函数值。

$m\left(x\right|{\mathrm{\π}}_k)={\∫}_0^{1}f\left(X\right|p){\mathrm{\π}}_k\left(p\right)\mathrm{dp},k=\mathrm{1,2}---\left(4.0\right)$

$L\left(X\right|{\mathrm{\π}}_k)=\underset{i=1}{\overset{n}{\mathrm{\Π}}}m\left(x_i\right|{\mathrm{\π}}_k),k=\mathrm{1,2}---\left(4.1\right)$

其中：m(x|π_k)为先验分布的边缘分布，L(X |π_k)为似然函数。

步骤4.2：根据公式(4.2)获取两种分布的置信因子，并将置信因子作为超参数的权重。

${\mathrm{\ϵ}}_k=\frac{L\left(X\right|{\mathrm{\π}}_k)}{\underset{k=1}{\overset{2}{\mathrm{\Σ}}}L\left(X\right|{\mathrm{\π}}_k)},k=\mathrm{1,2}---\left(4.2\right)$

其中：ε₁为先验矩方法取得的分布的置信因子，ε₂为最大熵方法取得的分布的置信因子。

步骤五、按不同的权重将两种方法估计出的超参数值进行融合，确定最终先验分布的超参数。

a＝ε₁a₁+ε₂a₂

                                    (5.0)

b＝ε₁b₁+ε₂b₂

步骤六、确定可靠性验证测试无失效最小测试用例数；

步骤6.1：通过上述步骤，得到了未知失效概率p的先验分布，确定第一轮可靠性验证测试的无失效最小测试用例数n₁，将公式(5.0)求出的a，b值代入公式(6.0)，得到n₁：

$p(p<p_0)={\&Integral;}_0^{p_0}f\left(p\right|0,n_1,a,b)=\mathrm{dp}={\&Integral;}_0^{p_0}\frac{p^{a-1}{(1-p)}^{b+n_1-1}\mathrm{dp}}{B(a,b+n_1)}\&GreaterEqual;c---\left(6.0\right)$

其中：(p₀，c)为要求的可靠性指标已知，p₀为失效率指标，c为置信度指标，f(p|0，n，a，b)为执行了n个测试用例发生0个失效的概率分布。

步骤6.2：若根据步骤6.1第一轮得到的n₁个测试用例全部无失效通过测试，则说明该软件符合给定的可靠性要求，验证结束；否则，若第一轮执行到第t₁(t₁≤n₁)个测试用例发生失效，说明不符合验收指标，排除故障后需进行第二轮可靠性验证测试，转入步骤6.3。

步骤6.3：把第一轮可靠性验证测试执行通过的测试用例数(t₁-1)和失效的用例数1，作为先验信息融入到先验分布(式(1.0))中，得到新的概率分布f(p|1，t₁，a，b)，见公式(6.1)所示。再根据公式(6.2)确定第二轮测试所需最小的无失效用例数n₂。

$f\left(p\right|1,t_1,a,b)=\frac{p^a{(1-p)}^{b+t_1-2}}{B(a+1,b+t_1-1)}---\left(6.1\right)$

$p(p<p_0)={\&Integral;}_0^{p_0}f\left(p\right|1,t_1+n_2,a,b)\mathrm{dp}={\&Integral;}_0^{p_0}\frac{p^a{(1-p)}^{b+t_1+n_2-2}\mathrm{dp}}{B(a+1,b+t_1+n_2-1)}\&GreaterEqual;c---\left(6.2\right)$

其中：f(p|1，t₁+n₂，a，b)为在第1轮的基础上，无失效执行n₂个测试用例的概率分布，公式(6.2)求解较复杂，本发明通过MATLAB仿真计算出n₂的值。

若n₂个无失效用例全部无失效通过测试，则说明该软件符合给定的可靠性要求，验证结束；否则，若第二轮测试执行到第t₂个测试用例发生失效，说明不符合验收指标，排除故障后需进行第下一轮可靠性验证测试，转入步骤6.4。

步骤6.4：以此类推，若可靠性验证测试进行了i轮，各轮失效分别发生在第t₁，t₁+t₂，，...t，₁+t₂+t₃+...+_i个测试用例上，则由公式(6.3)获取下一轮测试需要的无失效的最少用例数n_i+1。

${\&Integral;}_0^{p_0}\frac{p^{a+i-1}{(1-p)}^{b+\underset{1}{\overset{i}{\mathrm{\&Sigma;}}}{t}_i-i+n_{i+1}-1}\mathrm{dp}}{B(a+i,b+\underset{1}{\overset{i}{\mathrm{\&Sigma;}}}{t}_i-i+n_{i+1})}\&GreaterEqual;c---\left(6.3\right)$

其中：

为i轮测试共执行的测试用例数，为i轮测试共通过的测试用例数，令

则公式(6.3)变为：

${\&Integral;}_0^{p_0}\frac{p^{a+i-1}{(1-p)}^{b+N_{i+1}-i-1}\mathrm{dp}}{B(a+i,b+N_{i+1}-i)}\&GreaterEqual;c---\left(6.4\right)$

利用公式(6.4)，获取累积的总测试用例量N_i+1，然后再根据

得到具体的第i+1轮可靠性验证测试所需要的无失效测试用例量n_i+1。根据第i+1轮可靠性验证测试所需要的无失效测试用例量n_i+1进行软件可靠性测试。

实施例：

下面本发明实施例中，结合MATLAB程序仿真实现本发明方法的验证。

通过MATLAB仿真程序由beta(1，22)产生20个随机数，分别为：

x1＝0.0587，x2＝0.0027，x3＝0.0774，x4＝0.0958，x5＝0.0134，x6＝0.0635，x7＝0.0921，x8＝0.0016，x9＝0.0036，x10＝0.0552，x11＝0.0221，x12＝0.0539，x13＝0.0614，x14＝0.0448，x15＝0.0882，x16＝0.0681，x17＝0.0852，x18＝0.0051，x19＝0.0293，x20＝0.0249。

设某软件在某次可靠性测试中得到的20组失效概率值为X＝(x1，x2，x3，x4，x5，x6，x7，x8，x9，x10，x11，x12，x13，x14，x15，x16，x17，x18，x19，x20)。

基于以上先验信息，利用先验矩方法根据公式(2.4)和(2.5)计算出a₁＝0.8，b₁＝24.1利用最大熵方法求解方程组(3.4)得出a₂＝1.2，b₂＝23.3，再根据本发明中给出的权重的确定方法，得到混合后的参数a＝1.0，b＝23.6。

图2中，Line1为beta(1，22)，产生数据的真实分布，Line2为本发明所述的先验矩和最大熵混合方法得出的先验分布beta(1.0，23.6)，Line3为最大熵方法得出的先验分布beta(1.2，23.3)，Line4为先验矩方法得出的先验分布beta(0.8，24.1)。从图2的对比中可以看出混合方法比任一种方法更接近真实分布。

在不同的可靠性指标约束下，使用无先验知识的贝叶斯方法时，a＝1，b＝1，根据公式(6.0)计算出容忍失效数r为0时所需的最少用例数，结果见表1。

表1不同可靠性指标下无先验知识所需要的最小无失效用例量

	p＝0.1	p＝0.08	p＝0.06	p＝0.04	p＝0.02	p＝0.001
							c＝0.90	22	27	37	56	113	2301
c＝0.92	24	30	40	61	125	2524
							c＝0.94	27	33	45	68	139	2812
c＝0.96	31	38	52	78	159	3217
							c＝0.98	38	46	63	95	193	3910
c＝0.99	44	55	74	112	227	4602
							c＝0.999	66	82	111	169	341	6904

同理，使用先验知识，依据本发明提出的混合方法得到的a＝1.0，b＝23.6，计算出容忍失效数为零时所需的最小用例数，结果见表2。

表2不同可靠性指标下利用先验知识所需要的最小无失效用例量

	p＝0.1	p＝0.08	p＝0.06	p＝0.04	p＝0.02	p＝0.001
							c＝0.90		4	14	33	90	2278
c＝0.92	1	7	17	38	102	2501
							c＝0.94	4	10	22	45	117	2789
c＝0.96	7	15	28	55	136	3194
							c＝0.98	15	23	40	72	170	3887
c＝0.99	21	32	51	89	204	4579
							c＝0.999	43	59	88	146	317	6881

通过表1和表2的结果可以看出，基于本发明所采用的参数融合方法，得到的贝叶斯先验分布，利用先验知识的可靠性验证测试所需的最少无失效用例数相对于无先验知识的最小无失效用例数要少。

上述实施例说明了本发明采用基于先验矩和最大熵混合方法确定先验分布的超参数，得到的先验分布更接近真实的分布，并且相对于无先验知识的情况，可以有效的减少测试用例量。

Claims (1)

1.一种基于混合贝叶斯先验分布的可靠性验证测试方法，其特征在于，具体包括以下几个步骤：

步骤一、采用共轭先验分布的方法，确定未知失效概率p的先验分布；

假定某软件在任一运行时刻、任意选择输入的失效概率为p，且每次输入操作执行是均满足贝努利实验的独立统计性，则在n次执行中失效r次的概率为：

$p\left(r\right)=C_n^r{p}^r{(1-p)}^{n-r}$

因为二项式分布的共轭先验分布服从beta分布，所以失效概率p的概率密度函数的先验分布为：

f(p)＝p^a-1(1-p)^b-1/B(a，b)        (1.0)

其中：a＞0，b＞0，为beta函数的超参数，且 $B(a,b)={\&Integral;}_0^1{p}^{a-1}{(1-p)}^{b-1}\mathrm{dp};$

步骤二、利用先验矩方法确定共轭先验分布中的超参数；

首先收集测试增长阶段的测试信息，假设选取最后的m组测试信息作为先验信息，每组测试中含有n个测试用例，其中m组测试信息中导致失效的测试用例数的个数分别为k₁，k₂，…，k_m，k₁，k₂，…，k_m组成样本x，则样本x的边缘分布为：

$h\left(x\right)={\&Integral;}_0^{1}f\left(p\right)\mathrm{\&pi;}\left(x\right|p)\mathrm{dp}---\left(2.0\right)$

其中：π(x|p)为样本x对p条件分布，

将公式(1.0)带入公式(2.0)得出公式(2.1)：

$h\left(x\right)={\&Integral;}_0^1\frac{1}{B(a,b)}{p}^{a-1}{(1-p)}^{b-1}{C}_n^x{p}^x{(1-p)}^{n-x}\mathrm{dp}={\&Integral;}_0^1\frac{1}{B(a,b)}{C}_n^x{p}^{a+x-1}{(1-p)}^{b+n-x-1}\mathrm{dp}$

$=\frac{C_n^x}{B(a,b)}{\&Integral;}_0^1{p}^{a+x-1}{(1-p)}^{b+n-x-1}\mathrm{dp}$

因为 $B(a,b)={\&Integral;}_0^1{p}^{a-1}{(1-p)}^{b-1}\mathrm{dp},$ 所以 ${\&Integral;}_0^1{p}^{a+x-1}{(1-p)}^{b+n-x-1}\mathrm{dp}=B(a+x,n+b-x),$ 所以得到：

$h\left(x\right)=C_n^{x}B(a+x,n+b-x)/B(a,b)---\left(2.1\right)$

则h(x)的一阶矩，二阶矩具体为：

$E\left(x\right)=\&Integral;\mathrm{xh}\left(x\right)\mathrm{dx}=\underset{x=0}{\overset{n}{\mathrm{\&Sigma;}}}\mathrm{xh}\left(x\right)={\&Integral;}_0^{1}f\left(p\right)\underset{x=0}{\overset{n}{\mathrm{\&Sigma;}}}{C}_n^x{\mathrm{xp}}^x{(1-p)}^{n-x}\mathrm{dp}=\mathrm{na}/(a+b)---\left(2.2\right)$

$E\left[(x-1)x\right]=\&Integral;x(x-1)h\left(x\right)\mathrm{dx}=\underset{x=0}{\overset{n}{\mathrm{\&Sigma;}}}x(x-1)h\left(x\right)$

$={\&Integral;}_0^{1}f\left(p\right)\underset{x=0}{\overset{n}{\mathrm{\&Sigma;}}}{C}_n^{x}x(x-1)p^x{(1-p)}^{n-x}\mathrm{dp}=\frac{n(n-1)a(a+1)}{(a+b)(a+b+1)}---\left(2.3\right)$

通过公式(2.2)(2.3)获取a，b的值如下：

a＝Ex(Ex²-nEx)/[(n-1)(Ex)²+n(Ex-Ex²)]    (2.4)

b＝(n-Ex)(Ex²-nEx)/[(n-1)(Ex)²+n(Ex-Ex²)]    (2.5)

其中：Ex为h(x)的一阶矩，用样本均值

来估计超参数a和b，Ex²为h(x)的二阶矩，用

来估计超参数a和b，最后，通过先验矩方法得到的超参数a和b的值分别记为a₁和b₁；

步骤三、利用最大熵方法确定共轭先验分布中的超参数；

将失效概率p的先验信息用f(p)如式(3.0)所述的约束方式表示：

$E\left[g_k\left(p\right)\right]={\&Integral;}_0^{1}f\left(p\right)g_k\left(p\right)\mathrm{dp}={\mathrm{\&mu;}}_k---\left(3.0\right)$

其中：E[g_k(p)]为f(p)的k阶原点矩，即g_k(p)＝p^k，将E[g_k(p)]记为μ_k；本发明中E[p^k]采用f(p)的二阶原点矩，因此得到共轭最大熵求解公式如公式(3.1)所示；

$\left\{\begin{array}{c}\max H\left(p\right)=-{\&Integral;}_0^{1}f\left(p\right)\ln \left(\frac{f\left(p\right)}{f_0\left(p\right)}\mathrm{dp}\right)\\ s.t.E\left[g_k\left(p\right)\right]={\&Integral;}_0^{1}f\left(p\right)g_k\left(p\right)\mathrm{dp}={\mathrm{\&mu;}}_k\end{array}\right\}---\left(3.1\right)$

其中：H(p)表示p的熵，f₀(p)为问题自然的不变的无信息先验分布，取值为1；则将公式(1.0)代入公式(3.1)，得到公式(3.2)：

$H\left(p\right)=-\frac{1}{B(a,b)}{\&Integral;}_0^1{p}^{a-1}{(1-p)}^{b-1}[\ln p^{a-1}{(1-p)}^{b-1}-\ln B(a,b)]\mathrm{dp}$

$=\frac{\ln B(a,b)}{B(a,b)}-\frac{1}{B(a,b)}{\&Integral;}_0^1{p}^{a-1}{(1-p)}^{b-1}\ln p^{a-1}{(1-p)}^{b-1}\mathrm{dp}---\left(3.2\right)$

取k＝2，则：

$E\left[g_k\left(p\right)\right]=\frac{1}{B(a,b)}{\&Integral;}_0^1{p}^2{p}^{a-1}{(1-p)}^{b-1}\mathrm{dp}$

$=\frac{B(a+2,b)}{B(a,b)}=\frac{a(a+1)}{(a+b)(a+b+1)}={\mathrm{\&mu;}}_2---\left(3.3\right)$

公式(3.1)为条件极值问题，利用拉格朗日乘数法将条件极值转换为无条件极值；具体为：

令：

获取

的无条件极值，则；

式中：λ是拉格朗日求极值过程的一个参数，F_a(a，b，λ)表示F(a，b，λ)对a求偏导，F_b(a，b，λ)表示F(a，b，λ)对b求偏导，F_λ(a，b，λ)表示F(a，b，λ)对λ求偏导，h_a表示对a求偏导、h_b表示对b求偏导；根据方程组(3.5)获取得到a和b的值，分别记为a₂和b₂；

步骤四、采用第二类极大似然方法确定步骤二和步骤三中得出的两种先验分布的权重；

具体步骤如下：

步骤4.1：将先验矩方法得出的分布函数记为π₁(p)， ${\mathrm{\&pi;}}_1\left(p\right)=p^{a_1-1}{(1-p)}^{b_1-1}/B(a_1,b_1),$ 最大熵方法得出的先验分布函数记为π₂(p)， ${\mathrm{\&pi;}}_2\left(p\right)=p^{a_2-1}{(1-p)}^{b_2-1}/B(a_2,b_2),$ 根据公式(4.0)和(4.1)分别取出确定这两种不同方法得出的先验分布的似然函数值；

$m\left(x\right|{\mathrm{\&pi;}}_k)={\&Integral;}_0^{1}f\left(X\right|p){\mathrm{\&pi;}}_k\left(p\right)\mathrm{dp},k=\mathrm{1,2}---\left(4.0\right)$

$L\left(X\right|{\mathrm{\&pi;}}_k)=\underset{i=1}{\overset{n}{\mathrm{\&Pi;}}}m\left(x_i\right|{\mathrm{\&pi;}}_k),k=\mathrm{1,2}---\left(4.1\right)$

其中：m(x|π_k)为先验分布的边缘分布，L(X|π_k)为似然函数；

步骤4.2：根据公式(4.2)获取两种分布的置信因子，并将置信因子作为超参数的权重；

${\mathrm{\&epsiv;}}_k=\frac{L\left(X\right|{\mathrm{\&pi;}}_k)}{\underset{k=1}{\overset{2}{\mathrm{\&Sigma;}}}L\left(X\right|{\mathrm{\&pi;}}_k)},k=\mathrm{1,2}---\left(4.2\right)$

其中：ε₁为先验矩方法取得的分布的置信因子，ε₂为最大熵方法取得的分布的置信因子；

步骤五、按不同的权重将两种方法估计出的超参数值进行融合，确定最终先验分布的超参数；

a＝ε₁a₁+ε₂a₂

                                   (5.0)

b＝ε₁b₁+ε₂b₂

步骤六、确定可靠性验证测试无失效最小测试用例数；

步骤6.1：通过上述步骤，得到了未知失效概率p的先验分布，确定第一轮可靠性验证测试的无失效最小测试用例数n₁，将公式(5.0)求出的a，b值代入公式(6.0)，得到n₁：

$p(p<p_0)={\&Integral;}_0^{p_0}f\left(p\right|0,n_1,a,b)=\mathrm{dp}={\&Integral;}_0^{p_0}\frac{p^{a-1}{(1-p)}^{b+n_1-1}\mathrm{dp}}{B(a,b+n_1)}\&GreaterEqual;c---\left(6.0\right)$

其中：(p₀，c)为要求的可靠性指标已知，p₀为失效率指标，c为置信度指标，f(p|0，n，a，b)为执行了n个测试用例发生0个失效的概率分布；

步骤6.2：若根据步骤6.1第一轮得到的n₁个测试用例全部无失效通过测试，则说明该软件符合给定的可靠性要求，验证结束；否则，若第一轮执行到第t₁个测试用例发生失效，说明不符合验收指标，排除故障后需进行第二轮可靠性验证测试，转入步骤6.3；

步骤6.3：把第一轮可靠性验证测试执行通过的测试用例数(t₁-1)和失效的用例数1，作为先验信息融入到先验分布中，得到新的概率分布f(p|1，t₁，a，b)，见公式(6.1)所示；再根据公式(6.2)确定第二轮测试所需最小的无失效用例数n₂；

$f\left(p\right|1,t_1,a,b)=\frac{p^a{(1-p)}^{b+t_1-2}}{B(a+1,b+t_1-1)}---\left(6.1\right)$

$p(p<p_0)={\&Integral;}_0^{p_0}f\left(p\right|1,t_1+n_2,a,b)\mathrm{dp}={\&Integral;}_0^{p_0}\frac{p^a{(1-p)}^{b+t_1+n_2-2}\mathrm{dp}}{B(a+1,b+t_1+n_2-1)}\&GreaterEqual;c---\left(6.2\right)$

其中：f(p|1，t₁+n₂，a，b)为在第1轮的基础上，无失效执行n₂个测试用例的概率分布，仿真计算出n₂的值；

若n₂个无失效用例全部无失效通过测试，则说明该软件符合给定的可靠性要求，验证结束；否则，若第二轮测试执行到第t₂个测试用例发生失效，说明不符合验收指标，排除故障后需进行第下一轮可靠性验证测试，转入步骤6.4；

步骤6.4：以此类推，若可靠性验证测试进行了i轮，各轮失效分别发生在第t₁，t₁+t₂，，...t，₁+t₂+t₃+...+_i个测试用例上，则由公式(6.3)获取下一轮测试需要的无失效的最少用例数n_i+1；

${\&Integral;}_0^{p_0}\frac{p^{a+i-1}{(1-p)}^{b+\underset{1}{\overset{i}{\mathrm{\&Sigma;}}}{t}_i-i+n_{i+1}-1}\mathrm{dp}}{B(a+i,b+\underset{1}{\overset{i}{\mathrm{\&Sigma;}}}{t}_i-i+n_{i+1})}\&GreaterEqual;c---\left(6.3\right)$

其中：

为i轮测试共执行的测试用例数，

为i轮测试共通过的测试用例数，令

则公式(6.3)变为：

${\&Integral;}_0^{p_0}\frac{p^{a+i-1}{(1-p)}^{b+N_{i+1}-i-1}\mathrm{dp}}{B(a+i,b+N_{i+1}-i)}\&GreaterEqual;c---\left(6.4\right)$

利用公式(6.4)，获取累积的总测试用例量N_i+1，然后再根据

得到具体的第i+1轮可靠性验证测试所需要的无失效测试用例量n_i+1；根据第i+1轮可靠性验证测试所需要的无失效测试用例量n_i+1进行软件可靠性测试。

Images