CR <- read.table("CaseRandom.csv",header=TRUE,sep=",")
table(CR$snp1,CR$snp3)

N <- c(837,995,816,2112)
H <- 880
hetero <- function(th,N,H)
{
    f <- N + c(th,1-th,1-th,th)*H
    t <- f[1]*f[4]/(f[1]*f[4] + f[2]*f[3])
    return(th - t)
}
th <- uniroot(hetero,c(0,1),N=N,H=H)$root
th

f <- N + c(th,1-th,1-th,th)*H
p <- f/sum(f)
c <- p[1]*p[4]-p[2]*p[3]
p1 <- p[1]+p[2]
p2 <- p[1]+p[3]
r <- c/sqrt(p1*(1-p1)*p2*(1-p2))
r