R version 3.2.4 (2016-03-10) -- "Very Secure Dishes" Copyright (C) 2016 The R Foundation for Statistical Computing Platform: x86_64-pc-linux-gnu (64-bit) R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. Natural language support but running in an English locale R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > ## code to reproduce Table 2 > set.seed(5003) > options(width=150) > sessionInfo() R version 3.2.4 (2016-03-10) Platform: x86_64-pc-linux-gnu (64-bit) Running under: Red Hat Enterprise Linux Server release 6.8 (Santiago) locale: [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 LC_MONETARY=en_US.UTF-8 [6] LC_MESSAGES=en_US.UTF-8 LC_PAPER=en_US.UTF-8 LC_NAME=C LC_ADDRESS=C LC_TELEPHONE=C [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C attached base packages: [1] stats graphics grDevices utils datasets methods base > > ################################################ > ## Loop over lambda to create raw data > ## break simulation into 5 chunks of 50,000,000 each to keep RAM manageable > ################################################ > nchunk <- 5 > n <- 50000000 > bigres <- vector(nchunk, mode="list") > > for(ch in 1:nchunk){ + lambdaseq <- c(1,2) + nlam <- length(lambdaseq) + res <- vector(nlam, mode="list") + + ## loop over lambda + for(wh in 1:nlam){ + lambda <- lambdaseq[wh] + + cat(paste("\n\n#########################",ch,lambda,"#########################\n")) + + ## generate X + d <- data.frame(x1 = rnorm(n)) + + ## generate Ws as conditionally normal given X + g1 <- function(x){ + 0.07*(x+1)^2 + } + ## check reliabilty + ## var(d$x) / (var(d$x) + mean(g1(d$x))) + + d$w1 <- d$x1 + rnorm(n, sd = sqrt(g1(d$x1))) + + ## corresponding A(W) function, closed form in this case + A1 <- function(w){ + (.07/(.07 + 1))*(w+1)^2 + } + + ## get other plug-ins + d$egx <- mean(g1(d$x1)) + + ## brute force condition on w by rounding it to nearest 0.1 + d$wcat <- as.factor(round(d$w1, digits=1)) + d$exgivenw <- ave(d$x1, d$wcat) + d$egxgivenw <- ave(g1(d$x1), d$wcat) + + ####################################################### + ## generate different values of W* + ####################################################### + ## 1) using true variance + d$w1star.gx <- d$w1 + rnorm(n, sd = sqrt(lambda * g1(d$x1))) + ## 2) fresh sample using right variance function + d$w1star.fresh <- d$x1 + rnorm(n, sd = sqrt((1 + lambda) * g1(d$x1))) + ## 3) average g(X) + d$w1star.avggx <- d$w1 + rnorm(n, sd = sqrt(lambda * d$egx)) + ## 4) average g(W) + d$w1star.avggw <- d$w1 + rnorm(n, sd = sqrt(lambda * mean(g1(d$w1)))) + ## 5) g(W) plugin + d$w1star.gw <- d$w1 + rnorm(n, sd = sqrt(lambda * g1(d$w1))) + ## 6) g(E[X|W]) + d$w1star.gexw <- d$w1 + rnorm(n, sd = sqrt(lambda * g1(d$exgivenw))) + ## 7) A(W) + d$w1star.aw <- d$w1 + rnorm(n, sd = sqrt(lambda * A1(d$w1))) + ## 8) E[g(X)|W] + d$w1star.rcal <- d$w1 + rnorm(n, sd = sqrt(lambda * d$egxgivenw)) + + + m1 <- c(mean(d$w1star.gx), mean(d$w1star.gx^2), mean(d$w1star.gx^3), mean(d$w1star.gx^4), mean(d$w1star.gx^5), mean(d$w1star.gx^6), mean(d$w1star.gx^7),mean(d$w1star.gx^8)) + m2 <- c(mean(d$w1star.fresh), mean(d$w1star.fresh^2), mean(d$w1star.fresh^3), mean(d$w1star.fresh^4), mean(d$w1star.fresh^5), mean(d$w1star.fresh^6), mean(d$w1star.fresh^7),mean(d$w1star.fresh^8)) + m3 <- c(mean(d$w1star.avggx), mean(d$w1star.avggx^2), mean(d$w1star.avggx^3), mean(d$w1star.avggx^4), mean(d$w1star.avggx^5), mean(d$w1star.avggx^6), mean(d$w1star.avggx^7), mean(d$w1star.avggx^8)) + m4 <- c(mean(d$w1star.avggw), mean(d$w1star.avggw^2), mean(d$w1star.avggw^3), mean(d$w1star.avggw^4), mean(d$w1star.avggw^5), mean(d$w1star.avggw^6), mean(d$w1star.avggw^7),mean(d$w1star.avggw^8)) + m5 <- c(mean(d$w1star.gw), mean(d$w1star.gw^2), mean(d$w1star.gw^3), mean(d$w1star.gw^4), mean(d$w1star.gw^5), mean(d$w1star.gw^6), mean(d$w1star.gw^7),mean(d$w1star.gw^8)) + m6 <- c(mean(d$w1star.gexw), mean(d$w1star.gexw^2), mean(d$w1star.gexw^3), mean(d$w1star.gexw^4), mean(d$w1star.gexw^5), mean(d$w1star.gexw^6), mean(d$w1star.gexw^7),mean(d$w1star.gexw^8)) + m7 <- c(mean(d$w1star.aw), mean(d$w1star.aw^2), mean(d$w1star.aw^3), mean(d$w1star.aw^4), mean(d$w1star.aw^5), mean(d$w1star.aw^6), mean(d$w1star.aw^7), mean(d$w1star.aw^8)) + m8 <- c(mean(d$w1star.rcal), mean(d$w1star.rcal^2), mean(d$w1star.rcal^3), mean(d$w1star.rcal^4), mean(d$w1star.rcal^5), mean(d$w1star.rcal^6), mean(d$w1star.rcal^7), mean(d$w1star.rcal^8)) + + m <- round(rbind(m1, m2, m3, m4, m5, m6, m7, m8), digits=3) + colnames(m) <- paste("mom",1:8,sep="") + print(tmp <- data.frame(lambda = lambda, method=c("gx","fresh","avggx","avggw","gw","gexw","aw","rcal"), m)) + + ## rearrange tmp to align with table order in paper and save + res[[wh]] <- tmp[c(1,2,7,3,8,5,4,6),] + rm(d) + gc() + } + + ## store results for this chunk + bigres[[ch]] <- res + rm(res) + gc() + } ######################### 1 1 ######################### lambda method mom1 mom2 mom3 mom4 mom5 mom6 mom7 mom8 m1 1 gx 0 1.280 0.840 6.951 13.113 85.913 281.842 1855.344 m2 1 fresh 0 1.280 0.839 6.948 13.078 85.843 281.500 1859.818 m3 1 avggx 0 1.280 0.418 5.841 5.946 52.848 102.978 779.586 m4 1 avggw 0 1.290 0.419 5.918 5.996 53.719 104.149 792.977 m5 1 gw 0 1.290 0.985 7.881 18.433 122.888 514.188 3659.496 m6 1 gexw 0 1.273 0.826 6.853 12.736 83.455 268.485 1756.905 m7 1 aw 0 1.280 0.949 7.684 17.490 116.914 478.527 3411.777 m8 1 rcal 0 1.280 0.838 6.924 12.934 84.737 273.839 1802.559 ######################### 1 2 ######################### lambda method mom1 mom2 mom3 mom4 mom5 mom6 mom7 mom8 m1 2 gx 0 1.421 1.264 9.374 23.257 150.121 621.885 4318.304 m2 2 fresh 0 1.420 1.261 9.364 23.244 150.144 623.364 4316.698 m3 2 avggx 0 1.420 0.422 6.983 6.570 66.505 121.995 1020.633 m4 2 avggw 0 1.440 0.422 7.149 6.672 68.595 125.589 1061.236 m5 2 gw 0 1.441 1.561 11.422 36.586 255.116 1403.192 11400.655 m6 2 gexw 0 1.406 1.239 9.131 22.244 142.523 577.372 3945.949 m7 2 aw 0 1.421 1.482 10.915 33.801 234.008 1250.175 9978.260 m8 2 rcal 0 1.421 1.263 9.290 22.780 146.054 592.459 4050.075 ######################### 2 1 ######################### lambda method mom1 mom2 mom3 mom4 mom5 mom6 mom7 mom8 m1 1 gx 0 1.281 0.844 6.960 13.175 86.226 283.605 1867.153 m2 1 fresh 0 1.280 0.840 6.945 13.072 85.664 278.976 1834.865 m3 1 avggx 0 1.280 0.421 5.849 5.966 52.958 103.056 780.485 m4 1 avggw 0 1.290 0.422 5.924 6.023 53.867 104.696 797.869 m5 1 gw 0 1.291 0.992 7.918 18.652 124.428 524.770 3748.795 m6 1 gexw 0 1.273 0.830 6.866 12.830 84.101 273.698 1804.221 m7 1 aw 0 1.280 0.952 7.686 17.530 116.861 478.119 3379.345 m8 1 rcal 0 1.280 0.842 6.936 13.027 85.326 277.832 1831.358 ######################### 2 2 ######################### lambda method mom1 mom2 mom3 mom4 mom5 mom6 mom7 mom8 m1 2 gx 0 1.420 1.258 9.364 23.129 149.673 615.189 4260.932 m2 2 fresh 0 1.420 1.261 9.371 23.188 149.764 615.822 4239.415 m3 2 avggx 0 1.420 0.420 6.982 6.562 66.609 122.971 1029.817 m4 2 avggw 0 1.440 0.421 7.149 6.649 68.651 125.309 1064.117 m5 2 gw 0 1.440 1.554 11.399 36.391 254.405 1395.237 11352.065 m6 2 gexw 0 1.406 1.237 9.132 22.232 142.738 578.456 3964.286 m7 2 aw 0 1.420 1.479 10.908 33.789 234.217 1260.153 10123.492 m8 2 rcal 0 1.420 1.261 9.287 22.771 146.151 593.889 4065.909 ######################### 3 1 ######################### lambda method mom1 mom2 mom3 mom4 mom5 mom6 mom7 mom8 m1 1 gx 0 1.280 0.839 6.943 13.061 85.610 279.050 1835.934 m2 1 fresh 0 1.280 0.840 6.955 13.143 86.259 283.151 1869.396 m3 1 avggx 0 1.280 0.418 5.842 5.945 52.951 103.369 786.781 m4 1 avggw 0 1.290 0.418 5.915 5.979 53.775 104.452 800.376 m5 1 gw 0 1.289 0.987 7.901 18.612 124.665 530.362 3822.958 m6 1 gexw 0 1.273 0.826 6.850 12.704 83.348 267.167 1754.227 m7 1 aw 0 1.280 0.950 7.688 17.551 117.507 484.110 3464.903 m8 1 rcal 0 1.280 0.839 6.921 12.952 84.761 274.472 1802.315 ######################### 3 2 ######################### lambda method mom1 mom2 mom3 mom4 mom5 mom6 mom7 mom8 m1 2 gx 0 1.420 1.261 9.370 23.233 150.067 619.894 4281.107 m2 2 fresh 0 1.420 1.259 9.359 23.166 149.855 617.398 4274.820 m3 2 avggx 0 1.420 0.422 6.986 6.574 66.602 122.045 1023.454 m4 2 avggw 0 1.440 0.421 7.153 6.654 68.657 124.814 1061.273 m5 2 gw 0 1.440 1.557 11.401 36.422 253.890 1387.416 11250.230 m6 2 gexw 0 1.407 1.240 9.145 22.293 143.009 579.014 3955.562 m7 2 aw 0 1.421 1.483 10.933 33.854 234.249 1244.685 9841.000 m8 2 rcal 0 1.420 1.261 9.280 22.652 145.414 585.406 4009.775 ######################### 4 1 ######################### lambda method mom1 mom2 mom3 mom4 mom5 mom6 mom7 mom8 m1 1 gx 0 1.279 0.840 6.939 13.066 85.532 278.535 1828.451 m2 1 fresh 0 1.280 0.839 6.949 13.105 86.113 281.877 1869.507 m3 1 avggx 0 1.279 0.419 5.836 5.942 52.800 102.569 779.557 m4 1 avggw 0 1.289 0.420 5.913 5.991 53.682 103.917 794.191 m5 1 gw 0 1.289 0.987 7.897 18.556 124.284 525.527 3785.652 m6 1 gexw 0 1.272 0.827 6.848 12.738 83.383 268.308 1757.559 m7 1 aw 0 1.280 0.949 7.680 17.535 117.446 488.661 3556.413 m8 1 rcal 0 1.279 0.838 6.915 12.914 84.653 272.441 1799.649 ######################### 4 2 ######################### lambda method mom1 mom2 mom3 mom4 mom5 mom6 mom7 mom8 m1 2 gx 0 1.420 1.259 9.362 23.171 149.852 619.460 4300.750 m2 2 fresh 0 1.420 1.261 9.365 23.154 149.501 611.903 4218.553 m3 2 avggx 0 1.420 0.419 6.976 6.535 66.294 120.955 1012.454 m4 2 avggw 0 1.439 0.419 7.139 6.616 68.302 123.867 1048.665 m5 2 gw 0 1.440 1.553 11.391 36.269 252.564 1367.244 10984.393 m6 2 gexw 0 1.406 1.235 9.115 22.123 141.656 569.445 3883.532 m7 2 aw 0 1.420 1.478 10.898 33.619 232.037 1225.952 9673.318 m8 2 rcal 0 1.419 1.257 9.268 22.669 145.440 589.742 4028.323 ######################### 5 1 ######################### lambda method mom1 mom2 mom3 mom4 mom5 mom6 mom7 mom8 m1 1 gx 0 1.280 0.841 6.956 13.144 86.256 283.623 1870.018 m2 1 fresh 0 1.280 0.841 6.948 13.108 85.853 279.976 1838.431 m3 1 avggx 0 1.280 0.421 5.846 5.976 53.062 103.714 786.272 m4 1 avggw 0 1.290 0.421 5.925 6.036 54.038 105.784 807.470 m5 1 gw 0 1.290 0.986 7.893 18.444 123.381 513.643 3683.386 m6 1 gexw 0 1.273 0.829 6.869 12.806 84.055 271.625 1786.522 m7 1 aw 0 1.280 0.951 7.689 17.575 117.524 484.530 3447.542 m8 1 rcal 0 1.280 0.840 6.926 12.989 85.060 276.681 1823.427 ######################### 5 2 ######################### lambda method mom1 mom2 mom3 mom4 mom5 mom6 mom7 mom8 m1 2 gx 0 1.420 1.259 9.369 23.177 150.205 620.408 4319.487 m2 2 fresh 0 1.420 1.258 9.362 23.164 149.737 617.056 4257.517 m3 2 avggx 0 1.420 0.418 6.974 6.515 66.336 120.870 1013.967 m4 2 avggw 0 1.439 0.418 7.142 6.593 68.398 123.716 1052.036 m5 2 gw 0 1.439 1.550 11.372 36.201 252.266 1375.435 11074.974 m6 2 gexw 0 1.406 1.234 9.121 22.149 142.335 574.383 3936.537 m7 2 aw 0 1.420 1.478 10.910 33.757 233.857 1246.729 9903.431 m8 2 rcal 0 1.420 1.259 9.284 22.701 145.991 593.139 4084.999 > > ################################################ > ## Create table > ################################################ > res <- vector(nchunk, mode="list") > for(ch in 1:nchunk){ + res[[ch]] <- rbind(bigres[[ch]][[1]], bigres[[ch]][[2]]) + res[[ch]]$sim <- ch + } > > res <- do.call("rbind", res) > res$mom1 <- NULL > moms <- names(res)[grep("mom",names(res))] > res <- subset(res, method != "fresh") > > ## average moments across simulation replications > for(v in moms){ + res[,v] <- ave(res[,v], list(res$lambda, res$method)) + } > res <- unique(res[,c("lambda","method",moms)]) > res$index <- 1:nrow(res) > > ## get benchmark moments from gx cases and merge them onto res > b <- subset(res, method=="gx") > names(b)[which(names(b) %in% moms)] <- paste(names(b)[which(names(b) %in% moms)],"bench",sep="") > b$index <- NULL > b$method <- NULL > res <- merge(res, b, all.x=TRUE) > res <- res[order(res$index),] > > ## calculate percent bias > for(v in moms){ + res[,v] <- round(100 * (res[,v] - res[,paste(v,"bench",sep="")]) / res[,paste(v,"bench",sep="")], digits=0) + } > > ## reorganize table, get rid of irrelevant columns > res <- subset(res, method != "gx", select=c("method","lambda",moms)) > > ## resort to put the different lambdas adjacent for the same method > nest <- nrow(res)/2 > print(res[c(t(cbind(1:nest, (nest+1):nrow(res)))),]) method lambda mom2 mom3 mom4 mom5 mom6 mom7 mom8 2 aw 1 0 13 11 34 36 72 86 9 aw 2 0 17 16 46 56 101 131 3 avggx 1 0 -50 -16 -55 -38 -63 -58 10 avggx 2 0 -67 -25 -72 -56 -80 -76 4 rcal 1 0 0 0 -1 -1 -2 -2 11 rcal 2 0 0 -1 -2 -3 -5 -6 5 gw 1 1 17 14 41 44 85 102 12 gw 2 1 23 22 57 69 124 161 6 avggw 1 1 -50 -15 -54 -37 -63 -57 13 avggw 2 1 -67 -24 -71 -54 -80 -75 7 gexw 1 -1 -2 -1 -3 -3 -4 -4 14 gexw 2 -1 -2 -3 -4 -5 -7 -8 > > proc.time() user system elapsed 3384.991 703.986 4093.041