estimator <- function(X, htype, de, ae, ge, ive, itmax = 20000, k, 
                       ckpts = 2000,pl,pc1,pc2,pq)
 {
 #####
 ##
 ## The following function estimates weighted optimality criteria
 ##
 ## Function inputs are:
 ## X = k-variable design matrix;
 ## htype = 0 for weak heredity, 1 for strong heredity;
 ## Choose one or more weighted optimality criteria;
 ## Enter 1 for de if you want weighted D optimality, 0 otherwise;
 ## Enter 1 for ae if you want weighted A optimality, 0 otherwise;
 ## Enter 1 for ge if you want weighted G optimality, 0 otherwise;
 ## Enter 1 for ive if you want weighted IV optimality, 0 otherwise;
 ## itmax = number of estimated efficiency values;
 ## k = number of design variables;
 ## ckpts = number of random evaluation points in the hypercube;
 ## pl = probability x_i is in the model;
 ## pc1 = probability x_i x_j is in the model | x_i or x_j in model;
 ## pc2 = probability x_i x_j is in the model | x_i and x_j in model;
 ## pq = probability x_i^2 is in the model | x_i in model;
 ##
 #####
 
 n =  dim(X)[1] ## Obtain number of design points
 
 ## Generate a matrix of random values between -1 and 1 that will be used
 ## to create the random evaluation set for weighted IV efficiencies.
 
 if(ive == 1){xran <- 2*(matrix(runif(ckpts*k), ckpts, k) 
                      - 0.5*matrix(1, ckpts, k))}
 
 ## Generate a matrix of 3^k points with each coordinate in [-1, 1].
 ## This will be used as the evaluation set for weighted G efficiencies.
 
 if(ge == 1)
 {
 pts <- list()
 for(i in 1:k){pts[[i]] <- c(-1, 0, 1)}
 m <- expand.grid(pts)
 }
 pr <- matrix(NA, nrow = itmax, ncol = 4)
 
 for(iter in 1:itmax)
 {
 xvec <- matrix(0, k)
 xqvec <- matrix(0, k)
 xcmat <- matrix(0, k - 1, k - 1)
 spv <- NULL
 gspv <- NULL
 
 prl <- matrix(runif(k), k) ## Random vector of U(0,1) deviates for x_1 terms
 prc <- matrix(runif((k - 1)^2), k - 1, k - 1) 
                            ## Random matrix of U(0,1)deviates for x_i x_j
terms
 prq <- matrix(runif(k), k) ## Random vector of U(0,1) deviates for x_1^2
terms
 prm <- 1 ## Initial number of terms
 
 if(ive == 1){xiv <- matrix(1, ckpts)}
 if(ge == 1){xg <- matrix(1, 3^k)}
 xd <- matrix(1, n) ## Form intercept column for the design matrix
 
 ## Initialize efficiencies
 deff <- 1
 aeff <- 1
 geff <- 1
 iveff <- 1
 
 for(i in 1:k)
 {
 if(prl[i] <= pl) ## Determine if x_i is in the model
 {
 xvec[i] <- 1
 xd <- cbind(xd, X[,i])
 prm <- prm + 1
 if(ive == 1){xiv <- cbind(xiv, xran[,i])} ## For IV only
 if(ge == 1){xg <- cbind(xg, m[,i])} ## For G only
 if(prq[i] <= pq) ## Determine if x_i^2 is in the model
 {
 xqvec[i] <- 1
 prm <- prm + 1
 }
 }
 }
 for(i in 1:(k - 1)) ## Determine if x_i x_j is in the model
 {
 for(j in i:(k - 1))
 {
 if(xvec[i] == 1 & xvec[j + 1] == 1)
 {
 if(prc[i, j] <= pc2)
 {
 xcmat[i, j] <- 1
 prm <- prm + 1
 }
 }
 if(htype == 0) ## For weak heredity only
 {
 if(xvec[i] == 1 & xvec[j + 1] == 0)
 {
 if(prc[i, j] <= pc1)
 {
 xcmat[i, j] <- 1
 prm <- prm + 1
 }
 }
 if(xvec[i] == 0 & xvec[j + 1] == 1)
 {
 if(prc[i, j] <= pc1)
 {
 xcmat[i, j] <- 1
 prm <- prm + 1
 }
 }
 }
 }
 }
 for(i in 1:(k - 1))
 {
 for(j in i:(k - 1))
 {
 if(xcmat[i, j] == 1)
 {
 ctmp <- X[,i]*X[, j + 1] ## x_i x_j column for model matrix
 if(ive == 1) ## For IV only
 {
 xctmp <- xran[, i]*xran[, j + 1] ## x_i x_j column for IV evaluation set
 }
 if(ge == 1) ## For G only
 {
 gctmp <- m[, i]*m[, j + 1] ## x_i x_j column for G evaluation set
 }
 xd <- cbind(xd, ctmp) ## Add x_i x_j column to model matrix
 if(ive == 1) ## For IV only
 {
 xiv <- cbind(xiv, xctmp) ## Add x_i x_j column to IV evaluation set
 }
 if(ge == 1) ## For G only
 {
 xg <- cbind(xg, gctmp) ## Add x_i x_j column to G evaluation set
 }
 }
 }
 }
 for(i in 1:k)
 {
 if(xqvec[i] == 1)
 {
 qtmp <- X[, i]*X[, i] ## x_i^2 column for model matrix
 xd <- cbind(xd, qtmp) ## Add x_i^2 column to model matrix
 if(ive == 1) ## For IV only
 {
 xqtmp <- xran[, i]*xran[, i] ## x_i^2 column for IV evaluation set
 xiv <- cbind(xiv, xqtmp) ## Add x_i^2 column to IV evaluation set
 }
 if(ge == 1) ## For G only
 {
 gqtmp <- m[, i]*m[, i] ## x_i^2 column for G evaluation set
 xg <- cbind(xg, gqtmp) ## Add x_i^2 column to G evaluation set
 }
 }
 }
 if(de == 1)
 {
 deff <- 100*(det(t(xd) %*% xd)^(1/prm))/n ## Calculate D-efficiency
 }
 if(ae + ge + ive > 0)
 {
 xx <- t(xd) %*% xd
 xxi <- solve(xx) ## Find inverse of X'X
 }
 if(ae == 1)
 {
 aeff <- 100*prm/(sum(diag(n*xxi))) ## Calculate A-efficiency
 }
 if(ive == 1)
 {
 for(i in 1:ckpts)
 {
 spv[i] <- n * t(matrix(xiv[i,])) %*% xxi %*% matrix(xiv[i,]) 
           ## Scaled prediction variance for the ith evaluation point
 }        
 iveff <-mean(spv) ## Estimate IV-criterion
 }
 if(ge == 1)
 {
 for(g in 1:3^k)
 {
 gspv[g] <- n * t(matrix(xg[g,])) %*% xxi %*% matrix(xg[g,]) 
            ## Scaled prediction variance for the ith evaluation point
 }                                                           
 geff <- 100*prm/max(gspv) ## Estimate G-efficiency
 }
 options(warn = -1)
 rm(ctmp, qtmp, xd, xx, xxi)
 pr[iter, 1:4] <- c(deff, aeff, geff, iveff)
 }
 adwgt <- apply(pr, 2, mean) ## Mean of estimated efficiencies
 sewgt <- apply(pr, 2, sd)/sqrt(itmax)  
                             ## Standard error of estimated efficiencies
 if(de == 1)
 {
 cat("Weighted D and SE ", c(adwgt[1], sewgt[1]), "\n")
 }
 if(ae == 1)
 {
 cat("Weighted A and SE ", c(adwgt[2], sewgt[2]), "\n")
 }
 if(ge == 1)
 {
 cat("Weighted G and SE ", c(adwgt[3], sewgt[3]), "\n")
 }
 if(ive == 1)
 {
 cat("Weighted IV and SE ", c(adwgt[4], sewgt[4]), "\n")
 }
 }
 ##### Begin entering design matrices
 
 ## For Box-Draper, N = 10
 dmat <- matrix(c(-1, -1, -1, 1, -1, -1, -1, 1, -1, -1, -1,  1, 
 -1, 0.1925, 0.1925, 0.1925, -1,  0.1925, 0.1925, 0.1925, -1, 
  -0.2912, 1, 1, 1, -0.2912, 1, 1, 1, -0.2912), nrow = 10, byrow = TRUE)
 
 ##### CHECK THE DESIGN!
 
 ## Weak 
 estimator(dmat, htype = 0, de = 1, ae = 1, ge = 1, ive = 1, itmax = 20000, 
           k =3, ckpts = 2000, pl = 0.95, pc1 = 0.40, pc2 = 0.65, pq = 0.75)
 
 ## Strong
 estimator(dmat, htype = 1, de = 1, ae = 1, ge = 1, ive = 1, itmax = 20000, 
           k =3, ckpts = 2000, pl = 0.80, pc1 = 0.00, pc2 = 0.35, pq = 0.35)