# はじめに

[Linear Ballistic Accumulator modelを用いた反応時間分析](https://ykunisato.github.io/lbaStan/lbaStan.html)では，１名の参加者の反応時間データをRstanのユーザー定義関数を活用しながら解析する方法について解説を行いました。その際に，試行数は500施行という現実的ではない値を用いていました。今回は，10名が100試行反応するという，心理学実験でありえそうな状況を想定してみます。なお，以下でStanコードを載せていますが，HTML上で表記がおかしくなることがあるので，自分の環境で動かしてみたい方は，[ここ](hlbaStan.Rmd)から適宜RmdファイルをDLください。

# 使用パッケージ

今回は，以下のパッケージを用います。

```{r, message=FALSE, warning=FALSE}
rm(list=ls())
library(DiagrammeR)
library(ggplot2)
library(plotly)
library(bayesplot)
library(rstan)
library(loo)
library(rtdists)
library(msm)
```

# 10名分のデータの作成

Annis et al. (2016)による [Bayesian inference with Stan: A tutorial on adding custom distributions](http://www.ncbi.nlm.nih.gov/pubmed/27287444)のコードを参考にしつつ，10名の参加者が３条件の２選択肢の課題をそれぞれ100試行行ったというデータを作ってみる。

## パタメータの準備
まずは，人数分のパラメータを準備します。その際に，vは条件によって異なる（変化するor条件の効果がある）という設定にします。

```{r, message=FALSE, warning=FALSE}
set.seed(123)
# 参加者数，試行数，条件数
numberSubjects <- 10
numberTrials <- 100
numberConditions <- 3

#　ハイパーパラメータの設定
A_mu = .5
b_mu = 1
t0_mu = .5
v1_mu_c1 = 2
v2_mu_c1 = 2.5
v1_mu_c2 = 3
v2_mu_c2 = 3
v1_mu_c3 = 3.5
v2_mu_c3 = 2.5
s = 1

#　ハイパーパラメータをもとにした個人のパラメータ（切断正規分布を利用）
A = rtnorm(numberSubjects,A_mu,.5,0,Inf)
b = rtnorm(numberSubjects,b_mu,.5,A,Inf)
t0 = rtnorm(numberSubjects,t0_mu,.5,0,1)
v1_c1 = rtnorm(numberSubjects,v1_mu_c1,1,0,Inf)
v2_c1 = rtnorm(numberSubjects,v2_mu_c1,1,0,Inf)
v1_c2 = rtnorm(numberSubjects,v1_mu_c2,1,0,Inf)
v2_c2 = rtnorm(numberSubjects,v2_mu_c2,1,0,Inf)
v1_c3 = rtnorm(numberSubjects,v1_mu_c3,1,0,Inf)
v2_c3 = rtnorm(numberSubjects,v2_mu_c3,1,0,Inf)

# パラメータとハイパーパラメータのデータフレーム化
hyperParameters <- data.frame(A_mu,b_mu,t0_mu,v1_mu_c1,v2_mu_c1,v1_mu_c2,v2_mu_c2,v1_mu_c3,v2_mu_c3,s)
print(hyperParameters)
parameters <- data.frame(A,b,t0,v1_c1,v2_c1,v1_c2,v2_c2,v1_c3,v2_c3)
print(parameters)
```

## 反応と反応時間データの作成

上記で設定したパラメータをもとに，rtdistsのrLBA関数を用いて反応時間と反応のデータを作ります。作ったデータは，ロング形式にしました。そして，stanLongDataにStanで使えるようにリスト形式でデータを保存しています。

```{r, message=FALSE, warning=FALSE}
longData <- NULL
for(i in 1:numberSubjects){
    id <- rep(i,numberTrials)
    # 条件１
    tempRt <- rLBA(numberTrials, A=parameters$A[i], b=parameters$b[i], t0=parameters$t0[i], mean_v=c(parameters$v1_c1[i],parameters$v2_c1[i]), sd_v=c(s,s))
    condition <- rep(1,numberTrials)
    indDataC1 <- data.frame(id,condition,tempRt)
    
    # 条件2
    tempRt <- rLBA(numberTrials, A=parameters$A[i], b=parameters$b[i], t0=parameters$t0[i], mean_v=c(parameters$v1_c2[i],parameters$v2_c2[i]), sd_v=c(s,s))
    condition <- rep(2,numberTrials)
    indDataC2 <- data.frame(id,condition,tempRt)
    
    # 条件3
    tempRt <- rLBA(numberTrials, A=parameters$A[i], b=parameters$b[i], t0=parameters$t0[i], mean_v=c(parameters$v1_c3[i],parameters$v2_c3[i]), sd_v=c(s,s))
    condition <- rep(3,numberTrials)
    indDataC3 <- data.frame(id,condition,tempRt)
    indData <- rbind(indDataC1,indDataC2,indDataC3)
    
    # ロングデータ
    longData <- rbind(longData,indData)
}

numberAllTrials <- length(longData$rt)
# Stan用データの作成
stanLongData <- list(rt=longData$rt,response=longData$response,id=longData$id,condition = longData$condition, numberConditions=numberConditions,numberSubjects=numberSubjects,numberChoices=2,numberAllTrials=numberAllTrials)
```

## データの確認

longDataは，id，condition,rt,responseからなります。idは10名，conditionは3条件です。
```{r, message=FALSE, warning=FALSE}
head(longData)
```

全員の結果を確認するわけにもいかないので，適当に，ID=1の条件１，２，３の結果を見てみることにします。

```{r, message=FALSE, warning=FALSE}
# 参加者IDが1の人の各条件のデータを作る。
dataSub1 <- subset(longData,longData$id==1)
dataSub1Con1 <- subset(dataSub1,dataSub1$condition==1)
dataSub1Con2 <- subset(dataSub1,dataSub1$condition==2)
dataSub1Con3 <- subset(dataSub1,dataSub1$condition==3)
```

各条件の選択肢別の反応数

```{r, message=FALSE, warning=FALSE}
# 条件１
table(dataSub1Con1$response)
# 条件２
table(dataSub1Con2$response)
# 条件３
table(dataSub1Con3$response)
```

条件１の選択肢ごとの反応時間のヒストグラム

```{r, message=FALSE, warning=FALSE}
# 選択肢1の反応時間
dataSub1Con1Choice1 <- subset(dataSub1Con1,dataSub1Con1$response==1)
p1 <- ggplot(dataSub1Con1Choice1,aes(x = rt)) + geom_histogram()+coord_cartesian(xlim = c(0, 3), ylim = c(0, 10))

# 選択肢2の反応時間
dataSub1Con1Choice2 <- subset(dataSub1Con1,dataSub1Con1$response==2)
p2 <- ggplot(dataSub1Con1Choice2,aes(x = rt)) + geom_histogram()+coord_cartesian(xlim = c(0, 3), ylim = c(0, 10))
subplot(p1, p2)
```


条件２の選択肢ごとの反応時間のヒストグラム

```{r, message=FALSE, warning=FALSE}
# 選択肢1の反応時間
dataSub1Con2Choice1 <- subset(dataSub1Con2,dataSub1Con2$response==1)
p1 <- ggplot(dataSub1Con2Choice1,aes(x = rt)) + geom_histogram()+coord_cartesian(xlim = c(0, 3), ylim = c(0, 10))

# 選択肢2の反応時間
dataSub1Con2Choice2 <- subset(dataSub1Con2,dataSub1Con2$response==2)
p2 <- ggplot(dataSub1Con2Choice2,aes(x = rt)) + geom_histogram()+coord_cartesian(xlim = c(0, 3), ylim = c(0, 10))
subplot(p1, p2)
```


条件３の選択肢ごとの反応時間のヒストグラム

```{r, message=FALSE, warning=FALSE}
# 選択肢1の反応時間
dataSub1Con3Choice1 <- subset(dataSub1Con3,dataSub1Con3$response==1)
p1 <- ggplot(dataSub1Con3Choice1,aes(x = rt)) + geom_histogram()+coord_cartesian(xlim = c(0, 3), ylim = c(0, 10))

# 選択肢2の反応時間
dataSub1Con3Choice2 <- subset(dataSub1Con3,dataSub1Con3$response==2)
p2 <- ggplot(dataSub1Con3Choice2,aes(x = rt)) + geom_histogram()+coord_cartesian(xlim = c(0, 3), ylim = c(0, 10))
subplot(p1, p2)
```



# 階層ベイズモデル

　今回は，10名の参加者が３つの条件にたいしてそれぞれ100試行反応するというデータを作成しました。これを階層ベイズを使って推定します。参加者ごとにLinear Ballistic Accumulator Modelのパラメータがありますが，それらが切断正規分布に従うという事前分布をおいて，階層ベイズによる推定を行います。また，同じ参加者が３つの条件を経験しているので，こちらも階層化しますが，条件の効果がパラメータにかかる場合とかからない場合でモデルを組みます。
　条件の効果は，vのみにかかることにします。おそらく，無判断時間や閾値や開始点の上限は，条件によって変わらないのではないかという前提のものとで，vに条件の効果がかかるかどうかを検討します（実際は，もう少し複数のモデルを検討しても良いかと思います）。

## 条件による効果なしモデル

まず，k,A,v,tauは参加者ごとに異なるが，条件によってパラメータは異ならないという前提で解析を行います。グラフィカルモデルを描くと以下のようになる（Conditionは不要ですが，比較のために描いています）。

### グラフィカルモデル

```{r}
grViz("
  digraph dot {
    graph [splines = line,compound = true, nodesep = .5, ranksep = .25,
           color = black, label='Hierarchical Bayes Model of Linear Ballistic Accumulator']

      node [shape = circle,style = filled,fillcolor = white,color = black,label = 'k@_{n}'] k
      node [label = '&mu;@_{k}'] mk
      node [label = '&sigma;@_{k}'] sk
      node [label = 'A@_{n}'] A
      node [label = '&mu;@_{A}'] mA
      node [label = '&sigma;@_{A}'] sA
      node [label = '&tau;@_{n}'] t
      node [label = '&mu;@_{&tau;}'] mt
      node [label = '&sigma;@_{&tau;}'] st
      node [label = 'v[1]@_{n}'] v1
      node [label = '&mu;@_{v[1]}'] mv1
      node [label = '&sigma;@_{v[1]}'] sv1
      node [label = 'v[2]@_{n}'] v2
      node [label = '&mu;@_{v[2]}'] mv2
      node [label = '&sigma;@_{v[2]}'] sv2
      node [label = 'S'] s
      node [label = 'LBA@_{n,c,t}'] lba
      node [fillcolor = grey,label = 'RT@_{n,c,t}'] rt

        subgraph cluster1 {
          labelloc=b
          label = 'Participants n = 1...N'
            subgraph cluster2 {
              labelloc=b
              label = 'Condition c = 1...C'
                subgraph cluster3 {
                  labelloc=b
                  label = 'Trials t = 1...T'
                    edge [color = black]
                      lba -> rt
                }
              edge [color = black]
            }
          edge [color = black]
            k -> lba
            t -> lba
            A -> lba
            v1 -> lba
            v2 -> lba
        }
      edge [color = black]
        mk -> k
        sk -> k
        mA -> A
        sA -> A
        mt -> t
        st -> t
        mv1 -> v1
        sv1 -> v1
        mv2 -> v2
        sv2 -> v2 [tailport=s,headport=e]
        s -> lba
  }",engine = "dot")
```

事前分布はAnnisらを参考に以下にしました(ハイパーパラメータの分散の事前分布のみGamma(1,1)からstudent T(4,0,1変更）。
$$ 
k_{n} \sim Normal(\mu_{k},\sigma_{k})T(0,\infty)\\
A_{n} \sim Normal(\mu_{A},\sigma_{A})T(0,\infty)\\
\tau_{n} \sim Normal(\mu_{\tau},\sigma_{\tau})T(0,\infty)\\
v[1]_{n} \sim Normal(\mu_{v[1]},\sigma_{v[1]})T(0,\infty)\\
v[2]_{n} \sim Normal(\mu_{v[2]},\sigma_{v[2]})T(0,\infty)\\
\mu_{A},\mu_{k} \sim Normal(0.5,1)T(0,\infty)\\
\mu_{\tau} \sim Normal(0.5,0.5)T(0,\infty)\\
\mu_{v[1]},\mu_{v[2]} \sim Normal(2,1)T(0,\infty)\\
\sigma_{A},\sigma_{k}, \sigma_{v[1]},\sigma_{v[2]},\sigma_{\tau} \sim student\,t(4,0,1)T(0,\infty)\\
S = 1 \hspace{30pt}※Sは定数\\
RT[n,c, t] \sim LBA(k_{n},A_{n},\tau_{n},v[1]_{n},s) [選択肢１の時]\\
RT[n,c, t] \sim LBA(k_{n},A_{n},\tau_{n},v[2]_{n},s) [選択肢2の時]\\
$$

### Stanコード

基本的には，[「RStanによる反応時間の解析：Linear Ballistic Accumulator modelを用いて」](https://ykunisato.github.io/lbaStan/lbaStan.html)で解説したユーザー定義関数を用います。

```{stan, output.var="lbaModel1"}
functions{
     
     real lba_pdf(real t, real b, real A, real v, real s){
          //PDF of the LBA model
          
          real b_A_tv_ts;
          real b_tv_ts;
          real term_1;
          real term_2;
          real term_3;
          real term_4;
          real pdf;
          
          b_A_tv_ts = (b - A - t*v)/(t*s);
          b_tv_ts = (b - t*v)/(t*s);
          term_1 = v*Phi(b_A_tv_ts);
          term_2 = s*exp(normal_lpdf(b_A_tv_ts|0,1));
          term_3 = v*Phi(b_tv_ts);
          term_4 = s*exp(normal_lpdf(b_tv_ts|0,1)); 
          pdf = (1 / A)*(-term_1 + term_2 + term_3 - term_4);
          
          return pdf;
     }
     
     real lba_cdf(real t, real b, real A, real v, real s){
          //CDF of the LBA model
          
          real b_A_tv;
          real b_tv;
          real ts;
          real term_1;
          real term_2;
          real term_3;
          real term_4;
          real cdf;	
          
          b_A_tv = b - A - t*v;
          b_tv = b - t*v;
          ts = t*s;
          term_1 = b_A_tv/A * Phi(b_A_tv/ts);	
          term_2 = b_tv/A   * Phi(b_tv/ts);
          term_3 = ts/A     * exp(normal_lpdf(b_A_tv/ts|0,1)); 
          term_4 = ts/A     * exp(normal_lpdf(b_tv/ts|0,1)); 
          cdf = 1 + term_1 - term_2 + term_3 - term_4;
          
          return cdf;
          
     }
     
     real lba_lpdf(real rt, real res, real k, real A, vector v, real s, real tau){
          
          real t;
          real b;
          real cdf;
          real pdf;		
          real prob;
          real out;
          real prob_neg;

          b = A + k;
          t = rt - tau;
          
          if(t > 0){			
                cdf = 1;
                for(j in 1 : num_elements(v)){
                      if(res == j){
                            pdf = lba_pdf(t, b, A, v[j], s);
                      }else{	
                            cdf = (1 - lba_cdf(t, b, A, v[j], s)) * cdf;
                      }
                }
                
                prob_neg = 1;
                for(j in 1 : num_elements(v)){
                      prob_neg = Phi(-v[j]/s) * prob_neg;    
                }
                prob = pdf*cdf;		
                prob = prob/(1 - prob_neg);	
                if(prob < 1e-10){
                      prob = 1e-10;				
                }
          }else{
              prob = 1e-10;			
          }
          out = log(prob);
          return out;		
     }
     
    vector lba_rng(real k, real A, vector v, real s, real tau){
          
          int get_pos_drift;	
          int no_pos_drift;
          int get_first_pos;
          vector[num_elements(v)] drift;
          int max_iter;
          int iter;
          real start[num_elements(v)];
          real ttf[num_elements(v)];
          int resp[num_elements(v)];
          real rt;
          vector[2] pred;
          real b;
          
          //try to get a positive drift rate
          get_pos_drift = 1;
          no_pos_drift = 0;
          max_iter = 1000;
          iter = 0;
          while(get_pos_drift){
               for(j in 1 : num_elements(v)){
                    drift[j] = normal_rng(v[j],s);
                    if(drift[j] > 0){
                         get_pos_drift = 0;
                    }
               }
               iter = iter + 1;
               if(iter > max_iter){
                    get_pos_drift = 0;
                    no_pos_drift = 1;
               }	
          }
          //if both drift rates are <= 0
          //return an infinite response time
          if(no_pos_drift){
               pred[1] = -1;
               pred[2] = -1;
          }else{
               b = A + k;
               for(i in 1 : num_elements(v)){
                    //start time of each accumulator	
                    start[i] = uniform_rng(0, A);
                    //finish times
                    ttf[i] = (b-start[i])/drift[i];
               }
               //rt is the fastest accumulator finish time	
               //if one is negative get the positive drift
               resp = sort_indices_asc(ttf);
               ttf = sort_asc(ttf);
               get_first_pos = 1;
               iter = 1;
               while(get_first_pos){
                    if(ttf[iter] > 0){
                         pred[1] = ttf[iter] + tau;
                         pred[2] = resp[iter]; 
                         get_first_pos = 0;
                    }
                    iter = iter + 1;
               }
          }
          return pred;	
     }
}

data{
     int numberConditions;
     int numberSubjects;
     int numberChoices;
     int numberAllTrials;
     vector[numberAllTrials] rt;
     vector[numberAllTrials] response;
     int condition[numberAllTrials];
     int id[numberAllTrials];
}

parameters {
     real < lower = 0 > k[numberSubjects];
     real < lower = 0 > A[numberSubjects];
     real < lower = 0 > psi[numberSubjects];
     vector < lower = 0 > [numberChoices] v [numberSubjects];
     real < lower = 0 > k_sigma;
     real < lower = 0 > A_sigma;
     real < lower = 0 > psi_sigma;
     vector < lower = 0 > [numberChoices] v_sigma;
     real < lower = 0 > k_mu;
     real < lower = 0 > A_mu;
     real < lower = 0 > psi_mu;
     vector < lower=0> [numberChoices] v_mu;
}

transformed parameters {
     real s;
     s = 1;
}

model {
     int subIndex;
     // ハイパーパラメータの事前分布の設定
     k_mu ~ normal(.5, 1)T[0, ];
     A_mu ~ normal(.5, 1)T[0, ];
     psi_mu ~ normal(.5, .5)T[0, ];
     k_sigma ~ student_t(4, 0, 1)T[0, ];
     A_sigma ~ student_t(4, 0, 1)T[0, ];
     psi_sigma ~ student_t(4, 0, 1)T[0, ];
     
     for(n in 1 : numberChoices){
        v_mu[n] ~ normal(2, 1)T[0, ];
        v_sigma[n] ~ student_t(4, 0, 1)T[0, ];
      }
     
     // パラメータの事前分布の設定
     for(i in 1 : numberSubjects){		
          k[i] ~ normal(k_mu, k_sigma)T[0, ];
          A[i] ~ normal(A_mu, A_sigma)T[0, ];
          psi[i] ~ normal(psi_mu, psi_sigma)T[0, ];

          for(n in 1 : numberChoices){
              v[i,n] ~ normal(v_mu[n],v_sigma[n])T[0, ];
          }
     }
    
    // 対数事後確率の計算部分
    for(m in 1 : numberAllTrials){
           subIndex = id[m];
          rt[m] ~ lba(response[m], k[subIndex], A[subIndex], v[subIndex], s, psi[subIndex]);
    }
}

generated quantities {
     vector[2] pred[numberSubjects];
     vector[numberAllTrials] log_lik;
     int subIndex;
     for(i in 1 : numberSubjects){
          pred[i] = lba_rng(A[i], k[i], v[i], s, psi[i]);
     }
    for(m in 1 : numberAllTrials){
          subIndex = id[m];
          log_lik[m] = lba_lpdf(rt[m] | response[m], k[subIndex], A[subIndex], v[subIndex], s, psi[subIndex]);
    }
}
```

## 条件によるvへの効果ありモデル

次は，条件によるv（ドリフト率の平均）への効果を仮定したモデルになります。おそらく条件によって異なるという前提がある場合は，こちらのほうが自然に思います。

### グラフィカルモデル

```{r}
grViz("
  digraph dot {
    graph [splines = line,compound = true, nodesep = .5, ranksep = .25,
           color = black, label='Hierarchical Bayes Model of Linear Ballistic Accumulator']

      node [shape = circle,style = filled,fillcolor = white,color = black,label = 'k@_{n}'] k
      node [label = '&mu;@_{k}'] mk
      node [label = '&sigma;@_{k}'] sk
      node [label = 'A@_{n}'] A
      node [label = '&mu;@_{A}'] mA
      node [label = '&sigma;@_{A}'] sA
      node [label = '&tau;@_{n}'] t
      node [label = '&mu;@_{&tau;}'] mt
      node [label = '&sigma;@_{&tau;}'] st
      node [label = 'v[1]@_{n,c}'] v1
      node [label = '&mu;@_{v[1]}'] mv1
      node [label = '&sigma;@_{v[1]}'] sv1
      node [label = 'v[2]@_{n,c}'] v2
      node [label = '&mu;@_{v[2]}'] mv2
      node [label = '&sigma;@_{v[2]}'] sv2
      node [label = 'S'] s
      node [label = 'LBA@_{n,c,t}'] lba
      node [fillcolor = grey,label = 'RT@_{n,c,t}'] rt

        subgraph cluster1 {
          labelloc=b
          label = 'Participants n = 1...N'
            subgraph cluster2 {
              labelloc=b
              label = 'Condition c = 1...C'
                subgraph cluster3 {
                  labelloc=b
                  label = 'Trials t = 1...T'
                    edge [color = black]
                      lba -> rt
                }
              edge [color = black]
                v1 -> lba
                v2 -> lba
            }
          edge [color = black]
            A -> lba
            k -> lba
            t -> lba
        }
      edge [color = black]
        mk -> k
        sk -> k
        mA -> A
        sA -> A
        mt -> t
        st -> t
        mv1 -> v1
        sv1 -> v1
        mv2 -> v2
        sv2 -> v2 [tailport=s,headport=e]
        s -> lba
  }",engine = "dot")
```

事前分布はAnnisらを参考に以下にしました。
$$ 
k_{i} \sim Normal(\mu_{k},\sigma_{k})T(0,\infty)\\
A_{i} \sim Normal(\mu_{A},\sigma_{A})T(0,\infty)\\
\tau_{i} \sim Normal(\mu_{\tau},\sigma_{\tau})T(0,\infty)\\
v1_{i} \sim Normal(\mu_{v1},\sigma_{v1})T(0,\infty)\\
v2_{i} \sim| Normal(\mu_{v2},\sigma_{v2})T(0,\infty)\\
\mu^{A},\mu^{k} \sim Normal(0.5,1)T(0,\infty)\\
\mu^{\tau} \sim Normal(0.5,0.5)T(0,\infty)\\
\mu^{v1},\mu^{v2} \sim Normal(2,1)T(0,\infty)\\
\sigma^{A},\sigma^{k}, \sigma^{v1},\sigma^{v2},\sigma^{\tau} \sim student_t(4,0,1)T(0,\infty) [Gamma(1,1)から変更]\\ 
S = 1 [Sは定数]\\
RT[t, i] \sim LBA(k_{i},A_{i},\tau_{i},v1_{i,j},s) [選択肢１の時]\\
RT[t, i] \sim LBA(k_{i},A_{i},\tau_{i},v2_{i,j},s) [選択肢2の時]\\
k \sim Normal(.5,1)T[0,]\\
A \sim Normal(.5,1)T[0,]\\
\tau \sim Normal(.5,.5)T[0,]\\
v[1],v[2] \sim Normal(2,1)T[0,]\\
S = 1 \hspace{30pt}※Sは定数
$$


### Stanコード

```{stan, output.var="lbaModel2"}
functions{
     
     real lba_pdf(real t, real b, real A, real v, real s){
          //PDF of the LBA model
          
          real b_A_tv_ts;
          real b_tv_ts;
          real term_1;
          real term_2;
          real term_3;
          real term_4;
          real pdf;
          
          b_A_tv_ts = (b - A - t*v)/(t*s);
          b_tv_ts = (b - t*v)/(t*s);
          term_1 = v*Phi(b_A_tv_ts);
          term_2 = s*exp(normal_lpdf(b_A_tv_ts|0,1));
          term_3 = v*Phi(b_tv_ts);
          term_4 = s*exp(normal_lpdf(b_tv_ts|0,1)); 
          pdf = (1 / A)*(-term_1 + term_2 + term_3 - term_4);
          
          return pdf;
     }
     
     real lba_cdf(real t, real b, real A, real v, real s){
          //CDF of the LBA model
          
          real b_A_tv;
          real b_tv;
          real ts;
          real term_1;
          real term_2;
          real term_3;
          real term_4;
          real cdf;	
          
          b_A_tv = b - A - t*v;
          b_tv = b - t*v;
          ts = t*s;
          term_1 = b_A_tv/A * Phi(b_A_tv/ts);	
          term_2 = b_tv/A   * Phi(b_tv/ts);
          term_3 = ts/A     * exp(normal_lpdf(b_A_tv/ts|0,1)); 
          term_4 = ts/A     * exp(normal_lpdf(b_tv/ts|0,1)); 
          cdf = 1 + term_1 - term_2 + term_3 - term_4;
          
          return cdf;
          
     }
     
     real lba_lpdf(real rt, real res, real k, real A, vector v, real s, real tau){
          
          real t;
          real b;
          real cdf;
          real pdf;		
          real prob;
          real out;
          real prob_neg;

          b = A + k;
          t = rt - tau;
          
          if(t > 0){			
                cdf = 1;
                for(j in 1 : num_elements(v)){
                      if(res == j){
                            pdf = lba_pdf(t, b, A, v[j], s);
                      }else{	
                            cdf = (1 - lba_cdf(t, b, A, v[j], s)) * cdf;
                      }
                }
                
                prob_neg = 1;
                for(j in 1 : num_elements(v)){
                      prob_neg = Phi(-v[j]/s) * prob_neg;    
                }
                prob = pdf*cdf;		
                prob = prob/(1 - prob_neg);	
                if(prob < 1e-10){
                      prob = 1e-10;				
                }
          }else{
              prob = 1e-10;			
          }
          out = log(prob);
          return out;		
     }
     
    vector lba_rng(real k, real A, vector v, real s, real tau){
          
          int get_pos_drift;	
          int no_pos_drift;
          int get_first_pos;
          vector[num_elements(v)] drift;
          int max_iter;
          int iter;
          real start[num_elements(v)];
          real ttf[num_elements(v)];
          int resp[num_elements(v)];
          real rt;
          vector[2] pred;
          real b;
          
          //try to get a positive drift rate
          get_pos_drift = 1;
          no_pos_drift = 0;
          max_iter = 1000;
          iter = 0;
          while(get_pos_drift){
               for(j in 1 : num_elements(v)){
                    drift[j] = normal_rng(v[j],s);
                    if(drift[j] > 0){
                         get_pos_drift = 0;
                    }
               }
               iter = iter + 1;
               if(iter > max_iter){
                    get_pos_drift = 0;
                    no_pos_drift = 1;
               }	
          }
          //if both drift rates are <= 0
          //return an infinite response time
          if(no_pos_drift){
               pred[1] = -1;
               pred[2] = -1;
          }else{
               b = A + k;
               for(i in 1 : num_elements(v)){
                    //start time of each accumulator	
                    start[i] = uniform_rng(0, A);
                    //finish times
                    ttf[i] = (b-start[i])/drift[i];
               }
               //rt is the fastest accumulator finish time	
               //if one is negative get the positive drift
               resp = sort_indices_asc(ttf);
               ttf = sort_asc(ttf);
               get_first_pos = 1;
               iter = 1;
               while(get_first_pos){
                    if(ttf[iter] > 0){
                         pred[1] = ttf[iter] + tau;
                         pred[2] = resp[iter]; 
                         get_first_pos = 0;
                    }
                    iter = iter + 1;
               }
          }
          return pred;	
     }
}

data{
     int numberConditions;
     int numberSubjects;
     int numberChoices;
     int numberAllTrials;
     vector[numberAllTrials] rt;
     vector[numberAllTrials] response;
     int condition[numberAllTrials];
     int id[numberAllTrials];
}

parameters {
     real < lower = 0 > k[numberSubjects];
     real < lower = 0 > A [numberSubjects];
     real < lower = 0 > psi[numberSubjects];
     vector < lower = 0 > [numberChoices] v [numberSubjects, numberConditions];
     real < lower = 0 > k_sigma;
     real < lower = 0 > A_sigma;
     real < lower = 0 > psi_sigma;
     vector < lower = 0 > [numberChoices] v_sigma [numberConditions];
     real < lower = 0 > k_mu;
     real < lower = 0 > A_mu;
     real < lower = 0 > psi_mu;
     vector < lower=0> [numberChoices] v_mu [numberConditions];
}

transformed parameters {
     real s;
     s = 1;
}

model {
     int subIndex;
     int condIndex;
     // ハイパーパラメータの事前分布の設定
     k_mu ~ normal(.5, 1)T[0, ];
     psi_mu ~ normal(.5, .5)T[0, ];
     k_sigma ~ student_t(4, 0, 1)T[0, ];
     psi_sigma ~ student_t(4, 0, 1)T[0, ];
     A_sigma ~ student_t(4, 0, 1)T[0, ];
     A_mu ~ normal(.5, 1)T[0, ];
     
     for( j in 1 : numberConditions){
          for(n in 1 : numberChoices){
            v_mu[j, n] ~ normal(2, 1)T[0, ];
            v_sigma[j, n] ~ student_t(4, 0, 1)T[0, ];
          }
     }
      
     // パラメータの事前分布の設定
     for(i in 1 : numberSubjects){		
        k[i] ~ normal(k_mu, k_sigma)T[0, ];
        psi[i] ~ normal(psi_mu, psi_sigma)T[0, ];
        A[i] ~ normal(A_mu, A_sigma)T[0, ];
          for( j in 1 : numberConditions){
              for(n in 1 : numberChoices){
                v[i, j, n] ~ normal(v_mu[j, n],v_sigma[j, n])T[0, ];
              }
          }
     }
    
    for(m in 1 : numberAllTrials){
      subIndex = id[m];
      condIndex = condition[m];
      rt[m] ~ lba(response[m], k[subIndex], A[subIndex], v[subIndex, condIndex], s, psi[subIndex]);
    }
}

generated quantities {
     vector[2] pred[numberSubjects, numberConditions];
     vector[numberAllTrials] log_lik;
     int subIndex;
     int condIndex;
     
     for(i in 1 : numberSubjects){
          for(j in 1 : numberConditions){
               pred[i, j] = lba_rng(k[i], A[i], v[i, j], s, psi[i]);
          }
     }

    for(m in 1 : numberAllTrials){
          subIndex = id[m];
          condIndex = condition[m];
          log_lik[m] = lba_lpdf(rt[m] | response[m], k[subIndex], A[subIndex], v[subIndex, condIndex], s, psi[subIndex]);
    }
     
}
```

## Rstanによるサンプリング

さてさて，サンプリングしてみます。まずは，条件による効果なしモデルからサンプリングを行います。

```{r, message=FALSE, warning=FALSE}
rstan_options(auto_write=TRUE)
options(mc.cores = parallel::detectCores())
hmcIter = 1000
hmcChains = 4
hmcWarmup = 300
hmcThin = 2

fit1 <- sampling(lbaModel1, 
            seed = 1234, 
            data = stanLongData,
            warmup = hmcWarmup, 
            iter = hmcIter,
            chains = hmcChains,
            thin =hmcThin)
```
次に，条件による効果ありモデルのサンプリングを行います。
```{r, message=FALSE, warning=FALSE}
fit2 <- sampling(lbaModel2, 
            seed = 1234, 
            data = stanLongData,
            warmup = hmcWarmup, 
            iter = hmcIter,
            chains = hmcChains,
            thin =hmcThin)
```

## 収束判定(条件による効果なしモデル){.tabset}
収束してそうですね。収束判定には，bayesplotを活用しています。[このページ](https://cran.rstudio.com/web/packages/bayesplot/vignettes/MCMC-diagnostics.html)が参考になりました。なお，定数のSは除外しています。

### トレースプロット
```{r, message=FALSE, warning=FALSE}
posterior11 <- rstan::extract(fit1, inc_warmup = TRUE, permuted = FALSE)
# k[1:5]
color_scheme_set("mix-blue-pink")
pk1 <- mcmc_trace(posterior11,  pars = c("k[1]", "k[2]", "k[3]", "k[4]", "k[5]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pk1 <- pk1+ facet_text(size = 15)
plot(pk1)
# k[6:10]
pk2 <- mcmc_trace(posterior11,  pars = c("k[6]", "k[7]", "k[8]", "k[9]", "k[10]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pk2 <- pk2+ facet_text(size = 15)
plot(pk2)
# A[1:5]
pA1 <- mcmc_trace(posterior11,  pars = c("A[1]", "A[2]", "A[3]", "A[4]", "A[5]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pA1 <- pA1+ facet_text(size = 15)
plot(pA1)
# A[6:10]
pA2 <- mcmc_trace(posterior11,  pars = c( "A[6]", "A[7]", "A[8]", "A[9]", "A[10]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pA2 <- pA2+ facet_text(size = 15)
plot(pA2)
# tau[1:5]
pt <- mcmc_trace(posterior11,  pars = c("psi[1]", "psi[2]", "psi[3]", "psi[4]", "psi[5]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pt <- pt+ facet_text(size = 15)
plot(pt)
# tau[6:10]
pt <- mcmc_trace(posterior11,  pars = c("psi[6]", "psi[7]", "psi[8]", "psi[9]", "psi[10]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pt <- pt+ facet_text(size = 15)
plot(pt)
# v[1:3]
pv1 <- mcmc_trace(posterior11,  pars = c("v[1,1]", "v[1,2]", "v[2,1]", "v[2,2]", "v[3,1]", "v[3,2]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pv1 <- pv1+ facet_text(size = 15)
plot(pv1)
# v[4:6]
pv2 <- mcmc_trace(posterior11,  pars = c("v[4,1]", "v[4,2]", "v[5,1]", "v[5,2]", "v[6,1]", "v[6,2]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pv2 <- pv2+ facet_text(size = 15)
plot(pv2)
# v[7:8]
pv3 <- mcmc_trace(posterior11,  pars = c("v[7,1]", "v[7,2]", "v[8,1]", "v[8,2]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pv3 <- pv3+ facet_text(size = 15)
plot(pv3)
# v[9:10]
pv4 <- mcmc_trace(posterior11,  pars = c("v[9,1]", "v[9,2]", "v[10,1]", "v[10,2]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pv4 <- pv4+ facet_text(size = 15)
plot(pv4)
# hyperparameter 1
ph1 <- mcmc_trace(posterior11,  pars = c("k_mu", "A_mu", "psi_mu", "v_mu[1]", "v_mu[2]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
ph1 <- ph1+ facet_text(size = 15)
plot(ph1)
# hyperparameter 2
ph2 <- mcmc_trace(posterior11,  pars = c("k_sigma", "A_sigma", "psi_sigma", "v_sigma[1]", "v_sigma[2]","lp__"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
ph2 <- ph2+ facet_text(size = 15)
plot(ph2)
```

### $\hat{R}$の確認
bayesplotで確認をしても良いですが，ちょっとパラメータ多すぎるので，一括で検討します。$\hat{R}$が1.10以下ならTRUE，1.10以上のものがあればFALSEになります。
```{r, message=FALSE, warning=FALSE}
all(stan_rhat(fit1)$data <=1.10,na.rm=T)
```
    
### 有効サンプルサイズ
bayesplotで確認をしても良いですが，ちょっとパラメータ多すぎるので，一括で検討します。有効サンプルサイズが0.1以上ならTRUE，0.1以下のものがあればFALSEになります。
```{r, message=FALSE, warning=FALSE}
all(stan_ess(fit1)$data >= 0.1,na.rm=T)
```

### 自己相関
```{r, message=FALSE, warning=FALSE}
posterior12 <- as.matrix(fit1)
mcmc_acf(posterior12, par = c("k[1]","k[2]","k[3]","k[4]","k[5]"), lags = 10)
mcmc_acf(posterior12, par = c("k[6]","k[7]","k[8]","k[9]","k[10]"), lags = 10)
mcmc_acf(posterior12, par = c( "A[1]", "A[2]", "A[3]", "A[4]", "A[5]"), lags = 10)
mcmc_acf(posterior12, par = c("A[6]", "A[7]", "A[8]", "A[9]", "A[10]"), lags = 10)
mcmc_acf(posterior12, par = c("psi[1]", "psi[2]", "psi[3]", "psi[4]", "psi[5]"), lags = 10)
mcmc_acf(posterior12, par = c("psi[6]", "psi[7]", "psi[8]", "psi[9]", "psi[10]"), lags = 10)
mcmc_acf(posterior12, par = c("v[1,1]", "v[1,2]", "v[2,1]", "v[2,2]", "v[3,1]", "v[3,2]"), lags = 10)
mcmc_acf(posterior12, par = c("v[4,1]", "v[4,2]", "v[5,1]", "v[5,2]", "v[6,1]", "v[6,2]"), lags = 10)
mcmc_acf(posterior12, par = c("v[7,1]", "v[7,2]", "v[8,1]", "v[8,2]", "v[9,1]", "v[9,2]"), lags = 10)
mcmc_acf(posterior12, par = c("v[10,1]", "v[10,2]", "k_mu", "A_mu", "psi_mu", "v_mu[1]"), lags = 10)
mcmc_acf(posterior12, par = c("v_mu[2]", "k_sigma", "A_sigma", "psi_sigma", "v_sigma[1]", "v_sigma[2]","lp__"), lags = 10)
```

## 収束判定(条件による効果ありモデル){.tabset}
収束してそうですね。収束判定には，bayesplotを活用しています。[このページ](https://cran.rstudio.com/web/packages/bayesplot/vignettes/MCMC-diagnostics.html)が参考になりました。なお，定数のSは除外しています。

### トレースプロット
```{r, message=FALSE, warning=FALSE}
posterior21 <- rstan::extract(fit2, inc_warmup = TRUE, permuted = FALSE)
# k[1:5]
color_scheme_set("mix-blue-pink")
pk1 <- mcmc_trace(posterior21,  pars = c("k[1]", "k[2]", "k[3]", "k[4]", "k[5]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pk1 <- pk1+ facet_text(size = 15)
plot(pk1)
# k[6:10]
pk2 <- mcmc_trace(posterior21,  pars = c("k[6]", "k[7]", "k[8]", "k[9]", "k[10]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pk2 <- pk2+ facet_text(size = 15)
plot(pk2)
# A[1:5]
pA1 <- mcmc_trace(posterior21,  pars = c("A[1]", "A[2]", "A[3]", "A[4]", "A[5]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pA1 <- pA1+ facet_text(size = 15)
plot(pA1)
# A[6:10]
pA2 <- mcmc_trace(posterior21,  pars = c("A[1]", "A[2]", "A[3]", "A[4]", "A[5]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pA2 <- pA2+ facet_text(size = 15)
plot(pA2)
# tau[1:5]
pt <- mcmc_trace(posterior21,  pars = c("psi[1]", "psi[2]", "psi[3]", "psi[4]", "psi[5]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pt <- pt+ facet_text(size = 15)
plot(pt)
# tau[6:10]
pt2 <- mcmc_trace(posterior21,  pars = c("psi[6]", "psi[7]", "psi[8]", "psi[9]", "psi[10]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pt2 <- pt2+ facet_text(size = 15)
plot(pt2)
# v[1]
pv1 <- mcmc_trace(posterior21,  pars = c("v[1,1,1]", "v[1,1,2]", "v[1,2,1]", "v[1,2,2]", "v[1,3,1]", "v[1,3,2]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pv1 <- pv1+ facet_text(size = 15)
plot(pv1)
# v[2]
pv2 <- mcmc_trace(posterior21,  pars = c("v[2,1,1]", "v[2,1,2]", "v[2,2,1]", "v[2,2,2]", "v[2,3,1]", "v[2,3,2]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pv2 <- pv2+ facet_text(size = 15)
plot(pv2)
# v[3]
pv3 <- mcmc_trace(posterior21,  pars = c("v[3,1,1]", "v[3,1,2]", "v[3,2,1]", "v[3,2,2]", "v[3,3,1]", "v[3,3,2]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pv3 <- pv3+ facet_text(size = 15)
plot(pv3)
# v[4]
pv4 <- mcmc_trace(posterior21,  pars = c("v[4,1,1]", "v[4,1,2]", "v[4,2,1]", "v[4,2,2]", "v[4,3,1]", "v[4,3,2]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pv4 <- pv4+ facet_text(size = 15)
plot(pv4)
# v[5]
pv5 <- mcmc_trace(posterior21,  pars = c("v[5,1,1]", "v[5,1,2]", "v[5,2,1]", "v[5,2,2]", "v[5,3,1]", "v[5,3,2]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pv5 <- pv5+ facet_text(size = 15)
plot(pv5)
# v[6]
pv6 <- mcmc_trace(posterior21,  pars = c("v[6,1,1]", "v[6,1,2]", "v[6,2,1]", "v[6,2,2]", "v[6,3,1]", "v[6,3,2]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pv6 <- pv6+ facet_text(size = 15)
plot(pv6)
# v[7]
pv7 <- mcmc_trace(posterior21,  pars = c("v[7,1,1]", "v[7,1,2]", "v[7,2,1]", "v[7,2,2]", "v[7,3,1]", "v[7,3,2]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pv7 <- pv7+ facet_text(size = 15)
plot(pv7)
# v[8]
pv8 <- mcmc_trace(posterior21,  pars = c("v[8,1,1]", "v[8,1,2]", "v[8,2,1]", "v[8,2,2]", "v[8,3,1]", "v[8,3,2]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pv8 <- pv8+ facet_text(size = 15)
plot(pv8)
# v[9]
pv9 <- mcmc_trace(posterior21,  pars = c("v[9,1,1]", "v[9,1,2]", "v[9,2,1]", "v[9,2,2]", "v[9,3,1]", "v[9,3,2]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pv9 <- pv9+ facet_text(size = 15)
plot(pv9)
# v[10]
pv10 <- mcmc_trace(posterior21,  pars = c("v[10,1,1]", "v[10,1,2]", "v[10,2,1]", "v[10,2,2]", "v[10,3,1]", "v[10,3,2]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pv10 <- pv10+ facet_text(size = 15)
plot(pv10)
# hyperparameter 1
ph1 <- mcmc_trace(posterior21,  pars = c("k_mu","k_sigma", "A_mu","A_sigma","psi_mu","psi_sigma","lp__"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
ph1 <- ph1+ facet_text(size = 15)
plot(ph1)
# hyperparameter 2
ph2 <- mcmc_trace(posterior21,  pars = c("v_mu[1,1]", "v_mu[1,2]", "v_mu[2,1]", "v_mu[2,2]", "v_mu[3,1]", "v_mu[3,2]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
ph2 <- ph2+ facet_text(size = 15)
plot(ph2)
# hyperparameter 3
ph3 <- mcmc_trace(posterior21,  pars = c("v_sigma[1,1]", "v_sigma[1,2]", "v_sigma[2,1]", "v_sigma[2,2]", "v_sigma[3,1]", "v_mu[3,2]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
ph3 <- ph3+ facet_text(size = 15)
plot(ph3)
```

### $\hat{R}$の確認
bayesplotで確認をしても良いですが，ちょっとパラメータ多すぎるので，一括で検討します。$\hat{R}$が1.10以下ならTRUE，1.10以上のものがあればFALSEになります。
```{r, message=FALSE, warning=FALSE}
all(stan_rhat(fit2)$data <=1.10,na.rm=T)
```
    
### 有効サンプルサイズ
bayesplotで確認をしても良いですが，ちょっとパラメータ多すぎるので，一括で検討します。有効サンプルサイズが0.1以上ならTRUE，0.1以下のものがあればFALSEになります。
```{r, message=FALSE, warning=FALSE}
all(stan_ess(fit2)$data >= 0.1,na.rm=T)
```

### 自己相関
```{r, message=FALSE, warning=FALSE}
posterior22 <- as.matrix(fit2)
mcmc_acf(posterior22, par = c("k[1]", "k[2]", "k[3]", "k[4]", "k[5]"), lags = 10)
mcmc_acf(posterior22, par = c("k[6]", "k[7]", "k[8]", "k[9]", "k[10]"), lags = 10)
mcmc_acf(posterior22, par = c("A[1]","A[2]","A[3]","A[4]","A[5]"), lags = 10)
mcmc_acf(posterior22, par = c("A[6]","A[7]","A[8]","A[9]","A[10]"), lags = 10)
mcmc_acf(posterior22, par = c("psi[1]", "psi[2]", "psi[3]", "psi[4]", "psi[5]"), lags = 10)
mcmc_acf(posterior22, par = c("psi[6]", "psi[7]", "psi[8]", "psi[9]", "psi[10]"), lags = 10)
mcmc_acf(posterior22, par = c("v[1,1,1]", "v[1,1,2]", "v[1,2,1]", "v[1,2,2]", "v[1,3,1]", "v[1,3,2]"), lags = 10)
mcmc_acf(posterior22, par = c("v[2,1,1]", "v[2,1,2]", "v[2,2,1]", "v[2,2,2]", "v[2,3,1]", "v[2,3,2]"), lags = 10)
mcmc_acf(posterior22, par = c("v[3,1,1]", "v[3,1,2]", "v[3,2,1]", "v[3,2,2]", "v[3,3,1]", "v[3,3,2]"), lags = 10)
mcmc_acf(posterior22, par = c("v[4,1,1]", "v[4,1,2]", "v[4,2,1]", "v[4,2,2]", "v[4,3,1]", "v[4,3,2]"), lags = 10)
mcmc_acf(posterior22, par = c("v[5,1,1]", "v[5,1,2]", "v[5,2,1]", "v[5,2,2]", "v[5,3,1]", "v[5,3,2]"), lags = 10)
mcmc_acf(posterior22, par = c("v[6,1,1]", "v[6,1,2]", "v[6,2,1]", "v[6,2,2]", "v[6,3,1]", "v[6,3,2]"), lags = 10)
mcmc_acf(posterior22, par = c("v[7,1,1]", "v[7,1,2]", "v[7,2,1]", "v[7,2,2]", "v[7,3,1]", "v[7,3,2]"), lags = 10)
mcmc_acf(posterior22, par = c("v[8,1,1]", "v[8,1,2]", "v[8,2,1]", "v[8,2,2]", "v[8,3,1]", "v[8,3,2]"), lags = 10)
mcmc_acf(posterior22, par = c("v[9,1,1]", "v[9,1,2]", "v[9,2,1]", "v[9,2,2]", "v[9,3,1]", "v[9,3,2]"), lags = 10)
mcmc_acf(posterior22, par = c("v[10,1,1]", "v[10,1,2]", "v[10,2,1]", "v[10,2,2]", "v[10,3,1]", "v[10,3,2]"), lags = 10)
mcmc_acf(posterior22, par = c("k_mu", "k_sigma","A_mu","A_sigma","psi_mu","psi_sigma","lp__"), lags = 10)
mcmc_acf(posterior22, par = c("v_mu[1,1]", "v_mu[1,2]", "v_mu[2,1]", "v_mu[2,2]", "v_mu[3,1]", "v_mu[3,2]"), lags = 10)
mcmc_acf(posterior22, par = c("v_sigma[1,1]", "v_sigma[1,2]", "v_sigma[2,1]", "v_sigma[2,2]", "v_sigma[3,1]", "v_mu[3,2]"), lags = 10)
```

## 情報量基準(WAICとLOOCV)の算出
@berobero11氏の[Stan Advent 2016](http://statmodeling.hatenablog.com/entry/comparison-of-LOOCV-and-WAIC)を参考にしつつWAICとLOOCVの算出をしてみました。どちらの指標においても，条件の効果を考慮したモデルが当てはまりがう良さそうです。

```{r}
# WAICとLOOCVの計算
extractFit1 <- rstan::extract(fit1)
extractFit2 <- rstan::extract(fit2)
fitValue <- NULL
fitValue[1] <- loo::loo(extractFit1$log_lik)$looic/(2*numberAllTrials)
fitValue[2] <- loo::waic(extractFit1$log_lik)$waic/(2*numberAllTrials)
fitValue[3] <- loo::loo(extractFit2$log_lik)$looic/(2*numberAllTrials)
fitValue[4] <- loo::waic(extractFit2$log_lik)$waic/(2*numberAllTrials)

# プロットの準備
modelNames <- c("Without Condition Effect","Without Condition Effect","With Condition Effect","With Condition Effect")
valueNames <- c("LOOCV","WAIC","LOOCV","WAIC")
WAICData <- data.frame(modelNames,valueNames,fitValue)

# Bar plotの作成
bp <- ggplot(WAICData,aes(x=modelNames,y=fitValue,fill=valueNames)) + geom_bar(position="dodge",stat="identity")
bp <- bp + xlab("Model")+ylab("Values")+geom_text(aes(label = sprintf("%5.4f",fitValue)), vjust = 1.5, colour = "black", position = position_dodge(0.9),size = 5)+ylim(-1, NA)
plot(bp)
```


## 推定結果とパラメータリカバリのチェック

以降では，推定によって，元々設定したパラメータの値をリカバリできているのかチェックをしていきます（推定値の点推定値としてはMAPとEAPを算出しましたが，プロットはEAPのみにしています）。

### パラメータA
まあまあ，パラメターリカバリできているかな。

```{r, message=FALSE, warning=FALSE}
# MAP推定値計算関数を準備
mapEstimate <- function(z){
  density(z)$x[which.max(density(z)$y)]
}

fit1AEap <- apply(extractFit1$A,2,mean)
fit1AMap <- apply(extractFit1$A,2,mapEstimate)
fit2AEap <- apply(extractFit2$A,2,mean)
fit2AMap <- apply(extractFit2$A,2,mapEstimate)
parameterA <-data.frame(parameters$A,fit1AEap,fit1AMap,fit2AEap,fit2AMap)
print(parameterA)

Apfit1Eap <- ggplot(parameterA, aes(x=parameters.A, y=fit1AEap))+geom_point()+ xlab("Estimated parameter in Without Condition Effect Model")+ylab("Original value of A")
Apfit2Eap <- ggplot(parameterA, aes(x=parameters.A, y=fit2AEap))+geom_point()+ xlab("Estimated parameter in With Condition Effect Model")+ylab("Original value of A")
plot(Apfit1Eap)
plot(Apfit2Eap)
```

### パラメータk
まあまあ，パラメターリカバリできているかな。

```{r, message=FALSE, warning=FALSE}
# 元々はbからデータを作っていたので，kに変更(k=b-A)
originK <- parameters$b-parameters$A

fit1kEap <- apply(extractFit1$k,2,mean)
fit1kMap <- apply(extractFit1$k,2,mapEstimate)
fit2kEap <- apply(extractFit2$k,2,mean)
fit2kMap <- apply(extractFit2$k,2,mapEstimate)
parameterK <-data.frame(originK,fit1kEap,fit1kMap,fit2kEap,fit2kMap)
print(parameterK)

kpfit1Eap <- ggplot(parameterK, aes(x=originK, y=fit1kEap))+geom_point()+ xlab("Estimated parameter in Without Condition Effect Model")+ylab("Original value of k")
kpfit2Eap <- ggplot(parameterK, aes(x=originK, y=fit2kEap))+geom_point()+ xlab("Estimated parameter in With Condition Effect Model")+ylab("Original value of k")
plot(kpfit1Eap)
plot(kpfit2Eap)
```


### パラメータTau
まあまあ，パラメターリカバリできているかな。

```{r, message=FALSE, warning=FALSE}
fit1psiEap <- apply(extractFit1$psi,2,mean)
fit1psiMap <- apply(extractFit1$psi,2,mapEstimate)
fit2psiEap <- apply(extractFit2$psi,2,mean)
fit2psiMap <- apply(extractFit2$psi,2,mapEstimate)
parameterPsi <-data.frame(parameters$t0,fit1psiEap,fit1psiMap,fit2psiEap,fit2psiMap)
print(parameterPsi)

Psipfit1Eap <- ggplot(parameterPsi, aes(x=parameters.t0, y=fit1psiEap))+geom_point()+ xlab("Estimated parameter in Without Condition Effect Model")+ylab("Original value of Tau")
Psipfit2Eap <- ggplot(parameterPsi, aes(x=parameters.t0, y=fit2psiEap))+geom_point()+ xlab("Estimated parameter in With Condition Effect Model")+ylab("Original value of Tau")
plot(Psipfit1Eap)
plot(Psipfit2Eap)
```

### パラメータV(条件１の選択肢１)

当たり前だが，With Condition Effect Modelの方は，パラメータリカバリできている。

```{r, message=FALSE, warning=FALSE}
fit1v11Eap <- apply(extractFit1$v[,,1],2,mean)
fit1v11Map <- apply(extractFit1$v[,,1],2,mapEstimate)
fit2v11Eap <- apply(extractFit2$v[,,1,1],2,mean)
fit2v11Map <- apply(extractFit2$v[,,1,1],2,mapEstimate)
parameterV11 <-data.frame(parameters$v1_c1,fit1v11Eap,fit1v11Map,fit2v11Eap,fit2v11Map)
print(parameterV11)

V11pfit1Eap <- ggplot(parameterV11, aes(x=parameters.v1_c1, y=fit1v11Eap))+geom_point()+ xlab("Estimated parameter in Without Condition Effect Model")+ylab("Original value of v in condition 1 & choice1")
V11pfit2Eap <- ggplot(parameterV11, aes(x=parameters.v1_c1, y=fit2v11Eap))+geom_point()+ xlab("Estimated parameter in With Condition Effect Model")+ylab("Original value of v in condition 1 & choice1")
plot(V11pfit1Eap)
plot(V11pfit2Eap)
```

### パラメータV(条件１の選択肢2)

With Condition Effect Modelの方は，パラメータリカバリできている。

```{r, message=FALSE, warning=FALSE}
fit1v12Eap <- apply(extractFit1$v[,,2],2,mean)
fit1v12Map <- apply(extractFit1$v[,,2],2,mapEstimate)
fit2v12Eap <- apply(extractFit2$v[,,1,2],2,mean)
fit2v12Map <- apply(extractFit2$v[,,1,2],2,mapEstimate)
parameterV12 <-data.frame(parameters$v2_c1,fit1v12Eap,fit1v12Map,fit2v12Eap,fit2v12Map)
print(parameterV12)

V12pfit1Eap <- ggplot(parameterV12, aes(x=parameters.v2_c1, y=fit1v12Eap))+geom_point()+ xlab("Estimated parameter in With Condition Effect Model")+ylab("Original value of v in condition 1 & choice2")
V12pfit2Eap <- ggplot(parameterV12, aes(x=parameters.v2_c1, y=fit2v12Eap))+geom_point()+ xlab("Estimated parameter in Without Condition Effect Model")+ylab("Original value of v in condition 1 & choice2")
plot(V12pfit1Eap)
plot(V12pfit2Eap)
```

### パラメータV(条件２の選択肢１)

With Condition Effect Modelの方は，パラメータリカバリできている。

```{r, message=FALSE, warning=FALSE}
fit1v21Eap <- apply(extractFit1$v[,,1],2,mean)
fit1v21Map <- apply(extractFit1$v[,,1],2,mapEstimate)
fit2v21Eap <- apply(extractFit2$v[,,2,1],2,mean)
fit2v21Map <- apply(extractFit2$v[,,2,1],2,mapEstimate)
parameterV21 <-data.frame(parameters$v1_c2,fit1v21Eap,fit1v21Map,fit2v21Eap,fit2v21Map)
print(parameterV21)

V21pfit1Eap <- ggplot(parameterV21, aes(x=parameters.v1_c2, y=fit1v21Eap))+geom_point()+ xlab("Estimated parameter in Without Condition Effect Model")+ylab("Original value of v in condition 1 & choice1")
V21pfit2Eap <- ggplot(parameterV21, aes(x=parameters.v1_c2, y=fit2v21Eap))+geom_point()+ xlab("Estimated parameter in With Condition Effect Model")+ylab("Original value of v in condition 1 & choice1")
plot(V21pfit1Eap)
plot(V21pfit2Eap)
```

### パラメータV(条件２の選択肢2)

With Condition Effect Modelの方は，パラメータリカバリできている。

```{r, message=FALSE, warning=FALSE}
fit1v22Eap <- apply(extractFit1$v[,,2],2,mean)
fit1v22Map <- apply(extractFit1$v[,,2],2,mapEstimate)
fit2v22Eap <- apply(extractFit2$v[,,2,2],2,mean)
fit2v22Map <- apply(extractFit2$v[,,2,2],2,mapEstimate)
parameterV22 <-data.frame(parameters$v2_c2,fit1v22Eap,fit1v22Map,fit2v22Eap,fit2v22Map)
print(parameterV22)

V22pfit1Eap <- ggplot(parameterV22, aes(x=parameters.v2_c2, y=fit1v22Eap))+geom_point()+ xlab("Estimated parameter in Without Condition Effect Model")+ylab("Original value of v in condition 1 & choice2")
V22pfit2Eap <- ggplot(parameterV22, aes(x=parameters.v2_c2, y=fit2v22Eap))+geom_point()+ xlab("Estimated parameter in With Condition Effect Model")+ylab("Original value of v in condition 1 & choice2")
plot(V22pfit1Eap)
plot(V22pfit2Eap)
```

### パラメータV(条件３の選択肢１)

With Condition Effect Modelの方は，パラメータリカバリできている。

```{r, message=FALSE, warning=FALSE}
fit1v31Eap <- apply(extractFit1$v[,,1],2,mean)
fit1v31Map <- apply(extractFit1$v[,,1],2,mapEstimate)
fit2v31Eap <- apply(extractFit2$v[,,3,1],2,mean)
fit2v31Map <- apply(extractFit2$v[,,3,1],2,mapEstimate)
parameterV31 <-data.frame(parameters$v1_c3,fit1v31Eap,fit1v31Map,fit2v31Eap,fit2v31Map)
print(parameterV31)

V31pfit1Eap <- ggplot(parameterV31, aes(x=parameters.v1_c3, y=fit1v31Eap))+geom_point()+ xlab("Estimated parameter in Without Condition Effect Model")+ylab("Original value of v in condition 1 & choice1")
V31pfit2Eap <- ggplot(parameterV31, aes(x=parameters.v1_c3, y=fit2v31Eap))+geom_point()+ xlab("Estimated parameter in With Condition Effect Model")+ylab("Original value of v in condition 1 & choice1")
plot(V31pfit1Eap)
plot(V31pfit2Eap)
```

### パラメータV(条件３の選択肢2)

With Condition Effect Modelの方は，パラメータリカバリできている。

```{r, message=FALSE, warning=FALSE}
fit1v32Eap <- apply(extractFit1$v[,,2],2,mean)
fit1v32Map <- apply(extractFit1$v[,,2],2,mapEstimate)
fit2v32Eap <- apply(extractFit2$v[,,3,2],2,mean)
fit2v32Map <- apply(extractFit2$v[,,3,2],2,mapEstimate)
parameterV32 <-data.frame(parameters$v2_c3,fit1v32Eap,fit1v32Map,fit2v32Eap,fit2v32Map)
print(parameterV32)

V32pfit1Eap <- ggplot(parameterV32, aes(x=parameters.v2_c3, y=fit1v32Eap))+geom_point()+ xlab("Estimated parameter in Without Condition Effect Model")+ylab("Original value of v in condition 1 & choice2")
V32pfit2Eap <- ggplot(parameterV32, aes(x=parameters.v2_c3, y=fit2v32Eap))+geom_point()+ xlab("Estimated parameter in With Condition Effect Model")+ylab("Original value of v in condition 1 & choice2")
plot(V32pfit1Eap)
plot(V32pfit2Eap)
```

おおよそ，パラメータリカバリができていました。上記では，vのみが条件の効果を受けるという設定でしたが，適宜モデルを変えつつデータにあったモデルを設定できると良いかと思います。また，事後予測チェックはなかなか煩雑になるので省略しましたが，適宜行う必要があります。

**Enjoy Stan!**
