海量过程数据的 CPK 与 PPK 计算

海量过程数据的 CPK 与 PPK 分布式 SQL 计算

正态分布与 6 sigma

正态分布

正态分布，Normal Distribution，常见的分布，基本上能描述所有常见的事物和现象。比如正常人群的身高、体重、考试成绩、家庭收入等等。这里的描述是什么意思呢？就是说这些指标背后的数据都会呈现一种中间密集、两边稀疏的特征。以身高为例，服从正态分布意味着大部分人的身高都会在人群的平均身高上下波动，特别矮和特别高的都比较少见。

这是为什么？中心极限定理中指出，在自然界与生产中，一些现象受到许多相互独立的随机因素的影响，如果每个因素所产生的影响都很微小时，总的影响可以看作是服从正态分布的。
正态分布更准确的描述是如下的公式：
$$f(x;\mu ,\sigma )=\frac{1}{\sigma \sqrt[]{2\pi }}exp(-\frac{(x-\mu {)}^2}{2{\sigma }^2})$$
公式很难懂，但以图表形式表示，就比较直观了，我们只需记住公式中的字母 μ 代表均值，σ 代表方差，也就是 sigma。被评估的指标如果服从正态分布，它的取值会有 68% 的概率落在距离均值一个sigma的范围内，有99.7%的概率落在距离均值 3 个 sigma 内。

6 sigma

如果把范围扩大到6个sigma，则概率会扩大到99.9999998%。6 sigma理论就源自于此，其研究的就是如何控制对结果产生影响的因素，缩小其最终 sigma 的值，将其实际的取值范围控制在允许范围内。比如飞机事故的概率等。

工业过程控制

工业过程控制是极其复杂的一门理论，可以通过以下的表述简单理解。首先，是工序的概念，工业生产中会抽象出很多个工序。工序的划分类似于软件程序，工程师会根据指定的输入输出来将整个生产过程模块化，通过评估每个工序的输入与输出来评估其可控性。但对于工序输出的评估不同于软件程序，在一个生产工艺中，包含很多影响产品质量的因素：操作者，机器，原材料，生产方法，测试方法，生产环境。产品质量就是这些因素的综合表现。因此，对工序的评估需通过正态分布理论，并结合上下限与合格率来进行评估。如果一个工序的输出符合预期的正态分布，则其是可控的；如果一个工序的输入是来自可控的上游工序，但其自身输出给下游的指标不达预期，则其是不可控的。
对于工序的评估举一个具体的例子，一个工序的输出有是一个产品零件，其目标尺寸是200mm，对其的上下限要求是 200±2mm，合格率要求是99%。经过测量得出实际输出的产品零件尺寸均值为 199mm，方差为 0.5mm。则通过上图，我们可以得出其超出下限（198mm）为:(100%-95.44%)/2 = 2.28%。事实上还有可能是超出上限，但是概率为 正态分布数值位于6个σ外的概率，相比于超出下限的概率可忽略不计。因此其合格率约为：97.72%，不符合工序要求。
对于一个固定的生产工艺，产品质量总是分散在一个固定值的周围，通常呈正态分布。从前面的例子我们可以看出，工业过程控制的核心就是评估输出指标的分散区间。如果分散越大，说明生产工艺越不稳定，生产工艺能力越低；相反，如果分散越小，说明生产工艺能力越高。

CP、CPK与PP、PPK

CP(process capability index)、PP(process performance index) 就是通过统计这个分散程度来计算出一个值，来决定和衡量生产工艺能力的。
如果过程中只包含随机因素产生的波动，这个时候使用Cp(Cpk)进行过程能力评估，Cp(Cpk)被称为短期能力指数。如果过程波动由随机因素和特殊因素共同影响而产生，那么使用Pp(Ppk)进行过程能力评估，Pp(Ppk)被称为长期能力指数。Cpk之所以被称为是短期过程能力，是因为当过程受控的情况下，在短期内，过程只会受到随机因素的影响。而Ppk之所以被称为是长期过程能力，是因为在长期情况下(比如一周，一个月等），过程会受到诸如刀具变更，设备换型，员工变更，治具调整，原材料变化等多种因素的影响，过程波动由随机因素和特殊因素共同影响而产生。CP、CPK与PP、PPK的计算公式如下：
$$C_p/P_p=\frac{USL-LSL}{6\sigma }$$
C p k / P p k = m i n { U S L − μ 3 σ , μ − L S L 3 σ } C_pk/P_pk=min\{\frac{USL-\mu}{3\sigma},\frac{\mu-LSL}{3\sigma}\} Cp​k/Pp​k=min{3σUSL−μ​,3σμ−LSL​}
其中：
$LSL$：Lower spec Limit 生产工艺规定的最低值
$USL$：Upper spec Limit 生产工艺规定的最大值
Cp、Pp 与 Cpk、Ppk 的差异在于，前者不关注样本均值与理论均值的差异影响，后者关注更容易超限的一边，也就是均值偏向的一边。
Cp、Cpk 与 Pp、Ppk 的差异在于，当过程受控时，也就是一个工序过程的输出是已知状态时，也就是说其输出指标的均值 μ 与 方差 σ 已知且可信时，使用 Cp、Cpk 指标来评估。当不确定工序是否受控的情况下，只能够使用Pp、Ppk 来进行能力评估。
这也是是两类过程能力的计算公式在写法上一模一样，但其中的 σ 计算公式不一样的原因。

Cp、Cpk 中的 σ 计算公式	Pp、Ppk 中的 σ 计算公式
$\sigma =\overline{R}/d_2$ 子组数量大于1 $\sigma =\overline{MR}/d_2$ 子组数量为1	$$\sigma =\sqrt[]{\frac{\sum (x_i-\overline{x}{)}^2}{n-1}}$$

公式中的 $\overline{R}$ 代表极差的均值，也就是最大值 - 最小值。对于子组为1的情况，取其 $MR$ 的均值，即每两点差值绝对值的均值。CPK的样本容量通常是30~50，而PPK的样本容量是大于100。

基于 YMatrix 的具体计算实现

Cp、Cpk 与 Pp、Ppk 的计算方式成熟，但在大型生产企业的实际工作中，工序复杂，生产的产品众多。如何对大量工序数据进行计算，快速得到结果并实时更新成为工程师的挑战。接下来，本文将演示基于 YMatrix 超融合数据基座实现Cp、Cpk 与 Pp、Ppk 的快速计算。

关于 YMatrix

YMatrix 传承自 PostgreSQL 与 Greenplum ，它们都是世界顶级的数据库项目。其中，PostgreSQL 在近年被Stackflow 评为最受资深程序员欢迎的数据库，没有之一，因为其强大、开放、完善的数据库功能以及良好的 SQL 标准兼容性。而 Greenplum 是基于 PostgreSQL 的 分布式数据库，在 2019 年被 Gartner 评为世界排名第三的数据仓库产品。YMatrix 基于 PostgreSQL 与 Greenplum 进行深度的内核优化，保证了生态丰富的同时，进一步提高数据库在时序、GIS 等多个领域的适应性。它内置4个高性能微内核数据引擎，完善支持标准SQL（从SQL92到SQL2016），支持多模态数据类型，包括关系数据、时序数据、GIS数据、JSON数据、文本数据、图片等；支持多模式操作，包括高并发低延迟增删改查、点查、明细查询、聚合查询、窗口查询、关联查询、多维查询、库内机器学习等。实现了一库多用，全面支持用户的数据管理需求。基于 YMatrix 的数据基座内部构造简洁、操作简便，极大的降低了客户使用数据库产品的门槛。可以让用户更关注业务的价值创造，而不是繁琐的数据性能本身。

AD 验证

Cp、Cpk 的计算前提是抽样数据的正态性检验，即确保样本符合正态分布，再通过公式进行计算。
首先是抽样数据的正态性检验，这里我们采用AD检验的方式进行正态性检验。公式如下：

$$AD=-n-1/n\sum _{i=1}^n(2i-1)[lnF(X_i)+ln(1-F(X_{n-i+1}))]$$

其中：n 为样本数量、F(X) 为指定分布的累积分布函数、i 为数据升序排列时的第i个样本。

1. 示例样本
   首先，我们定义样例数据为抽查精细面粉的装包重量，结果如下所示（单位：kg）：

20.21	19.91	19.99	20.05
19.95	19.99	20.16	20.27
20.15	20.08	20.09	19.96
20.07	20.16	19.97	20.06

将样例数据入库：

CREATE TABLE demo_normal(value float,group int);
INSERT INTO demo_normal VALUES
                            (20.21,1),(19.91,1),(19.99,1),(20.05,1),
                            (19.95,1),(19.99,1),(20.16,1),(20.27,1),
                            (20.15,1),(20.08,1),(20.09,1),(19.96,1),
                            (20.07,1),(20.16,1),(19.97,1),(20.06,1);

2. 计算累计分布
   计算公式中的概率分布函数，由于 Python 比较方便做概率密度的计算，此处引入 plpython 的自定义函数实现，具体如下：

-- 函数注册
DROP FUNCTION if exists 
normcdf(x double precision,loc double precision,scale double precision);
create or replace function 
normcdf(x double precision,loc double precision,scale double precision)
returns double precision as
$$
from scipy.stats import norm
return norm.cdf(x,loc=loc,scale=scale)
$$language plpython3u immutable strict;

-- 查询使用
with mean_and_std as 
(select avg(value) as mean,stddev(value) as std from demo_normal)
select
    value,
    normcdf(value,mean,std) as fxi
from demo_normal CROSS JOIN mean_and_std order by value;

结果如下：

3. 计算$1-F(X_{n-i+1})$
   其是$1-F(X_i)$的逆排序。我们需要一个逆排序加自关联操作来实现：

with
    mean_and_std as 
    (select avg(value) as mean,stddev(value) as std from demo_normal),
    norm_cdf as 
    (select value,
            normcdf(value,mean,std) as fxi,
            row_number() over (order by value) as rn,
            row_number() over (order by value desc ) as rnd
     from demo_normal CROSS JOIN mean_and_std order by value)
select
    nc1.value,
    nc1.fxi,
    nc2.mfxi as mfxid
from norm_cdf nc1 
left join 
(select rnd,1-fxi as mfxi from norm_cdf) nc2 on nc1.rn = nc2.rnd;

结果如下：

4. 计算求和前每一项值
   此处需要两个新的输入，分别是i和n，i即value的序号，n即value的数量，计算过程如下：

with
    mean_and_std as 
    (select avg(value) as mean,
            stddev(value) as std,
            count(value) as number_value 
     from demo_normal),
    norm_cdf as 
    (select value,
            number_value,
            normcdf(value,mean,std) as fxi,
            row_number() over (order by value) as rn,
            row_number() over (order by value desc ) as rnd
     from demo_normal CROSS JOIN mean_and_std order by value),
    fxi_mfxid as 
    (select nc1.value,
            nc1.number_value,
            nc1.rn,
            nc1.fxi,
            nc2.mfxi as mfxid
     from norm_cdf nc1 left join 
     (select rnd,1-fxi as mfxi from norm_cdf) nc2 
     on nc1.rn = nc2.rnd)
select value ,(2*rn-1)*(ln(fxi)+ln(mfxid)) as s_for_sum from fxi_mfxid;

结果如下：

5. AD值计算
最后，将s_for_sum求和，与前面的项计算，得到最终的AD值：

with
    mean_and_std as 
    (select avg(value) as mean,stddev(value) as std,count(value) as number_value from demo_normal),
    norm_cdf as 
    (select value,
            number_value,
            normcdf(value,mean,std) as fxi,
            row_number() over (order by value) as rn,
            row_number() over (order by value desc ) as rnd
     from demo_normal CROSS JOIN mean_and_std order by value),
    fxi_mfxid as 
    (select nc1.value,
            nc1.number_value,
            nc1.rn,
            nc1.fxi,
            nc2.mfxi as mfxid
     from norm_cdf nc1 left join (select rnd,1-fxi as mfxi from norm_cdf) nc2 on nc1.rn = nc2.rnd),
    for_sum as 
    (select (2*rn-1)*(ln(fxi)+ln(mfxid)) as s_for_sum from fxi_mfxid)
select -1*count(s_for_sum) - sum(s_for_sum)/count(s_for_sum) as AD_value  
from for_sum;

结果如下：

6. 置信度计算
通过 AD 值求最终样本是否符合正态分布的置信度 P，根据不同的AD值，有着不同的计算方程，如下所示：
$$p=exp(1.2937-5.709\ast AD+0.0186\ast A{D}^2)AD\ge 0.6$$
$$p=exp(0.9177-4.279\ast AD-1.38\ast A{D}^2)0.34<AD<0.6$$
$$p=1-exp(-8.318+42.796\ast AD-59.938\ast A{D}^2)0.2<AD\le 0.34$$
$$p=1-exp(-13.436+101.14\ast AD-223.73\ast A{D}^2)AD\le 0.2$$

代码实现如下：

with
    mean_and_std as 
    (select avg(value) as mean,stddev(value) as std,count(value) as number_value from demo_normal),
    norm_cdf as 
    (select value,
            number_value,
            normcdf(value,mean,std) as fxi,
            row_number() over (order by value) as rn,
            row_number() over (order by value desc ) as rnd
     from demo_normal CROSS JOIN mean_and_std order by value),
    fxi_mfxid as 
    (select nc1.value,
            nc1.number_value,
            nc1.rn,
            nc1.fxi,
            nc2.mfxi as mfxid
     from norm_cdf nc1 left join (select rnd,1-fxi as mfxi from norm_cdf) nc2 on nc1.rn = nc2.rnd),
    for_sum as 
    (select number_value ,(2*rn-1)*(ln(fxi)+ln(mfxid)) as s_for_sum 
     from fxi_mfxid),
    ad as 
    (select -1*count(s_for_sum) - sum(s_for_sum)/count(s_for_sum) as AD_value 
     from for_sum)
select case
    when AD_value >= 0.6 then exp(1.2937-5.709*AD_value+0.0186*AD_value*AD_value)
    when 0.34 < AD_value and AD_value < 0.6 then exp(0.9177-4.279*AD_value-1.38*AD_value*AD_value)
    when 0.2 < AD_value and AD_value <= 0.34 then 1-exp(-8.318+42.796*AD_value-59.938*AD_value*AD_value)
    when AD_value <= 0.2 then 1-exp(-13.436+101.14*AD_value-223.73*AD_value*AD_value)
    end as P_valu
from ad;

解得p值

由于p>0.05，所以无法拒绝原假设，数据符合正态分布。

CP、CPK、Pp、PpK 计算实现

1. 样本导入
   有一个产品的某个尺寸的公差为5.33±0.05，每20分钟从产线上拿取一个产品进行测量，收集了30个数据，求Cp和Pp。（每次拿取一个产品可视为子组数为1）

   No.	Data	No.	Data	No.	Data	No.	Data	No.	Data
   1	5.343	7	5.304	13	5.314	19	5.308	25	5.280
   2	5.326	8	5.322	14	5.287	20	5.279	26	5.347
   3	5.307	9	5.297	15	5.321	21	5.316	27	5.329
   4	5.325	10	5.326	16	5.289	22	5.361	28	5.301
   5	5.302	11	5.316	17	5.299	23	5.319	29	5.313
   6	5.313	12	5.311	18	5.314	24	5.283	30	5.366

2. 示例数据入库

DROP TABLE IF EXISTS demo_cp_pp;
CREATE TABLE demo_cp_pp(value float,order_id int);
INSERT INTO demo_cp_pp VALUES
    (5.343,1),(5.326,2),(5.307,3),(5.325,4),(5.302,5),
    (5.313,6),(5.304,7),(5.322,8),(5.297,9),(5.326,10),
    (5.316,11),(5.311,12),(5.314,13),(5.287,14),(5.321,15),
    (5.289,16),(5.299,17),(5.314,18),(5.308,19),(5.279,20),
    (5.316,21),(5.361,22),(5.319,23),(5.283,24),(5.280,25),
    (5.347,26),(5.329,27),(5.301,28),(5.313,29),(5.366,30);

3. AD验证
   通过前文方法，对目标表进行AD验证，注意将SQL语句中目标表改成demo_cp_pp，得到P值为0.3268451871113953，大于0.05，目标样本符合正态分布。
4. Cp 中 σ 计算
   让我们回顾Cp 中 σ 的计算公式：

$$\sigma =\overline{R}/d_2$$
       子组数量大于1
$$\sigma =\overline{MR}/d_2$$
       子组数量为1

       其中的d2来自于控制图系数表，d2=1.128：

本示例子组数量为1，取第二个公式，其中 MR 值为移动极差均值，即对于每一列后一项减去前一项的绝对值再去平均：

select avg(mr_value)/1.128 as sigma from (
select value,
       abs(lag(value) over (order by order_id) - demo_cp_pp.value) as mr_value
from demo_cp_pp)mr_cal;

结果为：

5. Pp 中 σ 计算：
Pp 中 σ 即样本标准差，前文已有实现，如下：

select std_dev(value) from demo_cp_pp;

结果为：

6. Cp 、Cpk、Pp、Ppk计算
回顾公式：
其中，
LSL：Lower spec Limit 生产工艺规定的最低值，本例为5.28
USL：Upper spec Limit 生产工艺规定的最大值，本例为5.38
μ为均值，整体实现如下：

select (5.38-5.28)/(6*cp_sigma) as cp,
       (5.38-5.28)/(6*pp_sigma) as pp,
       LEAST((5.38-miu)/(3*cp_sigma),(miu-5.28)/(3*cp_sigma)) as cpk,
       LEAST((5.38-miu)/(3*pp_sigma),(miu-5.28)/(3*pp_sigma)) as ppk
from (select avg(mr_value)/1.128 as cp_sigma,
       stddev(value) as pp_sigma,
       avg(value) as miu
      from (select value,
                   abs(lag(value) over (order by order_id) - 
                       demo_cp_pp.value) as mr_value
            from demo_cp_pp)mr_cal)sigma_cal;

结果如下：

从结果上来看，样例中的情况属于比较危险的情况，其工艺控制的上下限并没有在样本随机分布的较安全区域，通常Cp 与 Pp 大于1.33会比较安全。该问题在图表分析上也可以看出：

上图中，蓝色柱状图是其样本的实际分布，黄色是概率密度分布，两侧蓝色的线即USL与LSL。可见其比较容易发生不合格现象，且更易于发生在LSL这一侧。

持续监测

实际工作中，对上述指标的持续监测会使业务响应与反馈更快，YMatrix 支持持续聚集功能，可以在数据入库时即完成计算，以Ppk为例，实现如下：

1. 目标工序事实表：

CREATE table public.ppk_demo_target
(
    station character varying(128),
    slot character varying(128),
    low_limit double precision,
    hi_limit double precision,
    result_val double precision
);
ALTER table public.ppk_demo_target
    OWNER TO mxadmin;

其中 station、slot为工序标识符，用于实际工作中按工序进行指标统计，low_limit、hi_limit 即对应工序的上下限。result_val为具体值
2. 创建聚集视图来计算μ与σ：

drop VIEW if exists ppk_demo_view;
CREATE VIEW ppk_demo_view WITH (CONTINUOUS) AS
select station,
       slot,
       stddev_samp(result_val) as sigma,
       avg(result_val) as miu,
       min(low_limit) as LSL,
       max(hi_limit) as USL
from ppk_demo_target GROUP BY station,slot order by station,slot;

3. 查询时通过视图中的μ与σ来完成Pp、Ppk的计算。

select (USL-LSL)/(6*sigma) as pp,
       LEAST((USL-miu)/(3*sigma),(miu-LSL)/(3*sigma)) as ppk
select * from ppk_demo_view order by station,slot;

如此一来，涉及到海量数据计算的部分由持续聚集完成，计算实际发生在数据入库阶段，即不会因为上层应用并发查询而增加计算量，同时也保证了查询的响应效率。

引用

Cpk和ppk的本质差别是什么? - 知乎
正态性检验Anderson-Darling在Excel中的实现!
Cp(Cpk)与Pp(Ppk)的解说和区分