KLopt.lnorm {rodd}

Calculation of KL-optimal discriminating design for lognormal errors

Description

Calculates an approximation \xi^{**} of the KL-optimal design (in case of lognormal errors) \xi^* for discrimination between a given list of error densities \{f_i(x,\theta_i),\; i = 1,\dots,\nu\}. This procedure is based on the work [8]. This function mimics tpopt almost entirely. It is planed to combine tpopt and KLopt.lnorm in the future. See tpopt for the detailed description of the arguments marked with “-//-”.

Usage

KLopt.lnorm(    x, 
                w = rep(1, length(x)) / length(x), 
                eta, 
                sq.var,
                theta.fix, 
                theta.var = NULL, 
                p, 
                x.lb = min(x), 
                x.rb = max(x), 
                opt = list())

Arguments

x

-//-

w

-//-

eta

a list of means for the error densities \{f_i(x,\theta_i),\; i = 1,\dots,\nu\} between which proposed optimization should be performed. Every function from this list should be defined in the form of \eta_i(x,\theta_i), where x is one dimensional variable from \mathcal{X} and \theta_i is a vector of corresponding model parameters. We will refer to length of this list as \nu.

sq.var

a list of variances for the error densities \{f_i(x,\theta_i),\; i = 1,\dots,\nu\} between which proposed optimization should be performed. Every function from this list should be defined in the form of v^2_i(x,\theta_i). This list also has the length equal to \nu.

theta.fix

-//-

theta.var

-//-

p

-//-

x.lb

-//-

x.rb

-//-

opt

-//-

Value

Object of class “KLopt.lnorm” which contains the following fields:

x, w, efficiency, functional

-//-

eta

a list of means from the input.

sq.var

a list of variances from the input.

theta.fix, theta.var, p, x.lb, x.rb, max.iter, done.iter, des.eff, time

-//-

See Also

plot.KLopt.lnorm, summary.KLopt.lnorm, print.KLopt.lnorm

Examples

## Not run: 
### Examples from [8]
### Cases 1 and 3 are presented here; case 2 can be computed using the
### function tpopt (see the description of this function for exact example)
 
library(mvtnorm)

### Example 1 from [8]; EMAX vs MM

#List of models
eta.1 <- function(x, theta.1) 
    theta.1[1] * x + theta.1[2] * x / (x + theta.1[3])

eta.2 <- function(x, theta.2) 
    theta.2[1] * x / (x + theta.2[2])

eta <- list(eta.1, eta.2)

#List of fixed parameters
theta.1 <- c(1, 1, 1)
theta.2 <- c(1, 1)
theta.fix <- list(theta.1, theta.2)

#Comparison table
p <- matrix(
    c(
        0,1,
        0,0
    ), c(length(eta), length(eta)), byrow = TRUE)

### Case 1 

#List of variances
sq.var.1 <- function(x, theta.1)
    1
    
sq.var.2 <- function(x, theta.2)
    1

sq.var <- list(sq.var.1, sq.var.2)
    
#Case 1, method 1
res <- KLopt.lnorm(
    x = seq(0.1, 5, length.out = 10), 
    eta = eta, sq.var = sq.var, theta.fix = theta.fix, p = p,  
    opt = list(method = 1)
)
plot(res)
summary(res)

#Case 1, method 2
res <- KLopt.lnorm(
    x = seq(0.1, 5, length.out = 10), 
    eta = eta, sq.var = sq.var, theta.fix = theta.fix, p = p,  
    opt = list(method = 2)
)
plot(res)
summary(res)

### case 3
#List of variances
sq.var.1 <- function(x, theta.1)
    exp(eta.1(x, theta.1))
    
sq.var.2 <- function(x, theta.2)
    exp(eta.2(x, theta.2))

sq.var <- list(sq.var.1, sq.var.2)

#Case 3, method 1
res <- KLopt.lnorm(
    x = seq(0.1, 5, length.out = 10), 
    eta = eta, sq.var = sq.var, theta.fix = theta.fix, p = p,  
    opt = list(method = 1)
)
plot(res)
summary(res)

#Case 3, method 2
res <- KLopt.lnorm(
    x = seq(0.1, 5, length.out = 10), 
    eta = eta, sq.var = sq.var, theta.fix = theta.fix, p = p,  
    opt = list(method = 2)
)
plot(res)
summary(res)

### Example 2 from [8]; sigmoidal

#List of models
eta.1 = function(x, theta.1)
    theta.1[1] - theta.1[2] * exp(-theta.1[3] * x ^ theta.1[4])

eta.2 <- function(x, theta.2)
    theta.2[1] - theta.2[2] * exp(-theta.2[3] * x)

#List of fixed parameters
theta.1.mean <- c(2, 1, 0.8, 1.5)
sigma <- 0.3
theta.1.sigma <- matrix(
    c(
        sigma,0,
        0,sigma
    ), c(2, 2), byrow = TRUE)
grid <- expand.grid(
    theta.1.mean[1],
    theta.1.mean[2],
    seq(theta.1.mean[3] - sqrt(sigma), theta.1.mean[3] + sqrt(sigma), length.out = 5),
    seq(theta.1.mean[4] - sqrt(sigma), theta.1.mean[4] + sqrt(sigma), length.out = 5)
)

theta.2 <- c(2,1,1)

theta.fix <- list()
for(i in 1:length(grid[,1]))
    theta.fix[[length(theta.fix)+1]] <- as.numeric(grid[i,])
theta.fix[[length(theta.fix)+1]] <- theta.2

density.on.grid <- dmvnorm(grid[,3:4], mean = theta.1.mean[3:4], sigma = theta.1.sigma)
density.on.grid <- density.on.grid / sum(density.on.grid)
    
eta <- list()
for(i in 1:length(grid[,1]))
    eta <- c(eta, eta.1)
eta <- c(eta, eta.2)

#Comparison table
p <- rep(0,length(eta))
for(i in 1:length(grid[,1]))
    p <- rbind(p, c(rep(0,length(eta)-1), density.on.grid[i]))
p <- rbind(p, rep(0,length(eta)))
p <- p[-1,]

### Case 1

sq.var.1 <- function(x, theta.1)
    1

sq.var.2 <- function(x, theta.2)
    1

sq.var <- list()
for(i in 1:length(grid[,1]))
    sq.var <- c(sq.var, sq.var.1)
sq.var <- c(sq.var, sq.var.2)

#Case 1, method 1
res <- KLopt.lnorm(
    x = c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10), 
    eta = eta, sq.var = sq.var, theta.fix = theta.fix, p = p, 
    opt = list(method = 1)
)
plot(res)
summary(res)

#Case 1, method 2 
res <- KLopt.lnorm(
    x = c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10), 
    eta = eta, sq.var = sq.var, theta.fix = theta.fix, p = p, 
    opt = list(method = 2)
)
plot(res)
summary(res)

### Case 3

sq.var.1 <- function(x, theta.1)
    exp(eta.1(x, theta.1))

sq.var.2 <- function(x, theta.2)
    exp(eta.2(x, theta.2))

sq.var <- list()
for(i in 1:length(grid[,1]))
    sq.var <- c(sq.var, sq.var.1)
sq.var <- c(sq.var, sq.var.2)

#Case 3, method 1 
res <- KLopt.lnorm(
    x = c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10), 
    eta = eta, sq.var = sq.var, theta.fix = theta.fix, p = p, 
    opt = list(method = 1)
)
plot(res)
summary(res)

#Case 3, method 2 
res <- KLopt.lnorm(
    x = c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10), 
    eta = eta, sq.var = sq.var, theta.fix = theta.fix, p = p, 
    opt = list(method = 2)
)
plot(res)
summary(res)

### Example 3 from [8]; dose response

#List of models
eta.1 <- function(x, theta.1)
    theta.1[1] + theta.1[2] * x

eta.2 <- function(x, theta.2)
    theta.2[1] + theta.2[2] * x * (theta.2[3] - x)

eta.3 <- function(x, theta.3)
    theta.3[1] + theta.3[2] * x / (theta.3[3] + x)

eta.4 <- function(x, theta.4)
    theta.4[1] + theta.4[2] / (1 + exp((theta.4[3] - x) / theta.4[4]))

#List of fixed parameters
theta.1 <- c(60, 0.56)
theta.2 <- c(60, 7 / 2250, 600)
theta.3 <- c(60, 294, 25)
theta.4.mean <- c(49.62, 290.51, 150, 45.51)
a <- 45
b <- 20
grid <- expand.grid(
        c(theta.4.mean[1] - b, theta.4.mean[1], theta.4.mean[1] + a), 
        c(theta.4.mean[2] - b, theta.4.mean[2], theta.4.mean[2] + a), 
        c(theta.4.mean[3] - b, theta.4.mean[3], theta.4.mean[3] + a), 
        c(theta.4.mean[4] - b, theta.4.mean[4], theta.4.mean[4] + a)  
        )

eta <- list()
eta <- c(eta, eta.1, eta.2, eta.3)
for(i in 1:length(grid[,1]))
    eta <- c(eta, eta.4)

theta.fix <- list(theta.1, theta.2, theta.3)
for(i in 1:length(grid[,1]))
    theta.fix[[length(theta.fix) + 1]] <- as.numeric(grid[i,])

density.on.grid <- rep(1,length(grid[,1]))
density.on.grid <- density.on.grid / sum(density.on.grid)

#Comparison table
p <- rep(0, length(eta))
p <- rbind(p, c(1, rep(0, length(eta) - 1)))
p <- rbind(p, c(1, 1, rep(0,length(eta) - 2)))
for(i in 1:length(grid[,1]))
    p <- rbind(p, c(rep(density.on.grid[i], 3), rep(0, length(eta) - 3)))

### Case 1

#List of variances
sq.var.1 <- function(x, theta.1)
    1
    
sq.var.2 <- function(x, theta.2)
    1

sq.var.3 <- function(x, theta.3)
    1

sq.var.4 <- function(x, theta.4)
    1

sq.var <- list()
sq.var <- c(sq.var, sq.var.1, sq.var.2, sq.var.3)
for(i in 1:length(grid[,1]))
    sq.var <- c(sq.var, sq.var.4)

#Case 1, method 1

#Design estimation
res <- KLopt.lnorm(
    x = seq(0, 500, length.out = 10), 
    eta = eta, sq.var = sq.var, theta.fix = theta.fix, p = p, 
    opt = list(max.iter = 10)
)
plot(res)
summary(res)

#Case 1, method 2

#Design estimation
res <- KLopt.lnorm(
    x = seq(0, 500, length.out = 10), 
    eta = eta, sq.var = sq.var, theta.fix = theta.fix, p = p, 
    opt = list(
        method = 2, 
        max.iter = 10, 
        weights.evaluation.max.iter = 50, 
        support.epsilon = 1e-4
    )
)
plot(res)
summary(res)

### Case 3

#List of variances
sq.var.1 <- function(x, theta.1)
    exp(1e-2 * eta.1(x,theta.1))
    
sq.var.2 <- function(x, theta.2)
    exp(1e-2 * eta.2(x,theta.2))

sq.var.3 <- function(x, theta.3)
    exp(1e-2 * eta.3(x,theta.3))

sq.var.4 <- function(x, theta.4)
    exp(1e-2 * eta.4(x,theta.4))
    
sq.var <- list()
sq.var <- c(sq.var, sq.var.1, sq.var.2, sq.var.3)
for(i in 1:length(grid[,1]))
    sq.var <- c(sq.var, sq.var.4)

#Case 3, method 1
    
#Design estimation
res <- KLopt.lnorm(
    x = seq(0, 500, length.out = 10), 
    eta = eta, sq.var = sq.var, theta.fix = theta.fix, p = p, 
    opt = list(max.iter = 10)
)
plot(res)
summary(res)

#Case 3, method 2

eta.2 <- function(x, theta.2)
    theta.2[1] + theta.2[2] * x - theta.2[3] * x * x

theta.2 <- c(60, 7 * 600 / 2250, 7 / 2250)

eta <- list()
eta <- c(eta, eta.1, eta.2, eta.3)
for(i in 1:length(grid[,1]))
    eta <- c(eta, eta.4)

theta.fix <- list(theta.1, theta.2, theta.3)
for(i in 1:length(grid[,1]))
    theta.fix[[length(theta.fix) + 1]] <- as.numeric(grid[i,])

#Design estimation
res <- KLopt.lnorm(
    x = seq(0, 500, length.out = 10), 
    eta = eta, sq.var = sq.var, theta.fix = theta.fix, p = p, 
    opt = list(max.iter = 6, method = 2)
)
plot(res)
summary(res)

## End(Not run)