決断したら、「次、行ってみよう！」

• こちらの続き
• 昨日は２ｘ２分割表の２群に１標本ずつ加えて、総標本数を２ずつ増やしながら２ｘ２分割表を大きくするときの分割表の推移パターンについて考えた
• 今日は、そんな分割表に関して、標本数が大きくなるたびに検定をして、帰無仮説が棄却されたら、「結論」を出して、有用な方に決め打ちにしてしまう(それ以降は蓄積するデータに頓着しない)ような方針のときの、成否の推移を「正確確率」的においかけよう、という話
• ソースは汚いので、これは解析もつけずにメモするだけ

my.test <- function(X){
	m <- matrix(X,2,2)
	tmp1 <- apply(m,1,sum)
	tmp2 <- apply(m,2,sum)
	if(prod(c(tmp1,tmp2))==0)return(1)
	#return(chisq.test(m,correct=FALSE)$p.value)
	return(fisher.test(m)$p.value)
}


ps <- c(0.1,0.9)
mean.ps <- mean(ps)

n <- 100

thress <- seq(from=0,to=1,by=0.1)
thress[1] <- 0.01
Xs <- X.fixs <- list()
X1s <- X2s <- list()
for(i in 1:length(thress)){
	Xs[[i]] <- list()
	Xs[[i]][[1]] <- matrix(1,1,1)
	X1s[[i]] <- X2s[[i]] <- list()
	X1s[[i]][[1]] <- X2s[[i]][[1]] <- matrix(0,1,1)
	X.fixs[[i]] <- list()
	X.fixs[[i]] <- matrix(0,n,3)
	X.fixs[[i]][1,2] <- 1
}

P.hist[[1]] <- matrix(1,1,1)



thres <- 0.1
for(i in 2:n){
	Y <- cbind(0:(i-1),(i-1):0)
	e <- expand.grid(1:(i),1:(i))
	Z <- cbind(Y[e[,1],],Y[e[,2],])
	Ps <- apply(Z,1,my.test)
	
	Ps.m <- matrix(Ps,i,i)
	P.hist[[i]] <- Ps.m
	Ps.m.sign <- matrix(0,i,i)
	for(j in 1:length(thress)){
		Xs[[j]][[i]] <- matrix(0,i,i)
		Xs[[j]][[i]][1:(i-1),1:(i-1)] <- Xs[[j]][[i]][1:(i-1),1:(i-1)] + Xs[[j]][[i-1]]*ps[1]*ps[2]
		Xs[[j]][[i]][1:(i-1),2:i] <- Xs[[j]][[i]][1:(i-1),2:i] + Xs[[j]][[i-1]]*ps[1]*(1-ps[2])
		Xs[[j]][[i]][2:i,1:(i-1)] <- Xs[[j]][[i]][2:i,1:(i-1)] + Xs[[j]][[i-1]]*(1-ps[1])*ps[2] 
		Xs[[j]][[i]][2:i,2:i] <- Xs[[j]][[i]][2:i,2:i] + Xs[[j]][[i-1]]*(1-ps[1])*(1-ps[2])
		
		X1s[[j]][[i]] <- X2s[[j]][[i]] <- matrix(0,i,i)
		X1s[[j]][[i]][1:(i-1),1:(i-1)] <- X1s[[j]][[i]][1:(i-1),1:(i-1)] + X1s[[j]][[i-1]]*ps[1]*ps[1]
		X1s[[j]][[i]][1:(i-1),2:i] <- X1s[[j]][[i]][1:(i-1),2:i] + X1s[[j]][[i-1]]*ps[1]*(1-ps[1])
		X1s[[j]][[i]][2:i,1:(i-1)] <- X1s[[j]][[i]][2:i,1:(i-1)] + X1s[[j]][[i-1]]*(1-ps[1])*ps[1] 
		X1s[[j]][[i]][2:i,2:i] <- X1s[[j]][[i]][2:i,2:i] + X1s[[j]][[i-1]]*(1-ps[1])*(1-ps[1])
		
		X2s[[j]][[i]][1:(i-1),1:(i-1)] <- X2s[[j]][[i]][1:(i-1),1:(i-1)] + X2s[[j]][[i-1]]*ps[2]*ps[2]
		X2s[[j]][[i]][1:(i-1),2:i] <- X2s[[j]][[i]][1:(i-1),2:i] + X2s[[j]][[i-1]]*ps[2]*(1-ps[2])
		X2s[[j]][[i]][2:i,1:(i-1)] <- X2s[[j]][[i]][2:i,1:(i-1)] + X2s[[j]][[i-1]]*(1-ps[2])*ps[2] 
		X2s[[j]][[i]][2:i,2:i] <- X2s[[j]][[i]][2:i,2:i] + X2s[[j]][[i-1]]*(1-ps[2])*(1-ps[2])

		Ps.m.sign[which(Ps.m < thress[j])] <- 1
		Ps.m.sign.2 <- -(Ps.m.sign-1)
		tmp.Xi <- Xs[[j]][[i]] * Ps.m.sign
		print(tmp.Xi)
		X1s[[j]][[i]][upper.tri(tmp.Xi)] <- X1s[[j]][[i]][upper.tri(tmp.Xi)] + tmp.Xi[upper.tri(tmp.Xi)]
		X2s[[j]][[i]][lower.tri(tmp.Xi)] <- X2s[[j]][[i]][lower.tri(tmp.Xi)] + tmp.Xi[lower.tri(tmp.Xi)]
		incl.1 <- sum(tmp.Xi[upper.tri(Xs[[j]][[i]])])
		X.fixs[[j]][i,1] <- X.fixs[[j]][i-1,1] + incl.1
		X.fixs[[j]][i,3] <- X.fixs[[j]][i-1,3] + sum(tmp.Xi)-incl.1
		Xs[[j]][[i]] <- Xs[[j]][[i]] * Ps.m.sign.2
		X.fixs[[j]][i,2] <- sum(Xs[[j]][[i]])
	}
}
matplot(X.fixs[[2]],type="l")

X.sum <- list()
X.saved <- list()
X.saved.ratio <- list()
for(i in 1:length(thress)){
	X.sum[[i]] <- list()
	X.saved[[i]] <- list()
	X.saved.ratio[[i]] <- rep(0,n)
	for(j in 1:n){
		X.sum[[i]][[j]] <- Xs[[i]][[j]] + X1s[[i]][[j]] + X2s[[i]][[j]]
		X.saved[[i]][[j]] <- rep(0,length(X.sum[[i]][[j]][,1])*2)
		for(k1 in 1:length(X.sum[[i]][[j]][,1])){
			for(k2 in 1:length(X.sum[[i]][[j]][,1])){
				X.saved[[i]][[j]][k1+k2-2] <- X.saved[[i]][[j]][k1+k2-2] + X.sum[[i]][[j]][k1,k2]
			}
		}
		X.saved.ratio[[i]][j] <- sum((length(X.saved[[i]][[j]])-1):0 * X.saved[[i]][[j]])
	}
}
plot(X.saved.ratio[[2]]/(2*(1:n)))

plot(X.saved.ratio[[2]],X.saved.ratio[[10]])
abline(0,1,col=2)