Shen, X., Tsui, K., Woodall, W., Zou, C. (2015) Self-starting monitoring scheme for Poisson count data with varying population sizes. Technometrics, 58(4), 460-471.
仕事の都合で読み漁った、カウント時系列の監視の論文のひとつ。残念ながら提案手法が理解できず、途中で読むのをやめてしまった奴である。いちおう記録しておくが、うーん、残念。
時系列監視の分野では、平常時のパラメータが未知なのに監視を始めないといけないという問題を自動スタートself-startingというらしい。機械学習とかでいうところのcold startと似た意味だと思う。
1. イントロダクション
統計的プロセス管理で、好ましくない出来事の発生率を監視する際、ふつうは発生数を、標本サイズを所与とした独立なポアソン確率変数と捉える。標本サイズが定数の場合も多い。Lucas(1985)とか。標本サイズが時変する場合にも関心が集まっている。レビューとしてPurdy et al.(2015 Technometrics)をみよ。
先行研究は主に前向きサーベイランスのフェイズ2に関心を向けてきた[履歴からin-control(IC)のパラメータを学習する段階をフェイズ1、管理限界を決めてout-of-control(OC)に遷移しないかどうか監視する段階をフェイズ2と呼んでいるのだと思う]。つまりICパラメータが既知である場面である。しかし監視は履歴情報が限られている段階で始めないといけないことが多い。こういう場合には自動スタート法が必要になる。
自動スタート法の歴史をたどると、まずHawkins(1987)が単変量正規データで自動スタートなCUSUM図を、Quesenbery(1991)が自動スタートなシューハート図を提案した。多変量版についても提案がある。自動スタートのプロファイル監視[なにそれ]の提案もある。
ポアソン率の監視の場合、まずQuesenberry(1991), Hawkins & Olwell(1998)による管理図の提案がある。後述するけど、前者はポアソン件数の期待値が低いときにICのARLが良くない[異常がなくても誤警告があるということであろう]。後者はIC, OCともARLは良いんだけど、標本サイズ固定が前提になっている。標本サイズ時変に拡張するのは結構難しい。[ここまではたぶん自動スタートの話抜きだと思うけど自信がない]
自動スタートな図を設計するのは難しい。観察がポアソン分布のとき、pivotalな検定統計量なし制御限界を決めるのが難しいからである[pivotalとは…?]。いいかえると、ポアソン過程を監視する際、統計量のIC下の分布は、ふつう基底にあるポアソン過程から切り離せないので、IC下のARLを固定してOC下のARLを満足いくレベルにするような管理限界を決めるのが難しいのである。
[あ、なるほど、云っている意味が分かったような気がする。ガウス過程だったら分散は時変せず期待値だけ時変すると仮定できるけど、ポアソン過程ではそうできない、ってことかな]
というわけで本論文では、ポアソン過程のIC下でのパラメータが未知なのに好ましくないイベントの発生率を監視するための、自動スタート型の指数重みづけ移動平均(EWMA)管理スキームを提案する。EWMA型で提案するけど、同じアイデアはシューハート管理図やCUSUM管理図にも使えるよ。
2. ポアソン率監視のための自動スタート管理図
2.1 モデルの仮定
時点\(t\)で観察した件数\(X_t\)が、平均\(\theta n_t\)のポアソン分布に独立に従うと仮定する。ある未知の時点\(\tau\)で\(\theta\)が\(\theta_0\)から\(\theta_1 (\gt \theta_0)\)に変化するのを検出したいとしよう。ふつうは独立な\(m_0\)個の履歴 \( \{(X_i, n_i)\}^0_{i=-m_0+1} \)を得ますよね。で、次の変化点モデルを考える。独立に分布していることを\(\sim_{id}\)と書く。$$ X_i \sim_{id} \begin{cases} Poisson(\theta_0 n_i | n_i) & \mathrm{for} \ i = -m_0+1, \ldots, 0, \ldots, \tau-1 \\ Poisson(\theta_1 n_i | n_i) & \mathrm{for} \ i = \tau, \ldots \end{cases} $$
2.2 先行研究
- TR-1。$$ \hat{\theta}_0 = \frac{\sum_{j = -m_0 + 1}^0 X_j}{\sum_{j = -m_0 + 1}^0 n_j} $$ を求め、これを真とする。\(m_0\)が小さいとひどい目にあう。
- TR-2。Shen et al.(2013)。確率的な管理限界を考え、発生率\(\theta_0\)のポアソン変数を生成するシミュレーションをやって、時点ごとに管理限界\(h_t\)を決める。\(\theta_0\)がわかんない場合、恣意的に決める。[おいおい、これもダメに決まっとるだろうが]
- TR-3。Quesenberry(1991)。観察されたポアソン件数を、標準正規分布に従う変数に変換します。[さあ歯を食いしばれ！！] $$ b(x; s, p) = \begin{cases} {}_s C_x p^x (1-p)^{s-x} & x = 0, \ldots, n \\ 0 & \mathrm{otherwise} \end{cases} $$ と定義し、$$ B(x; s, p) = \begin{cases} 0 & x \lt 0 \\ 1 & x \geq s \\ \sum^{[x]}_{k=0} b(k; s, p) & 0 \leq x \lt n \end{cases} $$ と定義しておいて、$$ s_t = \sum_{k=1}^t X_k, \ N_t = \sum_{k=1}^t n_t $$ $$ u_t = B(X_t; s_t, \frac{n_t}{N_t}) $$ $$ Q_t = \Phi^{-1}(u_t)$$ とすると、\(Q_1, Q_2, \ldots\)は独立で近似的に正規分布する統計量になるので、標準正規変数のEWMA監視スキーマを直接に使える。[へええええ！ と感心したんだけど、よく考えてみたら自動スタート問題は全然解決してないっすね]
この方法の弱点は、やはり\(\theta_0 n_t\)に強く依存するという点。ARL性能は保証されない。 - TR-4。Hawkins & Olwell (1998 書籍)。IC下でのポアソン件数を\(X_t \sim Poisson(u_0)\)とみる。\(W_t = \sum_{i=1}^t X_i\)として、標本サイズを\(t\)として[← ???? ほんとにこう書いてある。なにかの間違いじゃないですかね]、IC下では$$ X_t \sim Binomial(W_t, \frac{1}{t}) $$ と書ける。[全然わからん。なにがどうなっているんだ？？？]
\(B_t \sim Binomial(W_t, \frac{1}{t}), \ A_t = P(B_t \leq X_t)\)とする。事前の観察の平均を使って\(u_0\)の推定値\(\hat{u}_0\)を得て、\(X_t\)をパラメータ\(\hat{u}_0\)のポアソン変量\(Y_t\)に変換する。[式が書いてあるんだけど全然理解できないのでメモ省略]
で、\(Y_t\)をEWMAスキームで監視する。この方法の弱点は標本サイズが定数だという点である。
[いやー、これはさっぱり理解できなかった。完全にお手上げである。どこかに誤植があるんじゃないかしらん？]
2.3 提案手法
まず、各時点で過去を全部使って\(\hat{\theta}_t\)を得ます。たとえば$$ \hat{\theta}_t = \frac{\sum_{j=-m_0+1}^{t-1} X_j}{\sum_{j=-m_0+1}^{t-1} n_j} $$ とか。なんか直近に重い重みをつけてもよい。とにかく得ます。
次に、\(X_t\)を標準化します。$$ S_t = \frac{X_t – n_t \hat{\theta}_t}{\sqrt{n_t \hat{\theta}_t}} $$
次に、EWMA型の管理図統計量をつくる。\(\lambda \in (0, 1], Z_0 = 0\)として$$ Z_t = \mathrm{max}(0, (1-\lambda)Z_{t-1} + \lambda S_t) $$ \(Z_t\)が管理限界をこえたら警告する。
ここまではいいですね？
さあ、管理限界をどうするか。
特に自動スタート型管理図では、ARLを考えるだけでは不十分である。ICラン長分布は、もしそれが幾何分布に近かったら満足いくと考えられることが多い。指定されたIC ARLの下で、幾何分布と比べ、短いラン長で誤警告の確率が上がったら、そのチャートはふつう受け入れがたい。誤警告の数が多すぎることは、オペレータの警告への信頼を損なう。そこで、現時点より前に誤警告が出ていないことを所与としたとき、管理図統計量が現在の管理限界を超える条件付き確率が事前に特定した定数となるような管理限界が欲しい。[頭がこんがらがったのでほぼ逐語訳してみたけど、あんましたいしたことは云ってないわ… 下の式を言葉にしただけだ]
理論的には、管理限界\(h_t\)が以下の式を満たすようにしたい。\(\alpha\)を事前に決めた条件付き誤警告率として、$$ Pr(Z_1 \gt h_1(\alpha) | n_1) = \alpha $$ $$ Pr(Z_t \gt h_t(\alpha) | Z_i \lt h_i(\alpha), 1 \leq i \lt t, n_t) = \alpha \ \mathrm{for} \ t \gt 1$$ これは仮説検定で各時点での第一種過誤率を\(\alpha\)にするのと似たような話である。
伝統的には、動的管理限界はマルコフ連鎖近似かモンテカルロ・シミュレーションで得ていた。管理図統計量がふつうはpivotalな量だったからである[???]。しかし、上記の\(h_t(\alpha)\)の近似は難しい。\(Z_t|n_t\)のIC分布が\(\theta_0\)に依存するからである。そこで次のように提案する。パラメトリック・ブートストラップだ。各時点において、期待値\(n_t \hat{\theta}_t\)を持つ擬似的なポアソン観察を生成し、上の理論的要請を近似的に満たすような\(h_t\)を見つけるのである。
[はいはい！ アイデアはわかったぞ。時点ごとにブートストラップ信頼区間を出して管理限界をつくるのね？ なんでそんなややこしいことをするかというと、過去がいまいちあてにならないからであろう]
\(t=1\)から始めよう。\(Y_{i,1} \sim Poisson(n_1 \hat{\theta_1})\)を\(N\)個生成する。それぞれを\(S_{i,1}\)に変換し、さらに管理図統計量\(R_i(1) = \lambda S_{i,1}\)を求める。その長さ\(N\)のベクトルを\(\mathbf{R}(1)\)として、昇順に並び替え、\(h_1(\alpha)\)を近似する値\(R_{[H]}(1)\)を求める。\(Z_1\)と比べ、警告したりしなかったりする。
警告しなかったとしよう。\(t=2\)に進みます。\(\hat{\theta}_2\)を求めますわね。で、[わかりにくいので逐語訳] \(R_i(1)\)を上の理論的要請を満たすfeasibleな値へと制約する。具体的には、\(\mathbf{R}(1)\)の一部\(R_{[1]}(1), \ldots, R_{[H]}(1))\)を保持し、それを\(R_{[i]}(1)\)のfeasibleな値の空間として、この空間から\(R_{i}(1)\)の\(N\)変数をランダムにブートストラップし、更新された\(N\)次元ベクトル\(\mathbf{R}(1)\)をつくる。管理限界\(h_2\)は、\(Poisson(n_2 \hat{\theta}_2)\)確率変数のランダム数生成プロセスを繰り返すことによって同様に定義できる。
[何度も読み返したのだが理解できない… 毎時点で、過去を全部使って\(\hat{\theta}_t\)を作り、\(X_t\)を\(S_t\)に変換し\(Z_t\)に変換するよね？ こっちはブートストラップと関係ないよね？ 問題は管理限界\(h_t\)だよね？
えーと、\(t=1\)では、\(Y_{i,1} \sim Poisson(n_1 \hat{\theta_1})\)を\(N\)個生成し、それぞれを\(S_{i,1}\)に変換し、\(R_i(1)\)に変換し、ベクトル\(\mathbf{R}(1)\)とし、並び替えて\(R_{[H]}(1)\)をみつけて\(h_1\)にするよね。警告がなかったとして、\(t=2\)ではどうするのかがわからない。\(Y_{i,2} \sim Poisson(n_2 \hat{\theta_2})\)を\(N\)個生成し、それぞれを\(S_{i,2}\)に変換し、\(R_i(2)\)に変換し、並び替えて\(R_{[H]}(2)\)をみつけて\(h_2\)にするの？ 違うよね？ \(\mathbf{R}(1)\)のうち\(R_{[H]}(1)\)より下の奴だけを復元無作為抽出して\(\mathbf{R}(2)\)を作るんでしょう？ で、それをなににつかうの？ わからん！！]
[なにかを理解し損ねているようなので逐語訳]
この手続きはパラメトリック・ブートストラップ手法と捉えることができる。パラメトリックな仮定としてポアソンモデルを置き、リサンプリング法で管理図統計量の標本分布を近似する。ここで強調しておきたいのは、この手続きがMargavio et al.(1995)ほかの手続きとは著しく異なるという点である。我々の手法では、管理限界は観察とともにオンラインで決められるのであって、監視の前に決めるのではない。管理限界は時点\(t\)における\(\hat{\theta}_t, n_t\)を使っており、したがってデータ駆動的である。ポアソン分布の離散性により、条件付き誤警報確率は一般に\(\alpha\)と等しくはならないが、以下で示すように、近似は常に近い。
提案した管理図は、プロセスの初頭から実装できるという意味で自動スタート型スキームだが、\(m_0\)が小さすぎるのにテストを始めるのは良いアイデアではない。\(m_0\)が小さいと、短いランでの変化が起きたときにシビアな「マスキング効果」が生じてしまうだろう。実務家のみなさんには、フェイズ1である程度の数の観察を集め、また少なくともプロセスが実際に安定しているかどうかの初期検証が可能になる程度の先行知識を得てから管理図を開始することを勧める。[…]
[時間の都合でここで打ち止めにし、見出しのみメモする]
3. パフォーマンス研究
4. 応用例
5. 結論
——————-
うーん、肝心の部分を理解し損ね、残念至極である。著者様にお伺いのメールを危うく送りかけたほどだが、すんでのところで思いとどまった。Accepted manuscriptを読んでいるのがいかんのかもしれない(どっかに誤植が残っているとかで)。
おそらくアイデアとしては、統計量はポアソン件数をEWMAに変換したもので、管理限界のほうをいろいろ工夫して、毎時点の誤警報率をどうにかして\(\alpha\)に揃える、という話なのだろうと思う。で、いかに自動スタート型管理図といえど、それなりのフェイズ1は必要だ、ということだと思う。そらまあそうでしょうな。