Dear Professor Murdoch. Thank you for your help! 1. I believe c(0.5,0.3,0.5) satisfies the constrain because I did the following experiment ui=-1*ui ci=-1*ci constrOptim(c(0.5,0.3,0.5), f=fit.error, gr=fit.error.grr, ui=ui,ci=ci)
The same error message pops up. Any theta ( in this case, c(0.5,0.3,0.5)) cannot violate both ui%*%theta>=ci and -ui%*%theta>=-ci. 2. There is lambda1 available. The 0.3 in c(0.5,0.3,0.5) is lambda1. If you plug c(0.5,0.3,0.5) into fit.error and fit.error.grr by fit.error(0.5,0.3,0.5) fit.error.grr(0.5,0.3,0.5) It works. Best Wishes Yuchen Luo On 9/10/07, Duncan Murdoch <[EMAIL PROTECTED]> wrote: > > Yuchen Luo wrote: > > Dear Friends. > > I found something very puzzling with constOptim(). When I change the > > parameters for ConstrOptim, the error messages do not seem to be > > consistent with each other: > > > > > >> constrOptim(c(0.5,0.3,0.5), f=fit.error, gr=fit.error.grr, ui=ui,ci=ci) > >> > > Error in constrOptim(c(0.5, 0.3, 0.5), f = fit.error, gr = fit.error.grr > , : > > initial value not feasible > > > "Not feasible" means it doesn't satisfy the constraints. > >> constrOptim(c(0.5,0.9,0.5), f=fit.error, gr=fit.error.grr, ui=ui,ci=ci) > >> > > Error in constrOptim(c(0.5, 0.9, 0.5), f = fit.error, gr = fit.error.grr > , : > > initial value not feasible > > > >> constrOptim(c(0.3,0.5,0.5), f=fit.error, gr=fit.error.grr, ui=ui,ci=ci) > >> > > Error in f(theta, ...) : argument "lambda1" is missing, with no default > > > > This time your starting values satisfied the constraints, so your > objective function was called, but you didn't pass it a value for lambda1. > > I only changed the parameters, how come the lambda1 that is not > > missing in the first 2 cases suddently become missing? > > > > For your convenience, I put the complete code below: > > > > Best Wishes > > Yuchen Luo > > > > ######################################## > > rm(list = ls()) > > > > mat=5 > > > > rint=c(4.33,4.22,4.27,4.43,4.43,4.44,4.45,4.65,4.77,4.77) > > tot=rep(13319.17,10) > > sh=rep(1553656,10) > > sigmae=c(0.172239074,0.188209271,0.193703774,0.172659891,0.164427247, > 0.24602361,0.173555309,0.186701165,0.193150456, > > 0.1857315601) > > ss=c(56.49,56.39,56.55,57.49,57.37,55.02,56.02,54.35,54.09, 54.67) > > orange=rep(21.25,10) > > > > apple2=expression(rint*(1.0-rec)*(1.0- > (pnorm(-lambda/2.0+log(((ss+(tot/sh*1000.0)*lbar)/(tot/sh*1000.0 > )/lbar*exp(lambda*lambda)))/lambda)-((ss+(tot/sh*1000.0)*lbar)/(tot/sh* > 1000.0)/lbar*exp(lambda*lambda))*pnorm(-lambda/2.0-log(((ss+(tot/sh*1000.0 > )*lbar)/(tot/sh*1000.0 > )/lbar*exp(lambda*lambda)))/lambda))+(exp(rint*(lambda*lambda/(sigmae*ss/(ss+lbar*(tot/sh* > 1000.0)))/(sigmae*ss/(ss+lbar*(tot/sh*1000.0)))))*((((ss+(tot/sh*1000.0 > )*lbar)/(tot/sh*1000.0)/lbar*exp(lambda*lambda))^(sqrt(0.25+2.0*rint/ > (sigmae*ss/(ss+lbar*(tot/sh*1000.0)))/(sigmae*ss/(ss+lbar*(tot/sh*1000.0 > ))))+0.5)*pnorm(-log(((ss+(tot/sh*1000.0)*lbar)/(tot/sh*1000.0 > )/lbar*exp(lambda*lambda)))/((sigmae*ss/(ss+lbar*(tot/sh*1000.0 > )))*sqrt(mat+(lambda*lambda/(sigmae*ss/(ss+lbar*(tot/sh*1000.0 > )))/(sigmae*ss/(ss+lbar*(tot/sh*1000.0))))))-sqrt(0.25+2.0*rint/ > (sigmae*ss/(ss+lbar*(tot/sh*1000.0)))/(sigmae*ss/(ss+lbar*(tot/sh*1000.0 > ))))*(sigmae*ss/(ss+lbar*(tot/sh*1000.0 > )))*sqrt(mat+(lambda*lambda/(sigmae*ss/(ss+lbar*(tot/sh*! > > 1000.0)))/(sigmae*ss/(ss+lbar*(tot/sh*1000.0))))))+((ss+(tot/sh*1000.0 > )*lbar)/(tot/sh*1000.0)/lbar*exp(lambda*lambda))^(-sqrt(0.25+2.0*rint/ > (sigmae*ss/(ss+lbar*(tot/sh*1000.0)))/(sigmae*ss/(ss+lbar*(tot/sh*1000.0 > ))))+0.5)*pnorm(-log(((ss+(tot/sh*1000.0)*lbar)/(tot/sh*1000.0 > )/lbar*exp(lambda*lambda)))/((sigmae*ss/(ss+lbar*(tot/sh*1000.0 > )))*sqrt(mat+(lambda*lambda/(sigmae*ss/(ss+lbar*(tot/sh*1000.0 > )))/(sigmae*ss/(ss+lbar*(tot/sh*1000.0))))))+sqrt(0.25+2.0*rint/ > (sigmae*ss/(ss+lbar*(tot/sh*1000.0)))/(sigmae*ss/(ss+lbar*(tot/sh*1000.0 > ))))*(sigmae*ss/(ss+lbar*(tot/sh*1000.0 > )))*sqrt(mat+(lambda*lambda/(sigmae*ss/(ss+lbar*(tot/sh*1000.0 > )))/(sigmae*ss/(ss+lbar*(tot/sh*1000.0)))))))-(((ss+(tot/sh*1000.0 > )*lbar)/(tot/sh*1000.0)/lbar*exp(lambda*lambda))^(sqrt(0.25+2.0*rint/ > (sigmae*ss/(ss+lbar*(tot/sh*1000.0)))/(sigmae*ss/(ss+lbar*(tot/sh*1000.0 > ))))+0.5)*pnorm(-log(((ss+(tot/sh*1000.0)*lbar)/(tot/sh*1000.0 > )/lbar*exp(lambda*lambda)))/((sigmae*ss/(ss+lbar*(tot/sh*1000.0 > )))*sqrt((lambda*lam! > > bda/(sigmae*ss/(ss+lbar*(tot/sh*1000.0)))/(sigmae*ss/(ss+lbar*(tot/sh* > > 1000.0))))))-sqrt(0.25+2.0*rint/(sigmae*ss/(ss+lbar*(tot/sh*1000.0 > )))/(sigmae*ss/(ss+lbar*(tot/sh*1000.0))))*(sigmae*ss/(ss+lbar*(tot/sh* > 1000.0)))*sqrt((lambda*lambda/(sigmae*ss/(ss+lbar*(tot/sh*1000.0 > )))/(sigmae*ss/(ss+lbar*(tot/sh*1000.0))))))+((ss+(tot/sh*1000.0 > )*lbar)/(tot/sh*1000.0)/lbar*exp(lambda*lambda))^(-sqrt(0.25+2.0*rint/ > (sigmae*ss/(ss+lbar*(tot/sh*1000.0)))/(sigmae*ss/(ss+lbar*(tot/sh*1000.0 > ))))+0.5)*pnorm(-log(((ss+(tot/sh*1000.0)*lbar)/(tot/sh*1000.0 > )/lbar*exp(lambda*lambda)))/((sigmae*ss/(ss+lbar*(tot/sh*1000.0 > )))*sqrt((lambda*lambda/(sigmae*ss/(ss+lbar*(tot/sh*1000.0 > )))/(sigmae*ss/(ss+lbar*(tot/sh*1000.0))))))+sqrt(0.25+2.0*rint/ > (sigmae*ss/(ss+lbar*(tot/sh*1000.0)))/(sigmae*ss/(ss+lbar*(tot/sh*1000.0 > ))))*(sigmae*ss/(ss+lbar*(tot/sh*1000.0 > )))*sqrt((lambda*lambda/(sigmae*ss/(ss+lbar*(tot/sh*1000.0 > )))/(sigmae*ss/(ss+lbar*(tot/sh*1000.0 > ))))))))))/((pnorm(-lambda/2.0+log(((ss+(tot/sh*1000.0)*lbar)/(tot/sh* > 1000.0)/lbar*exp(lambda*lambda)))/lambda)-((ss+(tot/sh*100! > > 0.0)*lbar)/(tot/sh*1000.0 > )/lbar*exp(lambda*lambda))*pnorm(-lambda/2.0-log(((ss+(tot/sh*1000.0 > )*lbar)/(tot/sh*1000.0 > )/lbar*exp(lambda*lambda)))/lambda))-(pnorm(-sqrt((sigmae*ss/(ss+lbar*(tot/sh* > 1000.0)))*(sigmae*ss/(ss+lbar*(tot/sh*1000.0 > )))*mat+lambda*lambda)/2.0+log(((ss+(tot/sh*1000.0)*lbar)/(tot/sh*1000.0 > )/lbar*exp(lambda*lambda)))/sqrt((sigmae*ss/(ss+lbar*(tot/sh*1000.0 > )))*(sigmae*ss/(ss+lbar*(tot/sh*1000.0 > )))*mat+lambda*lambda))-((ss+(tot/sh*1000.0)*lbar)/(tot/sh*1000.0 > )/lbar*exp(lambda*lambda))*pnorm(-sqrt((sigmae*ss/(ss+lbar*(tot/sh*1000.0 > )))*(sigmae*ss/(ss+lbar*(tot/sh*1000.0 > )))*mat+lambda*lambda)/2.0-log(((ss+(tot/sh*1000.0)*lbar)/(tot/sh*1000.0 > )/lbar*exp(lambda*lambda)))/sqrt((sigmae*ss/(ss+lbar*(tot/sh*1000.0 > )))*(sigmae*ss/(ss+lbar*(tot/sh*1000.0 > )))*mat+lambda*lambda)))*exp(-rint*mat)-(exp(rint*(lambda*lambda/(sigmae*ss/(ss+lbar*(tot/sh* > 1000.0)))/(sigmae*ss/(ss+lbar*(tot/sh*1000.0)))))*((((ss+(tot/sh*1000.0 > )*lbar)/(tot/sh*1000.0)/lbar*exp(lambda*lambda))^(sqrt(0.! > > 25+2.0*rint/(sigmae*ss/(ss+lbar*(tot/sh*1000.0)))/(sigmae*ss/(ss+lbar* > > (tot/sh*1000.0))))+0.5)*pnorm(-log(((ss+(tot/sh*1000.0)*lbar)/(tot/sh* > 1000.0)/lbar*exp(lambda*lambda)))/((sigmae*ss/(ss+lbar*(tot/sh*1000.0 > )))*sqrt(mat+(lambda*lambda/(sigmae*ss/(ss+lbar*(tot/sh*1000.0 > )))/(sigmae*ss/(ss+lbar*(tot/sh*1000.0))))))-sqrt(0.25+2.0*rint/ > (sigmae*ss/(ss+lbar*(tot/sh*1000.0)))/(sigmae*ss/(ss+lbar*(tot/sh*1000.0 > ))))*(sigmae*ss/(ss+lbar*(tot/sh*1000.0 > )))*sqrt(mat+(lambda*lambda/(sigmae*ss/(ss+lbar*(tot/sh*1000.0 > )))/(sigmae*ss/(ss+lbar*(tot/sh*1000.0))))))+((ss+(tot/sh*1000.0 > )*lbar)/(tot/sh*1000.0)/lbar*exp(lambda*lambda))^(-sqrt(0.25+2.0*rint/ > (sigmae*ss/(ss+lbar*(tot/sh*1000.0)))/(sigmae*ss/(ss+lbar*(tot/sh*1000.0 > ))))+0.5)*pnorm(-log(((ss+(tot/sh*1000.0)*lbar)/(tot/sh*1000.0 > )/lbar*exp(lambda*lambda)))/((sigmae*ss/(ss+lbar*(tot/sh*1000.0 > )))*sqrt(mat+(lambda*lambda/(sigmae*ss/(ss+lbar*(tot/sh*1000.0 > )))/(sigmae*ss/(ss+lbar*(tot/sh*1000.0))))))+sqrt(0.25+2.0*rint/ > (sigmae*ss/(ss+lbar*(tot/sh*1000.0)))/(sigmae*ss/(ss+lbar*(tot/sh*1000.0 > ))))*(sigmae*ss/(ss+lb! > > ar*(tot/sh*1000.0 > )))*sqrt(mat+(lambda*lambda/(sigmae*ss/(ss+lbar*(tot/sh*1000.0 > )))/(sigmae*ss/(ss+lbar*(tot/sh*1000.0)))))))-(((ss+(tot/sh*1000.0 > )*lbar)/(tot/sh*1000.0)/lbar*exp(lambda*lambda))^(sqrt(0.25+2.0*rint/ > (sigmae*ss/(ss+lbar*(tot/sh*1000.0)))/(sigmae*ss/(ss+lbar*(tot/sh*1000.0 > ))))+0.5)*pnorm(-log(((ss+(tot/sh*1000.0)*lbar)/(tot/sh*1000.0 > )/lbar*exp(lambda*lambda)))/((sigmae*ss/(ss+lbar*(tot/sh*1000.0 > )))*sqrt((lambda*lambda/(sigmae*ss/(ss+lbar*(tot/sh*1000.0 > )))/(sigmae*ss/(ss+lbar*(tot/sh*1000.0))))))-sqrt(0.25+2.0*rint/ > (sigmae*ss/(ss+lbar*(tot/sh*1000.0)))/(sigmae*ss/(ss+lbar*(tot/sh*1000.0 > ))))*(sigmae*ss/(ss+lbar*(tot/sh*1000.0 > )))*sqrt((lambda*lambda/(sigmae*ss/(ss+lbar*(tot/sh*1000.0 > )))/(sigmae*ss/(ss+lbar*(tot/sh*1000.0))))))+((ss+(tot/sh*1000.0 > )*lbar)/(tot/sh*1000.0)/lbar*exp(lambda*lambda))^(-sqrt(0.25+2.0*rint/ > (sigmae*ss/(ss+lbar*(tot/sh*1000.0)))/(sigmae*ss/(ss+lbar*(tot/sh*1000.0 > ))))+0.5)*pnorm(-log(((ss+(tot/sh*1000.0)*lbar)/(tot/sh*1000.0 > )/lbar*exp(lambda*! > > lambda)))/((sigmae*ss/(ss+lbar*(tot/sh*1000.0)))*sqrt((lambda*lambda/( > > sigmae*ss/(ss+lbar*(tot/sh*1000.0)))/(sigmae*ss/(ss+lbar*(tot/sh*1000.0 > ))))))+sqrt(0.25+2.0*rint/(sigmae*ss/(ss+lbar*(tot/sh*1000.0 > )))/(sigmae*ss/(ss+lbar*(tot/sh*1000.0))))*(sigmae*ss/(ss+lbar*(tot/sh* > 1000.0)))*sqrt((lambda*lambda/(sigmae*ss/(ss+lbar*(tot/sh*1000.0 > )))/(sigmae*ss/(ss+lbar*(tot/sh*1000.0))))))))))) > > > > apple.ana= function(rec1,lambda1,lbar1) > > {rec=rec1 > > lambda=lambda1 > > lbar=lbar1 > > apple=eval(apple2) > > gradient=cbind(eval(D(apple2,'rec')),eval(D(apple2,'lambda')), > > eval(D(apple2,'lbar'))) > > attr(apple.ana,'gradient')=gradient > > apple > > } > > > > fit.error=function(rec1,lambda1,lbar1) > > {rec=rec1 > > lambda=lambda1 > > lbar=lbar1 > > sum((eval(apple2)*1000-orange)^2/(orange^2)) > > } > > > > This is still coded incorrectly. Objective functions optimize over the > first parameter only. See ?optim for the details. constrOptim is just > a wrapper for optim. > > Duncan Murdoch > > fit.error.grr=function(rec1,lambda1,lbar1) > > {rec=rec1 > > lambda=lambda1 > > lbar=lbar1 > > > > > drec=sum(20000*eval(D(apple2,'rec'))*(eval(apple2)*10000-orange)/(orange^2)) > > > dlambda=sum(20000*eval(D(apple2,'lambda'))*(eval(apple2)*10000-orange)/(orange^2)) > > > dlbar=sum(20000*eval(D(apple2,'lbar'))*(eval(apple2)*10000-orange)/(orange^2)) > > > > c(drec,dlambda,dlbar) > > } > > > > ui=matrix(c(1,-1,0,0,0,0,0,0,1,-1,0,0,0,0,0,0,1,-1),6,3) > > ci=c(0,-0.5,0,-2,0,-0.6) > > > > constrOptim(c(0.5,0.3,0.5), f=fit.error, gr=fit.error.grr, ui=ui,ci=ci) > > constrOptim(c(0.5,0.9,0.5), f=fit.error, gr=fit.error.grr, ui=ui,ci=ci) > > constrOptim(c(0.3,0.5,0.5), f=fit.error, gr=fit.error.grr, ui=ui,ci=ci) > > > > ########################################################### > > > > ______________________________________________ > > R-help@stat.math.ethz.ch mailing list > > https://stat.ethz.ch/mailman/listinfo/r-help > > PLEASE do read the posting guide > http://www.R-project.org/posting-guide.html > > and provide commented, minimal, self-contained, reproducible code. > > > > > [[alternative HTML version deleted]] ______________________________________________ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.