CN104182377A - 一种基于β似然函数的参数估计方法 - Google Patents

Abstract

一种基于β似然函数的参数估计方法，该方法有四大步骤：步骤一、故障数据的收集；步骤二、平均秩次的计算；步骤三、寿命分布的选取和β似然函数的构造；步骤四、寿命分布参数的求解。本发明是根据可靠度非参数估计中的β分布法，计算各试件故障时产品可靠度服从的β分布，进而使用该分布的概率密度来度量各试件故障时可靠度估计值的合理程度，以各试件故障时可靠度估计值的合理程度之积构造β似然函数，将使β似然函数取值最大的分布参数作为估计结果。本发明在可靠性数据分析领域里具有广泛的适用性。

Description

一种基于β似然函数的参数估计方法

技术领域

本发明涉及一种基于β似然函数的参数估计方法，它是一种涉及概率分布的参数估计的方法，属于数理统计领域，该方法适用但不局限于可靠性数据分析领域。

背景技术

概率分布的参数估计是根据从总体中抽取的样本估计总体分布中包含的未知参数的方法，常见的参数估计方法有矩估计、最小二乘估计、极大似然估计等方法。这些方法各有其优点，但也都存在不足之处：矩估计和最小二乘估计是否可用取决于所选分布的数学形式，如不能通过取对数等操作化为线性结构的分布是不能使用最小二乘估计的；而极大似然估计考察的是各试件发生处的概率密度，当某个试件发生处所选分布的概率密度可为无穷大时极大似然估计无效。

在可靠性数据分析中，通过获取产品故障数据、选择寿命分布类型、进行参数估计后，最终关注的是产品可靠度随时间的变化情况，而在医学或农业中的所关注的存活率也和可靠度一样，都是由故障(死亡)的累积情况决定。但极大似然估计“使当前样本各个体故障(死亡)时的概率密度估计值之积最大”与“使当前样本各个体故障(死亡)时的可靠度(存活率)估计值最合理”并不是等价的，即从理论上讲在可靠度或存活率的研究中使用极大似然估计是存在问题的。

发明内容

本发明的目的是提供一种基于β似然函数的参数估计方法，具体是根据可靠度非参数估计中的β分布法，计算各试件故障时产品可靠度服从的β分布，进而使用该分布的概率密度来度量各试件故障时可靠度估计值的合理程度，以各试件故障时可靠度估计值的合理程度之积构造β似然函数，将使β似然函数取值最大的分布参数作为估计结果。

本发明一种基于β似然函数的参数估计方法，其具体步骤为：

步骤一、故障数据的收集；

步骤二、平均秩次的计算；

步骤三、寿命分布的选取和β似然函数的构造；

步骤四、寿命分布参数的求解。

其中，在步骤一中所述的“故障数据”是指，在寿命试验中产生的故障数据，具体包括样本量n、发生故障的试件数m、各故障试件的工作时间(由短到长排序为t₁,t₂,...,t_i,...,t_m)、未故障而撤离试验的各试件的工作时间；

其中，在步骤二中所述的“平均秩次”是指，对于有未故障而撤离试件的寿命数据，由于无法得知中途撤离的试件将在何时故障，假设中途撤离的试件将同等可能的在其撤离后的试件故障之间发生故障，则第i个故障试件在所有试件中的故障序号期望值即平均秩次A_i；

其中，在步骤二中所述的“平均秩次的计算”，其递推计算公式如下：

$\left\{\begin{array}{c}{A}_0=0\\ {\mathrm{\ΔA}}_i=\frac{n+1-A_{i-1}}{n-j+2}\\ A_i=A_{i-1}+{\mathrm{\ΔA}}_i\end{array}\right.$

式中n为样本量、i为故障序号、j为第i个故障在故障时间和撤离时间中的共同排序号；

其中，在步骤三中所述的“寿命分布的选取”是指，根据产品类型选择合适的寿命分布，如电子产品可选择指数分布、机械产品可选择威布尔分布；

其中，在步骤三中所述的“β似然函数的构造”是指，假设所选取寿命分布的可靠度函数为R(t,θ)(其中θ为待估参数)，则各故障发生时的可靠度估计值为R(t_k,θ)，根据可靠度非参数估计中的β分布法，对n个寿命独立同分布的产品进行寿命试验，则第k个产品故障时产品的可靠度服从β分布β(n-A_k+1,A_k)。根据β分布的定义，该可靠度估计值为R(t_k,θ)在可靠度分布β(n-A_k+1,A_k)中的概率密度为：

$f\left(R(t_k,\mathrm{\θ})\right)=\frac{R{(t_k,\mathrm{\θ})}^{(n-A_k)}\·{(1-R(t_k,\mathrm{\θ}))}^{A_k-1}}{B(n-A_k+1,A_k)}---\left(1\right)$

式中B(n-A_k+1,A_k)为自变量取n-A_k+1和A_k的β函数值。以该概率密度的大小作为评价此时可靠度估计的合理程度，进而以各个f(R(t_k,θ))的乘积作为评价整个估计结果的合理程度，由此得到β似然函数为：

$L\left(\mathrm{\θ}\right)=\underset{k=1}{\overset{m}{\mathrm{\Π}}}f\left(R(t_k,\mathrm{\θ})\right)=\underset{k=1}{\overset{m}{\mathrm{\Π}}}\frac{R{(t_k,\mathrm{\θ})}^{(n-A_k)}\·{(1-R(t_k,\mathrm{\θ}))}^{A_k-1}}{B(n-A_k+1,A_k)}---\left(2\right)$

其中，在步骤四中所述的“寿命分布参数的求解”是指，将使似然函数L(θ)取值最大的参数作为寿命分布参数的估计结果，由于似然函数形式较为复杂而难以使用解析方法求解其极大值，因此可运用数值方法进行求解，如使用MATLAB中的fminsearch()函数寻找-L(θ)的极小值即为L(θ)的极大值，此时的θ取值即为参数估计结果。

本发明的优点在于：

(1)β似然函数实际意义明确，直接将可靠度(生存率)的合理程度作为参数估计的准则，适用于如可靠性研究、医学研究等关注于事件发生比例的参数估计场合；

(2)β似然估计方法对各种数学形式的分布都适用，只需保证各试件故障时所选分布对应的可靠度函数取值在0到1之间即可，而这是对一个分布最基本的要求；

(3)β似然估计方法不会有类似极大似然估计不存在的情况，因为对于任何可能出现的可靠度所服从的β分布，其概率密度都是在一个有限的范围内，也就不会出现似然函数中有无穷大项的情况，因此具有广泛的适用性。

附图说明

图1是本发明的流程图；

图2是本发明实施例的估计结果可靠度函数图；

图中符号说明如下：

t为产品的工作时间；R(t)为通过β似然估计得到的产品可靠度。

具体实施方式

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

本发明是一种基于β似然函数的参数估计方法，流程如图1所示，包括以下几个步骤：

步骤一、故障数据的收集；

所需故障数据为无维修和替换的产品寿命数据，需收集的故障数据包括样本量n、发生故障的试件数m、各故障试件的工作时间(工作时间由短到长排序为t₁,t₂,...,t_m)、未故障而撤离试验的各试件的工作时间。按工作时间由短到长将故障和撤离一起进行排序。

步骤二、平均秩次的计算；

对于有未故障而撤离试件的寿命数据，由于无法得知中途撤离的试件将在何时故障，需对失效的平均秩次A_k进行计算，平均秩次的递推公式为：

$\left\{\begin{array}{c}{A}_0=0\\ {\mathrm{\ΔA}}_k=\frac{n+1-A_{k-1}}{n-i+2}\\ A_k=A_{k-1}+{\mathrm{\ΔA}}_k\end{array}\right.---\left(3\right)$

式中i为故障和撤离的共同排列顺序号。

步骤三、寿命分布的选取和β似然函数的构造；

寿命分布的选取应当由产品类型决定，如电子产品可选择指数分布、机械产品可选择威布尔分布。

假设所选取寿命分布的可靠度函数为R(t,θ)(其中θ为待估参数)，则各故障发生时的可靠度估计值为R(t_k,θ)，根据可靠度非参数估计中的β分布法，对n个寿命独立同分布的产品进行寿命试验，则第k个产品故障时产品的可靠度服从β分布β(n-A_k+1,A_k)。根据β分布的定义，该可靠度估计值为R(t_k,θ)在可靠度分布β(n-A_k+1,A_k)中的概率密度为：

$f\left(R(t_k,\mathrm{\θ})\right)=\frac{R{(t_k,\mathrm{\θ})}^{(n-A_k)}\·{(1-R(t_k,\mathrm{\θ}))}^{A_k-1}}{B(n-A_k+1,A_k)}---\left(4\right)$

式中B(n-A_k+1,A_k)为自变量取n-A_k+1和A_k的β函数值。以该概率密度的大小作为评价此时可靠度估计的合理程度，进而以各个f(R(t_k,θ))的乘积作为评价整个估计结果的合理程度，由此得到β似然函数为：

$L\left(\mathrm{\θ}\right)=\underset{k=1}{\overset{m}{\mathrm{\Π}}}f\left(R(t_k,\mathrm{\θ})\right)=\underset{k=1}{\overset{m}{\mathrm{\Π}}}\frac{R{(t_k,\mathrm{\θ})}^{(n-A_k)}\·{(1-R(t_k,\mathrm{\θ}))}^{A_k-1}}{B(n-A_k+1,A_k)}---\left(5\right)$

步骤四、分布参数的求解。

将使似然函数L(θ)取值最大的参数作为参数的估计结果，由于似然函数形式较为复杂而难以使用解析方法求解其极大值，因此可运用数值方法进行求解，如使用MATLAB中的fminsearch()函数寻找-L(θ)的极小值即为L(θ)的极大值，此时的θ取值即为参数估计结果。

实施例：

对某电子产品10件进行1000h定时截尾的寿命试验，试验中记录了试件的故障情况和撤离情况如表1所示，其中3号试件为人为损坏，记为未故障撤离，为了估计该产品的可靠度函数，现使用本发明所述的方法进行参数估计。

表1试验记录表

试件序号	工作时间(h)	终止原因	试件序号	工作时间(h)	终止原因
						1	216	故障	6	460	故障
2	50	故障	7	1000	试验截尾
						3	514	撤离	8	183	故障
4	565	故障	9	131	故障
						5	298	故障	10	940	故障

步骤一、故障数据的收集；

根据试验记录表可知，试件数量n＝10，发生故障的试件数量m＝8，进行时间排序如表2所示。

表2事件排序表

步骤二、平均秩次的计算；

采用递推公式(3)求故障平均秩次，计算结果如表3所示。

表3平均秩次表

步骤三、寿命分布的选取和β似然函数的构造；

通常认为电子产品的寿命服从指数分布，故而选择指数分布来对产品的可靠度进行估计，指数分布对应的可靠度函数为：

R(t)＝Exp(-λt)  (6)

式中λ为待估参数。构造β似然函数为：

$L\left(\mathrm{\λ}\right)=\underset{k=1}{\overset{8}{\mathrm{\Π}}}f\left(\mathrm{Exp}(-{\mathrm{\λt}}_k)\right)=\underset{k=1}{\overset{8}{\mathrm{\Π}}}\frac{\mathrm{Exp}{(-{\mathrm{\λt}}_k)}^{(10-A_k)}\·{(1-\mathrm{Exp}(-{\mathrm{\λt}}_k))}^{A_k-1}}{B(11-A_k,A_k)}$

式中各A_k、t_k取值见表3。

步骤四、分布参数的求解。

由于难以使用解析方法找到似然函数的极大值，故使用数值方法求似然函数的极大值。得当λ取0.0017时L(λ)取到极值5811，即参数λ的β似然估计值为0.0017，估计得到的产品可靠度函数为:

R(t)＝Exp(-0.0017t)

图2为可靠度函数R(t)的β似然估计结果与近似中位秩公式计算的可靠度经验函数的对比。

Claims (5)

1.一种基于β似然函数的参数估计方法，其特征在于：该方法具体步骤如下：

步骤一、故障数据的收集；

步骤二、平均秩次的计算；

步骤三、寿命分布的选取和β似然函数的构造；

步骤四、寿命分布参数的求解。

2.根据权利要求1所述的一种基于β似然函数的参数估计方法，其特征在于：在步骤一中所述的“故障数据”是指在寿命试验中产生的故障数据，具体包括样本量n、发生故障的试件数m、各故障试件的工作时间即由短到长排序为t₁,t₂,...,t_i,...,t_m、未故障而撤离试验的各试件的工作时间。

3.根据权利要求1所述的一种基于β似然函数的参数估计方法，其特征在于：在步骤二中所述的“平均秩次”是指，对于有未故障而撤离试件的寿命数据，由于无法得知中途撤离的试件将在何时故障，假设中途撤离的试件将同等可能的在其撤离后的试件故障之间发生故障，则第i个故障试件在所有试件中的故障序号期望值即平均秩次A_i；在步骤二中所述的“平均秩次的计算”，其递推计算公式如下：

$\left\{\begin{array}{c}{A}_0=0\\ {\mathrm{\ΔA}}_i=\frac{n+1-A_{i-1}}{n-j+2}\\ A_i=A_{i-1}+{\mathrm{\ΔA}}_i\end{array}\right.$

式中，n为样本量、i为故障序号、j为第i个故障在故障时间和撤离时间中的共同排序号。

4.根据权利要求1所述的一种基于β似然函数的参数估计方法，其特征在于：在步骤三中所述的“寿命分布的选取”是指，根据产品类型选择合适的寿命分布，电子产品选择指数分布、机械产品选择威布尔分布；在步骤三中所述的“β似然函数的构造”是指，假设所选取寿命分布的可靠度函数为R(t,θ)，其中θ为待估参数，则各故障发生时的可靠度估计值为R(t_k,θ)，根据可靠度非参数估计中的β分布法，对n个寿命独立同分布的产品进行寿命试验，则第k个产品故障时产品的可靠度服从β分布β(n-A_k+1,A_k)；根据β分布的定义，该可靠度估计值为R(t_k,θ)在可靠度分布β(n-A_k+1,A_k)中的概率密度为：

$f\left(R(t_k,\mathrm{\θ})\right)=\frac{R{(t_k,\mathrm{\θ})}^{(n-A_k)}\·{(1-R(t_k,\mathrm{\θ}))}^{A_k-1}}{B(n-A_k+1,A_k)}---\left(1\right)$

式中，B(n-A_k+1,A_k)为自变量取n-A_k+1和A_k的β函数值；以该概率密度的大小作为评价此时可靠度估计的合理程度，进而以各个f(R(t_k,θ))的乘积作为评价整个估计结果的合理程度，由此得到β似然函数为：

$L\left(\mathrm{\θ}\right)=\underset{k=1}{\overset{m}{\mathrm{\Π}}}f\left(R(t_k,\mathrm{\θ})\right)=\underset{k=1}{\overset{m}{\mathrm{\Π}}}\frac{R{(t_k,\mathrm{\θ})}^{(n-A_k)}\·{(1-R(t_k,\mathrm{\θ}))}^{A_k-1}}{B(n-A_k+1,A_k)}.---\left(2\right)$

5.根据权利要求1所述的一种基于β似然函数的参数估计方法，其特征在于：在步骤四中所述的“寿命分布参数的求解”是指，将使似然函数L(θ)取值最大的参数作为寿命分布参数的估计结果，由于似然函数形式较为复杂而难以使用解析方法求解其极大值，因此运用数值方法进行求解，使用MATLAB中的fminsearch()函数寻找-L(θ)的极小值即为L(θ)的极大值，此时的θ取值即为参数估计结果。