my.2x2.dev <- function(p,q,r){
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]
Aa <- my.2x2.dev(pd,pa,ra)
Bb <- my.2x2.dev(pd,pb,rb)
AaBbDd <- my.2x2x2.dev(Aa,Bb,rd)
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)
}