Kosugitti's BLOG

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

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

R

最強の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ファイルをこちらからダウンロードしてください。

[crayon-5ba403284d497260643255/]

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

[crayon-5ba403284d4ae918938025/]

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

[crayon-5ba403284d4b7436359870/]

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

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

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

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

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

Who are king of kings ?

事後分布から

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

[crayon-5ba403284d622571249393/]

強いのはブラマヨ。ね。

事後予測分布から

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

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

[crayon-5ba403284d7ae370163335/]

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文で分岐させるようにしています。

[crayon-5ba4032855413709778426/]

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

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

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

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

[crayon-5ba4032855421878165952/]

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

 

[crayon-5ba4032855428646668059/]

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

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

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

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

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

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

[crayon-5ba403285542e659884007/]

これでやってみた結果がこちら。

[crayon-5ba4032855434178425860/]

大変結構。tauの50%HDIが324から338というところ。どうやらデブモードに入ったのは10月20日〜11月3日の頃だ!ということがわかったわけです。

3.じゃあもうinv_logitの関数でbetaをかけてやったらいいんじゃね?

で,ここまでやって思いましたけど,これ普通に変化日までのS字カーブに増加率betaをかけるモデルにした方がシンプルなんじゃないかな?と思ってやって見ました。

[crayon-5ba403285543a187430003/]

こちらで分析したのがこちら。

[crayon-5ba403285543f539732346/]

ちょっとtauのRhatが悪いので,微調整が必要かもですが,やはりXデーは330日目の付近なわけです。うーん,そうかぁ。。。

ということがわかったからといって,なんの解決にもならないことに,数時間自分の体重を見つめ続けたあとで気づいたような次第です。

来年はデブモードを脱出するぞ!なんならヘルスサイドに堕ちるぞ!

みなさんMerry Christmas and Enjoy Stan!



RStudioの俺の好きな機能TOP5

このエントリーはRStudio Advent Calendar 2017の12月4日の記事です。

RStudioもアドベントカレンダーやってるんだー,なんか盛り上げるために気軽な記事を書こう,とおもってエントリしました。
RStudioはRを使うためだけでなく,サイトを作ったり報告書作ったり,色々な拡張機能があるのは確かなんですが,そこまで使い込んでない私でもRStudioの高機能にはいつも助けられています。
どれぐらい助けられているかというと,RStudioを知ってからはMplusをはじめとする他の統計ソフトを使わなくなっちゃった,ということでお察しいただけないでしょうか。素のRだと味気なくて,おなじ味気なさならMplusの方がいいし,学生に教えるときはRcommanderとかその他GUIな市販ソフトに劣るわけです。RStudioがあるから,これ一本で行けるんですよ。これはすごいことだと思う。

という前提を踏まえて,あー,これは助かるなぁと思う私の好きなRStudio機能,独断と偏見でTop5を記事にします。異論反論はもちろんあると思うので,そういう人は同じテーマで記事書いてね。

第5位 パッケージ作成補助

Rでパッケージ作る人,すでに中級者以上のイメージかもしれませんが,RStudioが来るまではやろうにもどうやって手をつけていいかわからなかったです。
ところがRStudioには新しいプロジェクトを作る,という時に「パッケージを作る」というのがもう入ってる。

ビルドボタンもついていて,ここでエラーチェックやら何やらしながら進められるので,これでグーーーーッと敷居が下がったと思うなあ。

第4位 Rpubsとの連携

これは私が大学で講義するような身分だから,というのもあるんですが,ソースコードや結果をそのままネットで公開できるのはとても助かります。

Rpubsというサービスがそもそも素敵!とおもいますけど,それと密接に関連してて,作ったものをすぐアップできるという気軽さ!
学生さんにはこのURLにアクセスしてくれれば,今日の授業のコードがコピペでできるよー,というようにして教えています。是非多くの講師先生に使ってもらいたい機能ですね。

第3位 Plotのサイズ設定

論文を書くようなお仕事をしていると,図版を論文に挿入しようと思うんだけど,そういう時ってサイズ指定ができるのはとても助かる。
PlotタブのExportをクリックするだけで,イメージでセーブ,PDFでセーブ,クリップボードにコピーができる。しかもサイズ指定,ファイル形式指定もできる。

これって相当ありがたい機能なんですよねー。素のRだとここまで簡単にできなかったから。この機能があるだけでも,学生さんなんかには「RStudio一択」と思わせる,超すごい機能じゃないでしょうか。

第2位 Stanコードの文法チェック

これはちょっとマニアな方向に舵を切りましたが,ベイジアンモデリングが趣味な人にとっては,この機能は本当に日々感謝しながら使ってるんじゃないかしら。

コンパイル前にチェックできるありがたさ。特にとある秘密結社では「セミコロンを忘れた」という理由で背中を蹴られたりしますから。教育の助けにもなっていいですね。
私はStanが主な用途ですが,そのほかにもC++などの言語系でチェックができるようです。ハイライトも地味にありがたいですし。

ちなみにベイジアンモデリングといえばStan,StanといえばそちらもAdvent Calendarをやっていますから,興味ある人はどうぞ。

第1位 プロジェクト管理

これは何位なんだろう,とずっと考えてました。ごく基本的な,当たり前の,この機能のためのRStudioとでもいうもので,特筆しなくていいじゃないと思うかもしれないけれども!逆に考えると,結局これが一番なんじゃないか。
フォルダ毎にプロジェクトを管理し,データやコードが飛散しないようになる。RStudioを使っている人でこの機能を使っていない人なんて,いないんじゃないかしら?ということで一位に据えました。

学生に指導するときも,まずプロジェクトを作れと言います。授業のたびに,授業のプロジェクトファイルを開け,と。これだけで,ファイルパスの問題がなくなるんです。こんなに嬉しいことはない。

自分で仕事をするときも,複数のプロジェクトを開いていたい時がある。そんなときはOpen Project in New Sessionですよ。

StanでMCMCサンプリングの待ち時間の間,学生のレポートをチェックするとか,研究資料の図表作成とか,ちょっとしたことをやることって多いですよね。そんなとき,複数のRを起動できている,これだけで効率がどれほど上がることか。
そういうことを考えると,一番いい機能はこれなんじゃないかな,と思っています。

ランク外は?

なんだか当たり前すぎて面白くない,という人もいるかもしれませんが,優れた技術は往々にして当たり前のようなものだとも思います。ご容赦ください。
自分なりに,エントリーしてから10日ほど熟考してつけた順位でしたが,ここから漏れた機能として次のようなものがありました。

  • 最近できたターミナルとの連携!これは嬉しい
  • Clear Consoleも結構使うのよね,スクショの時に綺麗にしたいから
  • 複数行をComment/Uncommentできるのもありがたいんです
  • 初学者にありがたいのはパッケージタブかな。インストール,アップデート,ライブラリ・・・
  • githubとの連携もありがたい。新しいブランチを切ることも最近できるようになったんだっけ。

というあたりで悩んでました。もちろんこれは私の環境下での話で,データベース使う人とか,サイト構築する人とか,スライド・文書を作る人からしたら「いやいやまだまだ」って思うところもあるでしょう。もう一度言いますが,そういう人は是非このテーマでブログ書いてください。俺もいろいろ知りたいんで。

Enjoy RStudio!


  • Categories:

状態空間モデルにはどれほどデータが必要か

この記事はStan Advent Calendar2017のエントリー記事です。

Stan盛り上がってますかー。盛り上がってますね(断定)。
私は本当に,Stanのおかげで生き返ったように楽しくモデリング出来ています。

先日の日本心理学会でもベイジアンモデリングのチュートリアルが大にぎわいでした。
それともう一つ人気だったな,と思うのが「時系列データ分析」とか「経験サンプリング」の話だったような印象です。
最近は加速度計とかモーションキャプチャ,アイトラックレコーダなどで,大量のデータが取れるようになってきているので,心理学者的には「わぁい隅々まで行動を計測できる!」というあたりが嬉しかったりするわけです。

でもまあ,大量に,時系列的に取れたデータはそれにあった分析が必要。特に

  • 自己相関が高い;一時点,あるいはそれ以上前の数値に関係のある現時点のデータ
  • 誤差まみれのデータ;無意味な情報も含んでいる大量のデータから意味のある部分だけ取り出す必要がある

という二つの問題から,「適切なモデリング」が必要になってきます。
そうすると,データ生成メカニズムを考えているベイジアンモデリングとの相性がいい。特に状態空間モデルは,ベイジアンアプローチをした方が,素直にモデリングできるのでお勧めです。

ということで,「いよいよ心理学にも空間統計とかはいってくるようになるかぁ」なんて思ったりしますけど,ちょっと待って。
心理学業界の一般的なデータは,上記のような器具を使うのでなければ,2,3時点の銃弾的データがほとんどで,1,000や10,000点以上あるでーたなんて(今のところ)あんまり見かけない。
どれぐらいデータがあればちゃんとした推定ができるんですかね?とちょっと気になりました。

と言うことで,一次元の状態空間モデル(時系列的なデータ)の仮想データを作り,そのデータサイズがどれぐらいだったら推定値として使えるのかのチェックをしてみようと思い立ちました。感度分析っていうのかな?こういうのも。

状態空間モデルについては過去記事もご参照いただくとして,味噌は目に見えない「状態」をモデリングし,それに観測誤差がついてデータになる,と考えること。すなわち,状態をS,観測値をXとすると,

$$S^{t+1}=S^t + \sigma$$
$$X^{t} = S^t + \phi$$

とします。あとは適切な事前分布をおいてやればオッケー。Stanのコードも簡単ですよ。

[crayon-5ba4032855cdf450772599/]

ほとんどイメージのまま書けますし,なめら化とかもしなくていいので楽です(わざとです。念のため。)。
さてさて,どれほどのデータがいるのかを検証するために,データサイズLをいろいろ変えたサンプルデータを作り,このコードで推定させて見ましょう。シミュレーションなので各試行を100回やってその平均値を考えることにします。

データセットを作り,MED推定値を返してくる関数を作りました。データサイズNと変動sig,測定誤差ERを渡すと

[crayon-5ba4032855cef851798526/]

これをコンパイルして走らせます。データサイズは指数的に増えるようにコード化しました。ここでは,Nは3,4,5,7,9,12,16,22,30,40,55,74,99,134,181,245,330,446,602,812,1097と増えて行きます。

[crayon-5ba4032855cf6517799853/]

これでシャラーンと走らせた結果がこんな感じ。

状態の変動分,sigの真値は3で,観測誤差ERの真値は5です。実際のスコアで見るとこんな感じになります。

P key mean sd
3 sig 4.468877 1.6121974
4 sig 4.052790 1.3807442
5 sig 3.839792 1.2260976
7 sig 3.789735 1.3442228
9 sig 3.543699 1.3567688
12 sig 3.566722 1.3576163
16 sig 3.345746 1.4561492
22 sig 3.323947 1.3837509
30 sig 3.113396 1.1905917
40 sig 3.142804 0.9679240
55 sig 3.101363 0.9104340
74 sig 3.176283 0.7182003
99 sig 3.017480 0.5714553
134 sig 2.997657 0.5463738
181 sig 2.974278 0.5059406
245 sig 2.983619 0.3664567
330 sig 3.006930 0.2936964
446 sig 2.997985 0.2609966
602 sig 3.069371 0.2744286
812 sig 3.001113 0.1897300
1097 sig 3.022984 0.1875880

だいたい3桁(99か134)あたりで真値に辿り着いている感じ。40,55点ぐらいだと真値+0.1ぐらいで,SDも1.0ぐらいあるので当たってるけどちょっと精度悪いと言う感じでしょうか。30以下はちょっと外れすぎかな。もちろんそれぞれの領域においてどの程度の精度が必要か,と言うのは変わってくると思いますけどね。ちなみに測定誤差の方についてもほぼ同様。

P key mean sd
3 ER 4.214914 1.6879550
4 ER 4.045415 1.6517713
5 ER 3.827050 1.3597881
7 ER 4.267874 1.5982695
9 ER 4.256110 1.3975502
12 ER 4.378329 1.2951766
16 ER 4.406097 1.2772903
22 ER 4.751362 1.0473943
30 ER 4.736674 0.9558007
40 ER 4.907256 0.8717817
55 ER 4.811194 0.7367819
74 ER 4.859863 0.5464713
99 ER 4.949720 0.5134745
134 ER 4.942098 0.4262254
181 ER 5.004040 0.4054680
245 ER 4.991647 0.3249265
330 ER 4.955008 0.2834097
446 ER 4.970939 0.2508226
602 ER 4.961162 0.1931686
812 ER 4.991295 0.1552512
1097 ER 4.980203 0.1451442

それぞれの変動分を大きくしたり小さくしたりして,試して見ることができます。

いや,実は自分の持っているデータでなんとかしようと思ったのですが,測定が35点ぐらいしかなくて,うまくいってなかったので「じゃあ何点ならできるんだよ!」と思って作ったコードでした。
100点を超えるデータって,例えばアンケート調査とかではほぼ無理だろうと思います。心理屋さんに状態空間モデルが流行するのがいつになるでしょうか。

最後に手前味噌ながら,宣伝を。こちらの本「Doing Bayesian Data Analysis」は翻訳・監修で関わらせていただきました。StanのメカニズムやMCMCのアルゴリズムについて,おそらく日本で一番簡単な図入りで説明している本です。お子様のクリスマスプレゼントにいかがでしょうか。

犬4匹本サポートサイトでは正誤表や関係イベントなどの情報を発信しております。

もしちょっと専門性の高いませたお子様,お友達がいると言うのであれば,こちらの通称「コワい本」を贈りましょう。これを枕の下に入れて眠ると,赤くて怖い人たちがパラメータを推定せよと迫ってくる夢が見られますよ。素敵!

ツイッタのハッシュタグ,「#犬4匹本」「#コワい本」でも情報発信しておりますのでぜひご活用ください。



Bayesplotパッケージはいいよ

ベイズの話があちこちで盛り上がっていますね!

ところで,ベイジアン・モデリングで得られた結果を表現する切り口っていろいろあるのですが,これを一手に引き受けてくれる便利なパッケージ,bayespotの存在を知りました。

これまではshinystanでみると綺麗!だったのですが,ブラウザが立ち上がるのでちょっと面倒。bayesplotパッケージはplot領域に出力されるし,結果がggplot2オブジェクトなので,ggplot2のレイアを上書きして言ったりすることができるのも便利。

ということで,練習がてら挙動を確認してみる。

まずは簡単なモデルを書いてみることに。

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

parameters{
  real mu;
  real beta;
  real<lower=0> sigma;
}

model{
  Y ~ normal(mu+beta*X,sigma);
}

generated quantities{
  vector[N] pred;
  for(n in 1:N)
    pred[n] = normal_rng(mu+beta*X[n],sigma);
}
# ライブラリの読み込み
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
library(bayesplot)

# 擬似データ作成
N <- 100
x <- runif(N,1,10)
yhat <- 7 + 15 * x 
y <- yhat + rnorm(N,0,10)

# モデルフィット
fit <- sampling(sampleCode,data=list(N=N,X=x,Y=y))
print(fit,pars=c("mu","beta","sigma"))
## Inference for Stan model: 3c6ac53b227702d2d44ca516d9974187.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##        mean se_mean   sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
## mu     6.82    0.07 2.67  1.54  5.02  6.87  8.59 12.05  1506    1
## beta  14.94    0.01 0.43 14.09 14.65 14.93 15.23 15.78  1510    1
## sigma 11.67    0.02 0.86 10.15 11.06 11.62 12.22 13.46  2246    1
## 
## Samples were drawn using NUTS(diag_e) at Wed Nov 29 12:19:32 2017.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

で,これをbayesplotパッケージに与えて色々するんだけど,stanfitオブジェクトのままではうまくいかないので,MCMCサンプルをarray 型にしておく。これでもって,サンプリングのいろいろな結果を描いてみる。

fit.array <- as.array(fit)

# MCMCのトレースプロット
mcmc_trace(fit.array,pars=c("mu","beta","sigma"))