■

```library(combinat)

n <- 10
k1 <- k2 <- 2
k <- k1 + k2
knall <- nsimplex(k,0:n)
kn1 <- nsimplex(k1,0:n)
kn2 <- nsimplex(k2,n:0)
kn1 * kn2
knall

X <- list()
n <- 4
for(i in 0:n){
X[[i+1]] <- as.matrix(xsimplex(k1,i))
}

for(i in 1:length(X)){
for(j in 1:length(X[[i]][1,])){
for(j2 in 1:length(X[[n+1-i+1]][1,])){
case <- X[[i]][,j]
control <- X[[n+1-i+1]][,j2]
m2x2 <- matrix(c(case,control),byrow=TRUE,2,2)
print(m2x2)
}
}
}

ps <- c(0.3,0.7)
Ps <- list()
Ps[[1]] <- list()
Ps[[1]][[1]] <- matrix(1,1,1)
for(i in 1:n){
select1 <- 0.5
select2 <- 0.5
p11 <- select1 * ps[1]
p12 <- select1 * (1-ps[1])
p21 <- select2 * ps[2]
p22 <- select2 * (1-ps[2])
#p.matrix <- matrix(c(p11,p12,p21,p22),byrow=TRUE,2,2)
p.vec <- c(p11,p12,p21,p22)
Ps[[i+1]] <- list()
for(j in 1:(i+1)){
Ps[[i+1]][[j]] <- matrix(0,length(X[[j]][1,]),length(X[[i+1-j+1]][1,]))
}

for(j in 1:4){
if(j ==1){
xst <- 1
yst <- 1
list.st <- 0
}else if(j == 2){
xst <- 1
yst <- 2
list.st <- 0
}else if(j == 3){
xst <- 1
yst <- 1
list.st <- 1
}else if(j == 4){
xst <- 2
yst <- 1
list.st <- 1
}
for(j2 in 1:i){
dims <- dim(Ps[[i]][[j2]])
print(dims)
print(Ps[[i]][[j2]])
print(Ps[[i+1]][[j2+list.st]][xst:(xst+dims[1]-1),yst:(yst+dims[2]-1)])
print(Ps[[i]][[j2]])
Ps[[i+1]][[j2+list.st]][xst:(xst+dims[1]-1),yst:(yst+dims[2]-1)] <- Ps[[i+1]][[j2+list.st]][xst:(xst+dims[1]-1),yst:(yst+dims[2]-1)] + Ps[[i]][[j2]] * p.vec[j]
}

}

}
```
• 関数化してちょっときれいに
```library(combinat)

n <- 10
k1 <- k2 <- 2
k <- k1 + k2
knall <- nsimplex(k,0:n)
kn1 <- nsimplex(k1,0:n)
kn2 <- nsimplex(k2,n:0)
kn1 * kn2
knall

X <- list()
n <- 4
for(i in 0:n){
X[[i+1]] <- as.matrix(xsimplex(k1,i))
}

for(i in 1:length(X)){
for(j in 1:length(X[[i]][1,])){
for(j2 in 1:length(X[[n+1-i+1]][1,])){
case <- X[[i]][,j]
control <- X[[n+1-i+1]][,j2]
m2x2 <- matrix(c(case,control),byrow=TRUE,2,2)
print(m2x2)
}
}
}

inc1 <- function(Ps,X,ps,qs){
select1 <- qs[1]
select2 <- qs[2]
p11 <- select1 * ps[1]
p12 <- select1 * (1-ps[1])
p21 <- select2 * ps[2]
p22 <- select2 * (1-ps[2])
p.vec <- c(p11,p12,p21,p22)
Ret <- list()
n <- length(Ps)
for(j in 1:(n+1)){
Ret[[j]] <- matrix(0,length(X[[j]][1,]),length(X[[i+1-j+1]][1,]))
}
for(j in 1:4){
if(j ==1){
xst <- 1
yst <- 1
list.st <- 0
}else if(j == 2){
xst <- 1
yst <- 2
list.st <- 0
}else if(j == 3){
xst <- 1
yst <- 1
list.st <- 1
}else if(j == 4){
xst <- 2
yst <- 1
list.st <- 1
}
for(j2 in 1:n){
dims <- dim(Ps[[j2]])
Ret[[j2+list.st]][xst:(xst+dims[1]-1),yst:(yst+dims[2]-1)] <- Ret[[j2+list.st]][xst:(xst+dims[1]-1),yst:(yst+dims[2]-1)] + Ps[[j2]] * p.vec[j]
}

}
Ret
}
ps <- c(0.3,0.7)
Ps <- list()
Ps[[1]] <- list()
Ps[[1]][[1]] <- matrix(1,1,1)
for(i in 1:n){
select1 <- 0.5
select2 <- 0.5
p11 <- select1 * ps[1]
p12 <- select1 * (1-ps[1])
p21 <- select2 * ps[2]
p22 <- select2 * (1-ps[2])
#p.matrix <- matrix(c(p11,p12,p21,p22),byrow=TRUE,2,2)
p.vec <- c(p11,p12,p21,p22)
Ps[[i+1]] <- inc1(Ps[[i]],X,ps,c(0.5,0.5))
}
```