■

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