# a poor solution sim.new <- function(n,maf,nit,d=c(1,2)){ n0 <- n1 <- n/2 y <- c(rep(0,n0),rep(1,n1)) y.mean <- mean(y) p <- c(maf^2,2*maf*(1-maf),(1-maf)^2) z <- rep(0,nit) for(k in 1:nit){ g <- sample(0:2,n,p=p,replace=TRUE) x <- table(c(g,2),c(y,1))[-1,] x[2,2] <- x[2,2]-1 num <- (x[1,1]+2*x[2,1])*(0-y.mean)+(x[1,2]+2*x[2,2])*(1-y.mean) wbar <- table(c(g,2)) wbar[3] <- wbar[3]-1 wbar <- (wbar/n)[-1] den <- sqrt(n*y.mean*(1-y.mean)*(sum((d^2)*wbar*(1-wbar))-2*prod(d)*prod(wbar))) z[k] <- num/den } return(z) } # a much better solution sim.new2 <- function(n,maf,nit,d=c(1,2)){ p <- c(maf^2,2*maf*(1-maf),(1-maf)^2) p0 <- p[1] p12 <- p[-1]/sum(p[-1]) z <- rep(0,nit) for(k in 1:nit){ n0 <- n/2-rbinom(1,n/2,p[1]) n1 <- n/2-rbinom(1,n/2,p[1]) n01 <- rbinom(1,n0,p12[1]) n02 <- n0-n01 n11 <- rbinom(1,n1,p12[1]) n12 <- n1-n11 num <- -0.5*(n01+2*n02)+0.5*(n11+2*n12) w1 <- (n01+n11)/n w2 <- (n02+n12)/n den <- sqrt(n*0.25*(d[1]^2*w1*(1-w1)+d[2]^2*w2*(1-w2)-2*prod(d)*w1*w2)) z[k] <- num/den } return(z) } # simulations print(system.time(z <- sim.new(2000,0.95,1e3))) print(system.time(z <- sim.new2(2000,0.95,1e3)))