[R] Are the error messages of ConstrOptim() consisten with each other?

2007-09-10 Thread Yuchen Luo 
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
 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

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!

Re: [R] Are the error messages of ConstrOptim() consisten with each other?

2007-09-10 Thread Duncan Murdoch 
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/s!
 h*!
  
 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*l!
 am!
  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*1!
 00!

Re: [R] Are the error messages of ConstrOptim() consisten with each other?

2007-09-10 Thread Thomas Lumley 

The error message about the feasible region comes from constrOptim(), 
before your function is called.  The error message about missing lambda1 
comes from calling your function.

-thomas

On Sun, 9 Sep 2007, 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
 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

 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/s!
 h*!
 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*la!
 m!
 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*1!
 00!

Re: [R] Are the error messages of ConstrOptim() consisten with each other?

2007-09-10 Thread Yuchen Luo 
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/