メモ

ps <- c(0.2,0.8)

N <- 30
library(combinat)
calc.survival.2 <- function(X){
	ret <- 0
	for(i in 1:length(X)){
		tmp <- expand.grid(0:(length(X[[i]][,1])-1),0:(length(X[[i]][1,])-1))
		tmp.2 <- apply(tmp,1,sum)
		ret <- ret + sum(tmp.2 * X[[i]])
	}
	return(list(n=ret,f=ret/(length(X)-1)))
	
}
calc.survival.3 <- function(X){
	ret <- 0
	for(i in 1:length(X)){
		tmp <- expand.grid(0:(length(X[[i]][,1])-1),0:(length(X[[i]][1,])-1))
		tmp.2 <- apply(tmp,1,sum)
		#tmp.3 <- max(tmp.2)-tmp.2
		ret <- ret + sum(tmp.2 * X[[i]])
	}
	return(ret/(length(X)-1))
	
}
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)
}

BetaDecision <- function(ps,N=20,pl=TRUE,X0=rep(0,4)){
	k1 <- k2 <- 2
	k <- k1 + k2
	M <- list()
	for(i in 0:N){
		M[[i+1]] <- as.matrix(xsimplex(k1,i))
	}

	z <- rep(2,N)

	Xsm <- list()
	Xsm[[1]] <- list(matrix(1,1,1))

	for(i in 1:N){
		n <- i+1
		Xsm[[n]] <- list()
		for(j in 1:n){
			Xsm[[n]][[j]] <- matrix(0,length(M[[j]][1,]),length(M[[n-j+1]][1,]))
		}
		for(j2 in 1:i){
			Y1 <- t(M[[j2]])
			
			Y2 <- t(M[[n-j2]])
			e <- expand.grid(1:length(Y1[,1]),1:length(Y2[,1]))
			Z <- cbind(matrix(Y1[e[,1],],ncol=2),matrix(Y2[e[,2],],ncol=2))
			Z <- t(t(Z)+X0)
			Qs <- t(apply(Z,1,Decision_Dice_beta.2.3,2))
			#print(Qs)
			dims <- dim(Xsm[[n-1]][[j2]])
			# choice 1 & event 1
			xst <- 1;yst <- 2;list.st <- 0
			#print(matrix(Qs[,1],length(Y1[,1]),length(Y2[,1])))
			#print(Xsm[[n-1]][[j2]])
			Xsm[[n]][[j2+list.st]][xst:(xst+dims[1]-1),yst:(yst+dims[2]-1)] <- Xsm[[n]][[j2+list.st]][xst:(xst+dims[1]-1),yst:(yst+dims[2]-1)] + Xsm[[n-1]][[j2]] * ps[1] * matrix(Qs[,1],length(Y1[,1]),length(Y2[,1]))
			# choice 1 & event 2
			xst <- 1;yst <- 1;list.st <- 0
			Xsm[[n]][[j2+list.st]][xst:(xst+dims[1]-1),yst:(yst+dims[2]-1)] <- Xsm[[n]][[j2+list.st]][xst:(xst+dims[1]-1),yst:(yst+dims[2]-1)] + Xsm[[n-1]][[j2]] * (1-ps[1]) * matrix(Qs[,1],length(Y1[,1]),length(Y2[,1]))
			# choice 2 & event 1
			xst <- 2;yst <- 1;list.st <- 1
			Xsm[[n]][[j2+list.st]][xst:(xst+dims[1]-1),yst:(yst+dims[2]-1)] <- Xsm[[n]][[j2+list.st]][xst:(xst+dims[1]-1),yst:(yst+dims[2]-1)] + Xsm[[n-1]][[j2]] * ps[2] * matrix(Qs[,2],length(Y1[,1]),length(Y2[,1]))
			# choice 2 & event 2
			xst <- 1;yst <- 1;list.st <- 1
			Xsm[[n]][[j2+list.st]][xst:(xst+dims[1]-1),yst:(yst+dims[2]-1)] <- Xsm[[n]][[j2+list.st]][xst:(xst+dims[1]-1),yst:(yst+dims[2]-1)] + Xsm[[n-1]][[j2]] * (1-ps[2]) * matrix(Qs[,2],length(Y1[,1]),length(Y2[,1]))
		}
	}
	sv <- sapply(Xsm,calc.survival.2)
	if(pl){
		plot(unlist(sv[2,]))
	}
	

	return(list(X=Xsm,sv.n=sv[1,],sv.f=sv[2,]))
	


}

BetaDecision(c(0.2,0.6),20,X0=c(10,4,6,3))