StatModeling Memorandum

StatModeling Memorandum

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

『ベイズ統計モデリング ―R,JAGS, Stanによるチュートリアル 原著第2版―』 John Kruschke著、前田和寛・小杉考司監訳

タイトルの本を頂きました。ありがとうございます。僕は原著を少し読んだことがあり、こちらで非常に評判が高い本です。翻訳にもかかわらず原著とほぼ同じ値段で購入できます。

先にJAGSになじみのない方へ説明しておきますと、JAGSはRコアメンバーの一人でもあるMartyn Plummer氏によってC++で開発されたMCMCソフトウェアです。Rから使うのが多数派ですが、PythonからもPyJAGSによって使うことができます。

複雑なモデルでなければStanより収束が早く、離散値をとるパラメータも使えるため、プログラミングがそんなに得意でない人がベイズ統計モデリングをはじめるには一番向いていると思います。最近、再び活発に開発され始めたようで、先日JAGS 4.3.0がリリースされました。

RからのJAGSの使い方としては、hikaru1122さんの記事がよくまとまっていてオススメです。 hikaru1122.hatenadiary.jp

特長

一言で言えば「ベイズ統計モデリングチュートリアル」としてオールインワンとなっている本です。

方向性が類似している和書は、豊田秀樹(2016)『はじめての統計データ分析』です。すなわち、伝統的な統計分析(t検定・分散分析、ロジスティック回帰などの一般化線形モデル)をすべてベイズ統計で焼きなおしていきましょうという方向性です。理論だけでなく、きちんとソースコード(R,JAGS,Stan)の解説もついています。さらに階層ベイズモデルやモデル比較についても章を割いて丁寧に説明しています。

しかし、そういう発展的な話題だけでもないのです。本の前半では、確率とは何か、条件付き確率・同時確率・独立、ベイズ統計の初歩についても説明しているほか、Rの基本、MCMCの基本(数式が多く割と難しめ章ですが竹林先生の訳がいい感じです)、BUGS/JAGSの基本(これまでの和書にはほとんどない情報!)、Stanの基本と、「チュートリアルとしてはこの本さえ読めば他は読まなくていい」というぐらいのオールインワンとなっています。そりゃそこまでやれば770ページになるし、鈍器にもなります。『はじめての統計データ分析』は比較的コンパクトにまとまっているために、初学者にとって本当に「はじめて」と言えるのか疑問な箇所がいくつかありましたが、この本はそのページ数を見ても分かる通り、手取り足取り非常に丁寧で「チュートリアル」のサブタイトルに偽りありません。

ベイズ統計モデリング: R,JAGS, Stanによるチュートリアル 原著第2版

ベイズ統計モデリング: R,JAGS, Stanによるチュートリアル 原著第2版

  • 作者: John K. Kruschke,前田和寛,小杉考司,井関龍太,井上和哉,鬼田崇作,紀ノ定保礼,国里愛彦,坂本次郎,杣取恵太,高田菜美,竹林由武,徳岡大,難波修史,西田若葉,平川真,福屋いずみ,武藤杏里,山根嵩史,横山仁
  • 出版社/メーカー: 共立出版
  • 発売日: 2017/07/22
  • メディア: 大型本
  • この商品を含むブログを見る

一点、類書と異なり注意を要する点は、この本ではベイズ信頼区間(信用区間, 確信区間)の定義として、両端からα/2%ずつ切って定義する方法ではなく、ほぼすべてで最高密度区間(Highest Density Interval; HDI*1)を使っている点です。こだわりと思います。HDIの方が都合の良い場合もあって、本には例として、峰が離れた二峰の事後分布になった場合に、区間が2つになって直感的に解釈できるとあって、なるほどと思いました。

また、この本は図もとても豊富で充実しています。特に以下のようなモデルの図が頻繁にあって分かりやすいです(この図は原著者の以下のブログ記事からの引用です http://doingbayesiandataanalysis.blogspot.jp/2015/12/lessons-from-bayesian-disease-diagnosis_27.html)。

書籍にはこのような可視化のコードも解説があります。しかし、ソースコードGitHubで管理されているわけではなく、以下のサイトからダウンロードするようです。

Software, with programs for book - Doing Bayesian Data Analysis

読者の方が非公式に更新しているリポジトリはあるようです。

GitHub - boboppie/kruschke-doing_bayesian_data_analysis: John K. Kruschke's Doing Bayesian Data Analysis: A Tutorial with R and BUGS

このあたり、訳書の方にも演習の答え(日本語)や図を含めたリポジトリがあるいいなぁと思いますが、どうでしょうか。期待しましょう。

数少ない不満点は、教養のある原著者で洋書にありがちな、流れを妨げる(理解に時間を要する)オヤジギャグ部分がたびたびあること。僕が訳者なら、原著者に「この部分、通じない可能性が高いのでカットしていいですか?」とか問い合わせてカットを試みたいです。まぁ、オヤジギャグが好きな読者もいるので、好みの問題かもしれません。

翻訳について

最近この手の技術関連の翻訳本を読む際に真っ先に調べることがあります。それは訳者の背景・経歴・アウトプットです。いわゆるプロ訳者の方でも、技術の鍛錬を怠っていると、日本語としては読みやすくとも不思議なことに中身がピンとこない翻訳になるのです。普段からその技術を追っていない(or使っていない)ために、「パッケージの方向性」、「著者の熱意がある部分」、「現実の問題への使われ方」などをよく理解できていないからだと思います。結果的にシャープでない曖昧な訳になります。その点、この本の訳者は数多くいますが、実際に自分たちの研究論文を書くために「ベイズ統計モデリング」を理解し使っている人々なので、曖昧な訳で逃げているような箇所はいまのところほとんど見当たらなかったです(オヤジギャグ部分はちょっとつらそうですが…)。また、訳語やトーンもほぼ一貫しており、監訳者二人の厳しい努力が垣間見えます。

*1:Highest Posterior Density Interval; HPDIという用語もあります

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を超えると問題数が減ってくる、などが分かります。

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

横軸に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)にも詳しく載っています。