Birnbaum-Saunders疲劳寿命分布是概率物理方法中的一个重要失效模型，该模型是BIRNBANUM和SAUDERS^[1]于1969年在研究主因裂纹扩展导致材料失效过程中推导而来，主要应用于疲劳失效研究，它比常用寿命分布如Weibull分布等更适合描述某些由疲劳引起失效的产品寿命失效规律.

设T服从两参数Birnbaum-Saunders疲劳寿命分布BS(α, β)，其分布函数与密度函数分别为：

$ {F\left( t \right) = \mathit{\Phi }\left[ {\frac{1}{\alpha }\left( {\sqrt {\frac{t}{\beta }} - \sqrt {\frac{\beta }{t}} } \right)} \right],\;\;\;\;t > 0,} $

$ {f\left( t \right) = \frac{1}{{2\alpha \sqrt \beta }}\left( {\frac{1}{{\sqrt t }} + \frac{\beta }{{t\sqrt t }}} \right)\varphi \left[ {\frac{1}{\alpha }\left( {\sqrt {\frac{t}{\beta }} - \sqrt {\frac{\beta }{t}} } \right)} \right],\;\;\;\;t > 0,} $

其中，α＞0为形状参数，β＞0为刻度参数(或者称为尺度参数)，φ(x), Φ(x)分别为标准正态分布的密度函数与分布函数，即

$ \varphi \left( x \right) = \frac{1}{{\sqrt {2{\rm{ \mathsf{ π} }}} }}{{\rm{e}}^{ - \frac{{{x^2}}}{2}}},\;\;\;\;\mathit{\Phi }\left( x \right) = \int_{ - \infty }^x {\varphi \left( y \right){\rm{d}}y} . $

关于两参数Birnbaum-Saunders疲劳寿命分布BS(α, β)的统计分析方法已有众多研究^[1-29].需要指出的是，2014年BALAKRISHNAN等^[28]证明了在定数截尾和定时截尾下两参数BS分布参数的MLE存在且唯一，这一结果说明当形状参数α＞2时，文献[4]关于似然函数有多个极值点这一看法是不正确的.关于两参数BS分布刻度参数的区间估计通常采用文献[19-20]所提出的2种方法，徐晓岭等^[29]通过Monte-Carlo模拟说明这2种方法可能无法得到参数β的区间估计.

论文通过对数变换给出了求两参数BS疲劳寿命分布BS(α, β)在全样本场合下参数的对数矩估计.基于对数变换通过一阶泰勒展开，将两参数BS疲劳寿命分布BS(α, β)近似看作两参数对数正态分布，由此得到了2个参数α, β的近似区间估计，通过Monte-Carlo模拟发现，本文的近似方法比原有方法更精确.最后通过若干实例说明方法的可行性.

1 参数点估计的综合比较分析 1.1 方法1：参数点估计的新方法——对数矩估计

设T₁, T₂, …, T_n为来自Birnbaum-Saunders疲劳寿命分布总体T~BS(α, β)的一个容量为n的简单随机样本，其样本观察值为t₁, t₂, …, t_n.

记参数μ=lnβ，令Y=lnT, Y_i=lnT_i，i=1, 2, …, n，则Y₁, Y₂, …, Y_n为来自分布函数为${F_Y}\left( y \right) = \mathit{\Phi }\left\{ {\frac{1}{\alpha }\left[ {\exp \left( {\frac{{y - \mu }}{2}} \right) - \exp \left( { - \frac{{y - \mu }}{2}} \right)} \right]} \right\}$的一个容量为n的简单随机样本，其样本观察值记为y₁, y₂, …, y_n.

令$Z = \frac{{Y - \mu }}{2}$，则Z的分布函数F_Z(z)和密度函数f_Z(z)分别为：

$ {{F_Z}\left( z \right) = \mathit{\Phi }\left[ {\frac{1}{\alpha }\left( {{{\rm{e}}^z} - {{\rm{e}}^{ - z}}} \right)} \right],} $

$ {{f_Z}\left( z \right)\frac{1}{\alpha }\left( {{{\rm{e}}^z} - {{\rm{e}}^{ - z}}} \right)\varphi \left[ {\frac{1}{\alpha }\left( {{{\rm{e}}^z} - {{\rm{e}}^{ - z}}} \right)} \right], - \infty < z < + \infty .} $

若k为奇数时，E(Z^k)=0；若k为偶数时，

$ E\left( {{Z^k}} \right) = 2\int_0^{ + \infty } {{{\left( {\ln \frac{{\alpha t + \sqrt {{\alpha ^2}{t^2} + 4} }}{2}} \right)}^k}\varphi \left( t \right){\rm{d}}t} . $

Y的分布函数F_Y(y)和密度函数f_Y(y)分别为：对-∞＜y＜+∞，

$ {F_Y}\left( y \right) = \mathit{\Phi }\left\{ {\frac{1}{\alpha }\left[ {\exp \left( {\frac{{y - \mu }}{2}} \right) - \exp \left( { - \frac{{y - \mu }}{2}} \right)} \right]} \right\}, $

$ \begin{array}{l} {f_Y}\left( y \right) = \frac{1}{{2\alpha }}\left[ {\exp \left( {\frac{{y - \mu }}{2}} \right) - \exp \left( { - \frac{{y - \mu }}{2}} \right)} \right] \times \\ \;\;\;\;\;\;\;\;\;\;\;\;\varphi \left\{ {\frac{1}{\alpha }\left[ {\exp \left( {\frac{{y - \mu }}{2}} \right) - \exp \left( { - \frac{{y - \mu }}{2}} \right)} \right]} \right\}. \end{array} $

记$\bar Y = \frac{1}{n}\sum\limits_{i = 1}^n {{Y_i}} ,\overline {{Y^2}} = \frac{1}{n}\sum\limits_{i = 1}^n {Y_i^2} $，

$ E\left( Y \right) = \mu + 2E\left( Z \right) = \mu , $

$ \begin{array}{l} E\left( {{Y^2}} \right) = {\mu ^2} + 4E\left( {{Z^2}} \right) = \\ \;\;\;\;\;\;\;\;{\mu ^2} + 8\int_0^{ + \infty } {{{\left( {\ln \frac{{\alpha t + \sqrt {{\alpha ^2}{t^2} + 4} }}{2}} \right)}^2}\varphi \left( t \right){\rm{d}}t} . \end{array} $

由此，利用矩估计的思想可得参数μ的矩估计为$\hat \mu = \bar Y$，易见Y为μ的无偏估计，进而得参数β的点估计：${{\hat \beta }_1} = \exp \left( {\bar Y} \right) = {\left( {\prod\limits_{i = 1}^n {{T_i}} } \right)^{1/n}}$.

引理1^[30]  设X₁, X₂, …, X_n是来自总体X的容量为n的一个简单随机样本，记E(X)=μ，D(X)=σ²＜+∞，该总体X的4阶中心矩v₄=E(X-EX)⁴有限，若函数h(x)的四阶导数存在且有界，则有

$ E\left( {h\left( {\bar X} \right)} \right) = h\left( u \right) + \frac{1}{{2n}}h''\left( \mu \right){\sigma ^2} + O\left( {{n^{ - 2}}} \right), $

$ \begin{array}{*{20}{c}} {D\left[ {h\left( {\bar X} \right)} \right] = \frac{1}{n}{{\left[ {h'\left( \mu \right)} \right]}^2}{\sigma ^2} + \frac{1}{{{n^2}}}\left\{ {h'\left( \mu \right)h''\left( \mu \right){v_3} + } \right.}\\ {\left. {\frac{1}{2}{{\left[ {h''\left( \mu \right)} \right]}^2}{\sigma ^4} + h'\left( \mu \right)h''\left( \mu \right){\sigma ^4}} \right\} + O\left( {{n^{ - 3}}} \right).} \end{array} $

定理1  ${{\hat \beta }_1}$是β的渐近无偏估计和相合估计.

证明  Y=2Z+μ, E(Y)=μ, Y-E(Y)=2Z，则Y的一至四阶中心矩为：

$ {v_1} = E\left[ {Y - E\left( Y \right)} \right] = 0, $

$ \begin{array}{l} {v_2} = E{\left[ {Y - E\left( Y \right)} \right]^2} = 4E\left( {{Z^2}} \right) = \\ \;\;\;\;\;\;\;\;\;\;\;8\int_0^{ + \infty } {{{\left( {\ln \frac{{\alpha t + \sqrt {{\alpha ^2}{t^2} + 4} }}{2}} \right)}^2}\varphi \left( t \right){\rm{d}}t} , \end{array} $

$ {v_3} = E{\left[ {Y - E\left( Y \right)} \right]^3} = 8E\left( {{Z^3}} \right) = 0, $

$ \begin{array}{l} {v_4} = E{\left[ {Y - E\left( Y \right)} \right]^4} = 16E\left( {{Z^4}} \right) = \\ \;\;\;\;\;\;\;\;\;32\int_0^{ + \infty } {{{\left( {\ln \frac{{\alpha t + \sqrt {{\alpha ^2}{t^2} + 4} }}{2}} \right)}^2}\varphi \left( t \right){\rm{d}}t} . \end{array} $

令函数h(x)=e^x，h(x)的任意阶导数任为e^x，则

$ \begin{array}{l} E\left( {{{\mathop \beta \limits^ \wedge }_1}} \right) = E\left[ {\exp \left( {\bar Y} \right)} \right] = {{\rm{e}}^\mu } + \frac{1}{{2n}}{{\rm{e}}^\mu } \cdot 4E\left( {{Z^2}} \right) + O\left( {{n^{ - 2}}} \right) = \\ \;\;\;\;\;\;\;\;\;\;\;\beta + \frac{2}{n}\beta E\left( {{Z^2}} \right) + O\left( {{n^{ - 2}}} \right), \end{array} $

$ \begin{array}{l} D\left( {{{\mathop \beta \limits^ \wedge }_1}} \right) = D\left[ {\exp \left( {\bar Y} \right)} \right] = \frac{1}{{2n}}{{\rm{e}}^{2\mu }} \cdot 4E\left( {{Z^2}} \right) + \\ \frac{1}{{{n^2}}}\left\{ {\frac{1}{2}{{\rm{e}}^{2\mu }} \cdot 16{{\left[ {E\left( {{Z^2}} \right)} \right]}^2} + {{\rm{e}}^{2\mu }} \cdot 16{{\left[ {E\left( {{Z^2}} \right)} \right]}^2}} \right\} + \\ O\left( {{n^{ - 3}}} \right) = \frac{4}{n}{\beta ^2}E\left( {{Z^2}} \right) + \frac{{24}}{{{n^2}}}{\beta ^2}{\left[ {E\left( {{Z^2}} \right)} \right]^2} + O\left( {{n^{ - 3}}} \right). \end{array} $

易见，$\mathop {\lim }\limits_{n \to + \infty } E\left( {{{\hat \beta }_1}} \right) = \beta ,\mathop {\lim }\limits_{n \to + \infty } D\left( {{{\hat \beta }_1}} \right) = 0$.

则${{\hat \beta }_1}$是β的渐近无偏估计和相合估计.而参数α的矩估计${{\hat \alpha }_1}$为如下方程的根：

$ \int_0^{ + \infty } {{{\left( {\ln \frac{{\alpha t + \sqrt {{\alpha ^2}{t^2} + 4} }}{2}} \right)}^2}\varphi \left( t \right){\rm{d}}t} = \frac{{\overline {{Y^2}} - {{\bar Y}^2}}}{8}. $

易证上述方程有唯一正实根.

需要指出的是，${{\hat \alpha }_1}$是含有积分运算复杂的超越方程的根，不适合工程应用.所以在此将相应α的估计取为：

${{\overset{\wedge }{\mathop{\alpha }}\,}_{1}}=\sqrt{\frac{1}{n}\sum\limits_{i=1}^{n}{{{\left( \sqrt{\frac{{{t}_{i}}}{{{\overset{\wedge }{\mathop{\beta }}\,}_{1}}}}-\sqrt{\frac{{{\overset{\wedge }{\mathop{\beta }}\,}_{1}}}{{{t}_{i}}}} \right)}^{2}}}}.$

1.2 方法2：参数的极大似然估计^{[2, 14]}

似然函数为：

$ \begin{array}{l} L\left( {\alpha ,\beta } \right) = \\ \prod\limits_{i = 1}^n {\left\{ {\frac{1}{{2\alpha \sqrt \beta }}\left( {\frac{1}{{\sqrt {{t_i}} }} + \frac{\beta }{{{t_i}\sqrt {{t_i}} }}} \right)\varphi \left[ {\frac{1}{\alpha }\left( {\sqrt {\frac{{{t_i}}}{\beta }} - \sqrt {\frac{\beta }{{{t_i}}}} } \right)} \right]} \right\}} . \end{array} $

令$\frac{{\partial \ln L\left( {\alpha ,\beta } \right)}}{{\partial \alpha }} = 0,\frac{{\partial \ln L\left( {\alpha ,\beta } \right)}}{{\partial \alpha }} = 0$，得方程组：

$ \left\{ \begin{array}{l} - \frac{n}{\alpha } + \frac{1}{{{\alpha ^3}}}\sum\limits_{i = 1}^n {{{\left( {\sqrt {\frac{{{t_i}}}{\beta }} - \sqrt {\frac{\beta }{{{t_i}}}} } \right)}^2}} = 0\\ - \frac{n}{{2\beta }} + \sum\limits_{i = 1}^n {\frac{1}{{{t_i} + \beta }} - \frac{1}{{2{\alpha ^2}}}} \sum\limits_{i = 1}^n {\left( { - \frac{{{t_i}}}{{{\beta ^2}}} + \frac{1}{{{t_i}}}} \right) = 0} , \end{array} \right. $

化简得仅含参数β的超越方程：

$ \frac{{\sum\limits_{i = 1}^n {\left( { - {t_i} + \frac{{{\beta ^2}}}{{{t_i}}}} \right)} }}{{\sum\limits_{i = 1}^n {\left( {{t_i} + \frac{{{\beta ^2}}}{{{t_i}}} - 2\beta } \right)} }} - \frac{2}{n}\sum\limits_{t = 1}^n {\frac{\beta }{{{t_i} + \beta }} + 1 = 0.} $

从上述超越方程中可以解出${{\hat \beta }_2}$，进而可得参数α的极大似然估计${{\hat \alpha }_2}$，即

${{\overset{\wedge }{\mathop{\alpha }}\,}_{2}}=\sqrt{\frac{1}{n}\sum\limits_{i=1}^{n}{{{\left( \sqrt{\frac{{{t}_{i}}}{{{\overset{\wedge }{\mathop{\beta }}\,}_{2}}}}-\sqrt{\frac{{{\overset{\wedge }{\mathop{\beta }}\,}_{2}}}{{{t}_{i}}}} \right)}^{2}}}}.$

需要指出的是，极大似然估计需要解非线性方程组，且计算较为复杂.

1.3 方法3：矩估计

作简单运算，可得BS(α, β)分布的如下特征.

引理2^[14]  设随机变量T~BS(α, β)，则

(1) T^-1~BS(α, β^-1)；

(2) $E\left( T \right) = \frac{{{\beta ^2}}}{2}\left( {{\alpha ^2} + 2} \right)$,

$ D\left( T \right) = \frac{{{\beta ^2}}}{4}\left( {5{\alpha ^4} + 4{\alpha ^2}} \right); $

(3) ${\mathop{\rm cov}} \left( {T,{T^{ - 1}}} \right) = - \frac{{{\alpha ^2}\left( {{\alpha ^2} + 4} \right)}}{4}$；

(4) 变异系数为$\sqrt {\frac{{5{\alpha ^4} + 4{\alpha ^2}}}{{{{\left( {{\alpha ^2} + 2} \right)}^2}}}} $，偏度为$\frac{{16{\alpha ^2}\left( {11{\alpha ^2} + 6} \right)}}{{{{\left( {5{\alpha ^2} + 4} \right)}^3}}}$，峰度为$3 + \frac{{6{\alpha ^2}\left( {93{\alpha ^2} + 41} \right)}}{{{{\left( {5{\alpha ^2} + 4} \right)}^2}}}$.

记$\bar T = \frac{1}{n}\sum\limits_{i = 1}^n {{T_i}} ,\overline {{T^2}} = \frac{1}{n}\sum\limits_{i = 1}^n {T_i^2} $，则可建立如下矩方程：

$ \frac{\beta }{2}\left( {{\alpha ^2} + 2} \right) = \bar T,\frac{{{\beta ^2}}}{2}\left( {3{\alpha ^4} + 4{\alpha ^2} + 2} \right) = \overline {{T^2}} . $

化简得仅含参数α的超越方程：

$ 2\frac{{3{\alpha ^4} + 4{\alpha ^2} + 2}}{{{{\left( {{\alpha ^2} + 2} \right)}^2}}} = \frac{{\overline {{T^2}} }}{{{{\bar T}^2}}}. $

记$c = \frac{{\overline {{T^2}} }}{{\overline {{T^2}} }} ＞ 1$，此时上述方程化为

$ 6{\alpha ^4} + 8{\alpha ^2} + 4 = c{\alpha ^4} + 4c{\alpha ^2} + 4c, $

即

$ \left( {c - 6} \right){\alpha ^4} + 4\left( {c - 2} \right){\alpha ^2} + 4\left( {c - 1} \right) = 0, $

$ \begin{array}{l} \Delta = 16{\left( {c - 2} \right)^2} - 16\left( {c - 6} \right)\left( {c - 1} \right) = \\ \;\;\;\;\;\;16\left( {3c - 2} \right) > 0, \end{array} $

又

$ \begin{array}{l} {\left( {c - 2} \right)^2} - \left( {3c - 2} \right) = {c^2} - 7c + 6 = \\ \;\;\;\;\;\left( {c - 1} \right)\left( {c - 6} \right). \end{array} $

于是可知：当c=6时，方程无实根；当1＜c＜6时，方程存在唯一正实根：${{\hat \alpha }^2} = \frac{{2\left( {c - 2} \right) + 2\sqrt {3c - 2} }}{{6 - c}}$；而当c＞6时，方程无正实根.

由此，当1＜c＜6时，可得参数α的矩估计为${{\hat \alpha }_3} = \sqrt {\frac{{2\left( {c - 2} \right) + 2\sqrt {3c - 2} }}{{6 - c}}} $，进而得参数β的矩估计为${{\hat \beta }_3} = \frac{{2\bar T}}{{\hat \alpha _3^2 + 2}} = \frac{{6 - c}}{{4 + \sqrt {3c - 2} }}\bar T$.

要使方法3的矩估计存在，则样本数据应满足：1＜c＜6，否则点估计不存在.下表 1给定样本容量n=10, 20, 30, 40, 50以及参数真值β=1，α=0.1, 0.5, 1, 2, 5, 7, 10, 15，通过10 000次模拟统计了其数据满足1＜c＜6的次数，从中可以看到，在α较小时，方法3的矩估计基本是存在的.

1.4 方法4：矩估计

由于T^-1~BS(α, β^-1)，记$\overline {{{\bar T}^{ - 1}}} = \frac{1}{n}\sum\limits_{i = 1}^n {\frac{1}{{{T_i}}}} $，于是可建立如下矩方程：

$ \frac{\beta }{2}\left( {{\alpha ^2} + 2} \right) = \bar T,\;\;\;\frac{1}{{2\beta }}\left( {{\alpha ^2} + 2} \right) = \overline {{T^{ - 1}}} , $

从中可解得参数β, α的矩估计分别为：

${{\overset{\wedge }{\mathop{\beta }}\,}_{4}}=\sqrt{\frac{{\bar{T}}}{\overline{{{T}^{-1}}}}},\quad {{\overset{\wedge }{\mathop{\alpha }}\,}_{4}}=\sqrt{2\left( \sqrt{\bar{T}\overline{{{T}^{-1}}}}-1 \right)}.$

1.5 方法5：逆矩估计

文献[15]给出了参数的逆矩估计如下：

${{\overset{\wedge }{\mathop{\beta }}\,}_{5}}=\frac{\sum\limits_{i=1}^{n}{\sqrt{{{T}_{i}}}}}{\sum\limits_{i=1}^{n}{\sqrt{\frac{1}{{{T}_{i}}}}}}=\frac{\overline{\sqrt{T}}}{\overline{\sqrt{{{T}^{-1}}}}},$

${{\overset{\wedge }{\mathop{\alpha }}\,}_{5}}=\sqrt{\frac{1}{n}\sum\limits_{i=1}^{n}{{{\left( \sqrt{\frac{{{T}_{i}}}{{{\overset{\wedge }{\mathop{\beta }}\,}_{5}}}}-\sqrt{\frac{{{\overset{\wedge }{\mathop{\beta }}\,}_{5}}}{{{T}_{i}}}} \right)}^{2}}}},$

其中，$\overline {\sqrt T } = \frac{1}{n}\sum\limits_{i = 1}^n {\sqrt {{T_i}} } ,\overline {\sqrt {{T^{ - 1}}} } = \frac{1}{n}\sum\limits_{i = 1}^n {\sqrt {T_i^{ - 1}} } $.

1.6 方法6：分位数估计

设T₁, T₂, …, T_n为来自Birnbaum-Saunders疲劳寿命分布总体T~BS(α, β)的一个容量为n的简单随机样本，其次序统计量记为T₍₁₎≤T₍₂₎≤…≤T_(n).

由于$F\left( \beta \right) = \frac{1}{2}$，则参数β的点估计为

当n为奇数时，

${{\overset{\wedge }{\mathop{\beta }}\,}_{6}}={{T}_{\left( \left[ \left( n+1 \right)/2 \right] \right)}};$

当n为偶数时，

${{\overset{\wedge }{\mathop{\beta }}\,}_{6}}=\frac{{{T}_{\left( n/2 \right)}}+{{T}_{\left( n/2+1 \right)}}}{2}.$

进而参数α的点估计可取为

${{\overset{\wedge }{\mathop{\alpha }}\,}_{6}}=\sqrt{\frac{1}{n}\sum\limits_{i=1}^{n}{{{\left( \sqrt{\frac{{{T}_{i}}}{{{\overset{\wedge }{\mathop{\beta }}\,}_{6}}}}-\sqrt{\frac{{{\overset{\wedge }{\mathop{\beta }}\,}_{6}}}{{{T}_{i}}}} \right)}^{2}}}}.$

1.7 方法7：回归估计

文献[21]利用回归分析模型，给出了如下参数的最小二乘估计：

${{\overset{\wedge }{\mathop{\beta }}\,}_{7}}=\sqrt{\frac{\sum\limits_{i=1}^{n}{{{T}_{i}}}}{\sum\limits_{i=1}^{n}{\frac{1}{{{T}_{i}}}}}}=\sqrt{\frac{{\bar{T}}}{\overline{{{T}^{-1}}}}},$

${{\overset{\wedge }{\mathop{\alpha }}\,}_{7}}=\sqrt{\frac{2n}{n-1}\left( \sqrt{\bar{T}\overline{{{T}^{-1}}}}-1 \right)},$

注意到，${{\hat \beta }_7}$与${{\hat \beta }_4}$是一样的.

1.8 方法8：回归估计

由于

$ \frac{1}{\alpha }\left( {\sqrt {\frac{T}{\beta }} - \sqrt {\frac{\beta }{T}} } \right) \sim N\left( {0,1} \right), $

$ \frac{1}{{\alpha \sqrt \beta }}\left( {\sqrt T - \frac{\beta }{{\sqrt T }}} \right) \sim N\left( {0,1} \right), $

则

$ Z = \sqrt T - \frac{\beta }{{\sqrt T }} \sim N\left( {0,{\alpha ^2}\beta } \right). $

令$\sqrt T - \frac{\beta }{{\sqrt T }} = \varepsilon ,\varepsilon \sim N\left( {0,{\alpha ^2}\beta } \right)$，则令

$ Q = \sum\limits_{i = 1}^n {{{\left( {\sqrt {{T_i}} - \frac{\beta }{{\sqrt {{T_i}} }}} \right)}^2}} , $

$ \begin{array}{l} \frac{{{\rm{d}}Q}}{{{\rm{d}}\beta }} = - 2\sum\limits_{i = 1}^n {\left[ {\frac{1}{{\sqrt {{T_i}} }}\left( {\sqrt {{T_i}} - \frac{\beta }{{\sqrt {{T_i}} }}} \right)} \right]} = \\ \;\;\;\;\;\;\;\; - 2\sum\limits_{i = 1}^n {\left( {1 - \frac{\beta }{{{T_i}}}} \right)} . \end{array} $

令$\frac{{{\rm{d}}Q}}{{{\rm{d}}\beta }} = 0$，得$\sum\limits_{i = 1}^n {\left( {1 - \frac{\beta }{{{T_i}}}} \right)} = 0$，

解得

${{\hat \beta }_8} = \frac{n}{{\sum\limits_{i = 1}^n {\frac{1}{{{T_i}}}} }} = \frac{1}{{\overline {{T^{ - 1}}} }},$

进而，

$\begin{array}{*{35}{l}} \overset{\wedge }{\mathop{\alpha _{8}^{2}}}\,{{\overset{\wedge }{\mathop{\beta }}\,}_{8}}=\frac{1}{n-1}\sum\limits_{i=1}^{n}{{{\left( \sqrt{{{T}_{i}}}-\frac{{{\overset{\wedge }{\mathop{\beta }}\,}_{8}}}{\sqrt{{{T}_{i}}}} \right)}^{2}}}= \\ \ \ \ \ \ \ \ \ \ \frac{1}{n-1}\sum\limits_{i=1}^{n}{\left( {{T}_{i}}+\frac{\overset{\wedge }{\mathop{\beta _{8}^{2}}}\,}{{{T}_{i}}}-2{{\overset{\wedge }{\mathop{\beta }}\,}_{8}} \right)=} \\ \ \ \ \ \ \ \ \ \ \frac{1}{n-1}\left( \sum\limits_{i=1}^{n}{{{T}_{i}}+\overset{\wedge }{\mathop{\beta _{8}^{2}}}\,\sum\limits_{i=1}^{n}{\frac{1}{{{T}_{i}}}-2n{{\overset{\wedge }{\mathop{\beta }}\,}_{8}}}} \right)= \\ \ \ \ \ \ \ \ \ \ \frac{n}{n-1}\left( \bar{T}-\frac{1}{\overline{{{T}^{-1}}}} \right), \\ \end{array}$

${{\overset{\wedge }{\mathop{\alpha }}\,}_{8}}=\sqrt{\frac{n}{n-1}\left( \bar{T}\overline{{{T}^{-1}}}-1 \right)}.$

1.9 方法9：回归估计

由于

$ \frac{1}{\alpha }\left( {\sqrt {\frac{T}{\beta }} - \sqrt {\frac{\beta }{T}} } \right) \sim N\left( {0,1} \right), $

$ \frac{{\sqrt \beta }}{\alpha }\left( {\frac{{\sqrt T }}{\beta } - \frac{1}{{\sqrt T }}} \right) \sim N\left( {0,1} \right), $

则

$Z = \frac{{\sqrt T }}{\beta } - \frac{1}{{\sqrt T }} \sim N\left( {0,\frac{{{\alpha ^2}}}{\beta }} \right).$

令$\varepsilon = \frac{{\sqrt T }}{\beta } - \frac{1}{{\sqrt T }},\varepsilon \sim N\left( {0,\frac{{{\alpha ^2}}}{\beta }} \right)$，则令

$ Q = \sum\limits_{i = 1}^n {{{\left( {\frac{{\sqrt {{T_i}} }}{\beta } - \frac{1}{{\sqrt {{T_i}} }}} \right)}^2}} , $

$ \begin{array}{l} \frac{{{\rm{d}}Q}}{{{\rm{d}}\beta }} = - \frac{2}{{{\beta ^2}}}\sum\limits_{i = 1}^n {\left[ {\sqrt {{T_i}} \left( {\frac{{\sqrt {{T_i}} }}{\beta } - \frac{1}{{\sqrt {{T_i}} }}} \right)} \right]} = \\ \;\;\;\;\;\;\;\; - \frac{2}{{{\beta ^2}}}\sum\limits_{i = 1}^n {\left( {\frac{{{T_i}}}{\beta } - 1} \right)} , \end{array} $

令$\frac{{{\rm{d}}Q}}{{{\rm{d}}\beta }} = 0$，得$ - \frac{2}{{{\beta ^2}}}\sum\limits_{i = 1}^n {\left( {\frac{{{T_i}}}{\beta } - 1} \right)} = 0$，

解得${{\hat \beta }_9} = \bar T$，

$ \begin{array}{l} \frac{{\mathop {\alpha _9^2}\limits^ \wedge }}{{\mathop {{\beta _9}}\limits^ \wedge }} = \frac{1}{{n - 1}}\sum\limits_{i = 1}^n {{{\left( {\frac{{\sqrt {{T_i}} }}{{\mathop {{\beta _9}}\limits^ \wedge }} - \frac{1}{{\sqrt {{T_i}} }}} \right)}^2}} = \\ \;\;\;\;\;\;\;\;\frac{1}{{n - 1}}\sum\limits_{i = 1}^n {\left( {\frac{{{T_i}}}{{\mathop {\beta _9^2}\limits^ \wedge }} - \frac{2}{{\mathop {{\beta _9}}\limits^ \wedge }} + \frac{1}{{{T_i}}}} \right)} , \end{array} $

$ \begin{array}{l} \mathop {\alpha _9^2}\limits^ \wedge = \frac{1}{{n - 1}}\sum\limits_{i = 1}^n {\left( {\frac{{{T_i}}}{{\mathop {{\beta _9}}\limits^ \wedge }} - 2 + \frac{{\mathop {{\beta _9}}\limits^ \wedge }}{{{T_i}}}} \right)} = \\ \;\;\;\;\;\;\;\;\frac{n}{{n - 1}}\left( {\bar T\overline {{T^{ - 1}}} - 1} \right), \end{array} $

$\overset{\wedge }{\mathop{{{\alpha }_{9}}}}\,=\sqrt{\frac{n}{n-1}\left( \bar{T}\overline{{{T}^{-1}}}-1 \right)}.$

1.10 参数点估计的模拟及比较分析

比较上述9种点估计方法，方法2由于涉及复杂的超越方程求解，在此不推荐使用.

下面通过10 000次Monte-Carlo模拟比较方法1、方法3~方法9的β点估计的精度.给定样本容量n=10, 20, 30, 35，参数α=0.5, 1, 1.5，β=1，通过10 000次模拟得参数β的估计均值与均方差(其中方法3都满足1＜c＜6)，结果列于表 2.

(1) 固定n，随着α的增大，参数β估计的均值与真值的偏差增大，均方误差也增大；

(2) 固定参数α的真值，参数β的估计随着n的增加变精确；

(3) 方法1，4，5中参数β估计的无偏性较为明显，同时其均方误差也较小，相对而言，方法4的均方误差更小，方法5与方法1均方误差很接近.

下面通过10 000次Monte-Carlo模拟，比较方法1、方法3~方法8(方法9与方法8同)的α点估计的精度.给定样本容量n=10, 20, 30, 35，参数α=0.5, 1, 1.5，β=1，通过10 000次模拟得参数α的估计均值与均方差(其中方法3都满足1＜c＜6)，结果列于表 3.

(1) 固定α，随着n的增大，参数α估计的均方误差在减少，即估计变精确；

(2) 方法1、方法4~方法7的均方误差相差不大，考虑到其均值，使用方法6或方法7更合理.

2 参数区间估计的比较分析

首先，给出利用对数正态分布求参数的区间估计方法(记为方法2)，进而通过Monte-Carlo模拟与文献[20-21]的方法(记为方法1)比较分析(在方法1能得到刻度参数β的区间估计的前提下).

记参数μ=ln β，令Y=lnT, Y_i=lnT_i，i=1, 2, …, n，则Y₁, Y₂, …, Y_n为来自分布函数为

${{F}_{Y}}\left( y \right)=\mathit{\Phi }\left\{ \frac{1}{\alpha }\exp \left( \frac{y-\mu }{2} \right)-\exp \left( -\frac{y-\mu }{2} \right) \right\}$

的一个容量为n的简单随机样本，其次序观察值记为y₁, y₂, …, y_n.

将$\frac{1}{\alpha }\left[ {\exp \left( {\frac{{y - \mu }}{2}} \right) - \exp \left( { - \frac{{y - \mu }}{2}} \right)} \right]$在y=μ处一阶泰勒展开得

$\frac{1}{\alpha }\left[ {\exp \left( {\frac{{y - \mu }}{2}} \right) - \exp \left( { - \frac{{y - \mu }}{2}} \right)} \right] \approx \frac{{Y - \mu }}{\alpha }$, 此时，${{F}_{Y}}\left( y \right)=\mathsf{\Phi }\left\{ \frac{1}{\alpha }\exp \left( \frac{y-\mu }{2} \right)-\exp \left( -\frac{y-\mu }{2} \right) \right\}$近似为${F_Y}\left( y \right) \approx \mathit{\Phi }\left( {\frac{{y - \mu }}{\alpha }} \right)$，即可将Y=lnT近似看作正态分布，$Y\mathop \sim \limits^ \cdot N\left( {\mu ,{\alpha ^2}} \right)$，或者$T\mathop \sim \limits^ \cdot LN\left( {\mu ,{\alpha ^2}} \right)$，也即将T近似看作对数正态分布.记

$ \bar Y = \frac{1}{n}\sum\limits_{i = 1}^n {{Y_i}} = \frac{1}{n}\sum\limits_{i = 1}^n {\ln {T_i}} , $

$ S_Y^2 = \frac{1}{{n - 1}}\sum\limits_{i = 1}^n {{{\left( {{Y_i} - \bar Y} \right)}^2}} = \frac{n}{{n - 1}}\left( {\overline {{Y^2}} - {{\bar Y}^2}} \right). $

由于$\frac{{\bar Y - \mu }}{{{S_Y}/\sqrt n }}\mathop \sim \limits^ \cdot t\left( {n - 1} \right)$，给定显著性水平γ，参数μ的置信水平1-γ的近似区间估计为

$ \left[ {\bar Y - \frac{{{S_Y}}}{{\sqrt n }}{t_{\gamma /2}}\left( {n - 1} \right),\bar Y + \frac{{{S_Y}}}{{\sqrt n }}{t_{\gamma /2}}\left( {n - 1} \right)} \right]. $

进而参数β的置信水平1-γ的近似区间估计为

$ \left[ {\exp \left[ {\bar Y - \frac{{{S_Y}}}{{\sqrt n }}{t_{\gamma /2}}\left( {n - 1} \right)} \right],\exp \left[ {\bar Y + \frac{{{S_Y}}}{{\sqrt n }}{t_{\gamma /2}}\left( {n - 1} \right)} \right]} \right], $

又

$ \frac{{\left( {n - 1} \right)S_Y^2}}{{{\alpha ^2}}}\dot \sim {\chi ^2}\left( {n - 1} \right). $

参数α的置信水平1-γ的近似区间估计为

$ \left[ {\sqrt {\frac{{\left( {n - 1} \right)S_Y^2}}{{\chi _{\gamma /2}^2\left( {n - 1} \right)}}} ,\sqrt {\frac{{\left( {n - 1} \right)S_Y^2}}{{\chi _{1 - \gamma /2}^2\left( {n - 1} \right)}}} } \right]. $

下面通过10 000次Monte-Carlo模拟，比较方法1、方法2关于参数α, β的区间估计的优劣.给定样本容量n=5, 10, 15，参数α=0.5, 1, 1.5，β=1，置信水平1-γ=0.90，通过10 000次模拟(方法1中的参数β的区间估计都存在)得参数α, β区间估计的平均下限、平均上限、平均长度，以及区间估计包含参数真值的次数，结果列于表 4.可知方法1、方法2所得的区间估计包含参数的真值次数都大于9 000；方法2所得的区间估计的平均长度较方法1要短.可知方法2比方法1更优.

3 计算实例

例1^[14]  数据来自BJERKEDAL(1960), 也被GUPTAETAL(1997)分析过.数据表示几内亚猪注射不同剂量结核杆菌的生存时间.几内亚猪对结核杆菌的敏感性比人类更高，首先关注在相同笼子里用相同养殖法养殖的猪.对于文中的养殖方法，有如下72个观测值：

12，15，22，24，24，32，32，33，34，38，38，43，44，48，52，53，54，54，55，56，57，58，58，59，60，60，60，60，61，62，63，65，65，67，68，70，70，72，73，75，76，76，81，83，84，85，87，91，95，96，98，99，109，110，121，127，129，131，143，146，146，175，175，211，233，258，258，263，297，341，341，376.

经计算得：c=1.651 2, d=0.063 1.

取置信水平为0.90，

$ {t_{0.05}}\left( {71} \right) = 1.6666,\;\;\;\;\chi _{0.05}^2\left( {71} \right) = 91.6702, $

$ \chi _{0.95}^2\left( {71} \right) = 52.6003. $

利用区间估计方法1、方法2得参数β, α的置信水平为0.90的区间估计，如表 6如示.

例2^[7]  本实例中的数据集为6061-T6铝合金的疲劳寿命.铝合金的切取位置应平行于轧制方向，振荡频率为18Hz.该数据集包括101个观测值，最大应力为31 000Pa.数据如下：

70，90，96，97，99，100，103，104，104，105，107，108，108，108，109，109，112，112，113，114，114，114，116，119，120，120，120，121，121，123，124，124，124，124，124，128，128，129，129，130，130，130，131，131，131，131，131，132，132，132，133，134，134，134，134，134，136，136，137，138，138，138，139，139，141，141，142，142，142，142，142，142，144，144，145，146，148，148，149，151，151，152，155，156，157，157，157，157，158，159，162，163，163，164，166，166，168，170，174，196，212.

经计算得：c=1.027 7, d=0.003 6.

文献[7]还给出了其他几种估计：

取置信水平为0.90，

$ {t_{0.05}}\left( {100} \right) = 1.66023,\;\;\;\;\chi _{0.05}^2\left( {100} \right) = 124.342, $

$ \chi _{0.95}^2\left( {100} \right) = 77.9295. $

取置信水平为0.95，

$ {t_{0.025}}\left( {100} \right) = 1.98397,\;\;\;\;\chi _{0.025}^2\left( {100} \right) = 129.561, $

$ \chi _{0.975}^2\left( 9 \right) = 74.2219. $

利用区间估计方法1、方法2得参数β, α的置信水平为0.90, 0.95时的区间估计，如表 9如示.

例3^[31]  对于空气污染物浓度，假定数据不相关且相互独立，因此，一昼夜或循环趋势分析是无必要的.这一假设得到了很多学者(包括GOKHALE和KHARE等)的支持.例如，环境数据有时以均值作为指标，因此空间-时间依赖性不复存在.以下数据对应的是1973年5~9月纽约每日的臭氧浓度(由纽约州环境保护部提供)：41，36，12，18，28，23，19，8，7，16，11，14，18，14，34，6，30，1，11，4，32，23，45，115，37，29，71，39，23，21，37，20，12，13，49，32，64，40，77，97，97，85，10，27，7，48，35，61，79，63，16，108，20，52，82，50，64，59，39，9，16，78，35，66，122，89，110，44，65，22，59，23，31，44，21，9，45，168，73，76，118，84，85，96，78，91，47，32，20，23，21，24，44，21，28，9，13，46，18，13，24，16，23，36，7，14，30，14，18，20，11，135，80，28，73，13.

经计算得：c=1.607 8, d=0.092 8.

取置信水平为0.90，

$ {t_{0.05}}\left( {115} \right) = 1.65821,\;\;\;\;\chi _{0.05}^2\left( {115} \right) = 141.03, $

$ \chi _{0.95}^2\left( {115} \right) = 91.2422. $

利用区间估计方法1、方法2得参数β, α的置信水平为0.90时的区间估计，如表 11如示.

例4^[32]  (1)第1个实例是有关洪水峰值的超过数，含1958~1984年共72个超过数，四舍五入到小数点后1位，数据如下：

1.7，2.2，14.4，1.1，0.4，20.6，5.3，0.7，1.9，13.0，12.0，9.3，1.4，18.7，8.5，25.5，11.6，14.1，22.1，1.1，2.5，14.4，1.7，37.6，0.6，2.2，39.0，0.3，15.0，11.0，7.3，22.9，1.7，0.1，1.1，0.6，9.0，1.7，7.0，20.1，0.4，2.8，14.1，9.9，10.4，10.7，30.0，3.6，5.6，30.8，13.3，4.2，25.5，3.4，11.9，21.5，27.6，36.4，2.7，64.0，1.5，2.5，27.4，1.0，27.1，20.2，16.8，5.3，9.7，27.5，2.5，27.0.

经计算得：c=2.001 2, d=0.247 5.

取置信水平为0.90，

$ {t_{0.05}}\left( {71} \right) = 1.666,\;\;\;\;\chi _{0.05}^2\left( {71} \right) = 91.6702, $

$ \chi _{0.95}^2\left( {71} \right) = 52.6003. $

利用区间估计的方法1、方法2得参数β, α的置信水平为0.90时的区间估计，如表 13如示.

(2) 第2个实例是有关50个工业设备的使用期，时间为零时进行寿命测试，数据如下：

0.1，0.2，1，1，1，1，1，2，3，6，7，11，12，18，18，18，18，18，21，32，36，40，45，46，47，50，55，60，63，63，67，67，67，67，72，75，79，82，82，83，84，84，84，85，85，85，85，85，86，86.

经计算得：c=1.506 2, d=0.382.

取置信水平为0.90，

$ {t_{0.05}}\left( {49} \right) = 1.67655,\;\;\;\;\chi _{0.05}^2\left( {49} \right) = 66.3386, $

$ \chi _{0.95}^2\left( {49} \right) = 33.9303. $

利用区间估计方法1、方法2得参数β, α的置信水平为0.90时的区间估计，如表 15如示.