StatModeling Memorandum

StatModeling Memorandum

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

Osaka.Stan#5で「MCMCサンプルの使い方 ~見る・決める・探す・発生させる~」というタイトルで話しました

先日、以下のイベントで話しました。

『StanとRでベイズ統計モデリング』読書会(Osaka.Stan#5) : ATND

発表資料は以下です。

理論的には事後分布や予測分布の使い方というのが正しいですが、プログラミング言語との相性を考えてMCMCサンプルの使い方というタイトルにしました。自著ではモデリングのやり方の体得にフォーカスしていますが、事後分布や予測分布が得られるメリットについては分野や人によって異なるので詳細は省きました。いつか補おうと思っていたので良い機会でした。

読書会では、小杉先生の発表やLTもめちゃ面白く、東京のStan勉強会では見たことがない盛り上がりを見ました。ネット上でしか知らなかったベイジアンにたくさん会って話すことができてよかったです。調子に乗って3次会まで行きました。

清水先生、三浦先生をはじめ、みなさまありがとうございました!

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

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

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

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

データの分布の確認

difficultyの分布

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

f:id:StatModeling:20201106170105p:plain

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

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

f:id:StatModeling:20201106170108p:plain

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

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

f:id:StatModeling:20201106170113p:plain

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

統計モデリング

このような問題は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のユーザを示します。

f:id:StatModeling:20201106170126p:plain

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

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

q_skip vs. pfの比較

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

f:id:StatModeling:20201106170131p:plain

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

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

f:id:StatModeling:20201106170121p:plain

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

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

f:id:StatModeling:20201106170116p:plain

例えば問題番号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

  • 発売日: 2017/02/16
  • メディア: 単行本(ソフトカバー)

参考文献

  • [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項に記述があります。以下の多変量正規分布があるとします。

f:id:StatModeling:20201106171132p:plain

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

f:id:StatModeling:20201106171136p:plain

f:id:StatModeling:20201106171141p:plain

f:id:StatModeling:20201106171145p:plain

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

f:id:StatModeling:20201106171148p:plain

f:id:StatModeling:20201106171153p:plain

shrinkage factorの導出

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

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

f:id:StatModeling:20201106171157p:plain

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

f:id:StatModeling:20201106171202p:plain

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

f:id:StatModeling:20201106171205p:plain

 * * *

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

f:id:StatModeling:20201106171209p:plain

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

f:id:StatModeling:20201106171213p:plain

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

f:id:StatModeling:20201106171217p:plain

f:id:StatModeling:20201106171221p:plain

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

f:id:StatModeling:20201106171226p:plain

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

f:id:StatModeling:20201106171229p:plain

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

f:id:StatModeling:20201106171234p:plain

f:id:StatModeling:20201106171239p:plain

f:id:StatModeling:20201106171242p:plain

f:id:StatModeling:20201106171246p:plain

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

f:id:StatModeling:20201106171249p:plain

 * * *

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

f:id:StatModeling:20201106171252p:plain

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

f:id:StatModeling:20201106171256p:plain

f:id:StatModeling:20201106171300p:plain

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

f:id:StatModeling:20201106171304p:plain

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

f:id:StatModeling:20201106171307p:plain

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 f:id:StatModeling:20201106171311p:plain f:id:StatModeling:20201106171322p:plain
- Student-t*1 f:id:StatModeling:20201106171315p:plain f:id:StatModeling:20201106171326p:plain
horseshoe f:id:StatModeling:20201106171334p:plain f:id:StatModeling:20201106171319p:plain *2 f:id:StatModeling:20201106171330p:plain

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

f:id:StatModeling:20201106171128p:plain

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:正の部分だけ定義された半コーシー分布