行政院國家科學委員會專題研究計畫 成果報告
監控製程變異之 SPC 方法(2/2)
計畫類別： 個別型計畫 計畫編號： NSC92-2118-M-009-012- 執行期間： 92 年 08 月 01 日至 93 年 10 月 31 日 執行單位： 國立交通大學統計學研究所 計畫主持人： 洪志真 計畫參與人員： 顏家鈴、黃立芬、徐雅甄、翁崢珮 報告類型： 完整報告 報告附件： 國外研究心得報告 處理方式： 本計畫可公開查詢中 華 民 國 94 年 1 月 26 日
中文摘要 本研究主要探討如何監控多變量製程變異增加的問題。我們利用多變量單邊檢定建 構偵測製程變異增加的控制圖。考慮H0 :Σ=Σ0 vs. H1:Σ≥Σ0 且 Σ≠Σ0，其中Σ為 所監控品質特性的共變異矩陣及 為在控制狀態下的製程變異，並且分 已知或未知 兩種情形來討論。 0 Σ Σ0 本研究計畫為兩年期計畫，在第一年的研究中，我們導出Σ 已知或未知兩種情形0 下之概似比檢定統計量，並用靴環法得出管制界限。針對此控制圖之績效問題，我們經 統計模擬對幾種 的變化比較平均連串長度，並以一個實例和模擬例子，證實所提出的 單邊檢定方法對多變量製程變異增加的問題在偵測能力上有相當不錯的效率，且與允許 多變量製程變異性可增加或減少之雙邊檢定方法作比較，也有較佳的偵測效率。 Σ 在第二年的研究中，我們接續第一年的研究，對上述概似比統計量提出 EWMA 管 制圖，並模仿 Crowder (1987)所提出計算平均連串長度的方法進而找出管制界限。針對 不同的 EWMA 的權重λ，以模擬實驗來探討 EWMA 管制圖在多變量製程變異增加時之 平均連串長度表現。與 Shewhart 管制圖比較，發現 EWMA 管制圖在變異小至中幅度增 加下有較佳的偵測效率。 關鍵詞：多變量控制圖、製程變異、單邊概似比檢定、靴環法、平均連串長度、Shewhart 管制圖、EWMA 管制圖
英文摘要
In this project, we consider the problem of detecting variation increases in multivariate processes. Consider the one-sided hypotheses testing H0:Σ = Σ 0
, where is the covariance matrix of the multivariate quality characteristic of interest and is the corresponding in-control covariance matrix.
versus H1:Σ ≥ Σ0 and 0
Σ ≠ Σ Σ
0 Σ
In the first year of the project, we derive the likelihood ratio statistic and propose two Shewhart-type control charts for the cases of Σ known and unknown, respectively. A 0 comparative simulation study further shows that the proposed control charts outperform the control charts based on the two-sided likelihood ratio statistic in terms of the average run length. The applicability and effectiveness of the proposed control charts are demonstrated through a real example and two simulated examples.
In the second year of the project, we propose and study the EWMA version of the Shewhart-type charts derived from the first year project. The performances of the control charts are studied via the average run length behavior. The effect of the EWMA parameter λ is also studied. Simulation studies show that the EWMA control charts are more sensitive than their Shewhart-type counterpart for small to moderate increases on the process variability.
Keywords: Multivariate control charts, Process variation, One-sided likelihood ratio test,
Bootstrap, Average run length, Shewhart-type control charts, EWMA control charts
目 錄
一、前言………1 二、研究目的………1 三、文獻探討………2 四、研究方法與結果………3 五、結果與討論………10 參考文獻………11報告內容
一、前言 現今經濟不景氣，製造商對產品良率的要求十分的重視，每一個不良品對製造商來 說都是成本上的浪費。為了能以最低的成本創造最高的利潤，每個製造商都希望能在最 短的時間內找出製程問題的發生原因並排除之，以減少成本的浪費及提高良率，因此對 於品質管制十分重視。對品質管制而言，統計製程管制 (Statistical Process Control，簡稱 SPC) 一直是個 十分重要的課題，利用分析樣本資料來判斷製程是否處於穩定的狀態，是監控製程品質 非常重要的一個方法。 在生產過程中，總是會有一些變異，大致分成兩類：一般變異 (common causes) 及 特殊變異 (special causes) 。一般變異是不管設計得再好總是會存在的且不易消除，例 如：環境因素、機器本身特性等；而特殊變異是偶而才會出現的但對製程的影響卻很大， 例如：不良的材料、不當的機器操作等。SPC 的主要目的之一就是希望能迅速偵測出特 殊變異的發生及製程的改變，而目前管制圖 (control charts) 是最常用來監控製程的工 具，能盡快的偵測出特殊變異以便採取改善行動，藉以提昇產品之品質。所以，管制圖 的選擇及應用，成為 SPC 中非常重要且需要探討的課題之一。 在過去品質管制研究中有許多有關監控製程否失控狀態 (out of control)之研究，其 中大多是針對單一品質特性的討論。但在一個製程所需監控的品質特性有兩個或兩個以 上時，可能需要使用多變量的方法來監控。 如果個別對此 p ( )個品質特性作管制圖監控時，很容易會有誤導的情形。假設
針對每一品質特性管制的假警報率(false alram rate)為 2 p≥ α ，亦即資料顯示失控狀態而事實 上並無特殊變異的情況發生之機率為α ，如果分別對此 p 個品質特性作管制，而且假設 這 p 個品質特性是獨立的，那麼整體而言假警報率為 。此時雖然可以經由修 正 p ) 1 ( 1− −α α 大小得到想要的效果；但當此 p 個品質特性並不獨立時則正確的假警報率很難得 知，故應考慮用多變量的方法來監控一多變量製程。 二、研究目的 考慮有 p 個相關的品質特性之多變量製程，從線上作業中隨機抽取一組合理子群 (rational subgroups) 的 樣 本 ， 假 設 此 樣 本 的 每 一 觀 察 值 均 為 多 變 量 常 態 分 配 ) , (µ Σ p N ，其中平均值向量μ和共變異矩陣Σ是未知的。過去大多只討論平均數向量 的監控，但是失控時我們很難分辨出是平均數向量，還是共變異數矩陣發生改變。要改 善品質，往往偵測製程變異比偵測平均數來得重要，因而近年來對於共變異數矩陣的監 控也越來越重視，也有少數討論這方面的論文。然而變異有可能增加或減少，但在大部
分的情況下我們比較擔心變異增加，變異的增加會直接影響到產品的品質，因此本研究 的主要目的在於探討如何監控一多變量製程之變異是否在統計管制(in statistical control) 狀態，且針對變異是否增加的情形做討論。針對此問題，我們在第一年的研 究中，提出一個 Shewhat 式之管制圖並探討此控制圖之績效問題。 現今製程，追求的品質要求較高，對於製程小小的變動都十分在意。在偵測製程平 均值方面，傳統的 Shewhart 管制圖已有許多文獻指出其偵測效率在小變動較不靈敏，而 累和管制圖(Cumulative Sum Control Chart，簡稱 CUSUM 管制圖)或指數加權移動平均管 制圖 (Exponentially Weighted Moving Average Control Chart，簡稱 EWMA 管制圖)在小變 動時較敏感，因此在第二年的研究中，我們接續第一年的研究，對上述概似比統計量提 出 EWMA 管制圖來偵測製程的變動，看是否較 Shewhart 式之管制圖靈敏。 三、文獻探討 在以往大部分的多變量品質管制研究中，多是致力於製程平均數向量μ的控制，很 少把焦點放在監控共變異矩陣Σ，這是由於Σ所牽涉到的分布理論很複雜及監控Σ的重 要性較不受到重視，但近年來已有較多的研究在監控多變量製程變異上。 Alt (1985) 提 供 一 個 監 控 多 變 量 製 程 變 異 的 方 法 ， 重 點 放 在 0 1 0 0 :Σ=Σ vs. H :Σ≠Σ H 的雙邊檢定，其中Σ為所監控品質特性的共變異矩陣及Σ0
為在控制狀態下已知的製程變異。Alt and Bedewi (1986)提出兩個監控Σ之管制圖，一種 管制圖是以概似比準則為基礎，利用概似比檢定統計量；另一種管制圖是利用樣本的廣
義變異數來當成多變量製程變異性的測度，即是使用樣本共變異矩陣 的行列式S S ，稱
之為 S 管制圖。
Tang and Barnett (1996a,1996b)在考慮合理子群及管制值Σ 是已知或未知兩種情形0
下，提出對多變量製程變異的監控方法。在雙邊檢定中，對Σ 為已知或未知兩種情形0
分別提出各自的檢定統計量，並且以模擬方法去比較他們所提方法的檢定力比現存的其 他方法好。
Yeh, Huwang, and Wu (2004)引進一新的管制圖—指數加權移動概似比管制圖 (Exponentially Weighted Moving Likelihood Ratio Control Chart)來偵測多變量製程變異性 的變化。此法乃是利用檢定兩個獨立母體的共變異矩陣是否相等的概似比檢定統計量， 取自然對數後所建構的指數加權移動平均管制圖。文中經由實際例子和統計模擬證實所 提出的管制圖對多變量製程變異性的變化有相當不錯的偵測效率。 上述文獻均是利用雙邊檢定來偵測多變量製程變異性的變化，Σ的增加或減少均是 偵測的對象；然而在大部分的製程中，我們比較擔心的是Σ的增加而不是減少。在此情 形下共變異矩陣的單邊檢定對監控製程失控(out of control)狀態應該會比雙邊檢定來得 敏感。所以本研究探討如何來監控多變量製程之變異是否在統計管制(in statistical control) 狀態的問題，但我們的焦點是放在製程的變異是否有增加。
考慮單邊檢定 H0 :Σ=Σ0 vs. H1:Σ≥Σ0 且 Σ≠Σ0 ， (3.1) 其中Σ≥Σ0代表Σ−Σ0為半正定矩陣。以下是有關共變異矩陣Σ 的單邊檢定的文獻。 Calvin (1994) 討 論 Σ0 已 知 的 情 形 。 他 將 單 邊 檢 定 0 0 1 0 0 :Σ=Σ vs. H :Σ≥Σ 且 Σ≠Σ H ，分成下列兩個步驟的檢定： (1) H0 :Σ=Σ0 vs. H1:Σ≥Σ0 (2) H1:Σ≥Σ0 vs. H2 :Σ≠Σ0 此結果與(3.1)之單邊檢定是同義的。但對這種兩步驟的方法，我們必須做兩個管制圖來 監控；若第一步的管制圖及第二步的管制圖全部在控制外，我們才說製程變異增加，所 以 Calvin 的方法在製程監控上比只做一個管制圖來得複雜，因此我們不考慮使用此法。 Sakata (1987)中有討論管制值Σ 未知的情形，文中對兩個多維常態母體導出檢定0 (3.1)之概似比檢定統計量且在兩個母體樣本數相等的情況下，導出檢定統計量的近似分 布，但檢定統計量的分布與近似分布很複雜，所以在管制值Σ 為未知時我們亦不考慮此0 方法。 另外，要提四篇文獻，其結果與我們要檢定(3.1)所得到的概似比檢定統計量很相 似，主要差別在於他們所用的模型與我們不相同。Anderson, Anderson, and Olkin (1986) 中，是針對有隨機效應的多變量單因子模型且在有重覆(replication)的平衡(balanced)情形 中，當效應及誤差為多變量常態分布且效應的共變異矩陣Θ 為半正定及 Θ 具有最大秩 (rank)的條件下，找Θ 的最大概似估計量 ，假設Θˆ rank(Θ)=k，其假設檢定為 1 0 1 0 0 :k k vs. H :k k k H ≤ < ≤ ， (3.2) 其中 為給定的數。Anderson, Anderson, and Olkin (1986)導出(3.2)的概似比檢定統計 量 1 0, k k λ及Θˆ。Anderson (1989) 針對 (3.2)導出當p=2時T =−2logλ的近似分布Y。若檢 定是否有隨機效應，也就是當k0 =0和k1 = p時，即是檢定 0 : 0 : 1 0 Θ= versus H Θ≠ H 。 (3.3) Anderson (1989)也寫出(3.3)的概似比檢定統計量，與我們所要導出的概似比檢定統計量 表現型態很相似。Amemiya, Anderson, and Lewis (1990)中利用蒙第卡羅模擬方法對 找出(3.2)之檢定統計量的近似分布之分位數(quantile)估計值的表。而在 Kuriki (1993)中則對一般性的 導出(3.2)之檢定統計量之近似分布，並且給 到 之 近似分布分位數的表。 , 3 ≥ p , p p=2 p=14 四、研究方法與結果 當製程在管制中，我們考慮製程的平均數向量μ的管制值µ0為未知，而共變異矩 陣Σ的管制值為 。我們利用概似比檢定做為決策法則，分別導出當管制值 為已知 及未知情形時(3.1)的檢定統計量。 0 Σ Σ0 假設在控制之下之共變異矩陣為Σ 且為已知。我們在時間 時，從線上作業中隨機0 t
抽取的一組個數為 之合理子群，n {Xtj,j=1,",n}，為一組取自 p維常態分配Np(µ ,Σ) 的隨機樣本，其中μ、Σ為未知，我們利用此樣本來檢定製程變異是否有增加。首先定 義： ) )( ( , 1 1 1 ′ − − Σ = Σ = = = tj t tj t n j t tj n j t X B X X X X n X ，t=1,2,"，
則 Bt具 有 自 由 度 為n−1， 尺 度 矩 陣 (scale matrix) 為Σ 之 Wishart 分布，我們記為 。亦即 ) , 1 ( ~W n− Σ Bt p Xt和St =Bt n為時間 線上作業中 個觀察值的樣本平均數和樣 本共變異矩陣。此處 除以 而不是除以 t n t S n n−1的原因是要避免證明的複雜性。為簡單 計，假設St為正定矩陣。我們證明此檢定的概似比統計量為
[
]
⎪⎩ ⎪ ⎨ ⎧ = > − − Π = = 0 , 1 0 , } ) 1 ( exp{ * * 2 1 * p p d d n i i p i 若 若 λ 其中d1 ≥"≥dp >0為 St − dΣ0 =0的根， 為 的個數。令檢定統計量為 * p di >1 λ log 2 − = T[
]
⎪⎩ ⎪ ⎨ ⎧ = > − − Σ = = 0 , 0 0 , log ) 1 ( * * 1 * p p d d n i i p i 若 若 (4.1) 假設µ0和 為未知。因此我們必須在製程為管制中時，從製程中取得 m 組樣本大 小為 n 的訓練樣本(training sample)來估計 0 Σ 0 µ 和Σ 。假設訓練樣本0 ，i=1,2,…,m, j=1,2,…,n ， 是 m 組 取 自 常 態 分 配 ij X ) , (µ0 Σ0 p N 的 隨 機 樣 本 。 定 義 ： ) )( ( , 1 1 1 1 =Σ Σ − − ′ A Σ Σ = = = = = X X X X mn X X ij ij n j m i tj n j m i ，其中X 和 S0 = A (mn)為訓練樣本中 個 觀 察 值 的 樣 本 平 均 數 和 樣 本 共 變 異 矩 陣 。 則 mn A 具 有 Wishart 分 布 ， 。 除以mn的原因與前同。為簡單計，假設 為正定矩陣。我們 導出一個新的概似比檢定統計量 ) , 1 ( ~W mn− Σ0 A p S0 S0 ⎪ ⎩ ⎪ ⎨ ⎧ = > ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − + Π = + = 0 , 1 0 , ) 1 ( * * 2 1 * p p w w n mn i w i p i 若 若 β β λ 其中β1 ≥"≥βp >0為 St − Sβ 0 =0的根， 1 1 + = m w ，p*為 βi >1的個數。如前，令 檢定統計量λ log 2 − = T ⎪ ⎪ ⎩ ⎪⎪ ⎨ ⎧ ≤ > − − + Σ + = = 0 , 0 0 , ] log ) 1 [log( ) ( * * 1 * p p w w w n mn i i p i 若 若 β β (4.2) 現利用上述(概似比)檢定統計量來建構單邊檢定管制圖。因為管制統計量 T 是非負 的，所以我們只考慮管制上限(upper control limit, UCL)，令管制下限(lower control limit,
LCL)等於 0。因此若管制統計量 T 大於管制界限，即表示製程失控。至於如何求得管制 界限，因為管制統計量 T 的分布不好求，我們利用模擬的方法來求得管制界限。 若製程在穩定狀態下，控制變異的管制界限之產生是跟 及 的特徵值 (eigenvalues) 有關(或 及 的特徵值 n p, StΣ0−1 , 1 d ",dp p,n,m StS0−1 β1,",βp有關)，但我們可以證 明，不論多變量製程的共變異矩陣是為何，所產生的管制界限皆會相等，因此我們可藉 由Np(0, Ip)的多變量常態製程來產生控制變異的管制界限。證明如下： 爲簡單計，假設St,Σ0均為正定。令d1 ≥"≥dp >0為 St − dΣ0 =0的根。因為Σ 為0 對稱正定矩陣，所以唯一存在一個對稱正定矩陣 使得 ，我們將 記為 。令 2 / 1 0 Σ Σ0 =(Σ10/2)(Σ10/2) 1 2 / 1 0 ) (Σ − 01/2 − Σ Ztj =Σ0−12Xtj，則 ，j=1,2,…,n，是一組取自常態分配 的 隨 機 樣 本 ， 且 tj Z Np(0,Ip ) t tj n j tj n j t X X n Z n Z 01/2 01/2 1 1 ( ) 1 1 − − = = = Σ Σ =Σ Σ = ， 2 / 1 0 2 / 1 0 1 ) ( ) )( ( 1 − − = − − ′=Σ Σ Σ = n tj t tj t t j z t Z Z Z Z S n S 分別為 n 個觀察值的樣本平均數和樣本共變 異矩陣，則Xt 1/2Zt 0 Σ = 和 1/2。令 0 ) ( 2 / 1 0 Σ Σ = z t t S S λ1 ≥"≥λp >0為 0 ) ( − = p z t I S λ 的根， 而 由 1/2 0 2 / 1 0 2 / 1 0 ) ( 2 / 1 0 0 = Σ Σ − Σ Σ Σ −d S d St tz = Σ10/2(St(z)−dIp)Σ10/2 = Σ0 St(z)−dIp ， 可 知 0 0 = Σ − d St 與 St(z) −λIp =0有相同的根，亦即di = ，λi i=1,2,"p，也就是 的 p 個特徵值與 的 p 個特徵值是相同的。所以，製程在控制下，由 的特徵值與 的特徵值分別所產生的 T 分布是相同的，故 T 分布是不隨多變量製程共變異矩陣不同而 有所改變。 1 0 − Σ t S ) ( z t S StΣ0−1 St( z) 以下我們分別以Σ 為已知及未知情形來討論如何求得管制界限。令第一型誤差機0 率為α。 我們利用靴環法來找管制界限，其模擬的步驟如下： 1. 給定p, n,α 的值
2. 生成一組大小為 的新樣本，且這 個觀察值皆服從 ，並由(4.1)計算 統計量 。 n n Np(0,Ip) T 3. 重覆步驟 2 共N次，然後找出這N個 值的第(T 1−α)樣本分位數。 4. 重覆步驟 3 共 b 次，然後求此 個第(b 1−α)樣本分位數的平均值。此平均值即 定為在步驟 1 狀況下所模擬估計的管制界限CLp,n,α 。 對 p=2,3,4,5個品質特性及樣本大小n=5,10,15,20,25和α =0.05, 0.01,0.0027之 組合下，各模擬N=1,000,000 次及 b =100 次，算出在不同p, n,α 組合下的管制界限及 管制界限的標準差。從結果可以觀察到： 對相同的 、p α，當 n 愈大時，管制界限 愈大，這是因為檢定統計量 T 與 n 成正比。 α , ,n p CL 對相同的 n、α，當p愈大時，管制界限CLp,n,α愈大。 (1-α)愈大，管制界限的標準差會愈大。 當 為未知時，與前所述相同，在製程控制下，T之分布是不隨多變量製程共變異 矩陣的不同而有所改變，因此，在不失一般性下，我們可以藉由 的特徵值所 產生之 T 分布來決定某特定切點，當成管制界限。同樣的，我們可利用靴環法來找管制 界限，其模擬的步驟類似，只有在步驟 1 多給定 m 值、步驟 2 從 生成 組樣 本大小為 的 個訓練樣本，並由(4.2)計算統計量T，如此即可求得管制界限 。 0 Σ 1 ) ( 0 ) ( ) ( z − z t S S ) , 0 ( p p I N m n mn CLp,m,n,α 對 個品質特性、訓練樣本組數 m=25, 30, 35, 40, 45, 50,55, 60, 65, 70, 80, 90, 100、樣本大小 和 5 , 4 , 3 , 2 = p 5 = n α =0.05,0.01,0.0027之組合下，各模擬 =1,000,000 次及 =100 次，算出在不同 N b p,m,n,α 組合下的管制界限及管制界限的標準差。從結果可以 觀察到一些與Σ0已知情形之類似結果： 對相同的 、n、p α，當 m 愈大時，管制界限CLp,m,n,α愈小，會收斂至一個值。 對相同的 m、n、α，當p愈大時，管制界限CLp,m,n,α愈大。 (1-α)愈大，管制界限的標準差會愈大。 接著，我們利用模擬實驗來比較單邊與雙邊管制圖的平均連串長度表現。考慮一雙 變量(bivariate)常態製程。若Σ0已知，考慮觀察值個數n=5、10；若 未知，則考慮訓 練樣本組數 、50、樣本大小 0 Σ 25 = m n=5。兩者α 均為 ，也就是在製程穩定下 。我們重覆的次數 次。依前述，在不失一般性下，可假設 0027 . 0 370 0 ≈ ARL K = 50,000 Σ =0 Ip，
因此在失控狀態下的共變異矩陣為 ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ∆ ∆ ∆ ∆ ∆ ∆ = Σ 2 2 1 2 1 1 ρ ρ ，其中∆ ，i 為第 i 個 品質特性其變異數所增加的倍數及 2 , 1 = i ρ為兩個品質特性的相關係數。由模擬結果，我們觀 察到： 單邊的 都比雙邊的 小，表示單邊偵測製程變異增加的能力比雙邊偵測製 程變異增加的能力好。 1 ARL ARL1 或 愈大，偵測變異增加能力愈好。 n m 當 或 固定時，變異愈大愈容易被偵測出來。 n m 當∆1,∆2,n固定時，ρ增加，而 愈小，表示相關係數增加時偵測變異增加能力 愈好。 1 ARL 我們也應用本研究所提出的方法去監控一個晶圓資料實例及兩個模擬例子，表現結 果都證明我們的方法對偵測製程變異增加的能力有相當不錯的表現。 在第二年的研究中，我們研究 EWMA 單邊製成變異管制圖之建構與表現。一般而 言，EWMA 管制圖對製程變異性有小變動的情形較敏感。若應用 EWMA 的想法於本文 中之管制統計量 T，應該可以用來偵測製程變異的較小變動。 令 EWMA 統計量為Wt =λTt + −
(
1 λ)
Wt−1，t=1, 2,...，其中 為統計量((4.1)及(4.2)) 在時間 t 之觀察值。Crowder (1987) 針對計算單變量平均值的 EWMA 管制圖的平均連 串長度導出其迭代的式子： t T 1 (1 ( ) 1 h ( ) h x u L u L x f λ dx λ − λ − − ⎛ = + ⎜ ⎝ ⎠∫
) ⎞⎟ 其中L u( )代表 EWMA 統計量之起始值是u時的ARL、 f( )⋅ 是樣本的分配、h 是管制界 限上界、-h 是管制界限下界及λ代表 EWMA 的權重。若要利用 Crowder (1987)的方法來求得我們的 EWMA 管制圖的 ARL，首先我們需 知道上式中的機率密度函數 f
( )
⋅ ，亦即統計量 T ((4.1)及(4.2))的分配。在 已知的情況 下，假設製程在控制狀態下。當 和 時，對 ，模擬統計量 T 值 100,000 筆，我們觀察到其分配有下列特性： 0 Σ n = 5 n = 10 p=2, 3, 4, 5, 10, 15 ‧當品質特性 p 增加時，高峰有往右移的趨勢。 ‧ 當品質特性 p = 2, 3, 4, 5 時,在統計量 T=0 附近的機率密度估計長條圖的高度很 高。 ‧ 當品質特性 p = 2, 3, 4, 5 時, 統計量 T 的分配是一個混合(連續和離散)型的分 配。 在Σ0未知的情況下，當m=25和m=50時，對n=5、 模擬統計量 T 值 100,000 筆，我們觀察到，所得之統計量 T (4.1)和(4.2)之分配是非常相似的，亦有上 述 已知的情況下之特性。 = 2, 3, 4, 5, 10, 15 p 0 Σ或次數多邊形估計( frequency polygon estimation)，來估計 T 之機率密度函數。先模擬生 成N筆統計量 T 之值t t1, ,...,2 tN。本研究取N =1, 000, 000。 核密度估計量為： 1 1 ˆ ( ) , N i i t t f t K Nb = b − ⎛ ⎞ = ⎜ ⎟ ⎝ ⎠
∑
其中 b 為帶寬(bandwidth)，N代表樣本個數， 代表資料點，ti i=1,...,N， 為核函 數，此研究中我們採用 Gaussian kernel，( )
K ⋅( ) ( )
2 1 2 2 2 u K u = π − e− ，並取b = 1.059× ×σ N-15。用 核密度估計去估計分配的方式是很常用的方法，因此我們原先考慮使用核密度估計去估 計分配。 由前面的討論中得知，不管Σ0已知或未知，在品質特性p=2,3,4,5時，均是混合(連 續和離散)型分配，在接近 0 的附近高度特別高，因此先把 1,000,000 筆資料中為 0 的部 分挑出來，再將剩下的資料用核密度估計去估計分配。結果我們發現核密度估計量，在 邊界會有嚴重低估的情形，即使做了邊界修正還是十分的不理想。因此我們改考慮使用 次數多邊形估計去估計機率密函數。 次數多邊形估計的估計機率密函數估計式為：( )
2 1 1(
1)
1 ˆ j j j j j j f t n c n c n n t Nb ⎡ + + + ⎤ = ⎣ − + − ⎦，t∈ ⎣⎡c cj, j+1⎤⎦ ， j=1,...,k， 其中N代表樣本個數，k為區間個數， 1 2 j j j b b c = + + ， 0 1 2 b c = − ，b 1 1 2 k k b c + =b + + ，{
b1,...bk+1}
代表以b = 2.15× ×σ N−15為寬度所切的區間的端點，{
}
1, 2... k n n n 代表每個區間內資料點的 個數，{
c c1, ...2 ck}
代表每個區間的中點。我們發現這方式雖然較不平滑，但是可以解決 邊界問題。 因在 p=2,3,4,5時，是混合型(連續和離散)，而在p=10,15時，是連續型，因此我 們分成兩個情況來討論如何計算平均連串長度 ARL。 (1)當 ( )f ⋅ 是連續型機率密度函數時，ARL，L u( )，滿足下式： L u( )=( )
(
)
0 1 1 1 hL y f y λ u dy λ λ − − ⎛ ⎞ + ⎜ ⎝ ⎠∫
⎟ (4.3) (2)當 ( )f ⋅ 是混合型(連續和離散)機率密度函數，先令F = p F1⋅ d + −(1 p1)⋅ 代表 TFc 的累積機率密度函數，其中 代表離散的累積機率密度函數、 代表連續的累積機率 密度函數，而 d F Fc 1 p 為離散型的比例。因假設離散的分配只有在 T=0 的時候才有機率，所以 我們得到 L u( )= 1(
(
)
)
(
(
)
)
(
1)
( )
(
)
0 1 1 1 p L 1 λ u I 1 λ u h 1 p hL y f y λ u λ λ − − ⎛ + ⋅ − − ≤ + − ⋅ ⋅ ⎜ ⎝ ⎠∫
⎞⎟dy (4.4)其中I
( )
⋅ 代表指標函數 (indictor function)。我們可仿 Crowder (1987) 計算 ARL 的方法來解出(4.3)和(4.4)式之 ARL 值。管制界
限的選擇可利用製程在管制狀態下之ARL0來做決定。 當Σ0已知時之管制界限可由以下之演算法求得： 步驟一、先給定 p 、 ，指定 ，模擬統計量 T ((4.1)式) S 筆資料，利用次數多邊形 估計去估計統計量 T 的機率密度函數。 n ARL0 步驟二、利用模擬資料，找出一個 大於指定 的 和一個 小於指定 的 作為管制界限起始值 和 。 0
ARL ARL0 h1 ARL0 ARL0
2 h h1 h2 步驟三、利用上述的方法，計算當管制上界為h=(h1+h2) /2 ARL0 h 的 (h)。 步驟四、當ARL0(h)大於指定的ARL0則h1← ，否則就h2 ← ；重複步驟三、四直到找h 到一個ARL0(h)很接近指定的ARL0的 。 h 步驟五、重複步驟一到步驟四 B 次，這 B 次 值的平均為我們所需要的管制界限。 h 當Σ0未知時之管制界限之求得亦類似。 當Σ0已知時，對 p=2,10, n=5,10，和當Σ 未知時，p=2,10, n=5, m=25,50 各種情形，0 分別產生 S=1,000,000 筆的資料，並取 B=100，利用上述方法建構出 EWMA 管制圖之管 制界限。 我們觀察到，λ越小時管制界限越低；Σ 已知的情形下，0 n=5的管制界限較 的 管制界限低； 未知的情形下， 10 n= 0 Σ m=50的管制界限較m=25的管制界限低。 接著，我們探討由以上方法建構出來的 EWMA 管制圖在失控情況下之平均連串長 度的表現，來評估其績效。假設製程在失控的情況下其平均數向量不變，但其共變異矩 陣改變為 ，利用前述方法的管制界限，以模擬的方式來求得失控情況下之平均 連串長度 ，方法如下：
(
0 ∑ ≠ Σ)
1 ARL (1) 當Σ0已知時，失控下之平均連串長度演算法： 步驟一、給定品質特性個數 p 、子群樣本數n、假警報率α 、EWMA 參數λ、管制界限 , , , p n C α λ。 步驟二、令t=1。 步驟三、生成失控下的一組樣本數為 的合理子群，且這 個樣本皆服從 ，計 算其統計量T ((4.1)式)。 n n Np(0, )∑步驟四、代入 EWMA 的關係式Wt =λT+ −
(
1 λ)
Wt−1，其中初始值W0 =0。 步驟五、若Wt>Cp n, , ,α λ，則記錄連串長度 ；否則t t← +t 1，回到步驟三。 步驟六、重複步驟二至步驟五 S 次，計算這 S 次連串長度的平均值及標準誤差即為 和 。 ^ 1 ARL 1 SRL (2)當Σ0未知時，失控下之平均連串長度演算法： 步驟一、給定品質特性個數 p 、子群樣本數 、訓練樣本組數n m、假警報率α、EWMA 參數λ、管制界限Cm p n, , , ,α λ。 步驟二至步驟六如上，除了統計量 T 由(4.2)式計算。 在單變量製程監控平均數中，在 Crowder (1987)提到， EWMA 管制圖中，大部分 的情況下平均數改變時，λ越小偵測效率越好；但在比較λ=0.05與λ=0.1時偵測效率 時，有時λ=0.05的偵測效率比λ=0.1的偵測效率差；由 Crowder (1987)所附的表也觀 察出，隨著平均數改變量變大(約大於 2.5)時，λ=0.05的偵測效率比 Shewhart 管制圖的 偵測效率還差。因此並不是λ取越小就越好。我們也針對一些失控的情況，去觀察是否 有相同的現象。 今考慮品質特性個數 p=2的情況，令失控下的共變異矩陣 2 1 1 2 1 2 2 2 σ ρσ σ ρσ σ σ ⎡ ⎤ Σ = ⎢ ⎥ ⎣ ⎦， 在p=2時我們只需考慮相關係數ρ及變異數 2 1 σ 、 2 2 σ 增加的情況， 2 1 σ 代表第一個的品 質特性的變異數、 2 2 σ 代表第二個的品質特性的變異數。同時ARL0指定為 370.4。 針 對 Σ0 已 知 及 Σ0 未 知 兩 種 p=2 、 ρ=0 、 0.05, 0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1 λ= ， 的情況下，利用上述方 法取 S=100,000 筆計算 2 2 1 2 (σ σ, )=(1.05,1), (1.10,1) , (1.15,1) , (1.25,1) , (1.50,1) , (1.75,1) , (2.00,1) , (2.25,1) , (2.50,1) , (2.75,1) , (3.00,1) 1 ARL 的結果。 由結果我們觀察到： ‧在變異數改變較小的情況下隨著λ的值越小，ARL1也會越小；但變異數改變較大的情 況下，ARL1會先減小再增加。 ‧λ=0.05的偵測效率大部分都比λ=0.1時差，有時甚至比 Shewhart 管制圖的偵測效率 還差。 因此我們發現 EWMA 管制圖的偵測效率在小和中幅度變動下比 Shewhart 管制圖 好，而且在小變動下，較小的λ會有較好的偵測效率；但是λ=0.05和λ=0.1在較大變 動下的偵測效率較差。五、結果與討論 0 在本研究中，我們利用多變量單邊檢定發展偵測多變量製程變異增加的方法，即為 檢定H0:Σ = Σ0對H1:Σ ≥ Σ0 且 Σ ≠ Σ ，其中Σ 為所監控品質特性的共變異數矩陣及 為在管制狀態下的製程共變異矩陣，並且分成為 0 Σ Σ 已知或未知兩種情形來討論如何0 監控製程變異增加的問題。 在第一年的研究中，我們導出Σ 已知或未知兩種情形下之概似比檢定統計量，並0 用靴環法得出管制界限。針對此控制圖之績效問題，我們經統計模擬對幾種Σ 的變化比 較平均連串長度，並以一個實例和模擬例子，證實所提出的單邊檢定方法對多變量製程 變異增加的問題在偵測能力上有相當不錯的效率，且與允許多變量製程變異性可增加或 減少之雙邊檢定方法作比較，也有較佳的偵測效率。 在第二年的研究中，我們接續第一年的研究，對上述概似比統計量提出 EWMA 管 制圖，並模仿 Crowder (1987)所提出計算平均連串長度的方法進而找出管制界限。我們 發現概似比統計量在某些狀況下為混合型分配，因此 Crowder 的方法需加以修正。我們 同時討論了 EWMA 管制圖在不同的權重λ下的偵測效率。針對不同的 EWMA 的權重 λ，以模擬實驗來探討 EWMA 管制圖在多變量製程變異增加時之平均連串長度表現。 我們觀察到 EWMA 管制圖對多變量製程變異增加的問題，有不錯的偵測效率。與 Shewhart 管制圖比較，發現如所預期的，EWMA 管制圖在變異小至中幅度增加下有較 佳的偵測效率。 對多變量製程變異性作管制時，我們建議可考慮同時使用我們所提之 Shewhart-type 管制圖和其相對應之 EWMA 管制圖來偵測製程變異性增加，此監控技術可以涵蓋更寬 廣的多變量製程變異變動範圍，為一相當有價值的管制機制。 五、參考文獻
1. Alt, F. (1985). Multivariate Quality Control, Encyclopedia of Statistical Sciences 6, pp.110-122, S. Kotz and N. L. Johnson, eds. Wiley, New York.
2. Alt, F. B. and Bedewi, G. E. (1986). SPC of Dispersion for Multivariate Data. American
Society for Quality Control, pp.248-254.
3. Amemiya, Y., Anderson, T. W., and Lewis, P. A. (1990). Percentage Points for a Test of Rank in Multivariate Components of Variance. Biometrika 77, pp.637-641.
4. Anderson, B. M., Anderson, T. W., and Olkin, I. (1986). Maximum Likelihood Estimators and Likelihood Ratio Criteria in Multivariate Components of Variance. The
Annals of Statistics 14, pp.405-417.
5. Anderson, T. W. (1984). An Introduction to Multivariate Statistical Analysis ( Edition ). Wiley, New York.
nd
6. Anderson, T. W. (1989). The Asymptotic Distribution of the Likelihood Ratio Criterion for Testing Rank in Multivariate Components of Variance. Journal of multivariate
Analysis 30, pp.72-79.
7. Calvin, J. A. (1994). One-Sided Test of Covariance Matrix with a Known Null Value.
Communications in Statistics－Theory and Methods 23, pp. 3121-3140.
8. Crowder, S. V. (1987). Average Run Length of Exponentially Weighted Moving Average Control Charts. Journal of Quality Technology 19(3), pp. 161-164.
9. Kuriki, S. (1993). One-Sided Test for Equality of Two Covariance Matrices. The Annals
of Statistics 21, pp. 1379-1384.
10. Muirhead, R.J. (1982). Aspects of Multivariate Statistical Theory. Wiley, New York. 11. Sakata, T. (1987). Likelihood Ratio Test for One-Sided Hypothesis of Covariance
Matrices of Two Normal Populations. Communications in Statistics－Theory and Methods 16, pp. 3157-3168.
12. Shiau, J.-J. H. and Hsu, Y.-C. (2005). Robustness of the EWMA Control Charts to Non-normality for Autocorrelated Process. Quality Technology & Qualitative
Management 2(2).
13. Tang, P. F. and Barnett. N. S. (1996a). Dispersion Control for Multivariate Processes.Australian Journal of Statistics 38, pp. 235-251.
14. Tang, P. F. and Barnett. N. S. (1996b). Dispersion Control for Multivariate Processes－ Some Comparisons. Australian Journal of Statistics 38, pp. 253-273.
15. Yeh, A. B., Huwang, L., and Wu, Y.-F. (2004). A Likelihood Ratio based EWMA Control Chart for Monitoring Multivariate Process Variability. IIE Transactions in Quality and
Reliability Engineering 36, pp. 865-879.
16. Yen, C.-L. and Shiau, J.-J. H. (2004). A Multivariate Control Chart for Detecting Increases in Multivaraite Process Variation. Technical Report, Institute of Statistics, National Chiao Tung University, Hsinchu, Taiwan.
計劃成果自評
本計劃之執行相當順利，結果亦相當豐碩，三位碩士班學生和一位博士班學生在工 業統計上得到相當不錯的訓練。研究內容與原計畫相符程度極高，亦已達成預期目標。 研究成果預計將有兩篇論文可以在國際知名期刊發表，其中一篇初稿(Yen and Shaiu (2004))已完成，目前正在修改中，另一篇則正在撰寫中。本研究所提出之多變量製程變 異性控制圖在學術上和應用上均有貢獻。
另外，在計畫研究期間，我們同時也研究了自我相關製程 EWMA 控制圖對非常態 分配之穩健性，此研究結果(Shiau and Hsu (2005))已被國際期刊 Quality Technology & Qualitative Management 接受，並將於 2005 年 9 月發表。
