ベイズマックスへようこそ。
統計に関するあれやこれやを自由気ままに投稿します。
全記事一覧はこちらから。

2020年09月13日

Drift Diffusion Modelのパラメータから選択確率を計算する関数:RとStanによる実装

 Drift Diffusion Modelのパラメータから選択確率(一方の反応が生じる確率;choice probability)を計算する方法について調べたので,備忘録として,かつ,日本語で書かれた資料がWeb上にはなさそうなので,記事にすることにしました。

 Drift Diffusion Modelは反応時間分布を説明するときによく使われるモデルの1つです。2つのうちいずれの反応を実行するかについての情報(証拠)が時間とともに蓄積され,情報量が閾値を越えたときに反応が生じることを仮定します。反応時間分布と書きましたが,正確にいえば<いずれの反応が生じるか>と<反応時間>の同時分布です。このうち,<いずれの反応が生じるか>だけに注目すると,これは2値変数なので,その周辺分布は選択確率Pupperをパラメータとするベルヌーイ分布に従います。この選択確率Pupperを,Drift Diffusion Modelのパラメータから計算する方法についてこの記事で説明します。といっても導出は難しそうなので,計算式とRとStanで実装する方法だけを説明します。これを知っていれば,例えばStanで選択確率の事後予測分布を生成したいときなどに役立つと思います。なお,Drift Diffusion Modelについては初学者向けの優れた解説が既にあるので(e.g.,【RとIATとBayesian hierarchical diffusion model】拡散過程モデルとは?(西村友佳さんによる解説), Stanで階層diffusionモデル(小林穂波さんによる解説)),ここでは改めて解説はしません。モデル式は原著論文 (Ratcliff, 1978) に詳しく書かれていますが,実用上はWiener過程の分散パラメータを1に固定したり (RWienerパッケージやStanでは1に固定している;see Wabersich & Vandekerckhove, 2014),より効率的な計算法を用いたり(e.g., Blurton, Kesselmeier, & Gondan, 2012; Navarro & Fuss, 2009)します。

 Drift Diffusion Modelでは,情報が上側(正の方向)の閾値に到達するか下側(負の方向)の閾値に到達するかのいずれかが生じたときに反応が生じると考えるので,上側選択確率Pupperと下側選択確率Plowerの両方を考えることができます。もちろんPupper + Plower = 1なので,どちらか一方が分かれば十分です。Wiener過程の分散パラメータを1に固定すると,これらは以下の式で表せます。

ddf_p_formula.png

この式はWabersich & Vandekerckhove (2014)のAppendixに載っているものです。ここで,αは閾値を決めるパラメータ(≒慎重さ),βは初期値を決めるパラメータ(≒反応バイアス),δはドリフト率パラメータ(情報集積の速さ)を表します。非決定時間パラメータのτは選択率には影響しないのでこの式には含まれていません。αが大きいほどδの符号と同じ方向の選択確率が高くなること,βが大きいほど上側選択確率が高くなること,δが大きいほどその符号の方向の選択確率が高くなることがこの式から読み取れます。直感にも合いますね。なお,δ=0のときは元の式が0/0の不定形になってしまうので,そうならないようにifで場合分けされています。これらの式を使って選択確率をRで計算するための関数は以下の通りです。

選択確率を計算するR関数
## 関数の定義 ------
wiener_cp <- function(alpha,beta,delta,resp="upper",tau=NULL){
  if(resp!="upper" & resp!= "lower"){
    resp <- "upper"
    print("Warining: resp must be either 'upper' or 'lower.' The probability at the upper boundary was calculated.")
  }
  if(delta == 0){
    p <- 1-beta
  } else{
    tmp <- -2*delta*alpha*(1-beta)
    p <- (1-exp(tmp)) / (exp(2*delta*alpha*beta) - exp(tmp))
  }
  
  if(resp=="upper"){
    p <- 1-p
  }
  return(p)
}

## 関数を使ってみる ------
alpha <- 1
beta <- 0.6
tau <- 0.5
delta <- 1

wiener_cp(alpha,beta,delta, resp="upper") #[1] 0.8081812
wiener_cp(alpha,beta,delta, resp="lower") #[1] 0.1918188

## 合っているかを確認する ------
library(RWiener)
pwiener(100000,alpha,tau,beta,delta,resp="upper") #[1] 0.8081812
pwiener(100000,alpha,tau,beta,delta,resp="lower") #[1] 0.1918188

 Stanで選択確率の事後予測分布を生成するには次のように書けば良いです。関数は,使うことの多い上側選択確率を返すように設定しています。

選択確率の事後予測分布を生成するStanコード
functions{
  //上側選択確率を返す関数
  real wiener_cp_upper(real alpha, real beta, real delta){
    real p_lower;
    if(delta == 0){
      p_lower = 1 - beta;
    } else{
      real tmp = -2*delta*alpha*(1-beta);
      p_lower = (1-exp(tmp)) / (exp(2*delta*alpha*beta) - exp(tmp));
    }
    return(1-p_lower);
  }
}

data{
  int N; //データの数
  real<lower=0> RT[N]; //反応時間
  int resp[N]; //反応:0(下側) or 1(上側)
}

transformed data{
  real minRT = min(RT); //tauの上限
}

parameters{
  real<lower=0> alpha;
  rea<lower=0, upper=minRT> tau;
  real<lower=0,upper=1> beta;
  real delta;
}

model{
  for(i in 1:N){
    if(resp[i] == 1){ //上側反応のとき
      RT[i] ~ wiener(alpha,tau,beta,delta);
    } else{ //下側反応のとき
      RT[i] ~ wiener(alpha,tau,1-beta,-delta);
    }
  }
}

generated quantities{
  real p_upper = wiener_cp_upper(alpha,beta,delta); //上側選択確率
  real p_lower = 1- p_upper; //下側選択確率
  int resp_pred = bernoulli_rng(p_upper); //選択確率の事後予測分布
}

 以上です。ちょっとしたメモのつもりだったのになんやかんや書くのに2時間ぐらいかかってしまいました。Enjoy!

posted by mutopsy at 17:07 | 統計