Minitab在統(tǒng)計學(xué)專業(yè)實踐課程中的案例拓展

束 慧， 徐宏光

(南京郵電大學(xué) 理學(xué)院，南京210023)

1 引 言

統(tǒng)計學(xué)專業(yè)實踐課程是培養(yǎng)學(xué)生從“知識”到“應(yīng)用”的重要環(huán)節(jié)，在整個統(tǒng)計學(xué)專業(yè)培養(yǎng)方案中具有不可替代的地位.近年來，一些新的、側(cè)重理論學(xué)習(xí)與實踐操作相融合的教學(xué)理念在國內(nèi)高校中被廣泛推廣[1]，越來越多的高校統(tǒng)計學(xué)、數(shù)學(xué)方向的教師開始設(shè)計實踐教學(xué)的新模型、新方法[2-3].作為一款商業(yè)統(tǒng)計軟件，Minitab在統(tǒng)計學(xué)專業(yè)的實踐教學(xué)中具有不可替代的地位，尤其是涉及到6σ管理等質(zhì)量控制課程的實踐時，其地位更加凸顯，圍繞Minitab設(shè)計的實踐教學(xué)案例陸續(xù)涌現(xiàn)[4].

Minitab是由美國賓夕法尼亞州立大學(xué)開發(fā)的一款商用統(tǒng)計軟件，其特色功能為質(zhì)量控制模塊，被廣泛運用于工業(yè)企業(yè)的產(chǎn)品質(zhì)量控制過程[5].以Minitab17版本為例，該軟件的質(zhì)量控制模塊由控制圖、質(zhì)量控制兩個子模塊組成.在控制圖模塊中，整合了子組變量、單值變量、屬性、時間加權(quán)、多屬性變量以及稀有事件控制圖等各類控制圖的繪制功能.在運用Minitab繪制各類控制圖時，控制限的選擇是開放給用戶的，提供了“標(biāo)準(zhǔn)差的倍數(shù)”或者指定“控制上、下限”的權(quán)限[5].這樣的功能設(shè)計，具有操作上的簡便性，但同時也增加了用戶在繪制控制圖時的隨機(jī)性，尤其是繪制一些檢測生產(chǎn)過程中的中小變異的控制圖，如指數(shù)加權(quán)移動平均(Exponentially Weighted Moving Average，EWMA)[6-9]，累計和CUSUM (Cumulative Sum)[10-11]等控制圖時，控制限設(shè)置存在一定的優(yōu)化空間.

本文從統(tǒng)計學(xué)專業(yè)實踐教學(xué)案例拓展的角度出發(fā)，對控制圖的重要優(yōu)化指標(biāo)進(jìn)行以Minitab為基礎(chǔ)的算法設(shè)計和計算，并根據(jù)運算結(jié)果對控制圖的關(guān)鍵參數(shù)進(jìn)行優(yōu)化選擇.本文的研究既可以作為統(tǒng)計學(xué)專業(yè)實踐課程的案例拓展，也可以作為Minitab統(tǒng)計軟件質(zhì)量控制模塊功能優(yōu)化研究，可以為統(tǒng)計學(xué)專業(yè)實踐課程提供實踐內(nèi)容和實踐工具兩個方面的參考.

本文第二部分選取典型的控制圖，簡要介紹其優(yōu)化指標(biāo)的計算方法，以及基于指標(biāo)的控制圖優(yōu)化原理；以理論方法為基礎(chǔ)，設(shè)計控制圖優(yōu)化指標(biāo)在Minitab上的實現(xiàn)算法，并提供一個案例計算結(jié)果；第三部分是本文的分析總結(jié)，基于文章的結(jié)果提供一些Minitab軟件在高校統(tǒng)計學(xué)專業(yè)實踐教學(xué)環(huán)節(jié)中的應(yīng)用拓展建議.

2 實踐案例構(gòu)建——基于Minitab 的控制圖優(yōu)化

2.1 控制圖優(yōu)化理論模型

(1)

在進(jìn)行質(zhì)量檢驗過程中，如果抽樣的時間間隔及樣本容量n固定不變，則稱相應(yīng)的控制圖為固定抽樣率(Fixed Sampling Rate，F(xiàn)SR)控制圖.對于一個FSR控制圖，定義從檢測開始到它發(fā)出生產(chǎn)變異的報警標(biāo)識為止，抽取的平均樣本組數(shù)為平均運行長度(Average Run Length， ARL).當(dāng)過程可控時，質(zhì)量控制圖的報警就屬于誤報，可控的ARL越大越好；當(dāng)過程失控時，質(zhì)量控制圖應(yīng)盡早報警，失控的ARL越小越好.

2.2 控制圖優(yōu)化指標(biāo)ARL的計算方法

控制圖ARL的計算方法多樣，研究較豐富，運用有限狀態(tài)Markov鏈的平均首達(dá)時間計算ARL是比較常見的[6，7，12，13].這部分內(nèi)容結(jié)合統(tǒng)計學(xué)專業(yè)實踐環(huán)節(jié)教學(xué)要求，可以設(shè)置為對學(xué)生的統(tǒng)計學(xué)理論基礎(chǔ)的考察項目.

在上述前提下，控制上限可以表示為

UCL=L·σs=h，

(2)

相應(yīng)的控制下限

LCL=-h，

(3)

其中，L為控制限參數(shù)，是優(yōu)化控制圖時的重要選擇參數(shù).

ARL計算的基本依據(jù)是有限狀態(tài)Markov鏈從非吸收態(tài)首次轉(zhuǎn)移到吸收態(tài)的平均時間，因此，需要對控制圖界限進(jìn)行Markov鏈狀態(tài)空間轉(zhuǎn)換.將控制圖的上下控制限之間的區(qū)間分成k(奇數(shù))個子區(qū)間，小區(qū)間的寬度d

(4)

Markov鏈的狀態(tài)空間可以用{Ij}j=1，2，…，k表示，且設(shè)Ij為小區(qū)間的中心點

(5)

當(dāng)EWMA檢驗統(tǒng)計量St落在第j個區(qū)間內(nèi)時，即

就定義統(tǒng)計量St處于狀態(tài)Ij.EWMA控制圖是雙邊的，將大于上控制限h的區(qū)域與小于下控制限-h的區(qū)域這兩個吸收域合并為一個.在上述設(shè)置下，該Markov鏈的一步轉(zhuǎn)移概率pij為

pij=P{St+1∈Ij|St∈Ii}

phj=0，j=1，2，…，k;

phh=1.

(6)

檢測統(tǒng)計量St在m時刻的一步轉(zhuǎn)移概率陣為

(7)

矩陣P的最后一列表示從狀態(tài)Ii出發(fā)一步轉(zhuǎn)移到吸收狀態(tài)Ih的概率，因此可記

以fih表示從狀態(tài)Ii出發(fā)第一次轉(zhuǎn)移到吸收狀態(tài)Ih所需要的步數(shù)(首達(dá)時間)，fh=(f1h，…，f(h-1)h)T，則當(dāng)檢測統(tǒng)計量的初始態(tài)為Ii時的ARL值即為E(fh)的第i個分量，特別地，初始狀態(tài)是從中心點出發(fā)，則為中間位置的分量.根據(jù)有限狀態(tài)馬氏鏈的平均首達(dá)時間計算方法，有

ARL=E(fh)=(I-R)-11，

(8)

2.3 基于Minitab的ARL計算及控制圖優(yōu)化

根據(jù)上述理論模型，控制圖優(yōu)化的基本思路是：給定變異系數(shù)δ，狀態(tài)空間劃分k，優(yōu)化參數(shù)：平滑系數(shù)λ、控制限L，使得控制圖在受控態(tài)下(變異系數(shù)δ=0)的ARL較大，而失控狀態(tài)(變異系數(shù)δ≠0)下的ARL較小，同時分析狀態(tài)空間劃分對相應(yīng)的ARL優(yōu)化的影響.

借助Minitab的數(shù)據(jù)整理、數(shù)據(jù)(矩陣)計算功能，以Minitab基于列變量的運算規(guī)則為語句設(shè)計特點，在Minitab宏命令“GMACRO”功能框架下，定義ARL計算過程腳本，設(shè)計基于Minitab的控制圖優(yōu)化指標(biāo)ARL計算及反饋優(yōu)化機(jī)制(圖1).

圖1 基于Minitab的EWMA控制圖優(yōu)化指標(biāo)計算及反饋

Step1 設(shè)置不同的控制限L，偏移系數(shù)δ，狀態(tài)空間k，平滑系數(shù)λ值，計算相應(yīng)的LCL，d;

Step2 在計算轉(zhuǎn)移概率矩陣時，借助Minitab的“Mesh”命令，可以產(chǎn)生計算轉(zhuǎn)移概率矩陣的網(wǎng)格化指標(biāo)(i，j)，替代一般編程語言中的循環(huán)結(jié)構(gòu);

Step3 按列存儲(i，j，LCL，L，λ，d)，其中的(LCL，L，λ，d)利用Minitab的產(chǎn)生模板數(shù)據(jù)，和對應(yīng)的網(wǎng)格數(shù)據(jù)對一致， 以(i，j，LCL，L，λ，d)值為輸入，計算相應(yīng)的轉(zhuǎn)移概率值pij;

Step4 利用Minitab的拆分列“Unstack”、列轉(zhuǎn)置“Transpose”，以及數(shù)據(jù)復(fù)制“Copy”功能，將列變量存儲為轉(zhuǎn)移概率矩陣R；

Step5 計算ARL值，利用Minitab的“Diagonal”產(chǎn)生單位矩陣，“Subtract”矩陣運算，“Invert”矩陣求逆，“Multiply”矩陣乘積等宏命令；

Step6 (反饋機(jī)制)記錄不同的偏移系數(shù)δ、狀態(tài)空間劃分k，控制限L，平滑系數(shù)λ下的ARL計算結(jié)果，分組擬合ARL與k，λ值的關(guān)系曲線，從中篩選出符合要求的優(yōu)化參數(shù)設(shè)置.具體，根據(jù)EWMA控制圖檢測中小偏移為主的特征，設(shè)置偏移系數(shù)δ=0，δ=1，分別擬合出ARL與(k，λ)的等值線圖(圖2)； ARL與k的關(guān)系曲線圖(圖3)；ARL與λ的關(guān)系曲線圖(圖4).以上圖形都是用Minitab作圖功能繪制.

(k，λ)組合對ARL的影響是綜合性的：δ=0時，不同控制限制L的設(shè)置下，ARL取值隨著λ取值的變化，始終呈現(xiàn)“梯度”特征，且在平滑系數(shù)λ取值較低的階段，ARL均較長.這樣的結(jié)果說明了過程受控時，若平滑系數(shù)λ取值較低(0.1附近)，EWMA控制圖可以減少“誤報”(圖2).

δ=1時， ARL取值隨著λ取值的變化雖然仍呈現(xiàn)“梯度”特征，但在不同控制限制L的設(shè)置下，ARL與平滑系數(shù)λ取值的變化并沒有呈現(xiàn)相似的變化趨勢.L=3，平滑系數(shù)λ取值較低區(qū)間[0.1，0.4]，ARL較短；L=2，平滑系數(shù)λ取值落在[0.4，0.6]區(qū)間，ARL較短.這樣的結(jié)果說明了過程非受控時， EWMA控制圖可以盡早“報警”，但需要將控制限和平滑系數(shù)進(jìn)行動態(tài)的調(diào)整(圖2).

圖2 控制圖受控、失控時，ARL與(k，λ)的等值線圖

固定(k，λ)中的一個，ARL與另一參數(shù)間的關(guān)系則更加特征鮮明：

δ=0時，不同控制限制L的設(shè)置下，ARL受狀態(tài)空間k的影響較小，僅僅在平滑系數(shù)λ取值較低(0.1)時，出現(xiàn)一定的影響作用(圖3(a)); δ=1時，不同控制限制L的設(shè)置下，ARL受狀態(tài)空間劃分k的影響均較小(圖3(b)).

(a) 受控時，ARL與k的關(guān)系擬合圖 (b) 非受控時，ARL與k的關(guān)系擬合圖圖3 ARL與k的關(guān)系擬合圖

δ=0時，不同控制限制L的設(shè)置下，ARL受平滑系數(shù)λ設(shè)置的影響較大，和等值線圖顯示的結(jié)果一致，平滑系數(shù)λ取值較低(0.1)時，ARL均較長(圖4(a)).

(a) 受控時，ARL與λ的關(guān)系擬合圖 (b) 非受控時，ARL與λ的關(guān)系擬合圖圖4 ARL與λ的關(guān)系擬合圖

δ=1時，不同控制限制L的設(shè)置下，ARL隨著平滑系數(shù)λ的取值不同，呈現(xiàn)差異化的變化趨勢：控制限較大，ARL隨著λ的取值變大而變長；控制限中等，則ARL隨著λ的取值變大而先變短后又變長；控制限較小時，ARL隨著λ的取值變大而變短(圖4(b)).

綜合以上仿真分析結(jié)果，遵循控制圖優(yōu)化的基本原理，受控時，EWMA控制圖應(yīng)選取較小的平滑系數(shù)值；而當(dāng)過程非受控時，應(yīng)結(jié)合控制限L的寬度，相應(yīng)的平滑系數(shù)值也要作出變化：L寬，則λ小；L中等，則λ中等；L窄，則λ大.這樣的取值關(guān)系，可以使得受控時的ARL較長，非受控時ARL較短，實現(xiàn)控制圖優(yōu)化的目標(biāo).

3 結(jié) 論

本文以統(tǒng)計軟件Minitab為基礎(chǔ)，設(shè)計了一類統(tǒng)計學(xué)專業(yè)實踐教學(xué)中的拓展案例.該案例以質(zhì)量控制領(lǐng)域的控制圖優(yōu)化理論為基礎(chǔ)，結(jié)合Minitab的數(shù)據(jù)處理特征、功能模塊以及宏命令語言，實現(xiàn)了控制圖的仿真優(yōu)化.該案例完整、高效的實現(xiàn)過程說明，功能性統(tǒng)計軟件在高校統(tǒng)計學(xué)專業(yè)實踐教學(xué)中具有非常大的潛力，這類軟件的教學(xué)推廣，是對廣泛流行的SAS、SPSS、R以及Stata等主流統(tǒng)計軟件的有益補(bǔ)充.

致謝作者非常感謝相關(guān)文獻(xiàn)對本文的啟發(fā)以及審稿專家提出的寶貴意見.