Kosugitti's BLOG

アンドロイドは正規分布の夢を見るか

この画面は、簡易表示です

stan

メソ研初参戦

メソドロジー研究会@別府に参加。

ツイッターでゆるくつながってる、英語系研究者の研究会で、方法論が主体なので勉強になる。今回は特に、俺の好きな人がご自身の研究についてトークするので、勉強させてもらう身分で。

いや、楽しかった!やはり自分の専門に熱く語る人は、ジャンル関係なく楽しいな。一般発表ぶもんでも、「へえ、ほう」と思うことが多かったし。やはりちょくちょく出かけていかなあかんなと思った。

懇親会では、結婚の様々なスタイル、大学業界の昨今まで広く語り合い(というか愚痴を聞いてくれてありがとうございます)、結局創発が大事、金一元論に落ち着く大団円。

楽しかったです。

研究会や学会というのは、気の置けない間柄の人とやると楽しいし、それで内輪だけの閉じた世界になるんじゃなくて、オープンな感覚を持ってなきゃダメで、そのバランスがうまくいってる会だな、という印象。KSPもそれだよな。いまんとこ、この二事例ぐらいしかしらんけど、そういうのに参加させてもらうの、すごく楽しいですはい。



最強のM-1漫才師は誰だ

はじめに

この記事はStan Advent Calendar2017,12月11日のエントリー記事です。

Stanというかベイズ統計を勉強するには,犬4匹本もいいんですが,コワい本もまたいいですね。

コワい本の中に,7人の科学者というお題がありまして。7人がある実験でスコアをとったのだけど,どうも最初の二人は未熟者で変なスコアを出している。残りの5人の熟練者から考えるとスコアは10.00点ぐらいになるはずなんだけど・・・という話があります。これをモデリングするにあたって,真の値は一つなんだけど,評定者ごとに「測定誤差」を持っている,というモデルを立ててます。すなわち,科学者$$k$$によって得られた計測値を$$X_{k}$$とすると,$$X_{k} \sim N(\mu,\sigma_k)$$というわけです。

これを応用して,「評定者のくせを考慮したスコア」というのを算出することを考えました。 データはそう,M-1グランプリです。

データの出典はこちらです。 そして今年の分はビデオを見ながら手入力しました。

結論;ブラマヨが一番すごいんじゃないか

この結論に至る手続きを,以下解説いたします。

手続き

データを取り込みます。審査員は毎回変わりますから,欠損値が多いデータセットです。 元のファイルも置いておきますね。UTF-8のcsvファイルをこちらからダウンロードしてください。

データは次のようにして,縦長に整形しました。

モデルは先ほどと同様で,漫才師(演者)$$k$$の真のお笑いの実力$$\theta_k$$を持っていたとして,審査員(評定者)$$j$$がつけた得点$$X_{jk}$$は$$X_{jk}\sim N(\theta_k,\sigma_j)$$としています。

これで推定した結果をグラフにします。まずは各演者の「漫才力$$\theta_k$$」のプロット。 点がEAP推定値,横幅は50%と95%のHDIです。

かわいそうなDonDokoDon。っていうか皆さん覚えてます?ぐっさんと平畠のコンビなんですが。

ただ,このモデルの中に含まれていることは,「評定者のスコアは100点満点で絶対的な感覚を持って採点している」ということ。だんだんスコアがインフレしているかも・・・というのは仮定していません。初期のメンバーはスコアが低くなりがちですが,ひょっとしたら審査員もみんな慣れてきて,スコアが上がってきたとかってことがあるかもです。

次に審査員の$$\sigma_j$$の結果。これは歪んでいるのでMED推定値にしています。 

立川談志師匠はまあいいとして,意外と大吉先生は評価の幅が大きい。一番小さいのは哲夫(笑い飯)ですが。審査員歴が一番長い松本人志は5点弱のブレのようですね。

Who are king of kings ?

事後分布から

さて,これらのスコアを元に,もし歴代のM-1選手たちが一堂に会し,実力一発勝負をしたらどうなるでしょうか? MCMCサンプルはバーンイン期間を除いて10000点とりましたので,$$\theta_k$$の同時確率空間における順序で一位を取った回数をグラフにしました。それがこちら。

強いのはブラマヨ。ね。

事後予測分布から

違う角度から検討します。実力の通りではなく,評定者との関係というのもあるでしょうから,事後予測分布を使った予測をして見ます。

全てのプレイヤーが,全ての評定者に評価されたと考え,事後予測分布の平均スコアをもとに予測される得点順に並べたのがこちら!

1位はブラックマヨネーズ!

2位はパンクブーブーとミキが同率!

3位はアンタッチャブル,オードリー,サンドウィッチマン,かまいたち,とろサーモン,トレンディエンジェルが並ぶ!

という結果になりますた。 ちなみに,2017年の結果を入れなければ,ミキは入ってこないし(初挑戦だから当然),ジャルジャルは3位に入れてたんですよw

もちろん,順番の問題もあるでしょう。また複数回ファイナリストになった人の実力はひとつ(成長や変化がない)というモデルだという限界はあります。

なにより,これはあくまでも確率です。実現値は常に一つ。可能性は色々あるんだけど,結果は結果。そこに生まれる物語や感動について,ベイジアンは言葉を持たず,ただ感動に浸り,漫才という芸術に感じ入るしかないものだと思います。

Enjoy Bayes and OWARAI!

 

 

追伸 博多大吉先生の評価についてのコメント,ラジオクラウドで聴けます。期間限定かも。是非〜 > 博多大吉のM-1審査員をやって・・・ https://radiocloud.jp/archive/tama954/?content_id=24438



二階差分の状態空間モデル,それより俺はいつからデブキャラになったんだ

この記事はStan Advent Calendar 2017の12月8日の記事です。

Stanの盛り上がりに一役も二役も買った,「StanとRでベイズ統計モデリング」の,読書会があちこちで行われておりますが,去る11月18日に関西学院大学で行われたOsaka.stanも読書会。この本の11,12章でした。この最後の二つの章は,ちょっと複雑で難しいモデリングでもあって,上手に解説していただいたおかげで自分でもできるかも,という気になりますね。

さて,12章の状態空間モデル,これについては私も去年のアドカレ今年のアドカレでもやってるんですが,まあ手習い程度で複雑なモデリングはしてないのです。で,読書会でスピーカーの @masashikomori 先生が「まあ一階の状態空間による予測って,フエルトも減るとも言わずにわからねーっていう信用区間を示しますから,あんまり意味がないですよね。二階差分ぐらいしないとね」というわけです。ま,たしかにそうだと。式を見てもそんなに難しくないわけです。だからやってみようかな,と思ってこのエントリーもしてみました。

二階差分は「差分の差分」なので,「傾きを滑らかにする」ようにモデリングするから,方向性(上昇,下降)をもった予測ができるよ,と。三階になると加速度になるので現実的な空間モデルではよくわかんないとのこと。たしかにそうです。式としてはややこしく見えるんだけど,
$$ u^t – u^{t-1} = u^{t-1} – u^{t-2} + e $$
から
$$u^t = u^{t-1} + u^{t-1} – u^{t-2} + e$$
$$ u^t = 2u^{t-1} – u^{t-2} + e$$

となるだけで,簡単にできるのです。

さて,私は毎朝体重と体脂肪を測定するという趣味を持っていますから,そのデータを今回も使います。データそのものはかなり長くあるんだけど,これの今年度部分を使いたいと思います。2017年1月1日からの体重変化をモデリング!
(今回使ったデータも含めて公開中です。おヒマな方はどうぞ。)

コードはこちら。
観測が欠損値や予測の箇所は,データに9999というありえない数字を入れて,if文で分岐させるようにしています。

これで予測したのが次の図。

帯の濃い色は50%信用区間,薄い部分は95%信用区間。
なるほど,予測線が下を向いてます。これでできてるし,体重は下降傾向。よしよし。

 ・・・ちょっと待てよ。
薄々気づいてはいたんだけど,最近体重がちょっと増え気味じゃないかなぁ。ベスト体重は66kgですよ。そして結婚してからこっち,それを超えて76kgで安定しているのは知っていますよ。それにしても80kgの大台に乗るのはよろしくないなあ。これ,どこかで「デブモード(Dev Mode?)」に入ってしまったのではないか。これを考えるモデルはないだろうか。
ということで,色々考えて見ました。

1. デブモードモデル
読書会11章で読んだ知識を活用します。俺には「普通モード」と「デブモード」の二種類があって,デブモードのときは普通モードよりも体重が高い水準で推移するとします。普通モードになるか,デブモードになるかは確率$$\pi$$のベルヌーイ分布(0/1)になるとします。そんなモデルがこちら。

betaは下限0,つまりデブモードは増えている前提で,どれぐらい増えているかわからないので裾の思いコーシー分布を置いてみました。結果はこちら。

 

ちょっと収束が悪いところもありますが,betaの中央値が3.21kg,EAPが25kgですな。半コーシーだからだいぶん歪んでます。まあデブモードが平均25キロ増ってのは言い過ぎでしょうよ,このモードに入ったら3キロ太るのか。なるほど。で,このモードに入る確率は4%ぐらい。うーん,滅多にないことなんだけど,入ったら3キロ太るのかー,気をつけないと。

とは思ったんだけど,なんかシックリ来ないんですよねえ。

2. ただ切り替わるんじゃなくて,日付の関数だよね

このモデルだと,ただの確率で二つのモードがコロコロ切り替わるわけです。が,俺が考えたかったのは「あれ?ある日を基準にデブモードに入ってるんじゃない?」ということなのでした。これはどう考えたらいいか。
変化点モデルとか,隠れマルコフモデルとかが該当するんだろうか,とちょっと悩んだんですが,もう少し単純に表現できるのではないかと。
そこで考えたのが,ベルヌーイ分布の確率$$\pi$$に構造を入れること。イメージはIRTの1PLモデルで,要するに逆ロジット(ロジスティック)関数の中に困難度母数bを入れてS字カーブを右に左にずらす,あんな感じです。つまりある日を境に変化率$$\pi$$が1の方に上がっていくイメージ。

例えばXデーが300日目だったらこんなカーブになるわけです。

ということで,そんなモデルを書いてみました。ここで二階差分モデルを捨てたことに注意。