ベイズマックスへようこそ。
このブログでは,心理学と統計が大好きなmutopsyが,
統計に関するあれやこれやを自由気ままに投稿します。

2018年07月01日

認知心理学への実践:データ生成メカニズムのベイズモデリング[スライド紹介とWS感想]

 こんにちは。mutopsyです。昨日(2018年6月30日),専修大学にて,広島ベイズ塾第三回ワークショップ「心理学者のためのベイズ統計学:モデリングの実際と,モデル選択・評価」が開催されました。本ワークショップですが,告知後ほんの数日で200名以上の予約参加者が集まりまして,当日のYouTubeのライヴ配信もアクティブ人数が常に40〜60人程度となるほどの盛り上がりを見せました。会場にお越しの皆様,配信を視聴してくださった皆様,そして一緒にワークショップの運営をしてくださった皆様,本当にありがとうございました。ベイズ塾生として発表した僕にとってもただただ有益で幸せな時間でした。(他の皆様の発表もめちゃくちゃ勉強になりました!)

 午後からは,心理学の分野ごとにベイズモデリングの実践例を紹介するセッションが行われましたが,そのセッションで僕は「認知心理学への実践:データ生成メカニズムのベイズモデリング」というタイトルの発表を行いました。発表で使用したスライドの一部をSlideShareに公開しています(↓)。

 認知心理学の考え方とベイズモデリングの相性の良さについて冒頭で説明した後で,認知心理学におけるモデリングの実践例を3つ紹介し,そこから分かるベイズモデリングの利点を列挙していくという流れで発表をしました。実践例として,(1) 自由再生実験で得られる系列位置曲線の生成メカニズムをモデル化した記憶のSIMPLEモデルの紹介,(2) 心的回転に関する既存のモデルをベイズ流に拡張した事例の紹介,それから(3)心理学の概念モデルを確率モデルに翻訳し,直接的なモデル検証を試みた事例の紹介を行いました。それぞれのモデルの詳細には立ち入りませんでしたが,「ベイズモデリングを使えばこんなこともできるんだ!」という実感を少しでもお伝えできたのであれば幸いです。(なお,実践例の(2)と(3)は僕自身の未発表のアイディアとデータを含む内容でしたので,上記のスライドでは非公開にしております。)

 この発表で最も伝えたかったのは,実践例の中身ではなく,冒頭で説明した認知心理学とベイズモデリングの付き合い方に関する哲学じみた部分です。僕が本格的にベイズの勉強を始めたのは1年と数か月前ですが,ベイズを学び始めた当初はしばらくの間,僕の研究分野(認知心理学)においてベイズを用いる利点がいまいち見い出せずにおりました。ここ数年,ベイズは心理学でもずいぶん流行っているのですが,情報がかなり錯綜している印象で,ある種の「神話」じみた言説を未だに多く見かけます(「ベイズ万能論」のような風潮?)。僕自身,そうした「神話」がきっかけでベイズの扉を開くことができた訳ですが,ベイズをきちんと学ぶにつれてそういった神話に疑問を抱くようになり,どんどん幻想が崩れていくうちに,結局ベイズの何が良いのかが分からなくなってしまいました。けれど,ベイズ塾や学会・研究会・SNSを通じて出会った多くの方々とともに勉強を続け,Stanで遊び続けてきた結果,ようやく1つの答えにたどり着くことができました。それが冒頭の内容であり,この発表の芯の部分です。つまり今回の僕の発表は,ベイジアンとしての僕の(現時点での)集大成でもある訳です。これまでの学びを自分なりに整理できたという意味でも,このワークショップは僕にとって本当に意味のあるものでした。

title.png
posted by mutopsy at 22:32 | 統計

2017年12月30日

Stanでポン・デ・リングを作ってみた

 こんにちは。mutopsyです。この記事は,Stan Advent Calendar 2017の25日目のエントリー記事です(実際にこの記事を書いたのは12月30日ですが)。

 みなさんはドーナツは何が好きですか?僕はフレンチクルーラーとエンゼルクリームがお気に入りですが,たまーにポン・デ・リングなんかも食べたくなります。このブログを訪れてくださった皆様の中には「Stanでポン・デ・リングが作れたら嬉しいのになあ……」とお考えになられた経験のある方もおられるでしょう。本日,なんとなくTogetterまとめを見ていたら,こんなまとめを見つけました。

「おっぱい関数選手権」の名古屋代表の名大生さん、1つの数式だけでポン・デ・リングを表示することに成功「天才の所業」

 このまとめは,ポン・デ・リングの形状(三次元と二次元の両方)がたった1つの数式で表現可能であることを発見した@CHARTMANqさんのツイートをまとめたものです。このまとめを読んで甚く感動した僕は,さっそくStanで実装してみることにしました。Stanを乱数生成器として用いて,二次元のポン・デ・リング分布からのMCMCサンプルを発生させるという試みです。三次元でもできますが,とりあえず二次元ポン・デ・リング不等式を用いてみました。

 二次元ポン・デ・リング不等式は,極座標を用いて(r-2|cos4θ|-10)(r+|cos4θ|-6)≤0と表せます(上記まとめの@CHARTMANqさんのツイートをご参照ください)。極座標というのは,原点からの距離 (r) と角度 (θ) を使って空間上(ここでは平面上)の点の位置を表す座標系です。この不等式を満たすrとθの領域が 二次元のポン・デ・リングを表します。

 二次元ポン・デ・リング分布からの乱数(MCMCサンプル)を発生させるStanコードの例は以下の通りです。

Stanコード01: 二次元のポン・デ・リング分布から乱数を発生させるコード
parameters{
  real<lower=0> r;  //極座標のr
  real<lower=0,upper=2*pi()> theta;  //極座標のθ(0°〜360°の範囲)
}

transformed parameters{
  real<upper=0> D;  //上限を指定することで不等式を表現
  D = (r - 2*fabs(cos(4*theta)) - 10)*(r + fabs(cos(4*theta))-6); //二次元ポン・デ・リング(不)等式
}

generated quantities{
  real x;
  real y;
  x = cos(theta) * r;  //直交座標のxを計算
  y = sin(theta) * r;  //直交座標のyを計算
}

 今回はStanを乱数生成のためだけに用いているのでdataブロックは必要ありません。また,事前分布として一様分布を用いるため,modelブロックも不要です。Stanに慣れていると変な感じがするかもしれません。rとθの事前分布(一様分布)のうち,ポン・デ・リング不等式を満たすパラメータのみをMCMCサンプルとして生成することにより,二次元ポン・デ・リング一様分布からの乱数を発生させます。このモデルを実行するためのRコードは以下の通りです。

Rコード01: 上記モデルを走らせるためのRコード
library(rstan)
data <- list(Dummy = 1) #ダミーデータを作る(実際には使わない)
stanmodel <- stan_model(file="stancode/Ponde_2d.stan") #上記Stanコードを読み込む
fit <- sampling(
  stanmodel,data=data,seed=123,
  chains=4,
  iter=10000
)

 上記のコードではダミーデータを渡していますが,実際にはこのデータに意味はありません。Stanコードを走らせてみたところ,Rhatは全て1.01未満となり問題なく収束しました。各パラメータ (r, θ, x, y) の事後分布は以下のようになります。

ponde_marginal.png

 周辺分布だけを見ても,なんとなくポン・デ・リングっぽくありませんか。θの分布の凹凸やxとyの二峰分布具合はまさにポン・デ・リングですよね。ではいよいよxとyの同時分布を見てみましょう。

pon.png

 どう見てもポン・デ・リングです。本当にありがとうございました。ぜひ皆さんもStanでポン・デ・リングを作ってみてください。もちろん数式が分かればポン・デ・リング以外のものも作れます。お試しあれ。それでは皆様,よいお年を!

 (「別にStanじゃなくてもよくね?」とつっこんではいけません。)

posted by mutopsy at 21:15 | 雑記

2017年12月14日

実験心理学者のためのベイジアンモデリング──正答率データをどうやって分析するか──

 こんにちは。mutopsyです。この記事は,Stan Advent Calendar 2017の14日目のエントリー記事です。この記事では,実験心理学でよく扱われる正答率・誤答率(あるいは正答数・誤答数)データの分析を例に,よくある素朴な分析(平均値の点推定・区間推定およびt検定)の背後にあるモデルを可視化した上で,そのモデルを改善した別の3つのモデルを紹介します。また,情報量規準によるモデル比較も行います。そして最後に,実験心理学におけるベイジアンモデリングの有用性について述べたいと思います。

1. はじめに:正答率データを素朴に分析してみる

 改めまして,こんにちは。mutopsyこと武藤拓之です。大阪大学で認知心理学の研究をしています。この記事では,実験心理学でもそれ以外の分野でもお馴染みの,正答率(正答数)データを分析する方法について考えてみたいと思います。昨日はフォン・ミーゼス分布に関する若干マニアックな記事を投稿しましたが,この記事ではもう少しだけ需要のありそうな話をするつもりです。

 さて,いきなりですが,下のグラフ(図1)をご覧ください。このグラフは,筆者が過去に行った記憶に関する実験のデータを要約したものです1。この実験では,実験参加者さんに2枚の画像を順番に見てもらい,2枚目の画像に写っている対象が1枚目の画像と同じであったか異なっていたかを回答してもらいました。この試行を繰り返すと,実験参加者さん一人ひとりの正答数のデータを得ることができます。正答数を合計試行数で割れば正答率になるので,正答率のデータが実験参加者の人数分だけ手に入ると考えることもできます。図1のグラフの青い点は,実験参加者一人ひとりの正答率をプロットしたものです。この実験では難易度の異なる2つの条件が設定されており,条件1と条件2にはそれぞれ異なる実験参加者が割り当てられたと考えてください。なお,この記事の内容は,正答数と合計試行数が分かっているデータであればどんなデータにも当てはまる話なので,読者の皆様にとって分かりやすいデータをイメージして頂ければ良いかと思います。

ac_rawdata.png
図1. この記事で分析に用いる正答率データの平均値とその95%信頼区間。青い点は個人ごとの正答率。

 このデータのcsvファイルはこちらからダウンロードできます。→ accuracy_data.csv

 このデータは,各行が実験参加者1人分のデータに対応しており,以下の5つの変数を含んでいます。

  • ID: 実験参加者のID(1から24までの整数)
  • N.Corr: 正解した試行数(1から6までの整数)
  • N.Total: 合計試行数(6で固定)
  • Corr.Rate: 正答率(「正解した試行数」を「合計試行数」で割った値)
  • Group: 条件群(1または2)

 実験心理学の文脈では,このようなデータが得られたら,まずは平均値の点推定値と誤差範囲を計算するのが一般的だと思います。ここでは各条件における平均正答率に関心があるので,正答率の算術平均(=母平均の点推定値)と母平均の95%信頼区間を計算してみました。これらの値は先ほどの図1に折れ線グラフとエラーバーで示されています。赤色の四角形が各条件の平均正答率,エラーバーがその95%信頼区間です。信頼区間ではなく標準誤差をプロットする場合も多いですが,このようにして算術平均とエラーバーをグラフに描くのが一般的かと思います。データの傾向を確認したら,例えば2つの条件間で平均正答率に差があるか否かを検定したり(Welchのt検定),あるいはある条件の平均正答率がチャンスレベル(ここでは50%)よりも有意に高いと言えるかどうかを検定したり(一標本のt検定),研究仮説に合った分析を行います。実際にこのデータに対してWelchのt検定を行うと,p = .024 という結果が得られ,条件2の平均正答率は条件1の平均正答率よりも平均して23.6%高い,という結論が得られます。

 ただ,改めて図1のグラフを見ると,気になる点がいくつかあります。まず,条件2の平均正答率の95%信頼区間をよく見ると,区間の上側が1.0 (=100%) の値よりも上にあります。正答率は理論上0%から100%の範囲内に収まるはずなのに,平均正答率の母数が100%を超過しているかもしれないという奇妙な推定結果が得られてしまいました。また,実験参加者一人ひとりの正答率の散らばり方も気になります。例えば,条件1でも条件2でも正答率が100%の実験参加者がそれなりにいるように見えます。少なくとも正規分布に従っているようには見えません。たとえデータそのものが正規分布に従っていなくても,サンプルサイズが十分に大きければ平均値の分布は中心極限定理により正規分布に近付きますが,今回のように高々 n = 12 の場合には微妙な気もします。それに,やはり平均値の信頼区間が100%をまたいでしまうのは何となく気持ち悪いです。

 機械的に信頼区間を求めたりt検定を行ったり,あるいは分散分析したり回帰分析したり,ということを行ってしまうと,これらの分析が結局何をやっているのかが見えにくくなってしまいます。そこでまず,上述した分析の背後にあるモデルを可視化してみようと思います。そうすれば,先の分析の違和感の正体が掴めるはずです。

2. モデル1:制約のない正規分布モデル

 この節では,1節で説明したような素朴な分析が何をやっているのかを可視化するために,先の分析の背後にあるモデルをベイズ的に表現してみようと思います。具体的には平均値の区間推定の背後にあるモデルをStanで記述してみます(平均値差の推定については6節で触れます)。「頻度主義 vs. ベイジアン」という枠組みで考えてしまうと違和感を持たれるかもしれませんが,以下のモデルはいわゆる頻度主義の枠組みで想定されているモデルと(ほぼ)等価なモデルであり,後で確認するように,推定結果も(ほぼ)一致します。なお,この記事ではパラメータの解釈に関する頻度主義とベイズの違いについては触れません。前節で説明したようなよく使われる分析方法の背後にあるモデルは,以下に示す「モデル1:制約のない正規分布モデル」として記述することができます2。特に重要な行を赤色でハイライトしています。

Stanコード01: 【モデル1】制約のない正規分布モデル
data {
  int N;  //データの総数
  real CorrRate[N];  //正答率
  int Group[N];  //実験条件(群)
}

parameters {
  real mu[2];  //平均正答率 (μ)
  real<lower=0> sigma[2];  //標準偏差 (σ)
}

model {
  //モデル
  for (n in 1:N){
    CorrRate[n] ~ normal(mu[Group[n]],sigma[Group[n]]);
    //↑正答率データが正規分布に従う
  }
  
  //事前分布
  for (k in 1:2){
    mu[k] ~ uniform(-9999999,9999999); 
    sigma[k] ~ uniform(0,9999999);
  }
}

generated quantities{
  real diff_mu;
  diff_mu = mu[2] - mu[1];  //2群の平均正答率の差
}

 このモデルでは,条件1と条件2の正答率データがそれぞれ独立な正規分布から生成されると仮定しています。正規分布は平均値μと標準偏差σという2つのパラメータを持ちますが,今回の分析ではこの平均値μの分布が求めたい平均正答率の分布に相当します。以下,コードの解説です。

  • 1-5行目: dataブロックでは分析に用いるデータを宣言します。ここではデータの総数(今回はデータが24行なのでN=24),正答率,実験条件(群)を使用します。
  • 7-10行目: parametersプロックでは,このモデルで推定するパラメータを宣言します。ここでは正規分布のパラメータμとσを条件ごとに推定します。σは非負の値なので下限を0に指定しておきます(これを指定しないとうまく推定できません)。
  • 15行目: 正規分布から一人ひとりの正答率データが生成されることを表現しています。この部分がこのモデルの核となります。
  • 21-22行目: 事前分布を指定しています。この部分を省略した場合,Stanは十分に幅の広い一様分布を事前分布として採用します。なので本当はこの2行は不要なのですが,モデルの仮定を示すためにあえて明示的に書いています。伝統的な分析による推定結果は,-∞から+∞の範囲の一様分布3を事前分布としてベイズ推定した時の結果と一致します。ここで指定しているuniform(-9999999,9999999)は,-9999999から+9999999までの範囲の一様分布を表します。つまり,平均値がいくつになるか全く分からないという事前の信念を表現しています。本当はuniform(-Inf,Inf)とでも書きたいところですが,Stanではこのような書き方は(おそらく)できないので,あえて極端に範囲の広い事前分布を指定してみました。
  • 26-29行目: 条件間の平均値差を計算しています。単にmu[0]とmu[1]の引き算をしているだけです。

 このモデルに先ほどのデータを当てはめて,推定したμの事後モード(MAP推定値)と95%最高密度区間を求めると,以下の図2のような結果が得られます。個人のデータのプロットはローデータそのものなので,値は図1と厳密に同じです(描画の都合上,横位置は若干ずれています)。折れ線グラフで表された平均値(μ)の点推定値(MAP推定値)とその95%最高密度区間は上記のコードで推定したものです。推定上の誤差が多少あるとはいえ,図1とほとんど同じ結果になっていることが分かるかと思います。

ac_model1.png
図2. 制約のない正規分布モデル(モデル1)で推定した平均正答率。エラーバーは95%最高密度区間を表す。

 改めてこのモデルを見てみると,まず事前分布の箇所が気になるかと思います。正答率の平均値は理論上0から1の範囲に収まるはずなのに,事前知識としてほぼ-∞から+∞の範囲の一様分布を用いるのは明らかに不自然です。素朴な信頼区間の計算やt検定の背後にはこのようなモデルが想定されているのです。そこで次のモデルでは,「平均正答率(μ)は0から1までの値しか取らない」という事前知識をパラメータの事前分布として反映させてみようと思います。

3. モデル2:μの範囲を制限した正規分布モデル

 モデル2は,正規分布の平均値パラメータμについて「0から1までの値しか取らない」という事前知識を与えたモデルとなります。先ほどのモデル1を僅かに修正するだけでこのモデルを作ることができます。以下のコードの赤色でハイライトされた部分のみがモデル1と異なる部分です。

Stanコード02: 【モデル2】μの範囲を制限した正規分布モデル
data {
  int N;  //データの総数
  real CorrRate[N];  //正答率
  int Group[N];  //実験条件(群)
}

parameters {
  real<lower=0, upper=1> mu[2];  //平均正答率 (μ)
  real<lower=0> sigma[2];  //標準偏差 (σ)
}

model {
  //モデル
  for (n in 1:N){
    CorrRate[n] ~ normal(mu[Group[n]],sigma[Group[n]]);
    //↑正答率データが正規分布に従う
  }
  
  //事前分布
  for (k in 1:2){
    mu[k] ~ uniform(0,1); 
    sigma[k] ~ uniform(0,9999999);
  }
}

generated quantities{
  real diff_mu;
  diff_mu = mu[2] - mu[1];  //2群の平均正答率の差
  }
}

  • 8行目: parametersブロックで平均値パラメータを宣言する時に,あらかじめ下限0と上限1と指定しておきます。
  • 21行目: 平均値パラメータの事前分布として,0から1の範囲の一様分布を指定しています。parametersブロックで平均値パラメータの下限と上限を指定しているので,この行を省略しても自動的に範囲[0,1]の一様分布が事前分布として指定されます。しかし,この上限と下限の制約がパラメータの事前分布であるということを強調するために,あえてここでは明示的に書いています。

 このモデルによる推定結果は以下の図3のようになります。目論見通り,条件2の区間推定幅の上端が100%を越えなくなりました。なお,この記事の下部にある図6にこの記事で用いた全てのモデルの結果を示しているので,比較の際にはこちらをご覧いただくと分かりやすいかと思います。

ac_model2.png
図3. μの範囲を限定した正規分布モデル(モデル2)で推定した平均正答率。エラーバーは95%最高密度区間を表す。

 さて,平均値の範囲の問題は一応解決しましたが,そもそも正規分布を用いたモデルで本当に良いのだろうかという疑問は残ります。次の節では,正規分布ではなく二項分布を用いたモデルについて説明します。

4. モデル3:θを固定した二項分布モデル

 表が出る確率がθであるコインをn回投げた時に表が出る回数kは,pとnをパラメータに持つ二項分布で表現できます。もちろんコイン投げ以外のケースにも二項分布は使えます。例えば,記憶実験のある試行で正答する確率をθ,合計試行数をnと置けば,正答する回数kはやはり二項分布に従います(二項分布を使ったモベイズ推定については過去の記事でも紹介しています)。二項分布はnが大きくなるにつれて正規分布に近づきますが,今回分析に使用しているデータの合計試行数はたったの6なので,正規分布ではなく二項分布を使ったモデルのほうが適切であるように思われます(正規分布を用いるということは,二項分布モデルにおいてnが∞であると仮定することと一致します)。また,kとnは正の整数で,k > nとなることはなく,確率θは0から1の範囲に収まります。つまり,二項分布を用いる以上はモデル2のようにわざわざ範囲を制限する必要がありません(モデル上では制限をかけないと収束しませんが,この制限は事前の信念に基づく制限ではなく,モデルから自然に要求される制限であるという点で異なります)。

 それでは二項分布を使ったモデルについて考えてみましょう。上述したように,ある試行で正答する確率θ,合計試行数nの二項分布から正答回数が生じるモデルを想定しますが,まずは,θが条件のみによって決定するというモデルを考えます。つまり,条件1における正答確率をθ1と置くと,条件1の実験参加者は,どの参加者であっても任意の試行においてθ1の確率で正答し,1-θ1の確率で誤答すると考えます。このパラメータθは平均正答率ではなく各試行で共通の正答確率を表すという点でモデル1とモデル2のμとは厳密には異なります。ただし,このモデルでは同じ条件の参加者でθが共通しているという仮定を置いているため,各条件における正答確率θの参加者間平均は正答確率θと一致します(θ * 12 / 12 = θ)。そのため,ここではθをμと同じように解釈することにします(研究仮説において両者を区別する必要が特になければ問題ありません)。このモデルのStanコードは以下のようになります。推定するパラメータがmuからthetaに変わっていることと,データが正規分布ではなく二項分布に従うと仮定していることが先のモデルとの違いです。

Stanコード03: 【モデル3】θを固定した二項分布モデル
data {
  int N;  //データの総数
  int CorrResp[N];  //正答数
  int TotalResp[N]; //合計試行数
  int Group[N]; //実験条件(群)
}

parameters {
  real<lower=0, upper=1> theta[2]; //正答確率 (θ)
}

model {
  //モデル
  for (n in 1:N){
    CorrResp[n] ~ binomial(TotalResp[n],theta[Group[n]]);
    //↑正答数データが二項分布に従う
  }
  
  //事前分布
  for (k in 1:2){
    theta[k] ~ uniform(0,1);
  }
}

generated quantities{
  real diff_theta;
  diff_theta = theta[2] - theta[1];  //2群の正答確率の差
}

  • 3-4行目: 今度は正答率ではなく正答数と合計試行数をデータとして与えます。今回のデータでは合計試行数は6で固定ですが,この値が実験参加者によって違っていてもこのモデルは走ります。
  • 9行目: 正規分布の時とは推定するパラメータが違うので,名前もmuではなくthetaにしています。thetaの範囲は0から1に指定しておかないとうまく収束しません。
  • 15行目: 二項分布から正答数データが生成されるモデルを表現しています。

 このモデルを当てはめると以下の図4のような結果が得られます(他のモデルとの比較は図6参照)。縦軸が平均正答率(μ)ではなく正答確率(θ)になっていることに注意してください。パラメータのMAP推定値はモデル1・モデル2とあまり変わっていませんが,これまでのモデルよりも条件1における区間推定幅がかなり狭くなっていることが読み取れます。また,当然ながら区間推定幅は100%をまたぎません。

ac_model3.png
図4. θを固定した二項分布モデル(モデル3)で推定した正答確率。エラーバーは95%最高密度区間を表す。

 この節では二項分布を使ったモデルを紹介しましたが,「θが実験参加者にかかわらず一定である」という仮定はなんとなく直感に反するように思えます。記憶成績が良い人も悪い人もいる,というように,θにも個人差があると考える方が自然でしょう。そこで,次の節ではθが人によって異なるモデルの一例を示します。

5. モデル4:θが個人ごとに異なる二項分布モデル

 この節では,二項分布のパラメータθが人によって異なることを仮定したモデルの一例を紹介します。以下にそのStanコードを示します。赤色でハイライトされた行はモデル3と異なる箇所です。

Stanコード04: 【モデル4】θが個人ごとに異なる二項分布モデル
data {
  int N;  //データの総数
  int CorrResp[N];  //正答数
  int TotalResp[N]; //合計試行数
  int Group[N]; //実験条件(群)
}

parameters {
  real<lower=0, upper=1> theta[N]; //正答確率 (θ)
}

model {
  //モデル
  for (n in 1:N){
    CorrResp[n] ~ binomial(TotalResp[n],theta[n]);
    //↑正答数データが二項分布に従う
  }
  
  //事前分布
  for (n in 1:N){
    theta[n] ~ uniform(0,1);
  }
}

generated quantities{
  real diff_mu;
  real mu[2];
  mu[1] = 0;
  mu[2] = 0;
  for (n in 1:N) mu[Group[n]] = mu[Group[n]] + theta[n]/12;
  diff_mu = mu[2] - mu[1];
}

  • 9行目: 前節のモデル3ではパラメータthetaは条件によって異なると仮定していため,2つのthetaを推定すれば良かったのですが,今回のモデルでは実験参加者の人数分だけthetaを推定します。
  • 15行目: 条件ごとのthetaではなく実験参加者ごとのthetaを二項分布のパラメータに指定します。
  • 20-21行目: 事前分布も全てのthetaに対して指定します。
  • 27-30行目: 推定したthetaの平均値muを条件ごとに計算しています。このmuは正規分布のパラメータμではなく,ただの算術平均です(thetaが分布するのでmuも分布します)。

 このモデルでパラメータを推定した結果を以下の図5に示します(モデルの比較は図6を参照)。縦軸は正答確率 (θ) の参加者間平均です。モデル3に比べると,条件2のパラメータの推定値が低くなっています。

ac_model4.png
図5. θが個人ごとに異なる二項分布モデル(モデル4)で推定した正答確率の平均。エラーバーは95%最高密度区間を表す。

6. 4つのモデルを比較する

 さて,これまで4つのモデルを使って正答率データのモデリングを行ってきましたが,結局どのモデルを使うのが良いのでしょうか。図6は4つのモデルによる推定結果を並べたものですが,モデルによって結果が異なっている様子が見て取れます。

modelcomp.png
図6. 4つのモデルによる正答率(正答確率)の推定結果の比較。

 ここではモデル比較のために,情報量規準のWAICを用いることにします。WAICは予測の良さを表す指標で,値が小さいモデルほど真のモデルに近いことを意味します4。それぞれのモデルのWAICは図6に示した通りです5。この結果からは,モデルの良さは「モデル1 < モデル2 < モデル3 < モデル4」であると言えそうです。やはり,素朴な分析の背後にあるモデルと等価な正規分布モデルはあまり適切ではないようです。この中ではθが人によって異なる二項分布モデルが最も良いモデルであるといえそうです。

 これまでは触れてきませんでしたが,各モデルにおける条件間の(平均)正答率の差(diff_muまたはdiff_theta)の推定結果も見てみましょう。以下の図7に結果を示しました。これを見ると,どのモデルにおいても(平均)正答率の差の区間推定幅には0が含まれていませんが,幅の広さにはだいぶ違いが見られます。モデル1の95%最高密度区間はWelchのt検定の計算過程で用いられる95%信頼区間と一致します。どうやら二項分布モデルのほうが推定精度が良さそうです。

ac_diff.png
図7. 4つのモデルによる正答率の条件間差の推定結果。

 ただし,可能なモデルは他にもあるので(例えばベータ二項分布モデルのような,θに分布を仮定するモデルや,感度と反応バイアスを区別する信号検出理論モデル,項目困難度と能力を区別する項目反応理論モデルなど),今回のモデル4があらゆるモデルの中で最良であるとは限りません。研究仮説に応じて適切なモデルを選択することも必要でしょう。

7. まとめ

 この記事では,素朴な分析(平均値の点推定・区間推定およびt検定)の背後にあるモデルを可視化した上で,そのモデルを修正した別の3つのモデルを紹介しました。また,モデルごとの推定結果の違いを概観し,WAICを用いたモデル比較も行いました。

 この記事で最も主張したかったことは,分析の背後にあるモデルを意識することの重要性についてです。機械的に分析を実行することに慣れてしまうと,明らかにデータに合わないモデルを無意識に用いてしまう危険性があります。Stanを使ってベイズ的にモデルを記述すれば,そのモデルが何を仮定しているのかは簡単に分かりますし,この記事で実際にお見せしたように,モデルの修正も容易に行えます。この記事で紹介したモデル1-4のパラメータは,ベイズを用いないで最尤法で解くことも可能だとは思いますが,ベイズ推定で用いられるMCMCアルゴリズムは強力かつ汎用的で,またStanを使えばモデルの記述も修正も拡張も比較的簡単に行えるという点で,非常に一般性が高いです。共通の文法を用いてあらゆるモデルを記述しパラメータを効率的に推定できるベイジアンモデリングはとても魅力的であるように思います。ベイズというと,パラメータの解釈における頻度主義との違いが強調されることが多いですが,個人的にはモデリングの柔軟性こそがベイズの良さであるように感じています。あなたもベイジアンモデラーになってみませんか。

脚注

1本当はもっと多くの条件を設定していたのですが,ここでは2つの条件のデータのみを抽出しています。また,実際にはこの2条件は対応のあるデータなのですが,分かりやすさのために,ここでは対応のないデータとして扱うことにします。とはいえこの記事で書かれている内容は,もっと条件数の多い実験で得られたデータや対応のあるデータに対してももちろん適用可能です(Stanコードは適宜修正しなければなりませんが)。

2分かりやすさを優先し,各実験参加者の正答率が平均μ(=平均正答確率)の正規分布に従うモデル,すなわち平均正答確率を直接推定するモデルを書きました。しかし,正答率ではなく正答数が正規分布に従うモデルを書き,generated quantitiesブロックで平均正答数を合計試行数で割って平均正答率を計算するほうが,モデルの適用範囲は広がります。また,渡すデータが後の二項分布モデル(モデル3・モデル4)と同じになるのでそのままWAICを計算できるようになります。

3厳密に言えば,この関数は全区間で積分しても1にならないので確率密度関数ではない。

4情報量規準は万能ではなく,あくまでも1つの判断基準でしかないのでご注意ください。

5モデル比較のためにWAICを計算する際には,どのモデルにも同じデータを与える必要があります。上述のモデル1とモデル2では正答率データを渡していたのに対し,モデル3とモデル4では正答数と合計試行数のデータを渡していたので,そのままWAICを計算すると,モデル1とモデル2のWAICが極端に小さくなってしまいます。そのため,モデル1とモデル2のWAICを計算する際は,正答確率ではなく正答数が正規分布に従うモデルに修正しなければなりません。この時,修正したモデルには合計試行数のデータを渡す必要はありません(正規分布モデルは合計試行数が∞の二項分布モデルと等しいため)。このことは関西学院大学の清水先生に教えて頂きました!

2018年1月25日:csvデータへのリンク切れを修正しました。

posted by mutopsy at 01:45 | 統計

2017年12月13日

方向音痴の人は正しく「北」を指し示せるか──フォン・ミーゼス分布を使って角度データを分析する──

 こんにちは。mutopsyです。この記事は,Stan Advent Calendar 2017の13日目のエントリー記事です。この記事では,角度や色相,あるいは1日の時間のような,循環する変数のモデリングに使えるフォン・ミーゼス分布を用いた分析について,実験心理学っぽい例を挙げて説明していきます。実際に分析するためのRとStanのコードも紹介します。

1. はじめに

 という訳で,アドカレデビューしてみました。改めまして,mutopsyこと武藤拓之です。大阪大学で認知心理学の研究をしています。いろんなものに興味がありますが,主にはヒトの空間認知過程について実験心理学の手法を用いて研究しています。空間認知の研究では変数として角度(30度とか270度とか)を用いることがよくあります。角度は循環する変数なので,360度回転すると一周して元の角度に戻ってきます。このような角度データに正規分布を当てはめてしまうと,循環する性質をうまく表現できず,例えば「30度付近のデータが得られる確率は390度付近のデータが得られる確率よりも高い」といった奇妙なことが起こってしまいます。そこでフォン・ミーゼス分布の出番という訳です。

2. 架空の実験データの例

 分布の説明をする前に,まずはデータの具体例を挙げたいと思います。以下の図1は架空の実験場面を描いたものです。実験室を訪れた実験参加者さんに,「北の方角を指差してください」と教示します。この実験室がある学部棟は独特な形をしており,実験参加者の正面(0度)から時計回りに120度の方向が北になっています。つまり,120度付近を正しく指差せたら「正解」となります。実験参加者が実際に指を差した方向(正面を0度とした時計回りの角度)をデータとして記録します。

vm_experiment.png
図1. 実験場面の例。実験室を訪れた実験参加者は,「北」を指差すように求められる。

 この実験を200人の実験参加者に対して行ったと想定します。そのうちの100人は方向感覚が良い人たちで,残りの100人は方向感覚が悪い人たち,すなわち方向音痴な人たちであったことが事前に分かっていたとします。このようにして得られた200人分のデータをプロットしたら,以下のような結果が得られました。なお,このデータはこの分析のために作った架空データです。このデータを生成するRコードは後述します。

vm_dataplot.png
図2. 指差し実験から得られた架空の角度データの分布。左右のグラフは同じデータを異なるグラフで表している。

 データを見てみると,方向感覚が良い人たちのグループは,多少の誤差はあれども,ほぼ正確に120度の方向を指し示せていることが読み取れます。一方,方向感覚が悪い人たちのグループの回答はかなりバラついています。バラついてはいるものの,均等に散らばっているというよりは,0度(=正面)付近を指す人の割合が若干高そうです。

3. フォン・ミーゼス分布について

 循環する変数が従う分布としては,フォン・ミーゼス分布がよく用いられます。フォン・ミーゼス分布は単峰かつ対称な分布で,循環正規分布とも呼ばれます。平均値パラメータμと集中度パラメータκ(カッパ) の2つのパラメータで形が決まります。集中度パラメータのκは非負の実数で,値が大きいほどμの付近に値が集中し,かつ正規分布に近づきます。κが0の時には一様分布となります。また,κは第1種変形ベッセル関数を用いることで,円周分散Vと円周標準偏差ν(ニュー)に変換できます。この円周標準偏差νは,正規分布の標準偏差σと同様の解釈ができます (豊田, 2017)。以下の図3は,μとκを変えた時にフォン・ミーゼス分布がどのような形状になるのかを示しています。ここでは0度と360度を両端とする確率密度関数を描画していますが,実際には横軸が循環していることに注意してグラフを見る必要があります。例えば中央上のグラフと右上のグラフは一見すると二峰分布に見えてしまうかもしれませんが,いずれも値がμの時にピークとなる単峰分布です。

vms.png
図3. いろいろなフォン・ミーゼス分布。

4. シミュレーション用のデータセットを作る

 以下に,フォン・ミーゼス分布から乱数を発生させて架空データを作るRコードを示します。2節で紹介した架空データはこのコードを用いて作成しました。循環する変数を扱うにはcircularパッケージが便利です。このパッケージのrvonmises関数を使うことでフォン・ミーゼス分布からの乱数を発生させることができます。また,ここでは詳しく説明しませんが,plot.circularという関数を使えば循環する変数を円周上にプロットすることができます。図2の左のグラフはこの関数を使って作成しました(ちなみに右のグラフは皆大好きggplot2で作成しました)。

Rコード01: フォン・ミーゼス分布に従う架空データを2セット作る
##前準備
library(circular)
rad_trns <- pi / 180 #これを度数に掛けるとラジアンに変換できる

##架空データのパラメータを設定
n1 <- 100    #群1のサンプルサイズ
n2 <- 100    #群2のサンプルサイズ

m1 <- 120 * rad_trns  #群1の平均値
k1 <- 50  #群1の集中度

m2 <- 0 * rad_trns  #群2の平均値
k2 <- 0.5  #群2の集中度

##円周分散と円周標準偏差を計算する
A1 <- besselI(k1, 1, expon.scaled=F)/besselI(k1, 0, expon.scaled=F)
v1 <- 1 - A1  #群1の円周分散
s1 <- sqrt(-2 * log(A1))  #群1の円周標準偏差

A2 <- besselI(k2, 1, expon.scaled=F)/besselI(k2, 0, expon.scaled=F)
v2 <- 1 - A2  #群2の円周分散
s2 <- sqrt(-2 * log(A2))  #群2の円周標準偏差

##乱数を生成
d1 <- as.numeric(rvonmises(n1,circular(m1),k1))  #群1のデータ
d2 <- as.numeric(rvonmises(n2,circular(m2),k2))  #群2のデータ
  • 前準備(1-3行目):circularパッケージとrstanパッケージを読み込み,rad_trnsという変数を作る。この変数rad_trnsを度数(e.g., 90度)に掛けるとラジアン(e.g., π/2)になり,この変数でラジアンを割ると度数になる。角度データを扱う時には作っておくと便利。ちなみにcircularパッケージには度数をラジアンに変換する関数radがあるのでこれを使っても良い。
  • 架空データのパラメータを設定(5-13行目):群1と群2のパラメータの真値とサンプルサイズを設定する。ここの値をいろいろ変えてみながら遊ぶと楽しい。
  • 円周分散と円周標準偏差を計算する(15-22行目):上で設定したパラメータを変形して,各群の真の円周分散と円周標準偏差を得る。besselI(x, n, expon.scaled=F) はn次の第1種変形ベッセル関数。
  • 乱数を生成(24-26行目):平均m1 or m2,集中度k1 or k2のフォン・ミーゼス分布に従う乱数をn1 or n2個生成する。rvonmises関数で引数のfromとtoを指定すれば,生成する乱数の範囲を指定できる (例えばfrom = -pi, to = pi; デフォルトではfrom = 0, to = 2 * pi)。また,デフォルトでは弧度法(ラジアン)のデータを渡すようになっているが,control.circular = list(units = "degrees")と指定すれば度数法のデータを渡すこともできる。詳しくはヘルプを参照。

5. Stanで分析を実行

 それでは,先ほどR上で生成した架空データにフォン・ミーゼス分布を当てはめて,パラメータを推定してみましょう。さらに,(1) 方向感覚が良い人たちと悪い人たちとの平均値差(最短角距離)はどのくらいか,(2) 方向感覚の悪い人たちの円周標準偏差は方向感覚が良い人たち円周標準偏差の何倍か,という2つの問いに答えるための生成量も推定してみましょう。以下にそのためのStanコードを示します。

Stanコード01: フォン・ミーゼス分布のパラメータを推定する
data {
  int N1;  //群1(高群)のサンプルサイズ
  int N2;  //群2(低群)のサンプルサイズ
  real Angle1[N1];  //群1の角度データ(ラジアン)
  real Angle2[N2];  //群2の角度データ(ラジアン)
}

parameters {
   real<lower=-pi(), upper=pi()> mu1;  //群1の平均
   real<lower=0, upper=2*pi()> mu2;  //群2の平均
   real<lower=0, upper=200> kappa[2];  //群1と群2のそれぞれの集中度
}

model {
  for (n1 in 1:N1){
    Angle1[n1] ~ von_mises(mu1, kappa[1]);  //群1のモデル
  }
  for (n2 in 1:N2){
    Angle2[n2] ~ von_mises(mu2, kappa[2]);  //群2のモデル
  }
}

generated quantities{
  real mu1_trns;  //mu1(ラジアン)を度数に変換
  real mu2_trns;  //mu2(ラジアン)を度数に変換
  real a[2];  //平均合成ベクトル長(円周分散・円周標準偏差の計算に使用)
  real v[2];  //群1と群2の円周分散
  real nu[2];  //群1と群2の円周標準偏差
  real nu_trns[2];  //nu(ラジアン)を度数に変換
  real mu_trns_diff;  //mu1_trnsとmu2_trnsの差(最短距離)
  real ratio_nu; //nu1とnu2の比

  #上で宣言した生成量を計算
  for (k in 1:2){
    a[k] = modified_bessel_first_kind(1,kappa[k])/modified_bessel_first_kind(0,kappa[k]);
    v[k] = 1 - a[k];
    nu[k] = sqrt(-2 * log(a[k]));
    nu_trns[k] = nu[k] / pi() * 180; 
  }
  mu1_trns = mu1 / pi() * 180;
  mu2_trns = mu2 / pi() * 180;
  mu_trns_diff = mu2_trns - mu1_trns;
  ratio_nu = nu[2] / nu[1];
}

  • dataブロック(1-6行目):Rから渡すデータを記します。
  • parametersブロック(8-12行目):推定するパラメータを宣言します。平均パラメータmu1とmu2の範囲については,レンジが2π (=360度) になるように上限と下限を指定しておきます。また,上限と下限の境界周辺にパラメータが多く存在する場合には収束効率が著しく低下し推定がうまくいかなくなります。これは,境界部分が計算上非連続な点になってしまうことに起因すると思われます。これを防ぐために,あらかじめデータの分布を見て平均値の目星を付け,その平均値がパラメータの範囲の中心に位置するように上限と下限を設定するのが良いでしょう。例えば,群2の平均値の真値は0度ですが,mu2の範囲を<lower = 0, upper = 2 * pi()>としてしまうと,パラメータ空間が「0度〜」の部分と「〜360度(= 2π)」の部分の二か所に分断されてしまい,0度(= 360度)に近い値のパラメータをうまくサンプリングできなくなってしまいます。また,この状態で得られたMCMCサンプルは左右両極にピークを持つ二峰分布となるため,そのまま平均値を計算すると,EAP(mu2) ≒ π(= 180度)となり,真値(0度)と乖離してしまいます(こちらに関しては,事後分布の要約としてMAP推定値と最高密度区間を用いることで対処できます)。こういった問題を避けるために,上記のコードではmu2の範囲を-pi()からpi()までに指定してあります。興味のある方は,この範囲をあえて<lower = 0, upper = 2 * pi()>に変えたらどうなるのかをぜひ試してみてください。トレースプロットやn_effを見ると上記の説明が腑に落ちると思います。
  • modelブロック(14-21行目):互いに独立な2つのフォン・ミーゼス分布から群1と群2のデータが生成されるモデルを記述します。Stanでは,von_mises(mu, kappa)でフォン・ミーゼス分布からの乱数を発生させることができます。引数のmuには度数ではなくラジアンを与えます。
  • generated quantitiesブロック(23-44行目):modified_bessel_first_kind(n,x)はn次の第1種変形ベッセル関数です。RのbesselI関数とは引数の順番が異なるので気を付けましょう。このブロックで生成しているxxx_trnsという変数は,ラジアンを度数に変換したものです。パラメータを解釈する時には度数の方が分かりやすいのでここで変換しています。mu_trns_diffはmu1_trnsとmu2_trnsの差の生成量ですが,角度は循環する変数なので,モデルを書く時にはパラメータの範囲や引き算の順序に気を配る必要があります。

 データをStanに渡すRコードは以下の通りです。分かりやすさのために,MCMCのパラメータを明示的に指定しています。

Rコード02: 角度データをStanに渡す
library(rstan)

##データをStanに渡す
data <- list(N1 = n1, N2 = n2, Angle1 = d1, Angle2 = d2)
stanmodel_vm <- stan_model(file="Stancode/von_mises.stan")
fit <- sampling(
  stanmodel_vm,
  data=data,
  seed=123,
  chains=4,
  iter=5500,
  warmup=500,
  thin=1
)

 このコードを走らせて得られたMCMCサンプルをsummary関数で要約すると以下のような結果が得られます。

von_mises_output1.png
図4. パラメータと生成量の推定結果をsummary関数で要約したもの。

 フォン・ミーゼス分布の平均パラメータ(今回の場合はmu1とmu2)の結果をsummary関数から読み取る際には注意が必要です。この出力結果に表示されているmean (=EAP) と分位点の値はMCMCサンプルから計算された値ですが,このMCMCサンプル自体は指定したlowerとupperの範囲におけるパラメータの集まりで,循環する性質を持っていません。例えば,真のμが0度付近の場合にパラメータの範囲を<lower = 0, upper = 2*pi()>と指定してしまうとEAP推定値は180度付近になってしまいますし,分位点に基づく95%ベイズ信頼区間に0度(360度)が含まれなくなります。そのため,フォン・ミーゼス分布の平均パラメータを要約する際には最初に事後分布のヒストグラムを目視し,それからMAP推定値と最高密度区間を計算するという手順を踏むのが良いでしょう。ヒストグラムのピークがグラフの両端に分断されてしまっている場合には,パラメータの範囲を指定し直して推定をやり直すのが賢明です。今回の分析では,μの真値がパラメータ範囲の中央付近となるように範囲をあらかじめ指定しておいたので,図4の結果をそのまま用いてもそれなりに真値と一致します。とはいえ基本的にはMAP推定値と最高密度区間を用いるほうが安全です。

 推定したパラメータと生成量の真値・MAP推定値・95%最高密度区間を以下の表1にまとめました。推定はうまくいっているようですが,群1に比べると群2の最高密度区間の幅が若干広めで,推定精度が相対的に低くなっています。これはもともとの真のκが小さいことに由来すると思われます。結果をまとめると次のようになります:(1) 方向感覚の良い人たちは北をほぼ正確に指差すことができ,個人差も比較的小さい。(2) 方向音痴な人たちには正面(0度)付近を北であると誤って回答する傾向がある。(3) 方向音痴な人たちの回答は個人差が大きく,円周標準偏差で見ると,方向感覚が良い人たちの8〜13倍ばらついている。(4) 方向感覚の良い人たちと悪い人たちの回答の平均値差(最短角距離)は100〜147度程度。※ただし,これらの結果はあくまでも架空データの結果なので,実際にこのような結論が得られるとは限りません。

表1. パラメータと生成量の真値・MAP推定値・95%最高密度区間
vm_summary.png

6. 結び

 この記事では,フォン・ミーゼス分布に従うデータをRとStanで分析する方法について解説しました。今回は分析に架空データを用いましたが,将来的には実際に実験で集めたデータを使って分析をしてみたいものです。個人的には色相データ(実験参加者に色を選ばせる)にフォン・ミーゼス分布を適用したら楽しいことになる予感がします。フォン・ミーゼス分布の混合分布なんかも面白そうです。フォン・ミーゼス分布は実験心理学とは割に相性の良い分布だとは思うのですが,いざ実際に研究で使おうと思うとなかなか難しいものです。面白い適用例などがありましたらTwitterなどで教えて頂けると嬉しいです。あと,明日もStan Advent Calendar 2017の記事を投稿します。正答率データを分析する方法(モデル)をベイズで比較するという趣旨の記事です。そちらもご覧になってくださると嬉しいです!→実験心理学者のためのベイジアンモデリング──正答率データをどうやって分析するか──

引用・参考文献

posted by mutopsy at 00:00 | 統計

2017年12月11日

StanとRで折れ線回帰──空間的視点取得課題の反応時間データを説明する階層ベイズモデルを例に──[スライド紹介]

 こんにちは。mutopsyです。この記事は,『StanとRでベイズ統計モデリング』読書会最終回 (Osaka.Stan #6) のLT(ライトニングトーク)セッションで使用したスライドの紹介記事です。この記事では,このスライドで用いたStanコードの紹介と補足説明を行います。

 このスライドでは,空間的視点取得課題の反応時間データに対し,Stanを使って(階層)折れ線回帰モデル(分割点が1つの場合)を適用した例を紹介しています。その背景や詳しい方法については上記スライドをご覧ください。

 分析で用いたStanコードは以下の通りです。データとしては,各行がある参加者のある条件に対応したデータを用いています。使用している変数(列)は,「実験参加者ID」「角度(0度・40度・80度・120度・160度)」「平均反応時間」の3つです。したがって,行数は,角度の水準数×実験参加者の人数となります。もちろん,以下のコードは他のデータセットにも適用可能です(事前分布やパラメーターの範囲などを修正する必要はありますが)。

Stanコード01: 階層折れ線回帰モデルのStanコード例 (スライド12枚目)
data {
  //分析するデータ
  int N;  //合計試行数 (データの行数)
  int Np;  //実験参加者の人数
  real Angle[N];  //角度の配列
  real RT[N];  //平均反応時間の配列
  int Par[N];  //実験参加者IDの配列(整数)
  
  //以下は,平均値の事後分布を求めるために使用。
 //分析そのものとは無関係なので,不要なら消す。
  int N_new;  //事後分布用のN
  real Angle_new[N_new];  //事後分布用の角度の配列
}

parameters {
  //4つのパラメータの全体平均
  vector[4] mu0;  //BP, b0, b1, b2の全体平均を格納するベクトル

  //実験参加者ごとのパラメータ
  real<lower=1, upper=159> BP[Np];  //分割点:適切な範囲を指定しないとうまく収束しない。
  real<lower=0> b0[Np];  //切片:反応時間なので非負の実数。
  real b1[Np];  //分割点未満の角度における傾き
  real b2[Np];  //分割点以上の角度における傾き
  real<lower = 0> resid[Np];  //残差の標準偏差

  //パラメータ間の分散・共分散構造(階層化のために使用)
  cov_matrix[4] cov;  //パラメータ間の分散・共分散行列
}

transformed parameters{
  //階層化しやすいように,実験参加者ごとのパラメータベクトルを作る。
  vector[4] mu_p[Np];  //1人分のパラメータを格納するベクトル × 人数分
  mu_p[,1] = BP;
  mu_p[,2] = b0;
  mu_p[,3] = b1;
  mu_p[,4] = b2;
}

model {
  //階層化(こうすると短く書けて便利)
  for (np in 1:Np){
    mu_p[np,] ~ multi_normal(mu0, cov);
    //↑4つのパラメータの個人差が正規分布に従うと仮定。
    // 多変量正規分布を指定することで,パラメータ間の相関も推定できる。
  }

  //折れ線回帰モデル
  for (n in 1:N){
    if (Angle[n] < BP[Par[n]]){
      //角度が分割点未満の時のモデル
      RT[n] ~ normal(b0[Par[n]] + (b1[Par[n]] * Angle[n]),
                     resid[Par[n]]);
    } else {
      //角度が分割点以上の時のモデル
      RT[n] ~ normal(b0[Par[n]] + (b1[Par[n]] * BP[Par[n]]) 
                      + (b2[Par[n]] * (Angle[n] - BP[Par[n]])),
                     resid[Par[n]]);
    }
  }

  //事前分布
  for (np in 1:Np){
    resid[np] ~ cauchy(0,1000);  //残差の標準偏差の事前分布(半コーシー分布)
    BP[np] ~ normal(90,180);  //分割点の事前分布として弱情報を与えておく
  }
}

generated quantities{
  corr_matrix[4] rho;  //4つのパラメータの相関行列
  real RT_new[N_new];  //平均値の事後分布用

  //パラメータ間の相関係数を求める
  for (i in 1:4){
    for (j in 1:4){
      rho[i,j] = cov[i,j]/sqrt(cov[i,i]*cov[j,j]);
    }
  }

  //平均値の事後分布を求める(不要なら削除)
  for (n in 1:N_new){
    if (Angle_new[n] < mu0[1]){
      RT_new[n] = mu0[2] + (mu0[3] * Angle_new[n]);
    } else {
      RT_new[n] = mu0[2] + (mu0[3] * mu0[1]) + (mu0[4] * (Angle_new[n] - mu0[1]));
    }
  }
}


 なお,スライドでも触れている通り,if文を含むモデルはStanでは収束効率が悪くなります。モデルの書き方を変えることで収束効率が改善する可能性があるので,後日,また比較してみて記事にしようと思います。

 ちなみに余談ですが,この記事およびスライドのタイトル(「StanとRで〜」)は,アヒル本の書名(『StanとRでベイズ統計モデリング』)のオマージュです。Osaka.Stanでのアヒル本読書会最終回ということで,それにちなんだタイトルにしてみました。アヒル本読書会は終わりましたが,Osaka.Stan自体はまだまだ続くそうですので,乞うご期待。

posted by mutopsy at 04:00 | 統計