Borror, C.M., Shamp, C.W., Rigdon, S. (1998) Poisson EWMA Control Charts. Journal of Quality Technology, 30(4), 352-361.
カウント時系列の監視手法のひとつ、ポアソンEWMA管理図について知りたくて読んだ。初学者向けの解説論文である。ありがてえ、ありがてえ。
1. イントロダクション
反復的な生産プロセスで、ある生産単位における不適合数\(X\)を品質指標にしているとき、\(X_1, X_2, \ldots, \)をiidなポアソン確率変数と捉えることが多い。多くの管理図手続きでは、out-of-control状況が観察における非ランダムなパターンとして表現される。というわけで、ポアソン分布に基づくc管理図やu管理図が用いられることが多い[品質管理の用語だなあ。不適合数のことをc, 単位当たり不適合数のことをuというらしい]。
そのひとつがシューハートのc管理図である。生産プロセスの出力から一単位を取ってきて不適合数を数える。時点\(t\)における不適合数を\(X_t\)として、\(h_L \leq X_t \leq h_H\)ならin-controlと捉える。もし監視したいのが平均\(\mu\)のある瞬間のシフトで、品質指標がiidならこれでよい。なお、複数単位を抜きだせるんならu管理図が良い。
c管理図にせよu管理図にせよ、プロセスの状態は直近の標本のみで決めるわけである。もっと前の情報を含めようという提案がいろいろある。ルール付きシューハート図、移動平均図、CUSUM図、EWMA図など。
本論文ではEWMA図について論じる。なお、EWMA統計量は自己相関するので非ランダムなパターンを探すのには向いてない。
[…先行論文…]
管理図のARLとは、OC条件が発生するまでの間のステージ数の期待値として定義される。プロセスがICならそれは長くあってほしいし、シフトが起きているんなら短くあってほしい。詳しくはWoodall(1997 J.Qual.Tech.)をみよ。
2. ポアソンEWMA管理図
\(X_1, X_2, \ldots\)がiidなポアソン確率変数で、IC下で平均\(\mu_0\)とする。
EWMA管理図というのはこうである。$$ Z_0 = \mu_0 $$ $$ Z_t = \lambda X_t + (1-\lambda) Z_{t-1} $$ この統計量の分布は$$ E(Z_t) = \mu_0$$ $$ Var(Z_t) = \frac{\lambda}{2-\lambda} \left( 1 – (1-\lambda)^{2t} \right) \mu_0$$ \(t\)が大きければ$$ Var(Z_t) \approx \frac{\lambda \mu_0}{2-\lambda} = Var(Z_\infty)$$ となる。
管理限界をどう決めるか。\(Var(Z_t)\)の一本目の式の基づく方法もあるけど、ここでは二本目の式に基づくほうについて述べる。$$ h_L = \mu_0 -A_L \sqrt{ \frac{\lambda \mu_0}{2-\lambda} }$$ $$ h_L = \mu_0 -A_U \sqrt{ \frac{\lambda \mu_0}{2-\lambda} }$$ \(A_L = A_U\)とすることが多いがそうでなくてもよい。値の決め方は後で触れる。
\(X_1, X_2, \ldots\)がポアソン確率変数ならば、\(Z_t\)は非負だから、もし\(h_L \lt 0\)になったら0にしてよい。こういう場合には\(A_L neq A_U\)とするのがよい。下方シフトを見つけたいからね。
3. マルコフ連鎖アプローチによるARL計算
過程平均にシフトがあったときのARL[\(ARL^1\)のことかな]は、マルコフ連鎖アプローチで近似できる。
\((h_L, H_U)\)を\(N\)区間に均等に分割する。\(j\)番目の区間を\((L_j, U_j)\)とする。\(N+1\)番目はabsorbingな状態である(OC警告が出るから)。ARLはマルコフ過程がabsorptionするまでの期待時間となる。
区間\(i\)から\(j\)への遷移確率を$$ P_{ij} = P(L_j \lt Z_t \lt U_j | Z_{t-1} = m_i)$$ とする。\(m_i\)というのは区間\(i\)の中点である。この式を展開するとこうなる。$$ P_{ij} = P(h_L + \frac{h_U – h_L}{2N\lambda} (2(j-1)-(1-\lambda)(2i-1)) \lt X_t \lt h_L + \frac{h_U – h_L}{2N\lambda} (2j-(1-\lambda)(2i-1)) $$
4. \(A\)の値の選択
\(A\)をどう決めるか。
たとえば、IC平均が\(\mu_0 = 7\)で、\(\lambda = 0.2\)で、IC ARLを500にしたいなら、\(A\)は2.975にするとよい。[チャートが載っている]
プロセスが区間\(i\)から始まった(つまり\(Z_0 = m_i\))ときのARLを\(R_i\)とし、\(\mathbf{R} = \{R_1, R_2, \ldots, R_N\}\)とする。遷移行列を\(\mathbf{P}\)とし、そこから\(N+1\)行目と\(N+1\)列目をとった行列(つまりIC状態内部の遷移行列)を\(\mathbf{Q}\)とする。すると、$$ (\mathbf{I} – \mathbf{Q}) \mathbf{R} = \mathbf{1} $$ が成り立つことが知られている。[へえええええ！]
ここから\(\mathbf{R}\)が求まる。その中央にある値がARLとなる。(なので、\(N\)は奇数にしておくとよい。)
上の方法は、プロセスのシフトがシステムの初期段階で起きると仮定している(よって初期状態ARLという)。これとは別に、安定状態ARLというのも求めることができる。
まあとにかく、このやり方でARLを求めてくれるSAS IMLコードがあるから使うとよい。
5. 計算例
[略]
6. 他のチャートとの比較
Gan(1990)は、EWMA統計量の3つの求め方を提案している。どのタイミングでどういう風に\(Z_i\)を整数に丸めるかが異なる。[…後略…]
7. 結論
本論文のアプローチの利点は、ARLがシューハートのc管理図やGanの修正版よりも小さいという点である。さらに、下側の管理限界が正になるので下方変化も検出できる。
—————-
どの論文か忘れたけど、この数日間の間に読んだ論文のなかに、ARLの算出にマルコフ連鎖を使うという話が出てきて、なんのこっちゃ？ と疑問に思ったのだが、この論文を読んで意味がわかった。
それはともかく、品質管理の世界では、時点ごとに標本サイズが変動するという発想はあまりないようね… やはり疫学的サーベイランスの論文を探したほうがよいのかなあ。