放っておけば相乗的

my.2x2.dev <- function(p,q,r){
	#tmp <- r*sqrt(p*q*(1-p)*(1-q))
	x <- c(p*q,p*(1-q),(1-p)*q,(1-p)*(1-q))
	tmp <- r*min(x)
	matrix(x+c(tmp,-tmp,-tmp,tmp),byrow=TRUE,2,2)
}

my.2x2x2.dev <- function(p,q,r){
	tmp <- r*sqrt(prod(p)*prod(q)*r)
	tmp1 <- outer(p[1,],q[1,],"*")/sum(p[1,])
	tmp2 <- outer(p[2,],q[2,],"*")/sum(p[2,])
	tmp1 <- tmp1 + c(tmp,-tmp,-tmp,tmp)
	tmp2 <- tmp2 - c(tmp,-tmp,-tmp,tmp)
	array(c(tmp1,tmp2),rep(2,3))
}
n.iter <- 20

pas <- runif(n.iter)
pbs <- runif(n.iter)
pds <- runif(n.iter)
ras <- runif(n.iter)
rbs <- runif(n.iter)
rd <- 0
par(ask=TRUE)
for(i in 1:n.iter){
	pa <- pas[i]
	pb <- pbs[i]
	pd <- pds[i]
	ra <- ras[i]
	rb <- rbs[i]
	#rd <- 0
	Aa <- my.2x2.dev(pd,pa,ra)
	#apply(Aa,1,sum)
	Bb <- my.2x2.dev(pd,pb,rb)
	#apply(Bb,1,sum)
	#print(Aa[1,]/Aa[2,])
	#print(Bb[1,]/Bb[2,])
	AaBbDd <- my.2x2x2.dev(Aa,Bb,rd)
	#AaBbDd

	#sum(AaBbDd)
	RR <- AaBbDd[,,1]/AaBbDd[,,2]
	RR.st <- RR/RR[4]

	barplot(c(RR.st))
	RR.Aa <- Aa[1,]/Aa[2,]
	print(RR.Aa[1]/RR.Aa[2])
	RR.Bb <- Bb[1,]/Bb[2,]
	print(RR.Bb[1]/RR.Bb[2])
	print(RR.Aa[1]*RR.Bb[1]/(RR.Aa[2]*RR.Bb[2]))
	print(RR.st)
}