[R] logistic transformation using nlminb
John Pitchard
johnpitchard at googlemail.com
Thu May 15 14:41:01 CEST 2008
Dear all,
I want to find the optimal values of a vector, x (with 6 elements)
say, satisfying the following conditions:
1. for all x>=0
2. sum(x)=1
3. x[5]<=0.5 and x[6]<=0.5
For the minimisation I'm using nlminb and to satisfy the first 2
conditions the logistic transformation is used with box constraints
for condition 3. However, I don't seem to be able to get the values x
that i'm expecting from fitting a simpler model. For example,
set.seed(0)
# creating a correlation matrix
corr <- diag(5)
corr[lower.tri(corr)] <- 0.5
corr[upper.tri(corr)] <- 0.5
# Data for the minimisation
mat <- rmvnorm(10000, mean=c(3, 2.75, 2.8, 3, 2.9), sd=c(0.1,
0.1,0.12, 0.15, 0.10), cov=corr)
# here is the simple optimisation function that allows the 5th
element to be potentially negative
obj.funA <- function(opt, mat) {
opt <- c(opt, 1-sum(opt))
LinearComb <- mat%*%opt
obj <- -min(LinearComb)
obj
}
# I want to put an upper boundary constraint on the first variable -
not being greater than 0.35 and the rest being between 0 and 1
opt <- nlminb(rep(0,4), lower=rep(0,4), upper=c(0.35, 1, 1, 1) ,
obj.funA, mat=mat)
opt.x <- opt$parameters
opt.x <- c(opt.x, 1-sum(opt.x))
opt.x
# Result
[1] 0.34999902 0.06475651 0.00000000 0.16561760 0.41962686
However, I don't get the same result from the logistic transformation
obj.funB <- function(opt, mat) {
opt <- c(opt, 0)
opt <- exp(opt)/sum(exp(opt))
LinearComb <- mat%*%opt
obj <- -min(LinearComb)
obj
}
opt <- nlminb(rep(0,4), upper=c(log(0.35), NA, NA, NA), obj.funB, mat=mat)
opt.x <- opt$parameters
opt.x <- c(opt.x, 0)
opt.x <- exp(opt.x)/sum(exp(opt.x))
opt.x
# Result
[1] 2.355325e-001 6.339398e-009 1.202751e-004 9.139718e-002 6.729500e-001
I don't know how to obtain the same answer for both optimisations. In
reality, my own optimisation typically gives negative values for the standard
minimisation- so I have no choice but to use a more advanced method.
Also, there appears to be a dependency between the first and the last,
i.e. 2.355325e-001/6.729500e-001 =0.35
Does anyone know why the logistic doesn't give the same answer as the
simpler method?
Any help is much appreciated.
Regards,
John
More information about the R-help
mailing list