蟻本シリーズ 3 スライド最小値

今回は以下の問題を考えます。

長さNの数列x[1], x[2], ..., x[N]と数Kが与えられます。y[i] = min{x[i], x[i+1], ..., x[i+K-1]} (i = 1, ..., N-K+1)として定義される数列y[i]を計算しなさい。

この問題は両端キュー（デック, deque）を用いることで、なんとのオーダーで解けます。詳しくは蟻本の4-4節（p.300）を読んでください。

プログラミングコンテストチャレンジブック [第2版]　～問題解決のアルゴリズム活用力とコーディングテクニックを鍛える～

ここでは蟻本から例題を引用してN = 5, K = 3, x = {1, 3, 5, 4, 2}, y = {1, 3, 2}の場合について簡単に説明します。

• まずdeqは余裕をもって長さIの配列として定義し、両端キューの先頭を指すインデックスをs、末尾の次の要素を指すインデックスをtとします。そして、deqのsからt-1の部分を両端キューとします。
• 両端キューdeqはxのインデックスiを格納します。ただし、両端キューの中でインデックスは昇順になるように格納します（すなわちdeq[j] < deq[j+1]）。さらに各インデックスに対応するxの値も昇順になるように格納します（すなわちx[deq[j]] < x[deq[j+1]]）。この結果、この先最小値となりうるxのインデックスだけを格納することができます。

具体的な構成の仕方は以下の通りです。

• はじめにxのK個分のインデックスを両端キューの末尾に格納していきます。ただし、途中で末尾のインデックスに対応するxの値（x[deq[t-1]]）が、いま入れようとしているxの値（x[i]）以上の場合には、x[deq[t-1]] < x[i]を満たすようになるまで末尾tを縮めます。（★）
• K個分追加したら、K個のxの中の最小値を計算できます。両端キューの先頭の値deq[s]を使ったx[deq[s]]が最小値となります。したがってy[1] = x[deq[s]]です。そして、deq[s]が1ならば、次回から最小値を求める範囲の外になるため、sを1つ増やして両端キューを縮めます。（♪）
• 次にy[2]を求めるため、x[deq[t-1]]とx[K+1]の大きさを比べ、同じようにK+1を両端キューの末尾に追加します。以下この操作を繰り返します。

例題の場合のdeqの挙動を見ると以下になります。

---

以上は自分用の備忘録でした。以下では次のような問題を考えます。

上記のデータがあります。i列はインデックスです。X_true列は長さ100の数列です。X_obs列はそのうち部分的に与えられたものです。Y列はy[i] = min{X_true[i], X_true[i+1], ..., X_true[i+K-1]} (i = 1, ..., N-K+1)として定義される数列y[i]に、さらに測定ノイズを加えたものです。ここでK = 8とします。さてX_obsとYからX_trueを推定することを考えます。ただし、X_trueの欠測部分に関しては事前分布をnormal(0, 5)とします。

Stanでの実装例は以下になりました。

• 2～21行目： 両端キューを使ってスライド最小値を計算する関数です。
• 10～11行目： 上記で説明した構成の仕方の（★）に対応しています。
• 14～18行目：　同じく（♪）に対応しています。
• 46行目： 値の分からないx_misを混ぜた状態で、スライド最小値を求めてmuに代入しています。
• 52行目： muに測定誤差が加わって観測したデータYになっています。

実行するRコードは以下の通りです。

推定結果は以下になりました。

	mean	se_mean	sd	2.5%	25%	50%	75%	97.5%	n_eff	Rhat
x_mis[1]	0.89	0.07	4.20	-6.11	-2.20	0.60	3.66	9.99	3686	1.00
x_mis[2]	0.84	0.07	4.25	-5.99	-2.53	0.55	3.72	9.83	3287	1.00
x_mis[3]	0.89	0.07	4.38	-6.17	-2.51	0.58	3.81	10.15	3787	1.00
…	…	…	…	…	…	…	…	…	…	…
x_mis[68]	1.80	0.08	3.88	-3.53	-1.33	1.41	4.17	10.54	2310	1.00
x_mis[69]	1.68	0.08	3.86	-3.77	-1.46	1.26	4.19	10.08	2534	1.00
x_mis[70]	-5.34	0.01	0.39	-6.12	-5.60	-5.35	-5.09	-4.60	4000	1.00
s_y	0.37	0.00	0.03	0.32	0.35	0.37	0.39	0.45	1342	1.00
lp__	11.97	0.23	6.32	-1.35	7.82	12.36	16.47	23.28	777	1.00

推定結果の図を描くと以下になります。黒点はX_obs、ピンク点はY、黒線はX_true、青帯はMCMCサンプルから算出した、未知のX_trueの80%ベイズ信頼区間です。