StatModeling Memorandum

StatModeling Memorandum

StanとRとPythonでベイズ統計モデリングします. たまに書評.

蟻本シリーズ 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)を用いることで、なんと \mathcal{O} (N) のオーダーで解けます。詳しくは蟻本の4-4節(p.300)を読んでください。

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

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

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

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

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

  • はじめにxK個分のインデックスを両端キューの末尾に格納していきます。ただし、途中で末尾のインデックスに対応する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の挙動を見ると以下になります。

fig1.png


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

上記のデータがあります。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_obsYからX_trueを推定することを考えます。ただし、X_trueの欠測部分に関しては事前分布をnormal(0, 5)とします。

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

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

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

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

meanse_meansd2.5%25%50%75%97.5%n_effRhat
x_mis[1]0.890.074.20-6.11-2.200.603.669.9936861.00
x_mis[2]0.840.074.25-5.99-2.530.553.729.8332871.00
x_mis[3]0.890.074.38-6.17-2.510.583.8110.1537871.00
x_mis[68]1.800.083.88-3.53-1.331.414.1710.5423101.00
x_mis[69]1.680.083.86-3.77-1.461.264.1910.0825341.00
x_mis[70]-5.340.010.39-6.12-5.60-5.35-5.09-4.6040001.00
s_y0.370.000.030.320.350.370.390.4513421.00
lp__11.970.236.32-1.357.8212.3616.4723.287771.00

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

fig2.png