On Mon, 10 Sep 2007, Yuchen Luo wrote:
> 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.
Why not? ui%*%theta>=ci is a six-dimensional constraint in your example,
so a different element of the constraint can be violated.
-thomas
> 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)
>>>
>>> ###########################################################
>>>
>>> ______________________________________________
>>> [EMAIL PROTECTED] 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]]
>
> ______________________________________________
> [EMAIL PROTECTED] 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.
>
Thomas Lumley Assoc. Professor, Biostatistics
[EMAIL PROTECTED] University of Washington, Seattle
______________________________________________
R-help@r-project.org 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.
