Glickman, M.E., Jensen, S.T., (2005) Adaptive paired comparison design. Journal of Statistical Planning and Inference, 127, 279-293.
仕事の都合で一対比較法の実験計画の動的最適化について考えていて(なぜこんな柄にもないことを考えているのか…)、タイトルがそのものずばりだったのでめくってみた。Google様いわく被引用件数73件。この分野では多いのだろうか、少ないのだろうか。
本論文はあいにくトーナメント戦の設計を考えていて(大相撲の15日間の各日の取り組み表を決めるというような状況)、私が抱えているところの、調査における刺激の心理量の一対比較測定という問題とは、いっけん乖離が著しいのだが、でも本質的には同じ話なんだろうなと思う。
1. イントロダクション
一対比較の計画で、プレイヤーのmeritsについての事前知識が使える場面というのがある。本論文ではトーナメントのスケジューリングにおいて先行する情報を使う方法について述べる。
トーナメントのスケジュールリングに用いられている一対比較計画はふつう、meritが最高であるプレイヤーを推論するために計画する。よく使われるのがノックアウト・トーナメント(削除トーナメント)とラウンドロビン・トーナメントである。
さらによく使われている計画として、それぞれのプレイヤーを削除なしで複数回比較する(プレイヤーあたり比較数はプレイヤー数よりずっと少ない)、という計画がある。たとえば「スイス・システム」がそうである。このやり方は、トーナメント中の累積の結果が似ているプレイヤーをペアにするのが好ましいという、合理的だがアドホックな考え方に基づいている。スイス・システムの場合は、累積の結果が似ていて事前の強さができるだけ異なるようなペアを組む(最高のプレイヤーたちが最初にぶつかるのを避けるため)。
我々のアプローチはスイス・システムの代替案としてみることができる。ペア集合の決定を、meritの事前密度と事後密度のKL距離の期待値が最大化するように決める。我々はこの問題をベイジアン最適計画のひとつとして定式化し、この方法をマルチ・ラウンドのトーナメントで適応的に用いる方法を述べる。
以下ではチェスのトーナメントの用語を使って説明するけれど、いくつかの比較セットが系列的に生じる一対比較実験の全般に使える手法である。
2. 最適トーナメント計画の方法
以下で述べる方法は１ステップ先の計画を最適化する適応的手法なので、そこんところよろしく。[初日に千秋楽までの取り組みを決める気はさらさら無くて、本日の中入り後の勝敗を観た後で翌日の取り組みを決めるような状況を考えているってことね]
2.1 ペアリングの効用関数
\(N\) (奇数) プレイヤーのトーナメントのあるラウンドがあって、ゲーム\(k\)でプレイヤー\(i_k, j_k\)が試合するとしよう。\(i_k\)が勝ったら\(Y_{i_k, j_k} = 1\)とし負けたら0とする。引き分けはないとする。
$$ p_k = P(Y_{i_k, j_k} = 1|\theta_{i_k}, \theta_{j_k}) = F(\theta_{i_k} – \theta_{j_k})$$ と仮定する。\(F\)はなんらかの確率分布関数で単調増加する。Brandley-Terryモデルなら\(F\)は標準ロジスティック分布、Thurstone-Mostellerモデルなら\(F\)は標準正規分布関数である。説明の都合上、サポートは実数全体で、\(\theta_i\)が大きいと勝率が大きく、その値は無制約とする。
プレイヤーたちの強さを\(\Theta = (\theta_1, \ldots, \theta_N)\)と書く。\(\Theta\)は平均\(\mu\), 共分散行列\(\Sigma\)のMVNに従うとする。事前知識は多変量正規で、各プレイヤーについて独立なスカラー密度として表現できるとする[つまり事前分布もMVNだけど共分散はゼロってこと？？ いまいちぴんとこないが先に進もう]
\(\Theta\)の事前分布と事後分布の期待KL距離を最大化するペアリングを考える。つまりKL距離が効用関数である。
ペアリングは\(N/2\)ペア出来るわけだが、その特定の集合を\(s\)と書き、これを計画と呼ぶ。すべての計画の空間を\(\mathcal{T}\)とする。個々の\(s\)について、潜在的に可能なアウトカムの二値ベクトル\(\mathbf{y}\)の\(2^{N/2}\)本のコレクションを\(\mathcal{Y}_s\)とする。
準備が出来ました。ある計画の効用を次のように定義する。$$ U(s) = \int \sum_{\mathbf{y} \in \mathcal{Y}_s} \log \left( \frac{p(\Theta|\mathbf{y})}{p(\Theta)} \right) p(\mathbf{y}|\Theta) p(\Theta) d\Theta$$ で、すべての\(s\in \mathcal{T}\)について\(U(s^*) \geq U(s)\)であるような計画\(s^*\)をみつけるのがゴールである。
[あー、なるほど、云ってる意味がわかった。明日の結びの一番が終わった後で\(\Theta\)をベイズ更新するとき、事前と事後のKL距離がでかくなりそうな取り組みにしようぜ、ってことね。ってことは、実力伯仲のカードばかりでもだめだし、勝敗がほぼ見えているカードばかりでもだめなわけだ。うーん、これでいいんだろうか。目指すべきは、千秋楽に向けて事後分布の分散がどんどん小さくなることなんじゃないの？ これだと、何でもいいから事後分布が変わりそうな奴にすればいいってことにならない？ もっとも、期待KL距離を効用にするというのはベイジアン最適計画の一般的なフレームワークらしいので、私がなにか理解し損ねているんだと思うけど…]
これを書き換えると、\(q_k = 1-p_k\)とし、\(E(\cdot)\)を\(\Theta\)のMVN事前分布に関する期待値として、こうなる。詳細はAppendixをみてね。$$ U(s) = \sum_{k=1}^{N/2} \left( E(p_k \log p_k) + E(q_k \log q_k) – E(p_k) \log E(p_k) – E(q_k) \log E(q_k) \right) $$ このように、KLダイバージェンスは結局はゲームごとの貢献の和になる。
[狐につままれたような気分なので、Appendixで示された式変形を追いかけてみる。原文を勝手に補足している。
まず、\( p(\Theta|\mathbf{y}) p(\mathbf{y}) = p(\mathbf{y}|\Theta) p(\Theta) \) より \( \frac{p(\Theta|\mathbf{y})}{p(\Theta)} = \frac{p(\mathbf{y}|\Theta)}{p(\mathbf{y})} \)である。元の式に代入して $$ U(s) = \int \sum_{\mathbf{y} \in \mathcal{Y}_s} \log \left( \frac{p(\mathbf{y}|\Theta)}{p(\mathbf{y})} \right) p(\mathbf{y}|\Theta) p(\Theta) d\Theta$$ 対数の内側に分数が入っているのを対数の引き算に書き換えて$$ = \int \left[ \sum_{\mathbf{y} \in \mathcal{Y}_s} \log( p(\Theta|\mathbf{y}) )p(\mathbf{y}|\Theta) \right] p(\Theta) d\Theta – \int \left[ \sum_{\mathbf{y} \in \mathcal{Y}_s} \log( p(\mathbf{y})) p(\mathbf{y}|\Theta) \right] p(\Theta) d\Theta $$ $$ = \int \left[ \sum_{\mathbf{y} \in \mathcal{Y}_s} \log( p(\Theta|\mathbf{y}) )p(\mathbf{y}|\Theta) \right] p(\Theta) d\Theta – \sum_{\mathbf{y} \in \mathcal{Y}_s} \log( p(\mathbf{y})) p(\mathbf{y}|\Theta) $$
ここから、この式を第一項と第二項についてそれぞれ考えるぞ。
まずは第一項。積分記号の内側を周辺化する。\(\mathbf{y}\)とは長さ\(N/2\)の二値ベクトルであったから、$$ \sum_{\mathbf{y} \in \mathcal{Y}_s} \log( p(\Theta|\mathbf{y}) )p(\mathbf{y}|\Theta) $$ $$ = \sum_{y_1 = 0}^1 \sum_{y_2 = 0}^1 \cdots \sum_{y_{N/2}= 0}^1 p(y_1, \ldots, y_{N/2} | \Theta) \log p(y_1, \ldots, y_{N/2}|\Theta) $$ $$ = \sum_{y_1 = 0}^1 \sum_{y_2 = 0}^1 \cdots \sum_{y_{N/2}= 0}^1 \left( \prod_{k=1}^{N/2} p(y_k | \Theta) \sum_{k=1}^{N/2} \log p(y_k|\Theta) \right)$$ 上の式では\(y_1, \ldots, y_{N/2}\)をインデクスにした多重ループのなかで\(k\)をインデクスにしたループを2個回して掛けているけど、そのうち2個目の\(k\)のループを外側に出す。1個目のループのインデクスは\(m\)に書き換える。$$ = \sum_{k=1}^{N/2} \left( \sum \sum \cdots \sum \log p(y_k|\Theta) \prod_{m=1}^{N/2} p(y_m|\Theta) \right)$$ ここからがヤヤコシイので良く聞け。大括弧の内側に注目しよう。\(k\)はもう決まっていて、\(y_1, \ldots, y_{N/2}\)をインデクスにした多重ループのなかでごにょごにょしている。ごにょごにょというのは2つの項の掛け算だ。すなわち、\(\log p(y_k|\Theta)\)と\( \prod_{m} p(y_m|\Theta)\)である。いや、\(\log p(y_k|\Theta)\)と\( \prod_{m \neq k} p(y_m|\Theta)\) と\(p(y_k|\Theta)\)の3つを掛けているといってもいいですね。このうちひとつめとみっつめは、\(y_k\)だけで決まる。だから、多重ループの順番を変えて\(y_k\)を一番外側のループで回し、ひとつめとみっつめを一番外側のループに書けば、こうなるわけよ。$$ = \sum_{k=1}^{N/2} \left( \sum_{y_k = 0}^1 p(y_k|\Theta) \log p(y_k|\Theta) \left( \sum \cdots \sum \prod_{m \neq k} p(y_m|\Theta) \right) \right) $$ 内側の大括弧のなかにある\(\sum \cdots \sum\)には、\( \sum_{y_k = 0}^1 \)以外のすべてが書いてあるわけ。この総和ってなにを表しているかというと、たとえば\(k=3, N/2=4\)として、\(p(y_1=0,y_2=0,y_4=0|\Theta)\)と\(p(y_1=0,y_2=0,y_4=1|\Theta)\)と\(p(y_1=0,y_2=1,y_4=0|\Theta)\)と…\(p(y_1=1,y_2=1,y_4=1|\Theta)\)の和だ。つまり1だ。従って$$ = \sum_{k=1}^{N/2} \sum_{y_k=0}^1 p(y_k|\Theta) \log p(y_k|\Theta) $$ となる。
第二項も同じように処理できる。$$ \sum_{\mathbf{y} \in \mathcal{Y}_s} \log( p(\mathbf{y})) p(\mathbf{y}|\Theta) = \sum_{k=1}^{N/2} \sum_{y_k}^1 p(y_k) \log p(y_k)$$
従って、$$ U(s) = \int \left[ \sum_{k=1}^{N/2} \sum_{y_k=0}^1 p(y_k|\Theta) \log p(y_k|\Theta) \right] p(\Theta) d\Theta – \sum_{k=1}^{N/2} \sum_{y_k}^1 p(y_k) \log p(y_k)$$ $$ = \sum_{k=1}^{N/2} \sum_{y_k=0}^1 \left( \int p(y_k|\Theta) \log p(y_k|\Theta) p(\Theta) d\Theta – p(y_k) \log p(y_k) \right)$$ ひとつめの総和記号のなかには次の4つが入っていることになる。$$ \int Prob(y_k = 0|\Theta) \log Prob(y_k = 0|\Theta) p(\Theta) d\Theta = E(q_k \log q_k)$$ $$ – Prob(y_k = 0) \log Prob(y_k = 0) = E(q_k) \log E(q_k)$$ $$ \int Prob(y_k = 1|\Theta) \log Prob(y_k = 1|\Theta) p(\Theta) d\Theta = E(p_k \log p_k)$$ $$ – Prob(y_k = 1) \log Prob(y_k = 1) = E(p_k) \log E(p_k)$$ というわけで、$$ U(s) = \sum_{k=1}^{N/2} \left( E(p_k \log p_k) + E(q_k \log q_k) – E(p_k) \log E(p_k) – E(q_k) \log E(q_k) \right) $$
やれやれ、終わった…]
一対比較モデルの選び方によるけれど、ふつう総和記号の内側は解析的には求まらない。でもガウス・エルミート求積法で数値的に評価できる。[ええええ？ ここでかたくなに数値積分すんの？ てっきりMCMCかと思ってた… でも意外に楽なのかも]
2.2 最適計画の決定
というわけで、ある計画の効用関数は、個々の試合の貢献の和として評価できる。\((i, j)\) の試合の貢献は、\(i\)が勝つ確率を\(p_{ij}\), 負ける確率を\(q_{ij}\)として $$ C_{ij} = E(p_{ij} \log p_{ij}) + E(q_{ij} \log q_{ij}) – E(p_{ij}) \log E(p_{ij}) – E(q_{ij}) \log E(q_{ij}) $$ である。\(C_{ij} = C_{ji}\)である。
これを最大化するペア集合をみつけるという問題は、グラフ理論でいうところの最大ウェイト完全マッチング問題だ。つまり、個々のプレイヤーが頂点で、エッジがゲームごとの貢献で、それぞれの頂点がエッジを一本だけ持ち(完全マッチング)、かつエッジのウェイトの合計が最大になるようなグラフを求めればよい。そのためのアルゴリズムは、安心して下さい、たくさんありますよ。
もし\(N\)が奇数だったら… [パス]
2.3 ペア内の順序
[先手と後手を決めないといけない場合はどうすんのかという話。パス]
3. Bradley-Terryモデルの最適計画
みなさまご存じBTモデルってのは、プレイヤー\(i\)の強さを\(\theta_i\)として$$ P(Y_{ij} =1 | \theta_i, \theta_j) = \frac{\exp(\theta_i)}{\exp(\theta_i)+\exp(\theta_j)} = \frac{1}{1+\exp(-(\theta_i-\theta_j)} $$ ですね。よし、ではこれについて考えてみよう。
それまでのラウンドの結果を使ってMVN事前分布をつくる。これはLeonard(1973)によるベイジアンBTモデルでできる。ほかにもDavidson & Solomon(1973), Chen & Smith(1984) というのがあるけれど、Leonardのアプローチではこう考える。事後密度は、MVN事前密度とロジスティック確率のproduct-binomial尤度の積に比例する。だから、実際の事後分布の最初のふたつの積率を保持するMVNで近似できる。[わ・か・ら・ん。product-binomialってどういこうこと？ 学力不足だ…]
我々のセッティングでは、ひとたびゲームの結果が観察されたならば、近似された正規事後分布を次のラウンドの事前分布として使うことができる。
\(C_{ij}\)はガウス・エルミート求積で評価できる。\(\theta_i – \theta_j\)の事前平均と分散の関数となる。グリッド点が30もあれば十分。[あー、そうか… ラウンドにもプレイヤーペアにも関わらず、現時点での平均と分散を入れれば\(C_{ij}\)が出てくるような関数をつくれるのか。そりゃ楽だね]
観察するとわかるけど、差の事前平均が小さくてSEが大きいペアでC_{ij}が大きくなる。
4. デザイン実装の比較
他の方法と比較します。選手入場。
- 提案手法、無制約版。
- 提案手法、同じペアでは二度試合しない版。
- Olafsson(1990)による、スイス・システムの修正版。[これも適応的計画なのかなあ… たぶんそうなんだろうな]
[ここからちょっと関心があるので細かくメモ] プレイヤー数は50で、\(\theta_1=-2.5\)から\(\theta_{50}=2.5\)まで均等に並んでいるとする。
試合結果は次のように生成する。\(X_i, X_j\)を独立な指数確率変数としその平均を\(\exp(\theta_i), \exp(\theta_j\)\)とする。すると$$ Pr(X_i > X_j) = \frac{\exp(\theta_i)}{\exp(\theta_i)+\exp(\theta_j)}$$ がBTモデルの勝率になるので、\(X_i \gt X_j\)だったら\(i\)の勝ちとする。このアプローチは、それぞれのプレイヤーがゲームごとに「パフォーマンス」を示し、パフォーマンスの比較が目指す二値アウトカムを生んでいると理解できる。このアプローチの重要な特徴のひとつは、ゲームのアウトカムの生成が、ペアリングの方法とは独立に生成されており、プレイヤーは相手とは無関係にゲームのパフォーマンスを生むと仮定されているということである。
[えーっと、つまり指数分布に従う乱数を2つ生成してその大小を試合結果にするってことだよね。これって、\(\frac{\exp(\theta_i)}{\exp(\theta_i)+\exp(\theta_j)}\)を成功確率とする二値乱数を生成してそれを勝ち負けにする、っていってんのと同じ事だよね？ そうだよね？]
シミュレーション実験で動かす要因は、トーナメントの持続時間(4, 8, 16ラウンド)、色制約(奇数ラウンドでプレイヤーを黒か白に分け、次の試合では同じ色のプレイヤー同士は試合しない)、事前分布(平均は真値で分散が0.3で無相関; 平均がランダムで分散が4で無相関)。
評価は、トーナメント終了後の事後分散の行列式の対数(大きいときに不確実性が高い)、真の順位と事後平均の順位の偏差二乗和。
結果は…[面倒なのでとばしたけど、まあ当然優れてるんでしょうな]
5. 考察
本手法は推測の効率を高め、計算も楽で、組み合わせ最適化はグラフ理論の標準的問題に帰着する。
タイがある場合にも拡張できる[…中略…]
シミュレーションの結果は提案手法が事後分布の全分散の観点からみて優れていることを示した。
本手法では最初のラウンドから強さが似たペアの試合が組まれやすくなる(この点がスイス・システムと異なる)。事前分散が大きい場合には、十分なラウンド数が必要になる、という点に注意されたい(多くのトーナメントでは、たとえプレイヤー数が多くとも、せいぜい6～8ラウンドしか組まれない)。
——————
うううう… 備忘のため、わからんことを書き出しておく。恥をさらしているような気がするけれど。
要するにここでいっているのは、各日の結びの一番が終わったら、はね太鼓を聞きながら普通のBTモデルで各プレイヤーの\(\theta\)とそのSEを推定して(いったいどうやってんのか知らないけどソフトで出ますよね)、明日組みうる全ての取り組みペアについて\(\theta_i – \theta_j\)とそのSEを推定し、事前に用意した表を引いて\(C_{ij}\)を調べ、それをウェイトとして最大ウェイト完全マッチング問題を解け… ってことであってます？ ちがうの？
もしあってんなら、誰かその表をつくってくれ… 誰か… 頼む…
それはともかく、いろいろ勉強になりましたです。
ベイジアンBTモデルって70年代からあるのね、びっくり。これを読むまでうっかり忘れてたけど、要はRのbpcsパッケージのことだよね？