蟻本シリーズ 1 ナップサック問題

「プログラミングコンテストチャレンジブック [第2版]」（通称：蟻本）という本がとてもよかったので、これから3回にわたって統計モデルを絡めた感謝の記事を書こうと思います。

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

競技プログラミングはたまに実務で役に立たないとか言われていますが、僕はまったくそうは思いませんでした。プログラマーとしてもデータ解析者としても有益な考え方が満載で、ふと気づけば応用できる対象は至る所にあるように思いました。

---

今回は以下のナップサック問題と呼ばれる問題を考えます。

重さと価値がそれぞれWeight[i], Value[i]であるようなN個の品物があります。これらの品物から、重さの総和がMax_Wを超えないように選んだときの、価値の総和の最大値を求めなさい。

蟻本では2-3節（p.52）になります。検索しても解説がいくつか見つかります（例）。

この問題を速く解く方法は漸化式でdp（dynamic programming）テーブルを作る方法です。なんとのオーダーで解けます。ここでは詳しい解説は省略しますが、蟻本から例題を引用してN = 4,(Weight, Value) = {(2, 3), (1, 2), (3, 4), (2, 2)} (i = 1,2,3,4),Max_W = 5の場合について簡単に説明します。

• dp[i, w]を「i番目までの品物から重さの総和がw以下となるように選んだときの価値の総和の最大値」とします。dpは2次元配列でインデックスが0からはじまるとします。この設定がポイントです。
• dp[0, w] = 0 (w = 0, 1, ..., Max_W)になります。なぜならi = 0は「品物なし」の状態で価値の総和は常にゼロだからです。

ここからi = 1の場合、i = 2の場合、・・・のdp[i, w]を以下の漸化式で順番に決めていきます。

• もしw < Weight[i]ならばi番目の品物（重さWeight[i]）を入れることはできません。よって、その前のiまでの品物から選んだ時の価値の最大値に等しくなります。したがってdp[i, w] = dp[i-1, w]となります。
• もしw >= Weight[i]ならばi番目の品物を入れることができます。すると、入れない場合の最大値dp[i-1, w]と入れる場合の最大値dp[i-1, w-Weight[i]] + Value[i]を比べて大きい方の値がdp[i, w]となります。入れる場合の最大値は少々分かりにくいですが、dp[i-1, w-Weight[i]]はWeight[i]だけ余裕がある状況の価値の最大値であり、そこにi番目の品物（価値Value[i]）を入れるのでdp[i-1, w-Weight[i]] + Value[i]になります。

例題の数字を使って、dp[3, 4]を決めている状況が以下の図になります（蟻本のp.55から引用）。i = 3なので(Weight, Value) = (3, 4)です。よってdp[3, 4] = max(dp[2, 4], dp[2, 4-Weight[i]]) + Value[i] = max(5, 2+4) = 6となります。

このままdpテーブルを全て埋めてdp[N, Max_W]を求めると、それがナップサック問題の解となります。

---

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

上記のファイルのように重さはすべて既知である30個の品物があるとします。価値についてはそのうち25個しか分かっていません。ここで、あるブラックボックスの自然現象のルールがあって、30個の品物から重さの総和が40を超えないように選んで、価値を最大化すると考えられています。そして、最大化された価値を8回測定しました。結果は測定誤差を含めて、Y <- c(63.1, 66.8, 72.5, 58.1, 55.2, 67.5, 69.1, 70.2)でした。この結果から未知である5つの品物の価値の（事後）分布を推定することを考えます。ただし、価値の事前分布をnormal(3, 5)とします。

この問題はナップサック問題の逆問題のようになっています。Stanでの実装例は以下になりました。

• 2～15行目：冒頭で説明したアルゴリズムを実装しています。ただし、Stanの配列はインデックスが1からはじまるので、2次元配列dp[,]に関係するインデックスはすべて+1しています。
• 21, 23, 32行目：I_misは価値が未知の品物の個数、Idx_misはその品物のインデックス（ここでは26～30）、v_misはその品物のそれぞれの価値（パラメータ）です。
• 37～43行目：価値の分からない品物を混ぜた状態で、価値を最大化してmuに代入しています。
• 48行目：muに測定誤差が加わって観測したデータYになっています。

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

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

	mean	se_mean	sd	2.5%	25%	50%	75%	97.5%	n_eff	Rhat
v_mis[1]	5.88	0.09	4.21	0.38	2.63	5.07	8.17	16.27	2294	1.00
v_mis[2]	5.18	0.09	3.60	0.27	2.29	4.57	7.47	13.28	1628	1.00
v_mis[3]	9.87	0.13	4.54	0.85	6.97	10.16	13.02	18.26	1202	1.00
v_mis[4]	7.75	0.13	5.05	0.23	3.51	7.19	11.59	17.65	1454	1.00
v_mis[5]	10.67	0.10	4.08	2.65	7.83	10.73	13.40	18.61	1567	1.00
s_Y	8.07	0.08	3.22	4.30	5.88	7.31	9.35	16.55	1594	1.00
mu	62.69	0.06	3.17	55.12	61.06	63.14	64.83	67.67	2652	1.00
lp__	-13.70	0.07	2.19	-18.59	-14.95	-13.34	-12.06	-10.49	861	1.00

測定で得られたYの値はかなり高いので、28と30番目の品物（それぞれv_mis[3]とv_mis[5]）が価値最大化の品物リストに加わる可能性が高く、ともに事後平均が10ぐらいになっています。

ちなみにナップサック問題のバリエーションには以下のものがあるようです。解き方は一例です。

• 各品物の個数が無限個の場合、あるいは個数に制限がある場合。→漸化式を工夫して解く（蟻本のp.58, p.302）
• 重さが連続値で価値が整数値の場合。→dp[i, v]を「i番目までの品物から価値の総和がvとなるように選んだときの、重さの総和の最小値（そのような解が存在しない場合は十分大きな値INF）」と定義する（蟻本のp.60）。
• 各品物の重さと価値が連続量の場合。→僕が調べた限りでは、枝刈りとかで地道に探索するかないっぽい…。詳しい人いましたら教えてください。
• 各品物を液体とし、個数が離散値ではなく好きな量だけ入れることができる連続値の場合。→東北大の塩浦先生の資料が詳しかったです(pdf)。