Skip to content

Lower bounds on selfStart function not working

4 messages · Schatzi

#
I adapted a selfStart function and the lower bounds are not working. The
parameter "b" is negative, whereas I would like the lower bound to be zero.
Any ideas? Thanks.

Here is my code (I am still figuring out how to easily make replicable
examples):
A<-1.75
mu<-.2
l<-2
b<-0
x<-seq(0,18,.25)
create.y<-function(x){
y<-b+A/(1+exp(4*mu/A*(l-x)+2))
return(y)
}
ys<-create.y(x)
yvec<-(rep(ys,5))*(.95+runif(length(x)*5)/10)
Trt<-factor(c(rep("A1",length(x)),rep("A2",length(x)),rep("A3",length(x)),rep("A4",length(x)),rep("A5",length(x))))
Data<-data.frame(Trt,rep(x,5),yvec)
names(Data)<-c("Trt","x","y")

NewData<-groupedData(y~x|Trt,data=Data)

powrDpltInit <-
function(mCall, LHS, data) {
   xy     <- sortedXyData(mCall[["x"]],LHS,data)
   A.s  <- max(xy$y)-min(xy$y)
   mu.s  <- A.s/7.5
   l.s <- 0
   b.s    <- max(min(xy$y),0.00001)
   value  <- c(A.s, l.s, mu.s, b.s)

#function to optimize
func1 <- function(value) {
   A.s  <- value[1]
   mu.s  <- value[2]
   l.s <- value[3]
   b.s    <- value[4]
  y1<-rep(0,length(xy$x)) # generate vector for predicted y (y1) to evaluate
against observed y
  for(cnt in 1:length(xy$x)){
  y1[cnt]<- b.s+A.s/(1+exp(4*mu.s/A.s*(l.s-x[cnt])+2))}  #predicting y1 for
values of y
  evl<-sum((xy$y-y1)^2) #sum of squares is function to minimize
  return(evl)}

#optimizing
oppar<-optim(c(A.s , mu.s , l.s , b.s),func1,method="L-BFGS-B",
lower=c(0.0001,0.0,0.0,0.0),
control=list(maxit=2000,trace=TRUE))

#saving optimized parameters
value<-c(oppar$par[1L],oppar$par[2L],oppar$par[3L],oppar$par[4L])

   names(value) <- mCall[c("A","mu","l","b")]
   value

}

SSpowrDplt<-selfStart(~b+A/(1+exp(4*mu/A*(l-x)+2)),initial=powrDpltInit,
parameters=c("A","mu","l","b"))

test1<-nlsList(SSpowrDplt,NewData)
coef(test1)

-----
In theory, practice and theory are the same. In practice, they are not - Albert Einstein
--
View this message in context: http://r.789695.n4.nabble.com/Lower-bounds-on-selfStart-function-not-working-tp3999231p3999231.html
Sent from the R help mailing list archive at Nabble.com.
#
I tested the "optim" function and that is returning non-negative parameter
values (meaning they are bound by the lower limits), but I think those are
the starting estimates for the nlsList model which is then finding negative
values for the solution.

-----
In theory, practice and theory are the same. In practice, they are not - Albert Einstein
--
View this message in context: http://r.789695.n4.nabble.com/Lower-bounds-on-selfStart-function-not-working-tp3999231p4001986.html
Sent from the R help mailing list archive at Nabble.com.
#
I ran the code again and got an error saying that the "x" was unknown. I
don't know why I hadn't seen that error before. Anyway, I made the edits to
"func1" so instead of "x", it is "xy$x."
#function to optimize
func1 <- function(value) {
   A.s  <- value[1]
   mu.s  <- value[2]
   l.s <- value[3]
   b.s    <- value[4]
  y1<-rep(0,length(xy$x)) # generate vector for predicted y (y1) to evaluate
against observed y
  for(cnt in 1:length(xy$x)){
  y1[cnt]<- b.s+A.s/(1+exp(4*mu.s/A.s*(l.s-xy$x[cnt])+2))}  #predicting y1
for values of y
  evl<-sum((xy$y-y1)^2) #sum of squares is function to minimize
  return(evl)}


There is another place where there is an "x" in the selfStart function:
SSpowrDplt<-selfStart(~b+A/(1+exp(4*mu/A*(l-x)+2)),initial=powrDpltInit,
parameters=c("A","mu","l","b"))

I don't know why that is working fine or how it knows that my "x" is that
specific one. It seems that I am not fully understanding how this is
working.

-----
In theory, practice and theory are the same. In practice, they are not - Albert Einstein
--
View this message in context: http://r.789695.n4.nabble.com/Lower-bounds-on-selfStart-function-not-working-tp3999231p4012805.html
Sent from the R help mailing list archive at Nabble.com.
7 days later
#
I was able to solve this problem by going back to nls and obtaining the
initial parameter estimates through optim. When I used nlsList with my
dataset, it took 2 minutes to solve and was not limited by the bounds. Now I
have the bounds working and it takes 45 seconds to solve. Here is the new
code:

A<-1.75 
mu<-.2 
l<-2 
b<-0 
x<-seq(0,18,.25) 
create.y<-function(x){ 
y<-b+A/(1+exp(4*mu/A*(l-x)+2)) 
return(y) 
} 
ys<-create.y(x) 
yvec<-(rep(ys,5))*(.9+runif(length(x)*5)/5) 
Trt<-factor(c(rep("A1",length(x)),rep("A2",length(x)),rep("A3",length(x)),rep("A4",length(x)),rep("A5",length(x)))) 
Data<-data.frame(Trt,rep(x,5),yvec) 
names(Data)<-c("Trt","x","y") 

NewData<-groupedData(y~x|Trt,data=Data) 

ids<-levels(factor(xy$Trt))
output<-matrix(0,length(ids),4)

modeltest<- function(A,mu,l,b,x){
out<-vector(length=length(x))
for (i in 1:length(x)) {
out[i]<-b+A/(1+exp(4*mu/A*(l-x[i])+2))
}
return(out)
}

lower.bound<-list(A=.01,mu=0,l=0,b=0)

for (i in 1:length(ids)){
xy<-subset(NewData,Trt==ids[i])
y<-xy$y
x<-xy$x


   A.s  <- max(xy$y)-min(xy$y) 
   mu.s  <- A.s/7.5 
   l.s <- 0 
   b.s    <- max(min(xy$y),0.00001) 
   value  <- c(A.s, l.s, mu.s, b.s) 

#function to optimize 
func1 <- function(value) { 
   A.s  <- value[1] 
   mu.s  <- value[2] 
   l.s <- value[3] 
   b.s    <- value[4] 
  y1<-rep(0,length(xy$x)) # generate vector for predicted y (y1) to evaluate
against observed y 
  for(cnt in 1:length(xy$x)){ 
  y1[cnt]<- b.s+A.s/(1+exp(4*mu.s/A.s*(l.s-x[cnt])+2))}  #predicting y1 for
values of y 
  evl<-sum((xy$y-y1)^2) #sum of squares is function to minimize 
  return(evl)} 

#optimizing 
oppar<-optim(c(A.s , mu.s , l.s , b.s),func1,method="L-BFGS-B", 
lower=c(0.0001,0.0,0.0,0.0), 
control=list(maxit=2000)) 

#saving optimized parameters 
value<-c(oppar$par[1L],oppar$par[2L],oppar$par[3L],oppar$par[4L]) 

   names(value) <- c("A","mu","l","b")

try(nmodel<-nls(y~modeltest(A,mu,l,b,x),data=xy,
start=value,
lower=lower.bound,
algorithm="port")
)

coefv<-coef(nmodel)
output[i,]<-coefv
   
} 



-----
In theory, practice and theory are the same. In practice, they are not - Albert Einstein
--
View this message in context: http://r.789695.n4.nabble.com/Lower-bounds-on-selfStart-function-not-working-tp3999231p4073639.html
Sent from the R help mailing list archive at Nabble.com.