set.seed(1234) num.sim <- 100000 num.digits <- 53 num.pull <- 6 ys.sum <- rep(0,num.pull) ys.sum2 <- matrix(rep(0,num.pull*num.pull),ncol=num.pull) p.ord <- matrix(rep(0,num.digits*num.pull),ncol=num.pull) for (i in 1:num.sim) { y <- sample(1:num.digits,num.pull,replace=F) ys <- sort(y) for (j1 in 1:num.pull) { ys.sum[j1] <- ys.sum[j1] + ys[j1]/num.sim p.ord[ys[j1],j1] <- p.ord[ys[j1],j1] + 1/num.sim for (j2 in 1:num.pull) { ys.sum2[j1,j2] <- ys.sum2[j1,j2] + ys[j1]*ys[j2]/num.sim }} } ys.sum ys.sum2 p.ord cov.y <- matrix(rep(0,num.pull*num.pull),ncol=num.pull) for (j1 in 1:num.pull) { for (j2 in 1:num.pull) { cov.y[j1,j2] <- ys.sum2[j1,j2]-ys.sum[j1]*ys.sum[j2] }} cov.y corr.y <- matrix(rep(0,num.pull*num.pull),ncol=num.pull) for (j1 in 1:num.pull) { for (j2 in 1:num.pull) { corr.y[j1,j2] <- cov.y[j1,j2]/sqrt(cov.y[j1,j1]*cov.y[j2,j2]) }} corr.y