OptimalDesign (version 0.2)

od.m1: Optimal design for models with a one-dimensional parameter

Description

Computes an optimal approximate or exact design under general linear constraints for a model with a one-dimensional parameter.

Usage

od.m1(F, b, A=NULL, w0=NULL, type="exact", kappa=1e-9, 
        tab=NULL, graph=NULL, t.max=120)

Arguments

F

The n times 1 matrix of regressors corresponding to n design points and 1 model parameter.

b, A

The real vector of length k and the k times n matrix of reals numbers. The linear constraints A%*%w<=b, w0<=w define the set of permissible designs w (where w0 is a described below.) The argument A can also be NULL; in that case b must be a positive number and A is set to the 1 times n matrix of ones.

w0

The non-negative vector of length n representing the design to be augmented. This argument can also be NULL; in that case, w0 is set to the vector of zeros.

type

Specifies whether exact or approximate design is to be computed.

kappa

A small non-negative perturbation parameter.

tab

A vector determining the regressor components to be printed with the resulting design. This argument should be a subvector of 1:n, or a subvector of colnames(F), or it can be NULL. If tab=NULL, the design is not printed.

graph

A vector determining the regressor components to be plotted with the resulting design. This argument should be a subvector of 1:n, or a subvector of colnames(F), or it can be NULL. If graph=NULL, the resulting design is not visualized.

t.max

The time limit for the computation.

Value

A list with the following components:

method

The method used for computing the design w.best.

w.best

the best permissible design found, or NULL. The value of w.best will be NULL if the computation fails. This can happen, if no permissible solution is found within the time limit, no permissible solution exists, or the problem is unbounded; see the status variable for more details. Note that even if w.best is a permissible design, then it still can have a singular information matrix; cf. the Phi.best variable.

Phi.best

The value of the criterion of optimality of the design w.best. If w.best has a singular information matrix or if the computation fails, the value of Phi.best will be 0.

status

The status variable of the gurobi optimization procedure; see the gurobi solver documentation for details.

t.act

The actual time taken by the computation.

Details

For the case of m=1 parameter, the problem of optimal design is much simpler than for m>=2. First, for m=1 all standardized information criteria coincide, therefore the D-, A-, and IV-optimal designs are the same, and we can call them just "optimal". Second, the information matrix is a real number, therefore we can call it just "information". Under the most common size constraint, the optimal design is any design supported on the set of maxima of F[1,1]^2,...,F[n,1]^2, which is straightforward to find. Under a non-standard linear constraint, however, the problem becomes a less trivial knapsack problem, which is here solved by the integer linear programming solver of gurobi.

The model should be non-singular in the sense that there exists an exact design w satisfying the constraints 0<=w0<=w and A%*%w<=b, with a nonzero information. If this requirement is not satisfied, the computation may fail, or produce a deficient design.

If the criterion of IV-optimality is selected, the region R should be chosen such that the associated matrix L (in this case a real number; see the help page of the function od.crit) is non-zero. If this requirement is not satisfied, the computation may fail, or it may produce a deficient design.

The perturbation parameter kappa can be used to add n*m iid random numbers from the uniform distribution in [-kappa,kappa] to the elements of F before the optimization is executed. This can be helpful for increasing the numerical stability of the computation or for generating a random design from the potentially large set of optimal or nearly-optimal designs.

The performance depends on the problem and on the hardware used, but in most cases the function can compute an optimal exact design for a problem with a thousand design points within seconds of computing time.

Examples

# NOT RUN {
if(require("gurobi")){
# We will demonstrate the procedure on a simple randomly generated 
# knapsack problem. Here, the squares of elements of F correspond 
# to the values of n available items, the elements of the 1 times n 
# matrix A correspond to the weights of the n items, and the real 
# number b is the upper limit on the total weight of the items that 
# can  be put into the knapsack. 
# The resulting binary "optimal design" determines which of the items 
# should we put into the knapsack to steal the highest possible value.

n <- 200 # There are this many items to choose from.
F.square <- matrix(sample(1:10, n, replace=TRUE), ncol=1) 
# Generate random prices of items.
A <- matrix(sample(1:10, n, replace=TRUE), nrow=1) 
# Generate random weights of items in kgs.
A <- rbind(A, diag(n))  
# We assume that there is just one copy of each item.
b <- c(n / 4, rep(1,n)) 
# The capacity of the knapsack is n/4 kgs.

# Compute the optimal design, which in this case determines how many 
# (0 or 1) of each of the n items should we put into the knapsack.
od.m1(sqrt(F.square), b, A)

# Note: one can compare the result with a specialized function 
# as follows:
# library(adagio); knapsack(A[1,], F.square[,1], n / 4)

# However, od.m1 is more general than the standard knapsack functions.
# Suppose, for instance, that the uncle asks that we must be sure to 
# take the items number 1, 13 and 66. We will compute the most valuable 
# selection of items that fit into our knapsack and contain the 
# required items.
w0 <- rep(0, n); w0[c(1, 13, 66)] <- 1
od.m1(sqrt(F.square), b, A, w0, t.max=2)
}
# }