R/DiscrimOD.R
In PingYangChen/DiscrimOD: DiscrimOD: Finding Optimal Discrimination Designs by Hybridizing Particle Swarm and L-BFGS Algorithms
Defines functions
emptyModelList
equivalence
designCriterion
DiscrimOD
Documented in
designCriterion
DiscrimOD
emptyModelList
equivalence
#' Hybrid Algorithm of PSO and L-BFGS for Finding Optimal Discrimination Design
#'
#' Implement the hybridized PSO algorithm to find the optimal discrimination designs
#' for 2 or more than 2 competing models
#'
#' @param MODEL_INFO A list of information of competing models.
#' For details, run \code{emptyModelList()} and see the examples below.
#' @param DISTANCE The function of distance measure. See the example.
#' @param nSupp The integer number fo support points.
#' This number should be larger than or equal to 2.
#' @param dsLower The vector of finite lower bounds of the design space.
#' Its length should be equal to the dimension of design space.
#' @param dsUpper The vector of finite upper bounds of the design space.
#' Its length should be equal to the dimension of design space.
#' @param minWt A small numerical value. The weight of resulting design should
#' be larger than this value. This option is useful to avoid a degenerated resulting design.
#' Its default is zero.
#' @param crit_type The name of the case of the discrimination design problem. See details.
#' @param MaxMinStdVals The vector of values of demoninators in the design efficiency
#' calculation for finding max-min discrimination design.
#' @param PSO_INFO The list of PSO parameters generated by \code{getPSOInfo()}.
#' @param LBFGS_INFO The list of L-BFGS parameters generated by \code{getLBFGSInfo()}.
#' @param seed The random seed that controls initial swarm of PSO. The default is \code{NULL}.
#' @param verbose The logical value controls if PSO would reports the updating progress. The default is \code{TRUE}.
#' @return An List.
#' \itemize{
#' \item{BESTDESIGN}{ the resulting design in matrix form. Each row is a support point with its weight.
#' The last column is the weights of the support points.}
#' \item{BESTVAL}{ the design criterion value of the resulting design.}
#' \item{GBESTHIST}{ a vector of design citerion values of the global best particle in PSO search history.}
#' \item{CPUTIME}{ the computational time in seconds.}
#' }
#' @details
#' Currently the \code{DiscrimOD} accepts the option \code{crit_type} to be
#' \code{'pair_fixed_true'} or \code{'maxmin_fixed_true'}.
#' If specified \code{crit_type = 'pair_fixed_true'}, the \code{DiscrimOD}
#' searches for the optimal discrimination design for two competing models
#' such as \eqn{T}-optimal design (Atkinson and Fedorov, 1975a) and \eqn{KL}-optimal
#' design (Lopez-Fidalgo et al., 2007).
#' If specified \code{crit_type = 'maxmin_fixed_true'}, the \code{DiscrimOD}
#' searches for the max-min type of optimal discrimination design
#' which is able to simultaneously discriminate among many competing models.
#' For example, the max-min \eqn{KL}-optimal design in Tommasi et al. (2016).
#'
#' @examples
#' # Atkinson and Fedorov (1975a): T-optimal
#' # Two R functions of competing models are given by
#' af1 <- function(x, p) p[1] + p[2]*exp(x) + p[3]*exp(-x)
#' af2 <- function(x, p) p[1] + p[2]*x + p[3]*x^2
#' # Set the model information
#' # The nominla value in 'm1' is 4.5, -1.5, -2.0
#' # For 'm2', we set the parameter space to be [-10, 10]^3 and
#' # the initial guess (for LBFGS) of the rival model parameter is zero vector
#' AF_para_af1 <- c(4.5, -1.5, -2)
#' af_info_12 <- list(
#' # The first list should be the true model and the specified nominal values
#' list(model = af1, para = AF_para_af1),
#' # Then the rival models are listed accordingly. We also need to specify the model space.
#' list(model = af2, paraLower = rep(-10, 3), paraUpper = rep(10, 3))
#' )
#' # Define the R function for the distance measure in T-optimal criterion
#' # xt is the mean values of the true model
#' # xr is the mean values of the rival model
#' sq_diff <- function(xt, xr) (xt - xr)^2
#'
#' # Initialize PSO and BFGS options
#' PSO_INFO <- getPSOInfo(nSwarm = 32, maxIter = 100)
#' LBFGS_INFO <- getLBFGSInfo(LBFGS_RETRY = 2)
#'
#' # Find T-optimal design for models af1 and af2
#' af_res_12 <- DiscrimOD(MODEL_INFO = af_info_12, DISTANCE = sq_diff,
#' nSupp = 4, dsLower = -1, dsUpper = 1, crit_type = "pair_fixed_true",
#' PSO_INFO = PSO_INFO, LBFGS_INFO = LBFGS_INFO, seed = NULL, verbose = FALSE)
#'
#' round(af_res_12$BESTDESIGN, 3) # The resulting design
#' af_res_12$BESTVAL # The T-optimal criterion value
#' af_res_12$CPUTIME # CPU time
#'
#' # Test optimality by equivalence theorem
#' af_eqv_12 <- equivalence(ngrid = 100, PSO_RESULT = af_res_12, MODEL_INFO = af_info_12,
#' DISTANCE = sq_diff, dsLower = -1, dsUpper = 1, crit_type = "pair_fixed_true",
#' PSO_INFO = PSO_INFO, LBFGS_INFO = LBFGS_INFO)
#'
#' # Draw the directional derivative curve
#' plot(af_eqv_12$Grid_1, af_eqv_12$DirDeriv, type = "l", col = "blue",
#' main = "af_res_12", xlab = "x", ylab = "Directional Derivative"); abline(h = 0)
#' points(af_res_12$BESTDESIGN[,1], rep(0, nrow(af_res_12$BESTDESIGN)), pch = 16)
#'
#' # Following above, we add the 3rd model in Atkinson and Fedorov (1975b)
#' af3 <- function(x, p) p[1] + p[2]*sin(0.5*pi*x) + p[3]*cos(0.5*pi*x) + p[4]*sin(pi*x)
#' af_info_13 <- list(
#' list(model = af1, para = AF_para_af1),
#' list(model = af3, paraLower = rep(-10, 4), paraUpper = rep(10, 4))
#' )
#' # Find another T-optimal design for models af1 and af3
#' af_res_13 <- DiscrimOD(MODEL_INFO = af_info_13, DISTANCE = sq_diff,
#' nSupp = 5, dsLower = -1.0, dsUpper = 1.0, crit_type = "pair_fixed_true",
#' PSO_INFO = PSO_INFO, LBFGS_INFO = LBFGS_INFO, seed = NULL, verbose = FALSE)
#'
#' # Re-organize model list for finding the max-min T-optimal design
#' af_info_maxmin <- list(af_info_12[[1]], af_info_12[[2]], af_info_13[[2]])
#' # Define the vector of optimal criterion values for efficiency computations
#' af_vals_pair <- c(af_res_12$BESTVAL, af_res_13$BESTVAL)
#'
#' # Search for max-min T-optimal design for discriminating af1, af2 and af3
#' af_res_maxmin <- DiscrimOD(MODEL_INFO = af_info_maxmin, DISTANCE = sq_diff,
#' nSupp = 5, dsLower = -1, dsUpper = 1, crit_type = "maxmin_fixed_true",
#' MaxMinStdVals = af_vals_pair, PSO_INFO = PSO_INFO, LBFGS_INFO = LBFGS_INFO,
#' seed = NULL, verbose = FALSE)
#'
#' round(af_res_maxmin$BESTDESIGN, 3) # The resulting design
#' af_res_maxmin$BESTVAL # The T-optimal criterion value
#' af_res_maxmin$CPUTIME # CPU time
#'
#' # Test optimality by equivalence theorem
#' af_eqv_maxmin <- equivalence(ngrid = 100, PSO_RESULT = af_res_maxmin,
#' MODEL_INFO = af_info_maxmin, DISTANCE = sq_diff, dsLower = -1, dsUpper = 1,
#' crit_type = "maxmin_fixed_true", MaxMinStdVals = af_vals_pair,
#' PSO_INFO = PSO_INFO, LBFGS_INFO = LBFGS_INFO)
#'
#' af_eqv_maxmin$alpha # The weight of efficiency values
#'
#' # Draw the directional derivative curve
#' plot(af_eqv_maxmin$Grid_1, af_eqv_maxmin$DirDeriv, type = "l", col = "blue",
#' main = "af_res_maxmin", xlab = "x", ylab = "Directional Derivative"); abline(h = 0)
#' points(af_res_maxmin$BESTDESIGN[,1], rep(0, nrow(af_res_maxmin$BESTDESIGN)), pch = 16)
#'
#' # Not Run
#' # For other distance measure, here are a few
#'
#' # Heteroscedastic Normal model
#' # heter_norm <- function(xt, xr) {
#' # var_t <- xt^2
#' # var_r <- xr^2
#' # (var_t + (xt - xr)^2)/var_r - log(var_t/var_r)
#' # }
#'
#' # Logistic regression model
#' # logit_diff <- function(xt, xr) {
#' # exp_t <- exp(xt)
#' # exp_r <- exp(xr)
#' # mu_t <- exp_t/(1 + exp_t)
#' # mu_r <- exp_r/(1 + exp_r)
#' # mu_t*(log(mu_t) - log(mu_r)) + (1 - mu_t)*(log(1.0 - mu_t) - log(1.0 - mu_r))
#' # }
#'
#' # Gamma regression model
#' # gamma_diff <- function(xt, xr) log(xr/xt) + (xt - xr)/xr
#'
#' @references Atkinson, A. C. and Fedorov, V. V. (1975a). The design of experiments for discriminating between two rival models. Biometrika, 62(1):57-70.
#' @references Atkinson, A. C. and Fedorov, V. V. (1975b). Optimal design: experiments for discriminating between several models. Biometrika, 62(2):289-303.
#' @references Lopez-Fidalgo, J., Tommasi, C., and Trandafir, P. C. (2007). An optimal experimental design criterion for discriminating between non-normal models. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 69(2):231-242.
#' @references Tommasi, C., Martin-Martin, R., and Lopez-Fidalgo, J. (2016). Max-min optimal discriminating designs for several statistical models. Statistics and Computing, 26(6):1163-1172.
#' @name DiscrimOD
#' @rdname DiscrimOD
#' @export
DiscrimOD <- function(MODEL_INFO, DISTANCE, nSupp, dsLower, dsUpper, minWt = 0.0, crit_type = "pair_fixed_true", MaxMinStdVals = NULL,
PSO_INFO = NULL, LBFGS_INFO = NULL, seed = NULL, verbose = TRUE, environment, ...) {
stopifnot(nSupp >= 2L, all(is.finite(dsLower)), all(is.finite(dsUpper)),
length(dsLower) == length(dsUpper), all(dsUpper > dsLower),
all(names(PSO_INFO) == names(getPSOInfo())),
crit_type %in% c("pair_fixed_true", "maxmin_fixed_true"))
MODEL_LIST <- lapply(1:length(MODEL_INFO), function(k) MODEL_INFO[[k]]$model)
if (crit_type == "maxmin_fixed_true") {
if (is.null(MaxMinStdVals)) {
stop("Need optimal values for pairwise discrimination in efficiency computations.")
}
if (length(MaxMinStdVals) < (length(MODEL_INFO) - 1)) {
stop(paste0("Need ", length(MODEL_INFO) - 1, " values for pairwise discrimination in efficiency computations."))
}
} else {
if (is.null(MaxMinStdVals)) MaxMinStdVals <- 0
}
dSupp <- length(dsLower)
D_INFO <- getDesignInfo(D_TYPE = "approx", MODEL_INFO = MODEL_INFO, dist_func = DISTANCE,
crit_type = crit_type, MaxMinStdVals = MaxMinStdVals,
dSupp = length(dsLower), nSupp = nSupp, dsLower = dsLower, dsUpper = dsUpper,
minWt = minWt)
if (is.null(PSO_INFO)) {
PSO_INFO <- getPSOInfo(nSwarm = c(32, 32), maxIter = c(100, 100))
if (verbose) message(paste0("Use the default settings for PSO. See '?getPSOInfo'."))
}
if (is.null(LBFGS_INFO)) {
LBFGS_INFO <- getLBFGSInfo()
if (verbose) message(paste0("Use the default settings for LBFGS. See '?getLBFGSInfo'."))
} else if (is.numeric(LBFGS_INFO) && (LBFGS_INFO == 0)) {
LBFGS_INFO <- getLBFGSInfo(IF_INNER_LBFGS = FALSE)
if (verbose) message(paste0("Use NestedPSO, not PSO-QN."))
stopifnot(length(PSO_INFO$nSwarm) < 2)
}
if (!hasArg(environment)) environment <- new.env()
# Adjust PSO_INFO according to D_INFO
swarmSetting <- algInfoUpdate(D_INFO)
PSO_INFO$varUpper <- matrix(swarmSetting$UB, length(PSO_INFO$nSwarm), ncol(swarmSetting$UB), byrow = TRUE)
PSO_INFO$varLower <- matrix(swarmSetting$LB, length(PSO_INFO$nSwarm), ncol(swarmSetting$UB), byrow = TRUE)
PSO_INFO$dSwarm <- rep(ncol(swarmSetting$UB), length(PSO_INFO$nSwarm))
# Find Initial Guess for parameters of rival model
propGuess <- NULL
for (i in 1:dSupp) {
propGuess <- cbind(propGuess, seq(dsLower[i], dsUpper[i], length = nSupp))
}
propGuess <- cbind(propGuess, 1/nSupp)
UNIFDESIGN_M <- designV2M(propGuess, D_INFO)
iniGuess <- cppDesignCriterion(PSO_INFO, LBFGS_INFO, D_INFO, MODEL_LIST, 0, environment, UNIFDESIGN_M)
D_INFO$parasInit <- iniGuess$theta2
# Start
set.seed(seed)
cputime <- system.time(
psoOut <- cppPSO(0, PSO_INFO, LBFGS_INFO, D_INFO, MODEL_LIST, 0, environment, FALSE, verbose)
)[3]
if (verbose) message(paste0("CPU time: ", round(cputime, 2), " seconds."))
BESTDESIGN <- designM2V(psoOut$GBest, D_INFO)
list(BESTDESIGN = BESTDESIGN, BESTVAL = -psoOut$fGBest, GBESTHIST = -psoOut$fGBestHist,
CPUTIME = cputime)
}
#' Function of criterion value computing
#'
#' Computes the discrimination design criterion value
#'
#' @param DESIGN1 matrix. The approximate design.
#' @param MODEL_INFO list of information of competing models.
#' For details, run \code{emptyModelList()} and see the instruction below.
#' @param DISTANCE function. The R/C++ function of distance measure.
#' For T-optimal design, the function is the squared difference of two models means.
#' @param dsLower vector. The finite lower bounds of the design space.
#' Its length should be equal to the dimension of design space.
#' @param dsUpper vector. The finite upper bounds of the design space.
#' Its length should be equal to the dimension of design space.
#' @param crit_type string. The name of the case of the discrimination design problem.
#' The default is 'pair_fixed_true'.
#' @param MaxMinStdVals vector. The values of demoninators in the design efficiency calculation for finding max-min discrimination design.
#' @param PSO_INFO list. PSO and BFGS options.
#' @return An List.
#' \itemize{
#' \item{cri_val}{ the design criterion value of the resulting design in the PSO result or the given design.}
#' \item{theta2}{ a matrix of the parameters of each input model that result the design criterion value.}
#' }
#' @name designCriterion
#' @rdname designCriterion
#' @export
designCriterion <- function(DESIGN1, MODEL_INFO, DISTANCE, dsLower, dsUpper, crit_type = "pair_fixed_true", MaxMinStdVals = NULL,
PSO_INFO = NULL, LBFGS_INFO = NULL, environment, ...) {
stopifnot(all(is.finite(dsLower)), all(is.finite(dsUpper)),
length(dsLower) == length(dsUpper), all(dsUpper > dsLower))
nSupp <- nrow(DESIGN1)
dSupp <- ncol(DESIGN1) - 1
MODEL_LIST <- lapply(1:length(MODEL_INFO), function(k) MODEL_INFO[[k]]$model)
if (is.null(MaxMinStdVals)) MaxMinStdVals <- 0
D_INFO <- getDesignInfo(D_TYPE = "approx", MODEL_INFO = MODEL_INFO, dist_func = DISTANCE,
crit_type = crit_type, MaxMinStdVals = MaxMinStdVals, minWt = .0,
dSupp = length(dsLower), nSupp = nSupp, dsLower = dsLower, dsUpper = dsUpper)
if (is.null(PSO_INFO)) { PSO_INFO <- getPSOInfo(nSwarm = c(32, 32), maxIter = c(100, 100)) }
if (is.null(LBFGS_INFO)) {
LBFGS_INFO <- getLBFGSInfo()
} else if (is.numeric(LBFGS_INFO) && (LBFGS_INFO == 0)) {
LBFGS_INFO <- getLBFGSInfo(IF_INNER_LBFGS = FALSE)
}
if (!hasArg(environment)) environment <- new.env()
# Adjust PSO_INFO according to D_INFO
swarmSetting <- algInfoUpdate(D_INFO)
PSO_INFO$varUpper <- matrix(swarmSetting$UB, length(PSO_INFO$nSwarm), ncol(swarmSetting$UB), byrow = TRUE)
PSO_INFO$varLower <- matrix(swarmSetting$LB, length(PSO_INFO$nSwarm), ncol(swarmSetting$UB), byrow = TRUE)
PSO_INFO$dSwarm <- rep(ncol(swarmSetting$UB), length(PSO_INFO$nSwarm))
# Find Initial Guess for parameters of rival model
propGuess <- NULL
for (i in 1:dSupp) {
propGuess <- cbind(propGuess, seq(dsLower[i], dsUpper[i], length = nSupp))
}
propGuess <- cbind(propGuess, 1/nSupp)
UNIFDESIGN_M <- designV2M(propGuess, D_INFO)
iniGuess <- cppDesignCriterion(PSO_INFO, LBFGS_INFO, D_INFO, MODEL_LIST, 0, environment, UNIFDESIGN_M)
D_INFO$parasInit <- iniGuess$theta2
# Compute the criterion value
DESIGN1_M <- designV2M(DESIGN1, D_INFO)
cri_1 <- cppDesignCriterion(PSO_INFO, LBFGS_INFO, D_INFO, MODEL_LIST, 0, environment, DESIGN1_M)
rownames(cri_1$theta2) <- paste0("model_", 1:length(MODEL_INFO))
return(list(cri_val = -cri_1$val, theta2 = cri_1$theta2))
}
#' Equivalence theorem for discrimination design
#'
#' Computes the values of directional derivative function on an user-defined grid and
#' returns its maximal value.
#'
#' @param DESIGN1 matrix.
#' @param PSO_RESULT list.
#' @param ngrid integer. The default is 100.
#' @param IFPLOT logical.
#' @param MODEL_INFO list of information of competing models. For details, run \code{emptyModelList()} and see the instruction below.
#' @param DISTANCE function.
#' @param dsLower vector. The finite lower bounds of the design space. Its length should be equal to the dimension of design space.
#' @param dsUpper vector. The finite upper bounds of the design space. Its length should be equal to the dimension of design space.
#' @param crit_type string. The name of the case of the discrimination design problem. The default is 'pair_fixed_true'.
#' @param MaxMinStdVals vector. The values of demoninators in the design efficiency calculation for finding max-min discrimination design.
#' @param PSO_INFO list. PSO and BFGS options.
#' @return An List.
#' \itemize{
#' \item{Grid_1}{ a vector of grid points on the design space. For 2-D case, this vector is the grid of the first dimension.}
#' \item{Grid_2}{ a vector of grid points on the second dimension of the 2-D design space. For 1_d case, this output can be ignored.}
#' \item{DirDeriv}{ a vector (1-D case) or a matrix (2-D case) of the directional derivative function values on the grid points.}
#' \item{MAX_DD}{ the maximal value of the directional derivative function in \code{DirDeriv}.}
#' \item{alpha}{ a vector of weights of efficiency values that satisfies the equivalence theorem for max-min optimal discrimination design.}
#' }
#' @name equivalence
#' @rdname equivalence
#' @export
equivalence <- function(DESIGN = NULL, PSO_RESULT = NULL, ngrid = 100, IFPLOT = FALSE,
MODEL_INFO, DISTANCE, dsLower, dsUpper, crit_type = "pair_fixed_true",
MaxMinStdVals = NULL, PSO_INFO = NULL, LBFGS_INFO = NULL,
ALPHA_PSO_INFO = NULL, environment, ...) {
stopifnot(all(is.finite(dsLower)), all(is.finite(dsUpper)),
length(dsLower) == length(dsUpper), all(dsUpper > dsLower))
if (is.null(DESIGN)) DESIGN <- PSO_RESULT$BESTDESIGN
nSupp <- nrow(DESIGN)
dSupp <- ncol(DESIGN) - 1
MODEL_LIST <- lapply(1:length(MODEL_INFO), function(k) MODEL_INFO[[k]]$model)
if (!hasArg(environment)) environment <- new.env()
if (is.null(MaxMinStdVals)) MaxMinStdVals <- 0
D_INFO <- getDesignInfo(D_TYPE = "approx", MODEL_INFO = MODEL_INFO, dist_func = DISTANCE,
crit_type = crit_type, MaxMinStdVals = MaxMinStdVals, minWt = .0,
dSupp = length(dsLower), nSupp = nSupp, dsLower = dsLower, dsUpper = dsUpper)
if (is.null(PSO_INFO)) { PSO_INFO <- getPSOInfo(nSwarm = c(32, 32), maxIter = c(100, 100)) }
if (is.null(LBFGS_INFO)) {
LBFGS_INFO <- getLBFGSInfo()
} else if (is.numeric(LBFGS_INFO) && (LBFGS_INFO == 0)) {
LBFGS_INFO <- getLBFGSInfo(IF_INNER_LBFGS = FALSE)
}
# Adjust PSO_INFO according to D_INFO
swarmSetting <- algInfoUpdate(D_INFO)
PSO_INFO$varUpper <- matrix(swarmSetting$UB, length(PSO_INFO$nSwarm), ncol(swarmSetting$UB), byrow = TRUE)
PSO_INFO$varLower <- matrix(swarmSetting$LB, length(PSO_INFO$nSwarm), ncol(swarmSetting$UB), byrow = TRUE)
PSO_INFO$dSwarm <- rep(ncol(swarmSetting$UB), length(PSO_INFO$nSwarm))
# Find Initial Guess for parameters of rival model
propGuess <- NULL
for (i in 1:dSupp) {
propGuess <- cbind(propGuess, seq(dsLower[i], dsUpper[i], length = nSupp))
}
propGuess <- cbind(propGuess, 1/nSupp)
UNIFDESIGN_M <- designV2M(propGuess, D_INFO)
iniGuess <- cppDesignCriterion(PSO_INFO, LBFGS_INFO, D_INFO, MODEL_LIST, 0, environment, UNIFDESIGN_M)
D_INFO$parasInit <- iniGuess$theta2
# Compute for the equivalence theorem
DESIGN_M <- designV2M(DESIGN, D_INFO)
CRIT_VAL <- cppDesignCriterion(PSO_INFO, LBFGS_INFO, D_INFO, MODEL_LIST, 0, environment, DESIGN_M)
PARA_SET <- CRIT_VAL$theta2
ALPHA <- 0
if (crit_type == "maxmin_fixed_true") {
#message("Looking for best weight...")
# Find the weight vector first
ALPHA_INFO <- getDesignInfo(D_TYPE = "maxmin_eqv_wt", MODEL_INFO = MODEL_INFO, dist_func = DISTANCE,
crit_type = crit_type, MaxMinStdVals = MaxMinStdVals, minWt = .0,
dSupp = length(dsLower), nSupp = nSupp, dsLower = dsLower, dsUpper = dsUpper)
ALPHA_INFO$paras <- PARA_SET
if (is.null(ALPHA_PSO_INFO)) { ALPHA_PSO_INFO <- getPSOInfo(nSwarm = 64, maxIter = 200) }
swarmSetting <- algInfoUpdate(ALPHA_INFO)
ALPHA_PSO_INFO$varUpper <- matrix(swarmSetting$UB, length(ALPHA_PSO_INFO$nSwarm), ncol(swarmSetting$UB), byrow = TRUE)
ALPHA_PSO_INFO$varLower <- matrix(swarmSetting$LB, length(ALPHA_PSO_INFO$nSwarm), ncol(swarmSetting$UB), byrow = TRUE)
ALPHA_PSO_INFO$dSwarm <- rep(ncol(swarmSetting$UB), length(ALPHA_PSO_INFO$nSwarm))
dimnames(DESIGN) <- NULL
EXTERNAL_LIST <- list(DESIGN = as.matrix(DESIGN[,-ncol(DESIGN)], nSupp, dSupp), CRIT_VAL = -CRIT_VAL$val)
tmp <- getLBFGSInfo()
psoOut <- cppPSO(0, ALPHA_PSO_INFO, tmp, ALPHA_INFO, MODEL_LIST, EXTERNAL_LIST, environment, FALSE, FALSE)
ALPHA <- designM2V(psoOut$GBest, ALPHA_INFO)
}
equiv <- cppEquivalence(D_INFO, MODEL_LIST, -CRIT_VAL$val, PARA_SET, ALPHA, environment, ngrid)
if (crit_type == "maxmin_fixed_true") { equiv$alpha <- ALPHA }
return(equiv)
}
#' Create An Empty Model List
#'
#' Run \code{emptyModelList()} for instruction of model settings.
#'
#' @param N_model Integer. The numebr of competing models including the true model. The default value is 2.
#' @return An List.
#' \itemize{
#' \item{True}{ a list of informaiton of true model.}
#' \itemize{
#' \item{model}{ the model formula. Input the R function or xPtr function pointer generated by \code{cxxfunction} in the \code{inline} package.}
#' \item{para}{ the nominal values of parameeters. Input a vector of parameter values.}
#' }
#' \item{Rival\code{k}}{ a vector of grid points on the second dimension of the 2-D design space. For 1_d case, this output can be ignored.}
#' \itemize{
#' \item{model}{ the model formula. Input the R function or xPtr function pointer generated by \code{cxxfunction} in the \code{inline} package.}
#' \item{paraLower}{ the lower bound of teh parameters in the k-th rival model. Input a vector of lower bound values.}
#' \item{paraUpper}{ the upper bound of teh parameters in the k-th rival model. Input a vector of upper bound values.}
#' \item{paraInit}{ the initial guess of parameters for the LBFGS algorithm. Input a vector of parameter values.}
#' }
#' }
#' @examples
#' emptyModelList() # True model and one rival model
#' emptyModelList(3) # True model and two rival models
#' emptyModelList(5) # True model and four rival models
#' @name emptyModelList
#' @rdname emptyModelList
#' @export
emptyModelList <- function(N_model = 2) {
out <- lapply(1:N_model, function(k) {
if (k == 1) list(model = 'R or C++ Function for True Model', para = 'Nominal Values of Parameters in True Model')
else list(model = paste0('R or C++ Function for Rival ', k-1),
paraLower = paste0('Lower Bound of Parameter Space for Rival ', k-1),
paraUpper = paste0('Upper Bound of Parameter Space for Rival ', k-1))
})
names(out) <- c("True", paste0("Rival", 1:(N_model - 1)))
return(out)
}
PingYangChen/DiscrimOD documentation built on Jan. 30, 2022, 5:25 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions emptyModelList equivalence designCriterion DiscrimOD
Documented in designCriterion DiscrimOD emptyModelList equivalence
#' Hybrid Algorithm of PSO and L-BFGS for Finding Optimal Discrimination Design
#'
#' Implement the hybridized PSO algorithm to find the optimal discrimination designs
#' for 2 or more than 2 competing models
#'
#' @param MODEL_INFO A list of information of competing models.
#' For details, run \code{emptyModelList()} and see the examples below.
#' @param DISTANCE The function of distance measure. See the example.
#' @param nSupp The integer number fo support points.
#' This number should be larger than or equal to 2.
#' @param dsLower The vector of finite lower bounds of the design space.
#' Its length should be equal to the dimension of design space.
#' @param dsUpper The vector of finite upper bounds of the design space.
#' Its length should be equal to the dimension of design space.
#' @param minWt A small numerical value. The weight of resulting design should
#' be larger than this value. This option is useful to avoid a degenerated resulting design.
#' Its default is zero.
#' @param crit_type The name of the case of the discrimination design problem. See details.
#' @param MaxMinStdVals The vector of values of demoninators in the design efficiency
#' calculation for finding max-min discrimination design.
#' @param PSO_INFO The list of PSO parameters generated by \code{getPSOInfo()}.
#' @param LBFGS_INFO The list of L-BFGS parameters generated by \code{getLBFGSInfo()}.
#' @param seed The random seed that controls initial swarm of PSO. The default is \code{NULL}.
#' @param verbose The logical value controls if PSO would reports the updating progress. The default is \code{TRUE}.
#' @return An List.
#' \itemize{
#' \item{BESTDESIGN}{ the resulting design in matrix form. Each row is a support point with its weight.
#' The last column is the weights of the support points.}
#' \item{BESTVAL}{ the design criterion value of the resulting design.}
#' \item{GBESTHIST}{ a vector of design citerion values of the global best particle in PSO search history.}
#' \item{CPUTIME}{ the computational time in seconds.}
#' }
#' @details
#' Currently the \code{DiscrimOD} accepts the option \code{crit_type} to be
#' \code{'pair_fixed_true'} or \code{'maxmin_fixed_true'}.
#' If specified \code{crit_type = 'pair_fixed_true'}, the \code{DiscrimOD}
#' searches for the optimal discrimination design for two competing models
#' such as \eqn{T}-optimal design (Atkinson and Fedorov, 1975a) and \eqn{KL}-optimal
#' design (Lopez-Fidalgo et al., 2007).
#' If specified \code{crit_type = 'maxmin_fixed_true'}, the \code{DiscrimOD}
#' searches for the max-min type of optimal discrimination design
#' which is able to simultaneously discriminate among many competing models.
#' For example, the max-min \eqn{KL}-optimal design in Tommasi et al. (2016).
#'
#' @examples
#' # Atkinson and Fedorov (1975a): T-optimal
#' # Two R functions of competing models are given by
#' af1 <- function(x, p) p[1] + p[2]*exp(x) + p[3]*exp(-x)
#' af2 <- function(x, p) p[1] + p[2]*x + p[3]*x^2
#' # Set the model information
#' # The nominla value in 'm1' is 4.5, -1.5, -2.0
#' # For 'm2', we set the parameter space to be [-10, 10]^3 and
#' # the initial guess (for LBFGS) of the rival model parameter is zero vector
#' AF_para_af1 <- c(4.5, -1.5, -2)
#' af_info_12 <- list(
#' # The first list should be the true model and the specified nominal values
#' list(model = af1, para = AF_para_af1),
#' # Then the rival models are listed accordingly. We also need to specify the model space.
#' list(model = af2, paraLower = rep(-10, 3), paraUpper = rep(10, 3))
#' )
#' # Define the R function for the distance measure in T-optimal criterion
#' # xt is the mean values of the true model
#' # xr is the mean values of the rival model
#' sq_diff <- function(xt, xr) (xt - xr)^2
#'
#' # Initialize PSO and BFGS options
#' PSO_INFO <- getPSOInfo(nSwarm = 32, maxIter = 100)
#' LBFGS_INFO <- getLBFGSInfo(LBFGS_RETRY = 2)
#'
#' # Find T-optimal design for models af1 and af2
#' af_res_12 <- DiscrimOD(MODEL_INFO = af_info_12, DISTANCE = sq_diff,
#' nSupp = 4, dsLower = -1, dsUpper = 1, crit_type = "pair_fixed_true",
#' PSO_INFO = PSO_INFO, LBFGS_INFO = LBFGS_INFO, seed = NULL, verbose = FALSE)
#'
#' round(af_res_12$BESTDESIGN, 3) # The resulting design
#' af_res_12$BESTVAL # The T-optimal criterion value
#' af_res_12$CPUTIME # CPU time
#'
#' # Test optimality by equivalence theorem
#' af_eqv_12 <- equivalence(ngrid = 100, PSO_RESULT = af_res_12, MODEL_INFO = af_info_12,
#' DISTANCE = sq_diff, dsLower = -1, dsUpper = 1, crit_type = "pair_fixed_true",
#' PSO_INFO = PSO_INFO, LBFGS_INFO = LBFGS_INFO)
#'
#' # Draw the directional derivative curve
#' plot(af_eqv_12$Grid_1, af_eqv_12$DirDeriv, type = "l", col = "blue",
#' main = "af_res_12", xlab = "x", ylab = "Directional Derivative"); abline(h = 0)
#' points(af_res_12$BESTDESIGN[,1], rep(0, nrow(af_res_12$BESTDESIGN)), pch = 16)
#'
#' # Following above, we add the 3rd model in Atkinson and Fedorov (1975b)
#' af3 <- function(x, p) p[1] + p[2]*sin(0.5*pi*x) + p[3]*cos(0.5*pi*x) + p[4]*sin(pi*x)
#' af_info_13 <- list(
#' list(model = af1, para = AF_para_af1),
#' list(model = af3, paraLower = rep(-10, 4), paraUpper = rep(10, 4))
#' )
#' # Find another T-optimal design for models af1 and af3
#' af_res_13 <- DiscrimOD(MODEL_INFO = af_info_13, DISTANCE = sq_diff,
#' nSupp = 5, dsLower = -1.0, dsUpper = 1.0, crit_type = "pair_fixed_true",
#' PSO_INFO = PSO_INFO, LBFGS_INFO = LBFGS_INFO, seed = NULL, verbose = FALSE)
#'
#' # Re-organize model list for finding the max-min T-optimal design
#' af_info_maxmin <- list(af_info_12[[1]], af_info_12[[2]], af_info_13[[2]])
#' # Define the vector of optimal criterion values for efficiency computations
#' af_vals_pair <- c(af_res_12$BESTVAL, af_res_13$BESTVAL)
#'
#' # Search for max-min T-optimal design for discriminating af1, af2 and af3
#' af_res_maxmin <- DiscrimOD(MODEL_INFO = af_info_maxmin, DISTANCE = sq_diff,
#' nSupp = 5, dsLower = -1, dsUpper = 1, crit_type = "maxmin_fixed_true",
#' MaxMinStdVals = af_vals_pair, PSO_INFO = PSO_INFO, LBFGS_INFO = LBFGS_INFO,
#' seed = NULL, verbose = FALSE)
#'
#' round(af_res_maxmin$BESTDESIGN, 3) # The resulting design
#' af_res_maxmin$BESTVAL # The T-optimal criterion value
#' af_res_maxmin$CPUTIME # CPU time
#'
#' # Test optimality by equivalence theorem
#' af_eqv_maxmin <- equivalence(ngrid = 100, PSO_RESULT = af_res_maxmin,
#' MODEL_INFO = af_info_maxmin, DISTANCE = sq_diff, dsLower = -1, dsUpper = 1,
#' crit_type = "maxmin_fixed_true", MaxMinStdVals = af_vals_pair,
#' PSO_INFO = PSO_INFO, LBFGS_INFO = LBFGS_INFO)
#'
#' af_eqv_maxmin$alpha # The weight of efficiency values
#'
#' # Draw the directional derivative curve
#' plot(af_eqv_maxmin$Grid_1, af_eqv_maxmin$DirDeriv, type = "l", col = "blue",
#' main = "af_res_maxmin", xlab = "x", ylab = "Directional Derivative"); abline(h = 0)
#' points(af_res_maxmin$BESTDESIGN[,1], rep(0, nrow(af_res_maxmin$BESTDESIGN)), pch = 16)
#'
#' # Not Run
#' # For other distance measure, here are a few
#'
#' # Heteroscedastic Normal model
#' # heter_norm <- function(xt, xr) {
#' # var_t <- xt^2
#' # var_r <- xr^2
#' # (var_t + (xt - xr)^2)/var_r - log(var_t/var_r)
#' # }
#'
#' # Logistic regression model
#' # logit_diff <- function(xt, xr) {
#' # exp_t <- exp(xt)
#' # exp_r <- exp(xr)
#' # mu_t <- exp_t/(1 + exp_t)
#' # mu_r <- exp_r/(1 + exp_r)
#' # mu_t*(log(mu_t) - log(mu_r)) + (1 - mu_t)*(log(1.0 - mu_t) - log(1.0 - mu_r))
#' # }
#'
#' # Gamma regression model
#' # gamma_diff <- function(xt, xr) log(xr/xt) + (xt - xr)/xr
#'
#' @references Atkinson, A. C. and Fedorov, V. V. (1975a). The design of experiments for discriminating between two rival models. Biometrika, 62(1):57-70.
#' @references Atkinson, A. C. and Fedorov, V. V. (1975b). Optimal design: experiments for discriminating between several models. Biometrika, 62(2):289-303.
#' @references Lopez-Fidalgo, J., Tommasi, C., and Trandafir, P. C. (2007). An optimal experimental design criterion for discriminating between non-normal models. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 69(2):231-242.
#' @references Tommasi, C., Martin-Martin, R., and Lopez-Fidalgo, J. (2016). Max-min optimal discriminating designs for several statistical models. Statistics and Computing, 26(6):1163-1172.
#' @name DiscrimOD
#' @rdname DiscrimOD
#' @export
DiscrimOD <- function(MODEL_INFO, DISTANCE, nSupp, dsLower, dsUpper, minWt = 0.0, crit_type = "pair_fixed_true", MaxMinStdVals = NULL,
PSO_INFO = NULL, LBFGS_INFO = NULL, seed = NULL, verbose = TRUE, environment, ...) {
stopifnot(nSupp >= 2L, all(is.finite(dsLower)), all(is.finite(dsUpper)),
length(dsLower) == length(dsUpper), all(dsUpper > dsLower),
all(names(PSO_INFO) == names(getPSOInfo())),
crit_type %in% c("pair_fixed_true", "maxmin_fixed_true"))
MODEL_LIST <- lapply(1:length(MODEL_INFO), function(k) MODEL_INFO[[k]]$model)
if (crit_type == "maxmin_fixed_true") {
if (is.null(MaxMinStdVals)) {
stop("Need optimal values for pairwise discrimination in efficiency computations.")
}
if (length(MaxMinStdVals) < (length(MODEL_INFO) - 1)) {
stop(paste0("Need ", length(MODEL_INFO) - 1, " values for pairwise discrimination in efficiency computations."))
}
} else {
if (is.null(MaxMinStdVals)) MaxMinStdVals <- 0
}
dSupp <- length(dsLower)
D_INFO <- getDesignInfo(D_TYPE = "approx", MODEL_INFO = MODEL_INFO, dist_func = DISTANCE,
crit_type = crit_type, MaxMinStdVals = MaxMinStdVals,
dSupp = length(dsLower), nSupp = nSupp, dsLower = dsLower, dsUpper = dsUpper,
minWt = minWt)
if (is.null(PSO_INFO)) {
PSO_INFO <- getPSOInfo(nSwarm = c(32, 32), maxIter = c(100, 100))
if (verbose) message(paste0("Use the default settings for PSO. See '?getPSOInfo'."))
}
if (is.null(LBFGS_INFO)) {
LBFGS_INFO <- getLBFGSInfo()
if (verbose) message(paste0("Use the default settings for LBFGS. See '?getLBFGSInfo'."))
} else if (is.numeric(LBFGS_INFO) && (LBFGS_INFO == 0)) {
LBFGS_INFO <- getLBFGSInfo(IF_INNER_LBFGS = FALSE)
if (verbose) message(paste0("Use NestedPSO, not PSO-QN."))
stopifnot(length(PSO_INFO$nSwarm) < 2)
}
if (!hasArg(environment)) environment <- new.env()
# Adjust PSO_INFO according to D_INFO
swarmSetting <- algInfoUpdate(D_INFO)
PSO_INFO$varUpper <- matrix(swarmSetting$UB, length(PSO_INFO$nSwarm), ncol(swarmSetting$UB), byrow = TRUE)
PSO_INFO$varLower <- matrix(swarmSetting$LB, length(PSO_INFO$nSwarm), ncol(swarmSetting$UB), byrow = TRUE)
PSO_INFO$dSwarm <- rep(ncol(swarmSetting$UB), length(PSO_INFO$nSwarm))
# Find Initial Guess for parameters of rival model
propGuess <- NULL
for (i in 1:dSupp) {
propGuess <- cbind(propGuess, seq(dsLower[i], dsUpper[i], length = nSupp))
}
propGuess <- cbind(propGuess, 1/nSupp)
UNIFDESIGN_M <- designV2M(propGuess, D_INFO)
iniGuess <- cppDesignCriterion(PSO_INFO, LBFGS_INFO, D_INFO, MODEL_LIST, 0, environment, UNIFDESIGN_M)
D_INFO$parasInit <- iniGuess$theta2
# Start
set.seed(seed)
cputime <- system.time(
psoOut <- cppPSO(0, PSO_INFO, LBFGS_INFO, D_INFO, MODEL_LIST, 0, environment, FALSE, verbose)
)[3]
if (verbose) message(paste0("CPU time: ", round(cputime, 2), " seconds."))
BESTDESIGN <- designM2V(psoOut$GBest, D_INFO)
list(BESTDESIGN = BESTDESIGN, BESTVAL = -psoOut$fGBest, GBESTHIST = -psoOut$fGBestHist,
CPUTIME = cputime)
}
#' Function of criterion value computing
#'
#' Computes the discrimination design criterion value
#'
#' @param DESIGN1 matrix. The approximate design.
#' @param MODEL_INFO list of information of competing models.
#' For details, run \code{emptyModelList()} and see the instruction below.
#' @param DISTANCE function. The R/C++ function of distance measure.
#' For T-optimal design, the function is the squared difference of two models means.
#' @param dsLower vector. The finite lower bounds of the design space.
#' Its length should be equal to the dimension of design space.
#' @param dsUpper vector. The finite upper bounds of the design space.
#' Its length should be equal to the dimension of design space.
#' @param crit_type string. The name of the case of the discrimination design problem.
#' The default is 'pair_fixed_true'.
#' @param MaxMinStdVals vector. The values of demoninators in the design efficiency calculation for finding max-min discrimination design.
#' @param PSO_INFO list. PSO and BFGS options.
#' @return An List.
#' \itemize{
#' \item{cri_val}{ the design criterion value of the resulting design in the PSO result or the given design.}
#' \item{theta2}{ a matrix of the parameters of each input model that result the design criterion value.}
#' }
#' @name designCriterion
#' @rdname designCriterion
#' @export
designCriterion <- function(DESIGN1, MODEL_INFO, DISTANCE, dsLower, dsUpper, crit_type = "pair_fixed_true", MaxMinStdVals = NULL,
PSO_INFO = NULL, LBFGS_INFO = NULL, environment, ...) {
stopifnot(all(is.finite(dsLower)), all(is.finite(dsUpper)),
length(dsLower) == length(dsUpper), all(dsUpper > dsLower))
nSupp <- nrow(DESIGN1)
dSupp <- ncol(DESIGN1) - 1
MODEL_LIST <- lapply(1:length(MODEL_INFO), function(k) MODEL_INFO[[k]]$model)
if (is.null(MaxMinStdVals)) MaxMinStdVals <- 0
D_INFO <- getDesignInfo(D_TYPE = "approx", MODEL_INFO = MODEL_INFO, dist_func = DISTANCE,
crit_type = crit_type, MaxMinStdVals = MaxMinStdVals, minWt = .0,
dSupp = length(dsLower), nSupp = nSupp, dsLower = dsLower, dsUpper = dsUpper)
if (is.null(PSO_INFO)) { PSO_INFO <- getPSOInfo(nSwarm = c(32, 32), maxIter = c(100, 100)) }
if (is.null(LBFGS_INFO)) {
LBFGS_INFO <- getLBFGSInfo()
} else if (is.numeric(LBFGS_INFO) && (LBFGS_INFO == 0)) {
LBFGS_INFO <- getLBFGSInfo(IF_INNER_LBFGS = FALSE)
}
if (!hasArg(environment)) environment <- new.env()
# Adjust PSO_INFO according to D_INFO
swarmSetting <- algInfoUpdate(D_INFO)
PSO_INFO$varUpper <- matrix(swarmSetting$UB, length(PSO_INFO$nSwarm), ncol(swarmSetting$UB), byrow = TRUE)
PSO_INFO$varLower <- matrix(swarmSetting$LB, length(PSO_INFO$nSwarm), ncol(swarmSetting$UB), byrow = TRUE)
PSO_INFO$dSwarm <- rep(ncol(swarmSetting$UB), length(PSO_INFO$nSwarm))
# Find Initial Guess for parameters of rival model
propGuess <- NULL
for (i in 1:dSupp) {
propGuess <- cbind(propGuess, seq(dsLower[i], dsUpper[i], length = nSupp))
}
propGuess <- cbind(propGuess, 1/nSupp)
UNIFDESIGN_M <- designV2M(propGuess, D_INFO)
iniGuess <- cppDesignCriterion(PSO_INFO, LBFGS_INFO, D_INFO, MODEL_LIST, 0, environment, UNIFDESIGN_M)
D_INFO$parasInit <- iniGuess$theta2
# Compute the criterion value
DESIGN1_M <- designV2M(DESIGN1, D_INFO)
cri_1 <- cppDesignCriterion(PSO_INFO, LBFGS_INFO, D_INFO, MODEL_LIST, 0, environment, DESIGN1_M)
rownames(cri_1$theta2) <- paste0("model_", 1:length(MODEL_INFO))
return(list(cri_val = -cri_1$val, theta2 = cri_1$theta2))
}
#' Equivalence theorem for discrimination design
#'
#' Computes the values of directional derivative function on an user-defined grid and
#' returns its maximal value.
#'
#' @param DESIGN1 matrix.
#' @param PSO_RESULT list.
#' @param ngrid integer. The default is 100.
#' @param IFPLOT logical.
#' @param MODEL_INFO list of information of competing models. For details, run \code{emptyModelList()} and see the instruction below.
#' @param DISTANCE function.
#' @param dsLower vector. The finite lower bounds of the design space. Its length should be equal to the dimension of design space.
#' @param dsUpper vector. The finite upper bounds of the design space. Its length should be equal to the dimension of design space.
#' @param crit_type string. The name of the case of the discrimination design problem. The default is 'pair_fixed_true'.
#' @param MaxMinStdVals vector. The values of demoninators in the design efficiency calculation for finding max-min discrimination design.
#' @param PSO_INFO list. PSO and BFGS options.
#' @return An List.
#' \itemize{
#' \item{Grid_1}{ a vector of grid points on the design space. For 2-D case, this vector is the grid of the first dimension.}
#' \item{Grid_2}{ a vector of grid points on the second dimension of the 2-D design space. For 1_d case, this output can be ignored.}
#' \item{DirDeriv}{ a vector (1-D case) or a matrix (2-D case) of the directional derivative function values on the grid points.}
#' \item{MAX_DD}{ the maximal value of the directional derivative function in \code{DirDeriv}.}
#' \item{alpha}{ a vector of weights of efficiency values that satisfies the equivalence theorem for max-min optimal discrimination design.}
#' }
#' @name equivalence
#' @rdname equivalence
#' @export
equivalence <- function(DESIGN = NULL, PSO_RESULT = NULL, ngrid = 100, IFPLOT = FALSE,
MODEL_INFO, DISTANCE, dsLower, dsUpper, crit_type = "pair_fixed_true",
MaxMinStdVals = NULL, PSO_INFO = NULL, LBFGS_INFO = NULL,
ALPHA_PSO_INFO = NULL, environment, ...) {
stopifnot(all(is.finite(dsLower)), all(is.finite(dsUpper)),
length(dsLower) == length(dsUpper), all(dsUpper > dsLower))
if (is.null(DESIGN)) DESIGN <- PSO_RESULT$BESTDESIGN
nSupp <- nrow(DESIGN)
dSupp <- ncol(DESIGN) - 1
MODEL_LIST <- lapply(1:length(MODEL_INFO), function(k) MODEL_INFO[[k]]$model)
if (!hasArg(environment)) environment <- new.env()
if (is.null(MaxMinStdVals)) MaxMinStdVals <- 0
D_INFO <- getDesignInfo(D_TYPE = "approx", MODEL_INFO = MODEL_INFO, dist_func = DISTANCE,
crit_type = crit_type, MaxMinStdVals = MaxMinStdVals, minWt = .0,
dSupp = length(dsLower), nSupp = nSupp, dsLower = dsLower, dsUpper = dsUpper)
if (is.null(PSO_INFO)) { PSO_INFO <- getPSOInfo(nSwarm = c(32, 32), maxIter = c(100, 100)) }
if (is.null(LBFGS_INFO)) {
LBFGS_INFO <- getLBFGSInfo()
} else if (is.numeric(LBFGS_INFO) && (LBFGS_INFO == 0)) {
LBFGS_INFO <- getLBFGSInfo(IF_INNER_LBFGS = FALSE)
}
# Adjust PSO_INFO according to D_INFO
swarmSetting <- algInfoUpdate(D_INFO)
PSO_INFO$varUpper <- matrix(swarmSetting$UB, length(PSO_INFO$nSwarm), ncol(swarmSetting$UB), byrow = TRUE)
PSO_INFO$varLower <- matrix(swarmSetting$LB, length(PSO_INFO$nSwarm), ncol(swarmSetting$UB), byrow = TRUE)
PSO_INFO$dSwarm <- rep(ncol(swarmSetting$UB), length(PSO_INFO$nSwarm))
# Find Initial Guess for parameters of rival model
propGuess <- NULL
for (i in 1:dSupp) {
propGuess <- cbind(propGuess, seq(dsLower[i], dsUpper[i], length = nSupp))
}
propGuess <- cbind(propGuess, 1/nSupp)
UNIFDESIGN_M <- designV2M(propGuess, D_INFO)
iniGuess <- cppDesignCriterion(PSO_INFO, LBFGS_INFO, D_INFO, MODEL_LIST, 0, environment, UNIFDESIGN_M)
D_INFO$parasInit <- iniGuess$theta2
# Compute for the equivalence theorem
DESIGN_M <- designV2M(DESIGN, D_INFO)
CRIT_VAL <- cppDesignCriterion(PSO_INFO, LBFGS_INFO, D_INFO, MODEL_LIST, 0, environment, DESIGN_M)
PARA_SET <- CRIT_VAL$theta2
ALPHA <- 0
if (crit_type == "maxmin_fixed_true") {
#message("Looking for best weight...")
# Find the weight vector first
ALPHA_INFO <- getDesignInfo(D_TYPE = "maxmin_eqv_wt", MODEL_INFO = MODEL_INFO, dist_func = DISTANCE,
crit_type = crit_type, MaxMinStdVals = MaxMinStdVals, minWt = .0,
dSupp = length(dsLower), nSupp = nSupp, dsLower = dsLower, dsUpper = dsUpper)
ALPHA_INFO$paras <- PARA_SET
if (is.null(ALPHA_PSO_INFO)) { ALPHA_PSO_INFO <- getPSOInfo(nSwarm = 64, maxIter = 200) }
swarmSetting <- algInfoUpdate(ALPHA_INFO)
ALPHA_PSO_INFO$varUpper <- matrix(swarmSetting$UB, length(ALPHA_PSO_INFO$nSwarm), ncol(swarmSetting$UB), byrow = TRUE)
ALPHA_PSO_INFO$varLower <- matrix(swarmSetting$LB, length(ALPHA_PSO_INFO$nSwarm), ncol(swarmSetting$UB), byrow = TRUE)
ALPHA_PSO_INFO$dSwarm <- rep(ncol(swarmSetting$UB), length(ALPHA_PSO_INFO$nSwarm))
dimnames(DESIGN) <- NULL
EXTERNAL_LIST <- list(DESIGN = as.matrix(DESIGN[,-ncol(DESIGN)], nSupp, dSupp), CRIT_VAL = -CRIT_VAL$val)
tmp <- getLBFGSInfo()
psoOut <- cppPSO(0, ALPHA_PSO_INFO, tmp, ALPHA_INFO, MODEL_LIST, EXTERNAL_LIST, environment, FALSE, FALSE)
ALPHA <- designM2V(psoOut$GBest, ALPHA_INFO)
}
equiv <- cppEquivalence(D_INFO, MODEL_LIST, -CRIT_VAL$val, PARA_SET, ALPHA, environment, ngrid)
if (crit_type == "maxmin_fixed_true") { equiv$alpha <- ALPHA }
return(equiv)
}
#' Create An Empty Model List
#'
#' Run \code{emptyModelList()} for instruction of model settings.
#'
#' @param N_model Integer. The numebr of competing models including the true model. The default value is 2.
#' @return An List.
#' \itemize{
#' \item{True}{ a list of informaiton of true model.}
#' \itemize{
#' \item{model}{ the model formula. Input the R function or xPtr function pointer generated by \code{cxxfunction} in the \code{inline} package.}
#' \item{para}{ the nominal values of parameeters. Input a vector of parameter values.}
#' }
#' \item{Rival\code{k}}{ a vector of grid points on the second dimension of the 2-D design space. For 1_d case, this output can be ignored.}
#' \itemize{
#' \item{model}{ the model formula. Input the R function or xPtr function pointer generated by \code{cxxfunction} in the \code{inline} package.}
#' \item{paraLower}{ the lower bound of teh parameters in the k-th rival model. Input a vector of lower bound values.}
#' \item{paraUpper}{ the upper bound of teh parameters in the k-th rival model. Input a vector of upper bound values.}
#' \item{paraInit}{ the initial guess of parameters for the LBFGS algorithm. Input a vector of parameter values.}
#' }
#' }
#' @examples
#' emptyModelList() # True model and one rival model
#' emptyModelList(3) # True model and two rival models
#' emptyModelList(5) # True model and four rival models
#' @name emptyModelList
#' @rdname emptyModelList
#' @export
emptyModelList <- function(N_model = 2) {
out <- lapply(1:N_model, function(k) {
if (k == 1) list(model = 'R or C++ Function for True Model', para = 'Nominal Values of Parameters in True Model')
else list(model = paste0('R or C++ Function for Rival ', k-1),
paraLower = paste0('Lower Bound of Parameter Space for Rival ', k-1),
paraUpper = paste0('Upper Bound of Parameter Space for Rival ', k-1))
})
names(out) <- c("True", paste0("Rival", 1:(N_model - 1)))
return(out)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.