決断ヒストリーのシミュレーション

• 約１か月前、決断の話を始めた(こちら)
• ●は赤を出して勝ち、○は赤を出して負け、■は青を出して勝ち、□は青を出して負けとした

●●○○■○■■■□●●□●○○■■□□□●

• こんな●■○□の並びが、「決断過程」として妥当な感じがするかどうか、とか、妥当って何か、とかそんなところがそもそもの発端だった
• その後、
  • 「決断」とは何か、とか
  • さらに、「決断」とは個人のものか「集団」のものか、と言う議論(の端緒)や
  • 数式・数学を使わない、生物学的な「決断」機構との関連についてや
  • 決断条件の複雑化の取り扱いなど
• に話は飛んだり戻ったりし、
• データ構造や解析の道具としては
  • 共役事前分布やベータ分布・ディリクレ分布の復習があり、
  • その関連で正則ベータ関数が登場し
  • 複合ベータ分布とその積分のことなども扱い
  • 多次元格子や整数分割と、そのオブジェクト扱いのこと
• が登場した
• それらを使って、振り出しに戻って、『それらしい』●■○□の並び(「決断過程」)をシミュレーションで作ってみた
  • ●○は0.5:0.5、■□は0.6:0.4にしてみた

●■●●●■■□●○□■□●●●○●○●□○■○●●○●●●●●○□□●●●○○○●○○□○○○●●○●○○●○●○■○■■○□■○□■●●●●■○○○■□■○□○○■□●□■■■○■■□□□□■□

• この全100選択の、つどつどの２ｘ２分割表の４セルの値は以下の通り

> X.bd
       [,1] [,2] [,3] [,4]
  [1,]    0    0    0    0
  [2,]    1    0    0    0
  [3,]    1    0    1    0
  [4,]    2    0    1    0
  [5,]    3    0    1    0
  [6,]    4    0    1    0
  [7,]    4    0    2    0
  [8,]    4    0    3    0
  [9,]    4    0    3    1
 [10,]    5    0    3    1
 [11,]    5    1    3    1
 [12,]    5    1    3    2
 [13,]    5    1    4    2
 [14,]    5    1    4    3
 [15,]    6    1    4    3
 [16,]    7    1    4    3
 [17,]    8    1    4    3
 [18,]    8    2    4    3
 [19,]    9    2    4    3
 [20,]    9    3    4    3
 [21,]   10    3    4    3
 [22,]   10    3    4    4
 [23,]   10    4    4    4
 [24,]   10    4    5    4
 [25,]   10    5    5    4
 [26,]   11    5    5    4
 [27,]   12    5    5    4
 [28,]   12    6    5    4
 [29,]   13    6    5    4
 [30,]   14    6    5    4
 [31,]   15    6    5    4
 [32,]   16    6    5    4
 [33,]   17    6    5    4
 [34,]   17    7    5    4
 [35,]   17    7    5    5
 [36,]   17    7    5    6
 [37,]   18    7    5    6
 [38,]   19    7    5    6
 [39,]   20    7    5    6
 [40,]   20    8    5    6
 [41,]   20    9    5    6
 [42,]   20   10    5    6
 [43,]   21   10    5    6
 [44,]   21   11    5    6
 [45,]   21   12    5    6
 [46,]   21   12    5    7
 [47,]   21   13    5    7
 [48,]   21   14    5    7
 [49,]   21   15    5    7
 [50,]   22   15    5    7
 [51,]   23   15    5    7
 [52,]   23   16    5    7
 [53,]   24   16    5    7
 [54,]   24   17    5    7
 [55,]   24   18    5    7
 [56,]   25   18    5    7
 [57,]   25   19    5    7
 [58,]   26   19    5    7
 [59,]   26   20    5    7
 [60,]   26   20    6    7
 [61,]   26   21    6    7
 [62,]   26   21    7    7
 [63,]   26   21    8    7
 [64,]   26   22    8    7
 [65,]   26   22    8    8
 [66,]   26   22    9    8
 [67,]   26   23    9    8
 [68,]   26   23    9    9
 [69,]   26   23   10    9
 [70,]   27   23   10    9
 [71,]   28   23   10    9
 [72,]   29   23   10    9
 [73,]   30   23   10    9
 [74,]   30   23   11    9
 [75,]   30   24   11    9
 [76,]   30   25   11    9
 [77,]   30   26   11    9
 [78,]   30   26   12    9
 [79,]   30   26   12   10
 [80,]   30   26   13   10
 [81,]   30   27   13   10
 [82,]   30   27   13   11
 [83,]   30   28   13   11
 [84,]   30   29   13   11
 [85,]   30   29   14   11
 [86,]   30   29   14   12
 [87,]   31   29   14   12
 [88,]   31   29   14   13
 [89,]   31   29   15   13
 [90,]   31   29   16   13
 [91,]   31   29   17   13
 [92,]   31   30   17   13
 [93,]   31   30   18   13
 [94,]   31   30   19   13
 [95,]   31   30   19   14
 [96,]   31   30   19   15
 [97,]   31   30   19   16
 [98,]   31   30   19   17
 [99,]   31   30   20   17
[100,]   31   30   20   18

• ソース

Decision_Dice_beta.2.3 <- function(X,a){

	tmp <- Decision_beta.2(X+1)
	tmp2 <- c(tmp,1-tmp,0)
	if(min(tmp2) > a){
		tmp2 <- c(0.5,0.5,1)
	}
	return(tmp2)
}
Decision_beta.2 <- function(x){
	if(x[1] < x[4]){
		x <- x[4:1]
	}
	ret <- 0
	common <- -log(x[3]+x[4])+lgamma(x[1]+x[2])+lgamma(x[3]+x[4]+1)-lgamma(x[1])-lgamma(x[2])-lgamma(sum(x)-1)
	for(j in x[3]:(x[3]+x[4]-1)){
		tmp <- lgamma(x[1]+j)+lgamma(sum(x[2:4])-1-j)-lgamma(j+1)-lgamma(x[3]+x[4]-j)
		ret <- ret + exp(tmp+common)
	}
	return(ret)
}
my.Chisq.Fisher.test <- function(X,a,Fisher=FALSE){
	X.m <- matrix(X,2,2)
	rows <- apply(X.m,1,sum)
	cols <- apply(X.m,2,sum)
	if(prod(c(rows,cols))==0){
		return(c(0,0,1))
	}
	if(Fisher){
		tmp <- fisher.test(X.m)
	}else{
		tmp <- chisq.test(X.m,correct=FALSE)
	}
	
	tmp.p <- tmp$p.value
	rs <- X.m[,1]/rows
	tmp.2 <- c(1,0,0)
	if(rs[1] < rs[2]){
		tmp.2 <- c(0,1,0)
	}
	if(tmp.p>a){
		tmp.2 <- c(0,0,1)
	}
	return(tmp.2)
}
my.Prob <- function(v,prob=NULL){
	if(is.null(prob)){
		prob=rep(1/length(v),length(v))
	}
	s <- sum(v)
	exp(lgamma(s+1)-sum(lgamma(v+1))+sum(v*log(prob)))
}

N.decision <- 20
N <- 100

N.iter <- 2

X0 <- rep(0,4)

ps <- c(0.5,0.6)
a <- 0.05
Fisher <- FALSE
sv.bd <- matrix(0,N.iter,N)
sv.cd <- sv.bd
for(i in 1:N.iter){
	selected <- 0

	X.bd <- matrix(0,N,4)
	X.cd <- X.bd
	for(j in 2:N){
		X.bd[j,] <- X.bd[j-1,]
		X.cd[j,] <- X.cd[j-1,]
		tmp.bd <- Decision_Dice_beta.2.3(X.bd[j-1,],2)
		s.bd <- sample(1:2,1,prob=c(tmp.bd[1],1-tmp.bd[1]))
		h.bd <- sample(1:2,1,prob=c(ps[s.bd],1-ps[s.bd]))
		X.bd[j,(s.bd-1)*2+h.bd] <- X.bd[j,(s.bd-1)*2+h.bd] + 1

		if(j == N.decision){
			#print("yes")
			tmp.cd <- my.Chisq.Fisher.test(X.cd[j-1,],a,Fisher)
			s.cd <- which(tmp.cd==1)
			if(s.cd == 3){
				s.cd <- sample(1:2,1)
			}else{
				selected <- s.cd
			}
		}else{
			#print("no")
			#print(selected)
			if(selected==0){
				s.cd <- sample(1:2,1)
			}else{
				s.cd <- selected
			}
		}
		#print(s.cd)
		h.cd <- sample(1:2,1,prob=c(ps[s.cd],1-ps[s.cd]))
		X.cd[j,(s.cd-1)*2+h.cd] <- X.cd[j,(s.cd-1)*2+h.cd] + 1

	}
	sv.bd[i,] <- apply(X.bd[,c(1,3)],1,sum)
	sv.cd[i,] <- apply(X.cd[,c(1,3)],1,sum)
}

boxplot(sv.bd)

sv.prob.bd <- apply(t(apply(sv.bd,1,diff)),2,mean)
sv.prob.cd <- apply(t(apply(sv.cd,1,diff)),2,mean)

matplot(cbind(sv.prob.bd,sv.prob.cd),type="l")

m.bd <- apply(sv.bd,2,mean)
m.cd <- apply(sv.cd,2,mean)

matplot(cbind(m.bd,m.cd),type="l")


incl <- apply(X.bd,2,diff)
stst <- apply(incl,1,which.max)

kigou <- c("●","○","■","□")

stst.k <- kigou[stst]
stst.k