2016年7月29日 (金)
久しぶりに整数計画法を使う仕事があって、案の定、大変に往生した。根っからの数学嫌いなのに、なんでこんな目に。
朝から晩まで問題が頭から離れず、日常生活に困る始末だったのだが、そんな日々もようやく一段落しそうなので、読み漁った資料のうちいくつかをメモしておく。
藤江哲也 (2012) 整数計画法による定式化入門. オペレーションズ・リサーチ, 57(4), 190-197.
はじめて整数計画法を勉強したとき(2008年だっけ。嗚呼)、もっとも助けになったのは、フリーの整数計画ソフトウェアLPSolveのマニュアル、そして宮代・松井(2006)によるチュートリアル記事であった。今回はそれらに加えて、著者の宮代先生のwebページ、そして上記記事による定式化についての説明がなによりも助けになった。専門家による社会的発信の価値を痛感した。実にありがたいことであります。
それにしても... 今回のひとつの反省点は、なによりも定式化に注意しないといけないという点であった。整数計画法ってのはですね、全く同じ問題であっても定式化の仕方がいっぱいあり、それによって解の求めやすさががらっと変わってくるのですよ、奥さん。実物を見ると仰天しますよ。いやあ、正直、悪い夢を見ているようであった。
藤江哲也 (2011) 最近の混合整数計画ソルバーの進展について. オペレーションズ・リサーチ、56(5), 263-268.
宮代隆平 (2012) 整数計画ソルバー入門. オペレーションズ・リサーチ, 57(4), 183-189.
はじめて整数計画を勉強したときにこうした論文を読むことができていたら、どんなに助かったかしれない... SCIPやNEOS Serverの存在を知るまでに、かなりの回り道が必要だったのである。
Berthold, T., Gleizner, A.M., Heinz, S., Koch, T., 品野勇治 (2012) SCIP Optimization Suiteを利用した混合整数(線形/非線形)計画問題の解法. ZIB-Report 12-24. Zuse Institute Berlin.
非商用界最強を謳われる整数計画ソルバーSCIPの開発チームがお送りする、なんと日本語によるモノグラフ。日本人の方がいらっしゃるんですね。先生、貴方は日本の誇りです(と、急に国粋主義的になったりして...)
宮本裕一郎 (2012) はじめての列生成法. オペレーション・リサーチ, 57(4), 198-204.
仕事の役には立たなかったし、内容を全部理解できたわけでもないんだけど、頭の体操としてすごく面白かった。暇ができたらもういちど読み直してみたい。文章も、いかにも工学部的なユーモアにあふれていて、とても楽しい。
今回のもうひとつの反省点は、やっぱり基礎ができてないと困る、という点であった。ソルバーが問題をどうやって解いていようが知ったことか、とにかく解さえ手に入ればいいんだと割り切り、branch-and-cut？なにそれ？というスタンスで取り組んでいたんだけど、そういうど素人にはやはり限界がある。こんな私でも勉強できるものなのかなあと心が揺らぎ、教科書を手に入れてみたんだけど(Williams, 2013)、一朝一夕で読めるものではなさそうだ。悲しいけど、人生は有限だしなあ...
論文：データ解析(2015-) - 読了: 藤江(2012); 藤江(2011); 宮代(2012); Berthold, et al.(2012); 宮本(2012) 初夏の整数計画祭りを振り返る
2016年7月12日 (火)
Madigan, D., Gavrin, J., & Raftery, A.E. (1995) Enhancing the predictive performance of bayesian graphical models. Communications in Statistics - Theory and Methods. 24(9), 2271-2292.
都合により目を通した。原文は手に入んなかったので、draftのpsファイルで。
グラフィカル・モデリングによる予測性能を高める２つの手法を提案します、という論文。２つの手法とは、ベイジアンモデル平均(BMA)と、専門家から事前分布をうまく引き出す方法。
まずBMAのほうから。
いま、関心ある量$\delta$について推測するための$k$個のモデルがある。データ$D$の下での$\delta$の事後確率は
$P(\delta | D) = \sum_k^K P(\delta | M_k, D) P(M_k | D)$
右端はモデル$M_k$の事後確率で、
$P(M_k | D) = \frac{P(D|M_k) P(M_k)}{\sum_l^K P(D|M_l) P(M_l)}$
ここで$P(D|M_k)$ってのは、モデルのパラメータを$\theta$として、
$P(D|M_l) = \int P(D|\theta, M_k) P(\theta |M_k) d\theta$
である。[←推定されたモデルの尤度ではなく、パラメータを通した周辺尤度だってことね]
Madigan & Rafery (1994JASA)のアプローチでは、データからの支持があまりに低いモデルはモデル集合から取り除かれる。さらに、オッカムの剃刀のひそみにならい、サブモデルよりも複雑なわりにぱっとしない親モデルも取り除かれる(オッカムの窓)。[←フォーマルに説明してあるけど、メモは省略]
さらに、MCMC model composition (MC$^3$)という方法もあって...[←難しいのでパス。いずれ必要になったら読みます]
[話はちょっと逸れているんだけど、ここで面白いくだりがあったのでメモ。数値例の誤記らしきものを勝手に訂正したけど、もしかすると私の理解に誤りがあるかも]
モデルの予測性能を以下のようにして評価しよう。データ$D$を学習データ$D^S$とテストデータ$D^T$に折半する。で、対数スコアリング・ルールを用い、スコア
$- \sum_{d \in D^T} \log P(d | M, D^S)$
を求め、これを予測性能とする。BMAの場合には、採用したモデル集合を$A$として、
$\ \sum_{d \in D^T} \log \{ \sum_{M_A} P(d | M, D^S) P(M| D^S) \}$
を求めることになるわけだ。
著者らが示した例では、テストデータのサイズは149。BMA(オッカムの窓アプローチ)によるスコアは、最良の単一モデルのスコアと比べて13.8小さかった。$\exp(13.8/149)=1.097$、つまりBMAの性能は約10%優れていた、といえる。
さてここで、ベルヌーイ試行から真の確率$\pi$を推定するという問題について考えよう。仮にわけもわからず$\hat\pi=1/2$とすれば、スコアの期待値は$\pi \log \hat\pi + (1-\pi) \log (1-\hat\pi) = -\log 2$。いっぽう正確に$\hat\pi = \pi$を推測したならば、スコアの期待値は$\pi \log \pi + (1-\pi) \log (1-\pi)$。このスコアの差がモデルの改善である。これをbiased coin prediction scaleと呼ぼう。
このスコアの差がさきほどの13.8と一致するのは、$\pi=0.72$。これがbiased coin prediction scale上でのモデル改善である。つまり、BMAを使わないということは、出目の確率が72:28であるコインを、あたかもまともなコインとして扱っているようなものだ。
[... うーん。かえって話をわかりにくくしているような気がするけど？]
後半戦。
BMAでは最初に$P(M_k)$を与える必要がある。たいていは無情報分布と称して一様分布を使ってんだけど、よろしくないことだ、ぜひ既存知識を有効活用したい。
しかしこれ、実は結構難しい。Lichtenstein & Fischhoff(1980 OBHP)は専門家に観察不可能な量を直接推定してもらっているけど(グラフィカル・モデルの場合で言うとリンク)、それって難しかろう。
間接的にうまく聞く方法はないか。以下のような提案がある。
- Winkler(1967 JASA): ベルヌーイ過程について[←なんと、Winklerはすでに60年代にこううことをやってたのか...]
- Chaloner & Duncan(1983 The Statistician): ベルヌーイ過程について
- Kadane et al.(1980 JASA): 正規線形モデルについて[←あ、これ面白そう]
- Garthwaite & Dickey(1990 JRSS, 1991 J.Behav.DicisionMaking): 正規線形モデルについて
- Laskey & Black (1989 Proc.): ANOVAについて
- Laud et al.(1992 Tech.Rep.): 一般化線形モデルについて
- Chaloner et al.(1993 The Statistician): 生存分析について
お待ちかねの提案手法。ひとことでいうと、まずは一様分布を考え、専門家が提供する「イマジナリー・データ」でこれをベイズ更新する。類似のアプローチとして、線形モデルの文脈でIbrahim & Laud (1994 JASA), Laud et al(1992), 二項分布の文脈で Gavasakar(1998 MgmtSci)がある。
適用する事例はscrotal swelling(陰嚢腫大)である。診断に関連する二値変数が7つある(痛みがあるかとか)。著者のひとりは陰嚢腫大について16年間の実践経験を持っている(患者じゃなくて医者としてね、もちろん)。この人に簡単なコンピュータ・プログラムを操作してもらう。このプログラムは、(1)ランダムに変数を選び、その状態(真偽)をランダムに選ぶ。(2)変数をもうひとつランダムに選び、そのありそうな状態を入力するように求める。(3)全変数が出てくるまで(2)を繰り返す。これを2時間繰り返し、95個のイマジナリー・ケースをつくった。
ここからモデルの事前分布を求め、やおらデータを使ってグラフィカル・モデルをBMA推定したら、予想性能がちょっと上がった由。
[ああ、なるほどな... 全然不思議な話じゃないな。イマジナリー・データが実データを補完したのであろう。逆に言うと、実データが膨大であれば、こういう工夫はあんまり意味がなくなりそうだ]
後半、思ったより簡単な話なんだけど、勉強になった。専門家からの主観確率抽出については、Spiegelhalter et al.(1990 Stat.Sci.)というのを読むとよいらしい。忘れないようにここでメモしておくけど、南風原「続・心理統計学の基礎」はGarthwaite, Kadane, & O'Hagan (2005 JASA)というのを挙げていた。
論文：データ解析(2015-) - 読了：Madigan, Gavrin, & Raftery (1995) 専門家の知識を生かしたベイジアン・モデル平均