読者です 読者をやめる 読者になる 読者になる

StatModeling Memorandum

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

しょラーさんのブログ記事「StanでAizu Online Judgeの難易度・習熟度を推定したい」の追加解析

背景やデータはしょラーさんの以下のブログ記事を読んでください。 kujira16.hateblo.jp この記事ではAOJ-ICPCで付加された貴重な難易度の情報をフル活用して、問題の真の難易度の推定と、各ユーザの習熟度の推定を行います。

この問題の難しさは「解いていない問題が、スキップして取り組んでいないのか、解こうとしたけど解けなかったのか区別できない」という点にあります。そこで、元記事にもあったように問題をスキップする確率を導入してモデリングする必要があります。

とはいえ、まずはモデルのヒントになりそうなグラフを作成します。 以下では元記事にあわせて、難易度をdifficulty(StanコードではD)、習熟度をperformance(Stanコードではpf)と表現します。

データの分布の確認

difficultyの分布

横軸にdifficulty、縦軸に問題の数をとったヒストグラムは以下になります。山型の分布、途中から100刻みしかない、1000を超えると問題数が減ってくる、などが分かります。

問題の解かれた割合(%)

横軸にdifficulty、縦軸に今回データに含まれる全ユーザー(1000人)のうち何人がその問題を解いたかの割合(%)をとった散布図は以下になります。1つの点が1つの問題です。difficultyがあがると解かれる問題の割合がなんとなくですが指数関数的に落ちていきます。簡単な問題でも解いている人が少ない問題もあります。はじめはこの性質をモデルに取り込もうとしましたが、うまくいきませんでした。

解いた問題のdifficultyの平均と標準偏差

横軸にユーザが解いた問題のdifficultyの平均、縦軸に同じく標準偏差をとった散布図は以下になります。1つの点が1つのユーザです。簡単な問題をスキップしまくって難しい問題だけチャレンジする集団(図だと右下らへん)とかいるのかなと思ってグラフを作りましたが、思ったよりみなさん幅広く解いています。

このように背景知識から仮説を確かめていく過程でグラフを作ることは大変有効です。

統計モデリング

このような問題はIRT(Item Response Theory, 項目応答理論)というロジスティック回帰の一種で扱うことが一般的です*1。しかし、個人的にはこのような対戦ムード(問題 vs ユーザ)があるような現象に関してはプロビット回帰を使うのがよいと思っています。理由は「ユーザが対戦相手の力量を上回る(勝利となる)確率」が、累積正規分布で気持ちよく表現できるからです。また、ある説明変数を変えた場合にオッズの観点で議論になりにくいと考えているためです。そこで、この記事ではプロビット回帰(すなわちしょラーさんの2つ目の記事にあるモデルとほぼ同じ)を使います。実用上はロジスティック回帰とプロビット回帰はそこまで差がないと思うのでどちらを使ってもよいと思います。

Stanコードは以下になりました。

  • 2~3行目: データに含まれる問題数をQ(ここでは608)、ユーザ数をN(ここでは1000)で宣言しました。
  • 4行目: AOJ-ICPCで付加された難易度です。あとで1000で割ってスケーリングして渡します。
  • 5行目: 元記事ではGに対応します。問題数×ユーザ数の2次元配列です。あるユーザがある問題を解いている場合に1、その他の場合に0となっています。
  • 6行目: ユーザが解いたことがある問題のうち、最も難しかった問題の難易度です。後で使います。
  • 16行目: 問題をスキップする確率です。ざっとデータを見るとユーザごとに大きく異なりそうだったので、ユーザごとに宣言します。以降で事前分布は設定しないので一様分布に従います。
  • 20行目: 今回は最も難しい問題の難易度が1200で、それ以上のperformanceは正確に測定できない(とても大きな値を取るか分からない)はずなので、ユーザのperformanceは正規分布に従うとしました。この仮定があると、解いた問題数が少ないユーザがいても、mu_pfあたりに推定してくれて推定が安定します。なお、正規分布の代わりに student_t(6, mu_pf, s_pf) のような少しだけ裾が長い分布を試してもほぼ同様の結果でした。
  • 21~22行目: 一応階層モデルにしてあります。
  • 15, 23行目: sigmaはプロビット回帰で使用する累積正規分布標準偏差です。対戦を扱うプロビット回帰においては、問題のdifficultyとユーザのperformanceに差がある場合に、どれぐらい勝負のアヤがあるかを表していると解釈できます。sigmaが小さいと少しでも差があると強い方が順当に勝つことが多く、sigmaが大きいと差が少々あっても確率的に弱い方が勝つことがあるといった具合です。IRTにおける「識別パラメータ」に相当します。本来はsigmaは問題ごとに推定できるとよいのですが、今回は推定が厳しかったので全問題で共通のsigmaとしました。
  • 14, 24行目: d_trueは問題の真の難易度です。元記事のようにDをそのまま使うことも考えられますが、一般的に人がつけたものはキリのいい数字に偏りやすく、また、誤差を含んでいると考えた方がよいでしょう。そこで、24行目では平均D標準偏差0.1正規分布に従うとしています。標準偏差0.1は元の難易度が100ぐらいはブレるかなと考えていることに相当します。なお、標準偏差0.05でも実行してみましたが、そこまで大きな違いはありませんでした。経験ではこのような変数を導入することで、Dのままだとどうしても矛盾してしまうようなところがフニャリと解消されて推定が安定化することが多いです。
  • 28~29行目: 問題が解けた場合です。ユーザが問題をスキップしないで解いた(勝利した)確率になります。なお、log1m(x)log(1-x)です。
  • 31~32行目: 問題を解いていない場合です。ユーザが解いたことがある最も難しい問題の難易度よりD_range下回っている場合、簡単すぎてつまらないからスキップしているとみなします。
  • 33~34行目: 問題を解いていない場合です。ユーザが解いたことがある最も難しい問題の難易度よりD_range上回っている場合、難しすぎてチャレンジしても解けないとみなします。
  • 35~39行目: 問題を解いていない場合で、その間の難易度の問題はスキップしたかチャレンジして解けなかったの混合分布になります。

31~39行目の仮説はかなり大胆ですが、このようにデータを稼がないと、「解いていない問題が、スキップして取り組んでいないのか、解こうとしたけど解けなかったのか区別できない」という問題を打破できずにq_skipが不自然な値に収束してしまいます。結局ここでハマってトータル50個ぐらいモデルを試行錯誤しました。このようにある仮説に従ってデータを置き換えるのは統計モデリングでは常套手段で、「StanとRでベイズ統計モデリング (Wonderful R)」では5.3.3項、「予測にいかす統計モデリングの基本―ベイズ統計入門から応用まで (KS理工学専門書)」では8.3節で扱っています。なお、if文の中のDd_trueにすると、if文のなかにパラメータを含む推定が非常に厳しいモデルとなってしまい、今回だと推定できなくなります。

以下はキックするRコードです。

  • 20~21行目: DD_maxはスケーリングして渡します。D_rangeは仮説および推定結果のsigmaとの兼ね合いになりますが、0.20.25ぐらいでsigmaより十分に大きい値となりましたので0.2としました(元の難易度で200)。
  • 24行目: ここでは初期値を設定していませんが収束しました。ただし、chainによっては時間がかかる場合があったので初期値を init=function() { list(pf=max_difficulty/1000, d_true=d_ori$difficulty/1000) } のように定めた方がよいかもしれません。

結果

Solvedの要素数が60万を超えていることもあり、推定に要した時間はおよそ13時間でした。もう少しデータ増えたら、自動変分ベイズであるADVI使った方がよさそうです(もしくはK年後のGPU化Stanを待つ)。

推定されたsigma

推定されたsigmaの値は中央値が0.077、95%ベイズ信頼区間が[0.074, 0.079]でした。2×標準偏差で考えると、元の難易度のスコアにしておおよそ2×77≒150ぐらい差があると、解ける・解けないがはっきりするという解釈になります。与えたD_rangeはそれよりも50ほど大きな値になっています。

ユーザランキング Top 50

推定した習熟度(performance, pf)のMCMCサンプルの中央値Top 50のユーザを示します。

凡例は元記事にあわせてあります。すなわち、ヒゲが95%ベイズ信頼区間、箱が50%ベイズ信頼区間、真ん中の印が中央値です。問題を多く解いている人はベイズ信頼区間がせまくなっているのが分かります。

また、このモデルでは難易度重視(難しい問題が解けるか)でユーザランキングが決まります。例えば、下から13番目のasi1024さんは問題を非常に多く解いているため、現時点でAOJ-ICPCのランキングでは1位です。しかしながら、今回のモデルですと、解いている問題が多い=q_skipが小さい、それにもかかわらず1000より難易度が高い問題を比較的あまり解いていない=スキップではなく解けなかった確率が高い、と解釈されてpfが若干小さくなります。すなわち実力ある人が簡単な問題を多めに解くと損になります。「昔から真面目に解いてきたけど、忙しくなって最近出た難しい問題は(解けるにもかかわらず)着手できていない」場合が損になるのを避けたい場合には、問題が発表された日時や、そのユーザのアクティブ日時などをモデルに組み込むと改善する可能性があります。

q_skip vs. pfの比較

横軸にq_skipの中央値、縦軸にpfの中央値をとった散布図は以下になります。Top 20だけラベルを付けました。

AOJ-ICPCによる難易度と推定された難易度の比較

横軸にD、縦軸にd_trueの中央値をとった散布図は以下になります。1つの点が1つの問題です。ちゃんとy=xの直線に載っているのでそこまで大きくは変わらないことが分かります。

AOJ-ICPCによる難易度と推定された難易度が大きく異なるTop 30

横軸・縦軸はひとつ前のグラフと同じです。中央値だけでなく、95%ベイズ信頼区間をヒゲで、50%ベイズ信頼区間で箱を表現しました。

例えば問題番号2710BNFで数式を定義している問題で、AOJ-ICPCで付加された難易度は400ですが、このモデルでは500超えてるんじゃないかなと推定されています。同様に、問題番号1185のチョコレートの問題は1000という難易度が付加されていますが、このモデルではせいぜい800ぐらいかなと推定されています。

僕はドメイン知識がないので、有識者の解釈を聞いてみたいです。結果はsummaryだけgistにあります。

まとめ

  • StanでAizu Online Judgeの難易度・習熟度を推定しました。
  • ユーザの力量にあった難易度の問題を、easyモード・normalモード・hardモードなどで推薦できそう。
  • 有志によって付加された難易度の情報は非常に貴重で、モデリングの重要な足掛かりとなりました。
  • 問題をスキップしたか、チャレンジして解けなかったかをデータとして保持できれば、よりよいモデリングができそう。
  • AOJだけでなく、問題が公開されていて自由にチャレンジできるタイプのプログラミング教育サイトに適用できそう。

Enjoy!

*1:Stanのマニュアル(v2.15だと8.11. Item-Response Theory Models)にも詳しく載っています。

スパースモデルではshrinkage factorの分布を考慮しよう ~馬蹄事前分布(horseshoe prior)の紹介~

ベイズ統計の枠組みにおいて、回帰係数の事前分布に二重指数分布(ラプラス分布)を設定し回帰を実行してMAP推定値を求めると、lassoに対応した結果になります。また、回帰係数にt分布を設定する手法もあります。これらの手法は「shrinkage factorの分布」という観点から見ると見通しがよいです。さらに、その観点から見ると、馬蹄事前分布が魅力的な性質を持っていることが分かります。この記事ではそれらを簡単に説明します。

なお、lassoそのものに関しては触れません。岩波DS5がlassoを中心にスパースモデリングを多角的に捉えた良い書籍になっているので、ぜひそちらを参照してください。

岩波データサイエンス Vol.5

岩波データサイエンス Vol.5

参考文献

  • [1] C. Carvalho et al. (2008). The Horseshoe Estimator for Sparse Signals. Discussion Paper 2008-31. Duke University Department of Statistical Science. (pdf file)
  • [2] J. Piironen and A. Vehtari (2016). On the Hyperprior Choice for the Global Shrinkage Parameter in the Horseshoe Prior. arXiv:1610.05559. (url)
  • [3] C. Carvalho et al. (2009). Handling Sparsity via the Horseshoe. Journal of Machine Learning Research - Proceedings Track. 2009;5:73–80. (pdf file)
  • [4] Epistemology of the corral: regression and variable selection with Stan and the Horseshoe prior:PyStanの開発者による、PyStanでの実装例とlassoとの比較です。

horseshoe+ priorなんてのも提案されています。

理論編

分割された正規分布

まず分割された正規分布の復習から。例えば、PRMLの2.3.2項に記述があります。以下の多変量正規分布があるとします。

 

ここで、以下のように2つに分割します。

 

 

 

すると条件付き分布は以下になります。

 

 

shrinkage factorの導出

上記の[2]を参考にしました。丁寧な導出が見つからなかったので、以下は手計算で求めました。間違っていましたら連絡ください。

以下の線形モデルを考えます。

 

ここで \overrightarrow{\beta}が回帰係数のベクトルで \bf{X}が説明変数の行列です。多変量正規分布で表現すると以下になります。

 

ここで、 \bf{\Lambda}は以下の対角行列です。

 

 * * *

さて、ここで以下の \overrightarrow{z}を考えます。

 

同時分布の対数は以下になります。

 

これは \overrightarrow{z}の要素の2次関数なので、 p(\overrightarrow{z})正規分布になります。PRMLの2.3.3項とほぼ同じように、2次の項と1次の項にわけて考えることで、その正規分布の精度行列と平均ベクトルを求めることができます。まずは2次の項を整理すると、

 

 

となり、精度行列 \bf{T}が求まります。1次の項はないので、平均ベクトルは以下になります。

 

以上より、前述の「分割された正規分布」の公式より、 \overrightarrow{y}が与えられたもとでの \overrightarrow{\beta}の分布は以下になります。

 

ここで、平均ベクトルは以下のように計算できます。

 

 

 

 

ここで、 \overrightarrow{\hat{\beta}} \overrightarrow{\beta}の分布を考えない場合の最尤推定の解で以下です。

 

 * * *

さて、ここで各説明変数列が相関がなく、平均がゼロで、分散が1とします。すなわち、以下を仮定します。

 

すると、平均ベクトルは以下のように変形できます。

 

 

対角行列になるので要素ごと表すと、以下になります。

 

ここで、 \kappa_jは以下となり、shrinkage factorと呼ばれます。

 

shrinkage factorが0に近いと係数が最尤推定値に近くなり(shrinkageしてない)、1に近いと0に近くなります(shrinkageする)。

馬蹄事前分布(horseshoe prior)

上のモデルで各 \lambda_jがある確率分布に従うと仮定します。この確率分布によって、元の \beta_jにどんな確率分布を設定したのと等価になるのか、そしてshrinkage factorの分布がどのようになるのか手計算で求めることができます(参考: 確率変数の変数変換とヤコビアンに慣れる - StatModeling Memorandum)。以下に結果だけ対応表として載せておきます(手計算で求めましたが間違っていたらすみません)。

法名  \beta_jの分布  \lambda_jの分布  \kappa_jの分布
lasso double-exponential
- Student-t*1
horseshoe *2

shrinkage factor( \kappa_j)の分布を描くと以下になります。

horseshoeの場合の \kappa_jの分布形が馬蹄に似ているので馬蹄事前分布(horseshoe prior)と名づけられました。なお、 n \sigma^{-2} \tau^{2} = 1のときはベータ分布 Beta(0.5, 0.5) に一致します。

馬蹄事前分布のようにshrinkage factorが0か1を生成しやすいということは、「0につぶしやすい一方で、0につぶれない係数は0に近づくような振る舞いをしないで最尤推定値に近づく」ことになります。係数にメリハリがあるわけです。これに対し、lassoの場合は n \sigma^{-2} \tau^{2}が1もしくは10だとあまり0につぶす傾向はなく、またつぶれなかった係数は少しだけ最尤推定値に近づく感じになります。 n \sigma^{-2} \tau^{2} = 0.1だと0につぶす傾向はあるのですが、つぶれなかった係数も0に近づけてしまいます。

この結果、C. Carvalhoらはシミュレーション実験から「一部だけ影響のない説明変数があるようなシミュレーションデータに対しては馬蹄事前分布の方が予測の性能がよい」と主張しています([3])。さらに馬蹄事前分布は二重指数分布(ラプラス分布)とは異なり至るところで微分可能なので、偏微分を使うStanでの推定が安定であることから、Stanの開発メンバーはスパースモデリングには二重指数分布よりも馬蹄事前分布の使用をすすめています。

Stanによる実装例

[4]の実装に尽きています。すなわち以下です。簡単です。

data {
  int<lower=0> N;
  int<lower=0> D;
  matrix[N,D] X;
  vector[N] Y;
}

parameters {
  vector[D] beta;
  vector<lower=0>[D] lambda;
  real<lower=0> tau;
  real<lower=0> sigma;
}

model {
  lambda ~ cauchy(0, 1);
  tau ~ cauchy(0, 1);
  for (d in 1:D)
    beta[d] ~ normal(0, lambda[d] * tau);
  Y ~ normal(X * beta, sigma);
}

Enjoy!

*1:参考: http://stats.stackexchange.com/questions/52906/student-t-as-mixture-of-gaussian

*2:正の部分だけ定義された半コーシー分布

『Pythonで体験するベイズ推論 ―PyMCによるMCMC入門―』の書評

特長

Pythonユーザが待ちに待ったPythonによるMCMC本ではないでしょうか。原著タイトルが『Bayesian Methods for Hackers』だけあって、プログラマ・エンジニア向きだと思います。数式はびっくりするほど出てこない代わりに、Pythonコードは非常にたくさんでてきます。そしてPyMCの使い方が基礎から説明してあって丁寧です。自分でコーディングする際は原著のGitHubリポジトリを活用しましょう(なんとStarが10000個を超えてる!)。

Pythonで体験するベイズ推論 PyMCによるMCMC入門

Pythonで体験するベイズ推論 PyMCによるMCMC入門

  • 作者: キャメロン・デビッドソン=ピロン,玉木徹
  • 出版社/メーカー: 森北出版
  • 発売日: 2017/04/06
  • メディア: 単行本(ソフトカバー)
  • この商品を含むブログを見る
購入を迷っている人の一番の心配は、本書のPyMCのバージョンが1つ前のPyMC2であることだと思います。しかし、そこは気にせず、まずは安定でドキュメントも多いPyMC2で勉強したらいいと思います。なぜなら、MCMCを使ったモデリングに慣れていない人のボトルネックは、現象を数式やコードに落とすモデル化という営みそのものに慣れていないことだからです。モデル化に慣れてPyMC2に慣れたら、PyMC3に移行するのはエンジニアの方なら大して時間かからないでしょう。僕自身、BUGS言語(WinBUGS, JAGS)に慣れてからStanに移行しましたが、おおよそ1週間で移行できました。迷わずPyMC2でモデリングライフを開始したらいいと思います。なお、原著のGitHubリポジトリにはPyMC3のコードも含まれています。

ちなみに2015年10月刊行の岩波データサイエンス Vol.1では、渡辺さんの頑張りのおかげでPyMC3の解説(20ページ程度)になっております。ご参考までに。

本書で最高だった点

5章です。秀逸すぎる。この章では事後分布と損失関数を組み合わせて意思決定する、いわゆるベイズ決定を行います。著者の金融分野における経験が生かされていて、さまざまな損失関数が出てきて有益情報がいっぱいです。以下では例を挙げますが、この他にも多数あります。詳しくは本を読んでください。

政治評論家向きの損失関数

 L(\theta, \hat{\theta})=\frac{|\theta - \hat{\theta}|}{\theta (1 - \theta)}

ここで、 \theta \hat{\theta}はともに [0,1]の範囲です。真の値 \thetaが0,1に近い場合には、優柔不断(ハッキリ0/1で答えない)だと損失が大きくなるような損失関数になっているとのことです。

株価の予測における下振れリスクを表す損失関数

下振れリスクとは「符号の異なる間違った方向に大きく予測すること」で、上振れリスクは「正しい方向に大きく予測すること」。下振れリスクの方が避けたいので、符号をまたいだところから損失がぐっと増えるような非対称な関数形になっているとのことです。図5.5(p.158)参照。

その他にも良かった点

  • ベイズの定理や大数の法則などの基礎となる理論も独特の切り口で扱っていて勉強になりました。
  • 頻繁に出てくる例題が軽妙でよく練られていて面白いです。例を挙げます。
    • ベイズの定理の説明で、農村における性格がこまかい人が将来司書になるか農家になるか問題。
    • 導入の1章において、著者に送られてきたメッセージの数の日次データに対する変化点検出。
    • Kaggleのダークマターハローの座標を予測する問題。銀河のプロットも面白い!
  • 事前分布の選び方で、統計の知識が全くない人から信念を聞き出して分布にする方法(ルーレット法)。

少し残念だった点

  • これはPythonのデータ解析・機械学習の本にありがちですが、可視化のコードの分量がやや多い気がします。matplotlibのせいかもしれません。仕方ないのかもしれません。
  • ベイズ推定の本領発揮とも言える、階層モデルなどのやや高度なモデルは扱っていないのは少し残念です。
  • MCMCの収束に関してやや議論が甘いと思います。他のMCMCソフトウェアの収束判定では \hat{R}(アールハット)がよく使われており、かなりよい指標だと思いますが、PyMCはなぜかchainを複数流して \hat{R}を求めることはしません。そして、書籍からもそこは抜けており、いくつかの解析で収束が怪しそうなものがあります。