用Python实现SPC统计过程控制

SPC与六西格玛

SPC (Statistical Process Control) 统计过程控制，是六西格玛工业管理理论的其中一个重要模块。SPC的控制图 (control chart) 是数据可视化的一个重要手段。而控制图的选择应该根据实际需求来，这里不展开讲控制图，关于控制图的细节可以查找其他资料。（7 种控制图，8 个判异准则。）

简单介绍一下六西格玛，就是 6 sigma 的音译，sigma 是什么？接触过统计的人应该会有印象，$\sigma$ 这个符号就是 sigma，一般代表偏差。

先介绍一下生产质量控制中常用的一个概念 DPMO (Defects Per Million Opportunities)，就是在生产过程中每 100 万个机会中出现的缺陷数。比如在 100 万个焊点里，出现了 1 个缺陷，那么 DPMO 就等于 1。可见 DPMO 是越小越好，最好是为 0。当然在大量的生产过程中这种理想状态出现的几率是非常非常小的。DPMO 的取值会作为衡量生产质量的一个重要指标，下面是 DPMO 对应 sigma 水平的对照表。

sigma level	DPMO	yield
6	3.4	99.99966%
5	230	99.977%
4	6210	99.38%
3	66800	93.32%
2	308000	69.15%
1	690000	30.85%

如果 DPMO 为 3.4 ，对应的 sigma 水平是 6，也就是可以达到 99.99966% 的合格率。所以 sigma 水平就可以作为生产质量的一个评判标准，达到 6 sigma 水平的生产就是非常牛逼的水准了！

6 个 $\sigma$ 的介绍就到这里，不过六西格玛管理理论远不止于此，除了统计学的知识，大部分是关于精益生产管理。

U Chart 案例

这是当时工作中遇到的实际需求，要观察生产产品的单位缺陷数。拿到检测机器的测试记录，包括每天检测的产品编号、检测时间、残次件的维修时间、维修人员等等。冗余信息非常多，用 Excel 处理，只提取检测的产品数量和缺陷件数量，按月份统计整理。（现在看回头当时为什么会弄了两个表示时间的列也是百思不得其解。。。。）n 列表示当月检测的产品数量，c 列表示当月的缺陷数。

因为每个月的产品数量是不同的，所以采用用于可变样本量的 U-chart。如果每个月的产品数量相同，用 C-chart。

Defects per unit: $$u=\frac{c}{n}$$

Central Limit: $$CL=\bar{u}=\frac{\sum c_i}{\sum n_i}$$

Upper Central Limit: $$UCL=\bar{u}+3\sqrt{\frac{\bar{u}}{n_i}}$$

（开方里的是 $\bar{u}$，好像叠起来看不见了）本来还应该算下界的，但是因为 u 肯定是大于 0，且越小越好，所以下界不需要考虑。

用 pandas 导入整理好的数据命名为 defect_test。计算好3个需要展示的数据：u，CL，UCL。

time = defect_test['Month']
u = defect_test['c']/defect_test['n']
cl = sum(defect_test['c'])/sum(defect_test['n'])
ucl = cl + 3*np.sqrt(cl/defect_test['n'])

用 matplotlib 可视化

cl_s = np.ones(12)*cl
plt.style.use('ggplot')
plt.figure(figsize=(15,6))
plt.plot(time,u,'k',marker='o',markersize=10,lw=3,label='u=c/n')
plt.plot(time,cl_s,'r',label='CL')
plt.plot(time,ucl,'b--',label='UCL')
plt.legend()
plt.title('U chart');

从上图可以看出，缺陷数u一直在CL附近波动。大部分的u都在UCL以下，但是有两个月是在控制以外 (out of control)，二月和三月。二月份的u值是0.005277，三月份的是0.006198。对应的DPMO值5277和6198，对应了4 sigma的水平，对应Cpk值（工程能力，要求达到1.33以上）是1.33。总而来说，在二月和三月，缺陷数在控制以外，但Cpk还是达到了应有的期望值1.33。整一年生产过程是达到要求的，但是还是需要注意质量控制，需要保持所有的缺陷数在控制范围内。

对于超出上界限UCL的点判断为异常，在图表上显示出来：

def uchart_test(c,n,time):
    import numpy as np
    import matplotlib.pyplot as plt
    
    u = c/n
    cl = sum(c)/sum(n)
    ucl = cl + 3 * np.sqrt(cl/n)
    cl_s = np.ones(len(time))*cl
    
    # out of control points
    ofc = u[u >= ucl]
    ofc_ind = list(ofc.index)
    ofc_time = time[ofc_ind]
    print("Out of Control:")
    if len(ofc) == 0:
        print("All under control.")
    else:
        for i,j in zip(ofc_ind,range(len(ofc))):
            print(str(j+1) + " - " + str(ofc_time[i]) + ", " + str(ofc[i]))

    plt.style.use('ggplot')
    plt.figure(figsize=(15,6))

    # plot out of control points
    for i in ofc_ind:
        plt.scatter(ofc_time[i],ofc[i],c='r',s=200)

    # plot u chart
    plt.plot(time,u,'k',marker='o',markersize=8,lw=3,label='u=c/n')
    plt.plot(time,cl_s,'g',label='CL')
    plt.plot(time,ucl,'b--',label='UCL')
    plt.legend()
    plt.ylim((0,max