メモ

```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)
}

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,]))

}