在实际的生产制造过程中，过程能力指数(process capability indices，PCIs)用于衡量一个过程的生产能力，是反映过程能力的一个重要的量化指标。工厂质量管理者可用PCIs判断一个过程是否处于控制状态，以此指导生产。根据方向的不同，PCIs的研究主要分为2类：一是在各种情况下，通过构建不同的PCIs更好地描述过程；二是对已构造的PCIs，探讨这些指数及其估计量的性质。

PCIs的研究以产品质量特性服从参数为μ和σ的正态概率模型N(μ, σ²)为前提，下文用USL和LSL分别表示产品质量规格的上、下限，μ表示过程均值，σ表示过程标准差，d=(USL-LSL)/2代表规格区间长度的一半，m=(USL+LSL)/2代表规格上限与规格下限之间的中点，T表示客户期望的目标值。第1个PCI是Juran等^[1]提出的C_p，其定义如下

$ C_p=\frac{\mathrm{USL}-\mathrm{LSL}}{6 \sigma} $

(1)

之后PCIs得到长足的发展，Kane^[2]、Chan等^[3]和Pearn等^[4]分别推广出C_pk、C_pm和C_pmk，这3个指数的定义依次为

$ C_{p k}=\min \left\{\frac{\mathrm{USL}-\mu}{3 \sigma}, \frac{\mu-\mathrm{LSL}}{3 \sigma}\right\}=\frac{d-|\mu-m|}{3 \sigma}, $

(2)

$ C_{p m}=\frac{\text { USL - LSL }}{6 \sqrt{\sigma^2+(\mu-T)^2}}=\frac{d}{3 \sqrt{\sigma^2+(\mu-T)^2}}, $

(3)

$ \begin{aligned} C_{p m k} & =\min \left\{\frac{\mathrm{USL}-\mu}{3 \sqrt{\sigma^2+(\mu-T)^2}}, \frac{\mu-L S L}{3 \sqrt{\sigma^2+(\mu-T)^2}}\right\} \\ & =\frac{d-|\mu-m|}{3 \sqrt{\sigma^2+(\mu-T)^2}} . \end{aligned} $

(4)

随后，Vännman^[5]将上述4个PCIs推广得到C_p(u, v)，该指数依赖于2个非负参数u和v，具体形式如下所示

$ C_p(u, v)=\frac{d-u|\mu-m|}{3 \sqrt{\sigma^2+v(\mu-T)^2}} . $

当质量特性不服从正态分布时，Chen和Pearn^[6]给出C_p(u, v)的扩展形式

$ C_{N\rho}(u, v)=\frac{d-u|M-m|}{3 \sqrt{\left(\frac{F_{99.865}-F_{0.135}}{6}\right)^2+v(M-T)^2}}, $

(5)

其中：F_α表示质量特性服从分布的α分位数，M为分布的中位数。

对于PCIs的统计推断，许多学者采用bootstrap方法构造PCIs的置信区间。如Franklin和Wasserman^[7]首先将bootstrap方法应用于求解C_pk的置信区间；Tong和Chen^[8]研究非正态情况下的bootstrap方法，之后又给出了正态假设下2个PCIs差的bootstrap置信区间^[9]；Panichkitkosolkul^[10]分析过程服从半逻辑分布(half Logistic distribution)时2个过程能力指数差的bootstrap置信区间；Rao等^[11]提出在逆瑞利分布(inverse Rayleigh distribution, IRD)和对数逻辑分布(log-logistic distribution, LLD)下非正态过程能力指数C_Npk的bootstrap置信区间; 颜斌等^[12]在前人研究的基础上重新设计指数C_p和C_pm，并构造了新指数的标准、分位数和偏差校正百分位数3种bootstrap置信区间；徐锋等^[13]在过程数据分布未知且存在自相关的情况下利用bootstrap方法构建了过程能力指数C_pk的置信区间。

广义p值的概念由Tsui和Weerahandi^[14]在1989年引入，并由Weerahandi^[15]于1993年推广至区间估计。此后广义推断方法成功地求解了一些复杂的推断问题，很多文献的模拟结果表明广义推断方法通常具有良好的频率特性。由此出现了大量的文献研究PCIs的广义推断问题。Mathew等^[16]运用广义推断方法构造了C_pk的置信下限；Kurian等^[17]基于单因素随机效应模型给出了3个指数C_pk、C_pmk以及C″_pk的广义置信区间，这里C″_pk的定义为$C_{p k}^{\prime \prime}=\frac{d^*-A^*}{3 \sigma}$，其中，d^*=min{USL-T, T-LSL}，$A^*=\max \left\{\frac{d^*(\mu-T)}{\mathrm{USL}-T}, \frac{d^*(T-\mu)}{T-\mathrm{LSL}}\right\}$，μ是过程均值，σ表示过程标准差，m=(USL+LSL)/2，d=(USL-LSL)/2，T是预设的目标值；Ye等^[18]研究具有平衡数据的一般随机效应模型下3个过程能力指数C_pk、C_pmk和C″_pk的广义置信区间；Kanichukattu和Luke^[19]运用广义枢轴量方法推导出基于正态分布样本3个指数差C_pk1-C_pk2, C_pmk1-C_pmk2, C_pm1-C_pm2的置信区间; Yao和Jun^[20]讨论在单向随机模型下基于对数正态分布的C_p的广义置信区间; 贺加贝和李新民^[21]分析不平衡单因素随机效应模型下，过程无能力指数$C_{p p}=\left(\frac{\mu-T}{D}\right)^2+\left(\frac{\sigma}{D}\right)^2$的广义置信区间、Fiducial区间和bootstrap置信区间，这里μ是过程均值，σ表示过程标准差，T是目标值，D=(USL-LSL)/6。

IRD和LLD在可靠性领域、质量控制和生存分析等发挥着极其重要的作用^[22-29]。本文重点研究过程能力指数C_Npk的区间估计，该指数定义为

$ C_{N p k}=2 \min \left\{\frac{\mathrm{USL}-\mathrm{M}}{F_{99.865}-F_{0.135}}, \frac{M-\mathrm{LSL}}{F_{99.865}-F_{0.135}}\right\} . $

(6)

结合式(5)和式(6)，可以发现当(u, v)=(1, 0)时，C_Npk=C_Np(1, 0)。通过对文献中关于PCIs的深入分析，注意到在IRD和LLD 2种分布下关于PCIs的广义置信区间至今还未被探讨。因此，本文的研究重点是利用广义推断的思想构造这2种分布下过程能力指数C_Npk的广义枢轴量与广义置信区间，然后与Rao等^[11]提出的bootstrap置信区间进行比较。

1 广义置信区间

枢轴量法是构造置信区间最常用的方法之一。然而，当统计推断问题中存在讨厌参数时，这种方法往往难以使用。尽管大样本近似方法被广泛应用于求解此类问题，但当样本量较小或中等时，其表现往往不够理想。为了解决这些问题，文献[15]引入广义枢轴量的概念，并提出基于广义枢轴量的广义置信区间。下面先简要回顾这种方法，更多细节可参考文献[30]。

令随机变量X=(X₁, X₂, …, X_n)，变量X的样本观测值是x，η=(θ, δ), θ是未知的兴趣参数，δ是一个讨厌参数向量，R=R(X; x, η)是关于变量X, x和η的函数。如果R满足下面2个条件：

1) 对于给定的观测值x，R的概率分布与未知参数无关；

2) R=R(X; x, η)的观测值r(x)=R(X; x, η)不依赖讨厌参数δ。

则称R是一个广义枢轴量。

若给定广义枢轴量R=R(X; x, η)和置信水平γ(0 < γ < 1)，在R的样本空间寻找一个满足P(R(X; x, η)∈C_γ)=γ的子集C_γ。

取Θ_γ={θ|R(X; x, η)∈C_γ}，则称Θ_γ为参数θ的一个置信系数为γ的广义置信区间。

2 方法 2.1 逆瑞利分布下C_Npk的广义置信区间

逆瑞利分布被应用于多种领域。Mukherjee和Saran^[22]指出IRD是一类特殊的非单调失效率分布，在可靠性研究中具有很重要的应用。文献[26, 28]研究了IRD分布下的验收抽样设计。一些学者又考虑了IRD的统计推断研究：Maleki Jebely等^[31]分析简单随机抽样下IRD的分布函数估计问题；龙兵和张忠占^[32]基于左删失样本，分析在部分加速寿命试验下IRD参数的极大似然估计(maximum likelihood estimation, MLE)和置信区间；陈蒙等^[33]给出基于简单随机抽样和排序集抽样IRD对应样本所含参数的Fisher信息量及参数的MLE、修正MLE、最优线性无偏估计和最优线性同变估计。

若连续型随机变量T具有概率密度函数

$f(t ; \sigma)= \begin{cases}\frac{t}{\sigma^2} \mathrm{e}^{\frac{-t^2}{2\sigma^2}}, & t \geqslant 0, \\ 0, & \text { others, }\end{cases}$，则称随机变量T是服从尺度参数为σ的瑞利分布(Rayleigh distribution)，记作T: Rayleigh(σ)，其中σ是大于0的常数。

若随机变量T: Rayleigh(σ)，则称X=1/T服从尺度参数为σ的IRD，记作X: IRD(σ)。IRD的概率密度函数、累分布函数和α分位数函数依次为

$ \begin{gathered} f(x ; \sigma)= \begin{cases}\frac{2 \sigma^2}{x^3} \mathrm{e}^{\frac{-\sigma^2}{x^2}}, & x \geqslant 0, \\ 0, & x<0, \end{cases} \\ F(x ; \sigma)= \begin{cases}\mathrm{e}^{\frac{-\sigma^2}{x^2}}, & x \geqslant 0, \\ 0, & x<0 .\end{cases} \\ \xi_\alpha^{\mathrm{IRD}}=\sigma(-\ln \alpha)^{-\frac{1}{2}}, \sigma>0 . \end{gathered} $

(7)

假设X=(X₁, X₂, …, X_n)是来自尺度参数为σ(σ>0)的逆瑞利随机变量，利用IRD的性质可得到$\left(\frac{1}{X_1}, \frac{1}{X_2}, \cdots, \frac{1}{X_n}\right): \operatorname{Rayleigh}\left(\frac{1}{\sqrt{2} \sigma}\right)$，进而有$\sum_{i=1}^n \frac{1}{X_i^2} /\left(\frac{1}{\sqrt{2} \sigma}\right)^2 \sim \chi^2(2 n)$, 则可推导出尺度参数σ的广义枢轴量为

$ R_\sigma=\sqrt{W / 2 \sum\limits_{i=1}^n \frac{1}{x_i^2}}, $

(8)

其中$W \triangleq \sum_{i=1}^n \frac{1}{X_i^2} /\left(\frac{1}{\sqrt{2} \sigma}\right)^2: \chi^2(2 n)$, x_i表示X_i的观测值，χ²(m)表示自由度为m的卡方分布。因此，结合式(6)、式(7)和式(8)可构造出过程能力指数C_Npk的广义枢轴量为

$ R_{C_{N p k}}(\sigma)=2 \min \left\{\frac{\mathrm{USL}-R_\sigma\left(-\ln \left(p_2\right)\right)^{-\frac{1}{2}}}{R_\sigma\left(-\ln \left(p_3\right)\right)^{-\frac{1}{2}}-R_\sigma\left(-\ln \left(p_1\right)\right)^{-\frac{1}{2}}}, \frac{R_\sigma\left(-\ln \left(p_2\right)\right)^{-\frac{1}{2}}-\mathrm{LSL}}{R_\sigma\left(-\ln \left(p_3\right)\right)^{-\frac{1}{2}}-R_\sigma\left(-\ln \left(p_1\right)\right)^{-\frac{1}{2}}}\right\}, $

(9)

这里USL=29, LSL=1, p₁=0.998 65, p₂=0.5, p₃=0.001 35。

下面给出计算IRD下过程能力指数C_Npk的双边广义置信区间的算法: 设数据为X₁, X₂, …, X_n。

Step 1   产生自由度为2n的卡方分布随机变量；

Step 2   通过式(8)和式(9)计算参数σ的广义枢轴量R_σ和指数C_Npk的广义枢轴量R_CNpk(σ)；

Step 3   重复Step 1和Step 2 M次。将这M个值的α/2和1-(α/2)分位点作为C_Npk的1-α置信区间。

2.2 对数逻辑分布下C_Npk的广义置信区间

对数逻辑分布是一类重要的统计模型。Verhulst^[34]在1838年首次提出用LLD模拟人口增长；文献[24]探究了这一分布在生存分析领域的应用；文献[23, 25, 27, 29]研究了LLD下截断寿命试验的几种验收抽样方案。He等^[35]主要研究排序集抽样下LLD中尺度参数和形状参数的估计问题；Hassan和Doori^[36]基于次序统计量研究平方损失函数下LLD生存函数和风险函数的贝叶斯估计；Ibrahim和Kalaf^[37]开发了一种将下降单纯形法(simplex downhill algorithm)与矩量法(moment method)相结合的SMOM新方法来研究LLD生存函数的估计。

若随机变量Z的概率密度函数为

$f(z ; \lambda, \beta)= \begin{cases}\frac{\mathrm{e}^{-(z-\lambda / \beta)}}{\beta\left(1+\mathrm{e}^{-(z-\lambda / \beta)}\right)^2}, & z>0, \\ 0, & \text { others , }\end{cases}$，则称Z服从尺度参数为λ、形状参数为β的双参数逻辑分布(logistic distribution)，记作Z~logistic(λ, β)，且λ, β>0。称λ=0, β=1的逻辑分布logistic(0, 1)为标准逻辑分布。

若随机变量Z~logistic(λ, β)，令X=logZ，则称随机变量X服从尺度参数为λ、形状参数为β的双参数LLD，记作X~LLD(λ, β)，其概率密度函数、累积分布函数和α分位数函数分别为

$ \begin{aligned} & f(x ; \lambda, \beta)=\frac{\frac{\beta}{\lambda}\left(\frac{x}{\lambda}\right)^{\beta-1}}{\left(1+\left(\frac{x}{\lambda}\right)^\beta\right)^2}, \\ & x>0, \lambda>0, \beta>1, \\ & F(x ; \lambda, \beta)=\frac{\left(\frac{x}{\lambda}\right)^\beta}{1+\left(\frac{x}{\lambda}\right)^\beta}=\frac{x^\beta}{x^\beta+\lambda^\beta}, \\ & x>0, \lambda>0, \beta>1, \\ & \xi_\alpha^{\text {LLD}}=\lambda\left(\frac{\alpha}{1-\alpha}\right)^{\frac{1}{\beta}} . \end{aligned} $

假定X=(X₁, X₂, …, X_n)是来自双参数LLD的简单随机样本，令Y_i=-logx_i~$\operatorname{logistic}\left(-\log \lambda, \frac{1}{\beta}\right), i=1, 2, \cdots, n$。接下来将利用矩估计和MLE 2种参数估计法推导出过程能力指数C_Npk的广义枢轴量，该推导的思路是基于文献[38]提出的位置-尺度族中估计的不变性和等方差特性。首先，得到参数λ和β的矩估计是$\begin{aligned} & \hat{\lambda}_{\text {矩 }}=\exp \left(-\frac{1}{n} \sum_{i=1}^n\left(-\log x_i\right)\right) \text { 和 } \\ & \hat{\beta}_{\text {矩 }}=\sqrt{\frac{n-1}{\sum_{i=1}^n\left(-\log x_i-\frac{1}{n} \sum\limits_{i=1}^n\left(-\log x_i\right)\right)^2}}, \text { 又 } Z_i=\frac{-\log X_i-(\log \lambda)}{1 / \beta} \sim \operatorname{logistic}(0, 1), (i=1, 2, \cdots, n) \end{aligned} $，则参数λ和β基于矩估计的广义枢轴量分别为

$ R_\lambda=\exp \left(\frac{1}{R_\beta} \bar{Z}+\frac{1}{n} \sum\limits_{i=1}^n \log x_i\right), $

(10)

$ R_\beta=\sqrt{\frac{\sum\limits_{i=1}^n Z_i^2}{\sum\limits_{i=1}^n\left(-\log x_i-\frac{1}{n} \sum\limits_{i=1}^n \log x_i\right)^2}}, $

(11)

其中Z是Z_i的样本均值。则LLD下过程能力指数C_Npk的广义枢轴量为

$ R_{C_{N_{p l}}}(\lambda, \beta)=2 \min \left\{\frac{\operatorname{USL}-R_\lambda\left(\frac{p_2}{1-p_2}\right)^{\frac{1}{R_\beta}}}{R_\lambda\left(\frac{p_3}{1-p_3}\right)^{\frac{1}{R_\beta}}-R_\lambda\left(\frac{p_1}{1-p_1}\right)^{\frac{1}{R_\beta}}}, \frac{R_\lambda\left(\frac{p_2}{1-p_2}\right)^{\frac{1}{R_\beta}}-\mathrm{LSL}}{R_\lambda\left(\frac{p_3}{1-p_3}\right)^{\frac{1}{R_\beta}}-R_\lambda\left(\frac{p_1}{1-p_1}\right)^{\frac{1}{R_\beta}}}\right\}, $

(12)

这里USL=29, LSL=1, p₁=0.998 65, p₂=0.5, p₃=0.001 35。

下面介绍基于MLE的构造方法。给定X的n个样本观测值x₁, x₂, …, x_n，参数是对数似然函数为$\ln L(\lambda, \beta)=n \ln \beta-n \ln \lambda+(\beta-1) \sum_{i=1}^n \ln \left(\frac{x_i}{\lambda}\right)-\left.2 \sum_{i=1}^n \ln \left[1+\left(\frac{x_i}{\lambda}\right)^\beta\right]\right)$, 似然方程组是

$ \left\{\begin{aligned} \frac{\partial \ln L(\lambda, \beta)}{\partial \lambda}= & \frac{\beta}{\lambda}\left[2 \sum\nolimits_{i=1}^n \frac{\left(\frac{x_i}{\lambda}\right)^\beta}{1+\left(\frac{x_i}{\lambda}\right)^\beta}-n\right]=0, \\ \frac{\partial \ln L(\lambda, \beta)}{\partial \beta}= & \frac{n}{\beta}+\sum\nolimits_{i=1}^n \ln \left(\frac{x_i}{\lambda}\right)- \\ & 2 \sum\nolimits_{i=1}^n \ln \left(\frac{x_i}{\lambda}\right) F\left(\frac{x_i}{\lambda}\right)=0 . \end{aligned}\right. $

(13)

通过求解上述方程组(13)得到尺度参数λ和形状参数β的$\text { MLEs } \hat{\lambda}_{\text {MIE }} \text { 与 } \hat{\beta}_{\text {MIE }}$。因此，参数λ和β的广义枢轴量R_λ^MLE和R_β^MLE分别为$R_\lambda^{\mathrm{MLE}}=e^{\left(\log \hat{\lambda}_{\mathrm{MLE}}-\left(\frac{1}{R_\beta^\mathrm{MLE}}\right) \log \hat{\lambda}_{\mathrm{MLE}}^*\right)} \text { 和 } R_\beta^{\mathrm{MLE}}=\hat{\beta}_{\mathrm{MLE}} / \hat{\beta}_{\mathrm{MLE}}^*$, 这里$\hat{\lambda}_{\mathrm{MIE}}^* \text { 与 } \hat{\beta}_{\mathrm{MLE}}^*$是基于LLD(0, 1)的样本Z₁, Z₂, …, Z_n求得的参数λ和β的MLEs。

因此，得到基于MLE双参数LLD下过程能力指数C_Npk的广义枢轴量为

$ \begin{aligned} & \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;R_{C_{N pk}}^{\mathrm{MIE}}(\lambda, \beta)= \\ & 2 \min \left\{\frac{\operatorname{USL}-R_\lambda^{\mathrm{MLE}}\left(\frac{p_2}{1-p_2}\right)^{\frac{1}{R_\beta^{\mathrm{MIE}}}}}{R_\lambda^{\mathrm{MLE}}\left(\frac{p_3}{1-p_3}\right)^{\frac{1}{R_\beta^{\mathrm{MIE}}}}-R_\lambda^{\mathrm{MLE}}\left(\frac{p_1}{1-p_1}\right)^{\frac{1}{R_\beta^{\mathrm{MIE}}}}}, \right. \\ & \left.\frac{R_\lambda^{\mathrm{MLE}}\left(\frac{p_2}{1-p_2}\right)^{\frac{1}{R_\beta^{\mathrm{MLE}}}}-\mathrm{LSL}}{R_\lambda^{\mathrm{MLE}}\left(\frac{p_3}{1-p_3}\right)^{\frac{1}{R_\beta^{\mathrm{MLE}}}}-R_\lambda^{\mathrm{MLE}}\left(\frac{p_1}{1-p_1}\right)^{\frac{1}{R_\beta^{\mathrm{MLE}}}}}\right\} \end{aligned} $

与上一小节类似，使用Monte Carlo方法计算对数逻辑分布下过程能力指数C_Npk的双边广义置信区间。

设数据为X₁, X₂, …, X_n。

Step 1   产生标准逻辑分布LLD(0, 1)的随机变量；

Step 2   通过式(10)、式(11)与式(12)计算λ和β以及C_Npk的广义枢轴量；

Step 3   将Step 1至Step 2重复操作M次，把得到的M个广义枢轴量按照从小到大的顺序排列，取它的α/2和1-(α/2)分位点作为过程能力指数C_Npk的一个双边置信区间。

3 数值模拟

这一小节通过Monte Carlo数值模拟比较IRD和LLD 2种分布下过程能力指数C_Npk的广义置信区间和文献[11]提出的bootstrap置信区间。在模拟过程中，基于5 000个蒙特卡罗样本计算置信区间，通过2 000次的重复计算区间的覆盖率和平均长度。在样本容量n与参数的多种组合下的仿真结果如表 1和表 2所示。

从表中可以看出：广义置信区间的覆盖率接近名义水平，而当样本量n较小时bootstrap置信区间覆盖率较低。我们提出的广义置信区间方法比bootstrap方法的总体性能更优。对于LLD而言，基于MLE的广义置信区间在平均宽度方面优于基于矩估计的广义置信区间。

4 实例分析

本节使用Zimmer等^[39]讨论的Burr XII可靠性分析的真实数据集来验证所提出的方法。数据观测值记录了在大型制造工厂中用于内部运输和交付的20个电动推车的首次故障时间(以月为单位)：

{0.9, 1.5, 2.3, 3.2, 3.9, 5.0, 6.2, 7.5, 8.3, 10.4, 11.1, 12.6, 15.0, 16.3, 19.3, 22.6, 24.8, 31.5, 38.1, 53.0}

文献[11]利用Kolmogorov-Smirnov (K-S)检验验证了数据符合LLD的合理性。在对数逻辑分布下，应用2种广义置信区间和bootstrap置信区间这3种方法构造基于此数据集过程能力指数C_Npk的置信区间。结果如表 3所示，可以看出本文提出的广义置信区间的长度小于bootstrap方法的区间长度，且基于MLE的广义置信区间长度最短。

5 结论

本文构造了逆瑞利分布和对数逻辑分布下过程能力指数C_Npk的广义置信区间。数值模拟和实例均表明，特别是在小样本容量的情况下，我们的方法优于现有的bootstrap方法。这与许多文献的结论一致，即广义推断方法是解决小样本推断问题的一种优良选择。未来的一个有趣课题是将广义推断方法应用于不同分布下其他过程能力指数的研究中。