Experimetrics: 17.2 A Level-k Model for the Beauty Contest Game

今回は、実験経済学ではよく知られた美人投票ゲーム(Beauty Contest Game)の架空の被験者データにLevel-kモデルを当てはめ、推論の深さの程度ごとに被験者の割合を推定します。


美人投票ゲームの概要はこちら。

  • 各被験者は0から100までの整数を1つ選ぶ
  • 全被験者が選んだ数字の平均値×(2/3)に最も近い数字を選んだ被験者が勝者

繰り返し支配戦略を消していくことで、最終的には「各被験者が0を選ぶ」というこのゲームのナッシュ均衡にたどり着きます。最初の実験であるNagel (1995)では平均値が35、また多くの被験者が33や22といった数字を選びました。


次にLevel-kモデルの概要です。

  • 人々は様々なレベルの推論を行うと仮定する
  • Level-0(0次)の推論を行うプレイヤー:0から100までの数字を一様にランダムに選択するプレイヤー
  • Level-1(1次)の推論を行うプレイヤー:他のプレイヤーは全てLevel-0の推論を行うとの信念を持ち(つまり、自分以外のプレイヤーが選ぶ数字の平均値は50付近であるとの信念を持つ)、その結果33を最適な推論として選択するプレイヤー
  • Level-2(2次)の推論を行うプレイヤー:他のプレイヤーは全てLevel-1の推論を行うとの信念を持ち(つまり、自分以外のプレイヤーが選ぶ数字の平均値は33付近であるとの信念を持つ)、その結果22を最適な推論として選択するプレイヤー
  • Level-3(3次)の推論を行うプレイヤー:他のプレイヤーは全てLevel-2の推論を行うとの信念を持ち(つまり、自分以外のプレイヤーが選ぶ数字の平均値は22付近であるとの信念を持つ)、その結果15を最適な推論として選択するプレイヤー
  • Level-4(4次)の推論を行うプレイヤー:他のプレイヤーは全てLevel-3の推論を行うとの信念を持ち(つまり、自分以外のプレイヤーが選ぶ数字の平均値は15付近であるとの信念を持つ)、その結果10を最適な推論として選択するプレイヤー
  • 同じ要領でより高次の推論を行うプレイヤーを想定することができるが、今回はLevel-5を最高次とし、他のプレイヤーは全て数字0を選ぶとするナイーブなNashタイプのプレイヤーとする

このように、Level-kモデルによれば、推論がより高次のプレイヤーほど0に近い数字を選択することになるわけです。

下のグラフは、架空データ(n=500)の分布になります。このデータからは33や0、また22付近の数字や、15付近の数字の選択頻度が高いことがわかります。
f:id:hiroecon:20180610092728p:plain



この節では、全ての被験者はLevel-kモデルに基づき数字を選択しているものと仮定し、各レベルの被験者の割合を推定します。推定に際し、次のように仮定します。

  • Level-0の推論を行うプレイヤーは0から100の数字を等確率で選ぶ
  • Level-j(j>0)の推論を行うプレイヤーについては、選択する数字 yは各レベルの最適な推論の数字 y_{j}^{\ast}に平均0、分散 \sigma^{2}正規分布に従う誤差項を加えた数字とする

 y|\text{level}\;j=y_{j}^{\ast}+\epsilon,\;\;\;\;\epsilon\sim N(0,\sigma^{2})

なお、 y_{1}^{\ast}=33 y_{2}^{\ast}=22 y_{3}^{\ast}=15 y_{4}^{\ast}=10 y_{5}^{\ast}=0となります。


これらの仮定によって、各レベルの条件付き確率密度関数
 f(y|\text{level 0})=\frac{1}{100},\;\;0\leq y\leq 100
 f(y|\text{level}\;j)=\frac{1}{\sigma}\phi \Bigl(\frac{y-y_{j}^{\ast}}{\sigma}\Bigr),\;\;0\leq y\leq 100,\;\;j=1,2,3,4,5
となります( \phi(\cdot)正規分布確率密度関数)。


架空データの被験者500名は、Level-0からLevel-5までのプレイヤーで構成されているとし、 p_{j}をLevel-jの被験者の割合とします。これにより、対数尤度関数は次のような形になります。
\displaystyle LogL=\Sigma_{i=1}^{500}\ln \Bigl\{\frac{p_{0}}{100}+\Sigma_{j=1}^{5}p_{j}\frac{1}{\sigma}\phi \Bigl(\frac{y_{i}-y_{j}^{\ast}}{\sigma}\Bigr) \Bigr\}


Rプログラムは以下のようになります。

# read STATA dta file
library(foreign)
d<-read.dta("/Users/Hiro/Documents/R/moffatt/data_files/beauty_sim.dta")
y<-d$y
summary(y)
hist(y,
     breaks=seq(0,100,1),
     xlab="Guess",
     main="Distribution of Simulated Guesses")

negLogL<-function(p1,p2,p3,p4,p5,sigma){
  LogL<-sum(log((1-p1-p2-p3-p4-p5)*(1/100)
                +p1*(1/sigma)*dnorm((y-33)/sigma)
                +p2*(1/sigma)*dnorm((y-22)/sigma)
                +p3*(1/sigma)*dnorm((y-15)/sigma)
                +p4*(1/sigma)*dnorm((y-10)/sigma)
                +p5*(1/sigma)*dnorm((y-0)/sigma)))
  negLogL<--LogL
}

library(stats4)
fit<-mle(negLogL, list(p1=0.3,p2=0.4,p3=0.1,p4=0.1,p5=0.05,sigma=2))
summary(fit)

coef(fit) # estimated coefficients
vcov(fit) # variance-covariance matrix

p1_hat<-coef(fit)[1]
p2_hat<-coef(fit)[2]
p3_hat<-coef(fit)[3]
p4_hat<-coef(fit)[4]
p5_hat<-coef(fit)[5]
sigma_hat<-coef(fit)[6]

p0_hat<-1-p1_hat-p2_hat-p3_hat-p4_hat-p5_hat
p0_hat

# delta method
library(msm)
p0_std.error<-deltamethod(~1-x1-x2-x3-x4-x5,coef(fit),vcov(fit)) # The variables must be labelled x1, x2,...
p0_std.error

# posterior level probabilities
posterior<-function(x) {
  denom<-p0_hat*(1/100)+p1_hat*(1/sigma_hat)*dnorm((x-33)/sigma_hat)+p2_hat*(1/sigma_hat)*dnorm((x-22)/sigma_hat)+p3_hat*(1/sigma_hat)*dnorm((x-15)/sigma_hat)+p4_hat*(1/sigma_hat)*dnorm((x-10)/sigma_hat)+p5_hat*(1/sigma_hat)*dnorm((x-0)/sigma_hat)
  post_0<-p0_hat*(1/100)/denom
  post_1<-p1_hat*(1/sigma_hat)*dnorm((x-33)/sigma_hat)/denom
  post_2<-p2_hat*(1/sigma_hat)*dnorm((x-22)/sigma_hat)/denom
  post_3<-p3_hat*(1/sigma_hat)*dnorm((x-15)/sigma_hat)/denom
  post_4<-p4_hat*(1/sigma_hat)*dnorm((x-10)/sigma_hat)/denom
  post_5<-p5_hat*(1/sigma_hat)*dnorm((x-0)/sigma_hat)/denom
  return(list(post_0,post_1,post_2,post_3,post_4,post_5))
}

guess<-c(1:100)
posterior_type0<-posterior(guess)[[1]]
posterior_type1<-posterior(guess)[[2]]
posterior_type2<-posterior(guess)[[3]]
posterior_type3<-posterior(guess)[[4]]
posterior_type4<-posterior(guess)[[5]]
posterior_type5<-posterior(guess)[[6]]

plot(guess,posterior_type0,type='l',ann=F,xlim=c(0,100),ylim=c(0,1),lty=1)
par(new=T)
plot(guess,posterior_type1,type='l',ann=F,xlim=c(0,100),ylim=c(0,1),lty=2)
par(new=T)
plot(guess,posterior_type2,type='l',ann=F,xlim=c(0,100),ylim=c(0,1),lty=3)
par(new=T)
plot(guess,posterior_type3,type='l',ann=F,xlim=c(0,100),ylim=c(0,1),lty=4)
par(new=T)
plot(guess,posterior_type4,type='l',ann=F,xlim=c(0,100),ylim=c(0,1),lty=5)
par(new=T)
plot(guess,posterior_type5,type='l',xlab="Guess",ylab="Posterior Probability",xlim=c(0,100),ylim=c(0,1),lty=6)
legend_names<-c("posterior type 0","posterior type 1","posterior type 2","posterior type 3","posterior type 4","posterior type 5")
legend("topright",legend=legend_names,lty=c(1:6))


推定結果はこちらになります。

Coefficients:
        Estimate Std. Error
p1    0.39762922 0.02388147
p2    0.11127580 0.01644880
p3    0.08911294 0.01594946
p4    0.04588258 0.01357794
p5    0.04982359 0.01180616
sigma 1.96095481 0.10848212

-2 log L: 3981.022 

なお、Level-0の割合の推定値 \;\hat{p}_{0}は0.3062759、デルタ法でもとめた標準誤差は0.02925361でした。Experimetrics本の411ページの結果とほぼ同じ推定結果が得られました。


最後に、推定値とベイズの定理を用いて、各レベルの事後分布をもとめてみます。例えば、数字 yを観測した条件のもとでの各レベルの事後分布は次のように求めることができます。
\displaystyle f(\text{level}\;0|y)=\frac{  \frac{\hat{p}_{0}}{100}   }{   \frac{\hat{p}_{0}}{100} + \Sigma_{k=1}^{5}\;\frac{1}{\hat{\sigma}}\phi \Bigl(\frac{y_{i}-y_{j}^{\ast}}{\hat{\sigma}}\Bigr)\hat{p}_{k}   }

\displaystyle f(\text{level}\;j|y)=\frac{   \frac{1}{\hat{\sigma}}\phi \Bigl(\frac{y_{i}-y_{j}^{\ast}}{\hat{\sigma}}\Bigr)\hat{p}_{j}   }{   \frac{\hat{p}_{0}}{100} + \Sigma_{k=1}^{5}\;\frac{1}{\hat{\sigma}}\phi \Bigl(\frac{y_{i}-y_{j}^{\ast}}{\hat{\sigma}}\Bigr)\hat{p}_{k}   }\;\;j=1,2,3,4,5


各レベルの事後分布のグラフです。Experimentrics本の412ページのFigure 17.2を再現できました。もし yが40以上ならその被験者はほぼ100%の確率でLevel-0のプレイヤーであることを意味します。他のレベルの事後分布は一つのピークしか無く、レベルが上がるほどピークの yの位置が0に近くなっていくのが特徴です。
f:id:hiroecon:20180610114224p:plain