認知情報解析 カジノ体験

サイコロを3つなげ、その和が何であるかをかけるギャンブルがあるそうです。
和が、4~10ならsmall, 11~17ならlargeとおおよそ2択のギャンブルだそうです。(3をx、18をyと呼ぶことにします)
Kくんが体験したのは「過去50回の結果がlargeが70%と表示されていたので、次はsmallじゃないかと思いsmallに賭けて、実際にsmallが出た」ということでした。この出来事がどの程度よく起きるのかシミュレーションを用いて検証しまししょう。

# 非効率だけど、直感的に分かりやすいと思われる例

# 1試行の関数
die <- function() {
  die1<-sample(1:6,1)
  die2<-sample(1:6,1)
  die3<-sample(1:6,1)
  die.sum = die1 + die2 + die3
  if (die.sum==3) {
     result="X"
  } else { 
    if (die.sum==18) {
      result="Y"
   } else {
     if (die.sum<11){ 
       result="small"
       } else {
         result="large"
       }
     }
   }
  return(result)
}

# 51試行の結果を表す関数
SL51=function( ) {
  n.rep=51
  res=rep("N",51)
  for (i_rep in 1:n.rep) {
    res[i_rep]=die()
  }
  p.large=sum(res[1:50]=="large")/50
  return(list(p.large=p.large,last.result=res[51]))
}
  
# 51試行をN.REP回繰り返す関数
ck.kurimoto <- function(n.rep) {
  p.hist = rep(0,n.rep)
  lr.hist = rep("N",n.rep)
  for (i_rep in 1:n.rep) {
     result=SL51()
     p.hist[i_rep] = result$p.large
     lr.hist[i_rep] = result$last.result
  }
  return(list(p.hist=p.hist,lr.hist=lr.hist))
}

# 50試行の内、Lが出た割合が70%以上のものを選択
index70=which(m.result$p.hist>=0.7)
sum(m.result$lr.hist[index70]=="small")/length(index70)

# 上の3つの関数をまとめたもの
n.rep=1e6
die.mat = matrix(sample(c("X","Y","S","L"),
 n.rep*51,prob=c(1/16,1/16,7/16,7/16),replace=T),nrow=n.rep)
p.large=rowSums(die.mat[,1:50]=="L")/50
index70=which(p.large>=0.7)
ck51L70=die.mat[index70,51]
table(ck51L70)/length(index70)

認知情報解析演習 ギャンブル

# 直感的にはわかりやすいけど、非効率な例
bet.example<-function(n.trial,p.win,max.bet) {
  rew.history=rep(0,n.trial)
  bet.history=rep(0,n.trial)
  WL.history=rep("W",n.trial)
  loss.counter=0
  i_trial=0;
  WL="L"
  loss.counter.reset="F"
  while (i_trial < n.trial | WL=="L" ) {
    i_trial=i_trial+1
    bet=2^loss.counter
    if (bet > max.bet) {
      bet = max.bet
      loss.counter.reset = "T"
    }
    bet.history[i_trial]=bet
    WL=sample(c("W","L"),1,prob=c(p.win, 1-p.win))
    WL.history[i_trial]=WL
    if (WL=="L") {
      loss.counter=loss.counter+1
      reward=-bet
    } else {
      loss.counter=0
      reward=bet
    }
    if (loss.counter.reset=="T"){loss.counter=0}
    rew.history[i_trial]=reward
  }
  return(list(reward=sum(rew.history),money.required=max(bet.history)))
}

# より効率的な例 - max.betは未導入
bet.example02<-function(n.trial,p.win) {
  WL=sample(c("W","L"),n.trial,replace=T,prob=c(p.win, 1-p.win))
  if (WL[n.trial]!="W") {
    WL.extra=sample(c("W","L"),1,prob=c(p.win,1-p.win))
    while (tail(WL.extra,1)=="L") {
      WL.extra=c(WL.extra,sample(c("W","L"),1,prob=c(p.win,1-p.win)))
     }
    WL=c(WL,WL.extra)
  }
  reward=sum(WL=="W")
  result=rle(WL)
  loss.index=result$values=="L"
  money.required=2^max(result$lengths[loss.index])
  return(list(reward,money.required))
}

データ解析基礎論 分散分析の例

# 可視化
interaction.plot(dat$gender,
  dat$affil,
  dat$shoesize, 
  pch=c(20,20), 
  col=c("skyblue","orange"), 
  xlab="gender", ylab="shoesize", 
  lwd=3,type='b',cex=2,
  trace.label="Affiliation")

# 可視化2(非効率)
means<-tapply(dat$shoesize,list(dat$gender, dat$affil),mean)
Ns<-tapply(dat$shoesize,list(dat$gender, dat$affil),length)
sds<-tapply(dat$shoesize,list(dat$gender, dat$affil),sd)
sems<-sds/sqrt(Ns)

plot(c(0,1),means[,1],type='o',col='skyblue',
 ylim=c(min(means)*0.975,max(means)*1.025),
  xlim=c(-0.1,1.1),lwd=2,cex=2,pch=20,xlab="gender",ylab="shoesize")
lines(c(0,1),means[,2],type='o',col='orange',lwd=2,cex=2,pch=20)
legend("topleft",c("CogSci","PsySci"),col=c("skyblue","orange"),lwd=2)
lines(c(0,0),c(means[1,1]-sems[1,1],means[1,1]+sems[1,1]),col="skyblue",lwd=2.5)
lines(c(1,1),c(means[2,1]-sems[2,1],means[2,1]+sems[2,1]),col="skyblue",lwd=2.5)
lines(c(0,0),c(means[1,2]-sems[1,2],means[1,2]+sems[1,2]),col="orange",lwd=2.5)
lines(c(1,1),c(means[2,2]-sems[2,2],means[2,2]+sems[2,2]),col="orange",lwd=2.5)

# 分散分析
dat.aov=aov(shoesize~gender*affil, data=dat)
dat.aov.sum=summary(dat.aov)

# 単純主効果の検定
#  各所属講座における性別の効果
means<-tapply(dat$shoesize, list(dat$gender,dat$affil), mean)
SS_gen_CS<- 5*(means[2,1]^2 + means[1,1]^2 -0.5*sum(means[,1])^2) # SS_gender CS
SS_gen_PS<- 5*(means[2,2]^2 + means[1,2]^2 -0.5*sum(means[,2])^2) # SS_gender PS
dat.aov.sum=summary(dat.aov)   # ANOVA table
MSe=dat.aov.sum[[1]][4,3]      # MSE from ANOVA table or MSe=0.62
dfE=dat.aov.sum[[1]][4,1]      # DF for error or dfE=16
dfG=1                          # DF for gender
F_gen_CS=(SS_gen_CS/dfG)/MSe   # F-value for gender effect given CS
F_gen_PS=(SS_gen_PS/dfG)/MSe   # F-value for gender effect given PS
P_gen_CS=1-pf(F_gen_CS,1,dfE)  # p-value for gender effect given CS
P_gen_PS=1-pf(F_gen_PS,1,dfE)  # p-value for gender effect given PS

#  各性別における所属講座の効果
SS_affil_F<- 5*(means[1,1]^2+means[1,2]^2-0.5*sum(means[1,])^2) #SS_affil | F
SS_affil_M<- 5*(means[2,1]^2+means[2,2]^2-0.5*sum(means[2,])^2) #SS_affil | M
dfA=1				          # DF for affil
F_affil_F=SS_affil_F/dfA/MSe         # F-value for affiliation effect | F
F_affil_M=SS_affil_M/dfA/MSe         # F-value for affiliation effect | M
P_affil_F=1-pf(F_affil_F,1,dfE)      # p-value for affiliation effect | F
P_affil_M=1-pf(F_affil_M,1,dfE)      # p-value for affiliation effect | M

# 例2
interaction.plot(dat$duration,dat$method,dat$result, 
  pch=c(20,20), col=c("skyblue","orange"), ylab="score", 
  xlab="Duration",lwd=3,type='b',cex=2,trace.label="Method")
mod1=aov(result~method+duration,data=dat)
mod1.sum=print(summary(mod1))
mod2=aov(result~method*duration,data=dat)
mod2.sum=print(summary(mod2))
means<-tapply(dat$result,list(dat$method,dat$duration),mean)
ssM_1=5*(sum(means[,1]^2)-0.5*(sum(means[,1])^2))
ssM_2=5*(sum(means[,2]^2)-0.5*(sum(means[,2])^2))
ssM_3=5*(sum(means[,3]^2)-0.5*(sum(means[,3])^2))
ssM_4=5*(sum(means[,4]^2)-0.5*(sum(means[,4])^2))
MSe=mod2.sum[[1]][4,3]
DFe=mod2.sum[[1]][4,1]
DFm=1
fM_1=(ssM_1/DFm)/MSe
1-pf(fM_1,DFm,DFe)
fM_2=(ssM_2/DFm)/MSe
1-pf(fM_2,DFm,DFe)
fM_3=(ssM_3/DFm)/MSe
1-pf(fM_3,DFm,DFe)
fM_4=(ssM_4/DFm)/MSe
1-pf(fM_4,DFm,DFe)
ssD_X=5*(sum(means[1,]^2)-1/4*(sum(means[1,])^2))
ssD_Y=5*(sum(means[2,]^2)-1/4*(sum(means[2,])^2))
DFd=3
fD_X=(ssD_X/DFd)/MSe
fD_Y=(ssD_Y/DFd)/MSe
1-pf(fD_X,DFd,DFe)
1-pf(fD_X,DFd,DFe)
qv=qtukey(0.95,DFd+1,DFe)
hsd=qv*(sqrt(MSe/5))
print(diffM<-outer(means[1,],means[1,],"-"))           
abs(diffM)>hsd