[R] problem with coding for 'optim' in R

[Roger D. Peng]

It's important to remember that in R functions return whatever happens 
to be the last element of the function block, unless there is an 
explicit 'return' statement.  Your function 'f' in the second example is 
written incorrectly and will not work in 'optim'.  The last element in 
the function block is:

write.table (bioday, file = "perch.csv", append = FALSE, sep=",", na = 
NA, col.names = TRUE)

which I assume is *not* the value you want the function return.  Your 
function 'f' is returning whatever 'write.table' returns, which is 
nothing useful.  My guess is that you want your function 'f' to return 
the value 'f' defined in the function as

f <- (((Wt-Wtmod)^2 + (Hgt-Hgtmod)^2)/2)^2

So this statement should be the last line of your function.  

Also, your function 'f' (still from the second output) doesn't use the value 'q' at all, so I can't see how the optimizer can optimize a function that ignores its parameters.

-roger

Michael Rennie wrote:

[snip]

>OUTPUT 2- program that doesn't work, and gets screwed up in the daily 
>iterations- cqan't recognize starting conditions for q, even though it 
>worked fine before I placed it within the 'optim' function;
>
>R : Copyright 2001, The R Development Core Team
>Version 1.4.0  (2001-12-19)
>
>R is free software and comes with ABSOLUTELY NO WARRANTY.
>You are welcome to redistribute it under certain conditions.
>Type `license()' or `licence()' for distribution details.
>
>R is a collaborative project with many contributors.
>Type `contributors()' for more information.
>
>Type `demo()' for some demos, `help()' for on-line help, or
>`help.start()' for a HTML browser interface to help.
>Type `q()' to quit R.
>
> > ####################################
> > #    perch.R                       #
> > # Hewett and Johnson bioenergetics #
> > # model combined with              #
> > # Trudel MMBM to estimate          #
> > # Consumption in perch in R code   #
> > # Execute with                     #
> > # R --vanilla < perch.R > perch.out#
> > ####################################
> >
> > #USER INPUT BELOW
> >
> > #Weight at time 0
> > Wo<- 9.2
> >
> > #Hg concentration at time 0 (ugHg/g wet weight)
> > Hgo<- 0.08
> >
> > #Weight at time t
> > Wt<- 32.2
> >
> > #Hg concentration at time t (ugHg/g wet weight)
> > Hgt<- 0.110
> >
> > #Prey methylmercury concentration (as constant)
> > Hgp<- 0.033
> >
> > #Prey caloric value (as constant)
> > Pc<- 800
> >
> > #Energy density of fish (as constant, calories)
> > Ef <- 1000
> >
> > #Maturity status, 0=immature, 1=mature
> > Mat<- 0
> >
> > #Sex, 1=male, 2=female
> > Sex<- 1
> >
> > #USER INPUT ABOVE
> >
> > #Bioenergetics parameters for perch
> > CA <- 0.25
> > CB <- 0.73  #same as 1+(-0.27)- convert g/g/d to g/d * Pc to get cal/d
> > CQ <- 2.3
> > CTO <- 23
> > CTM <- 28
> > Zc<- (log(CQ))*(CTM-CTO)
> > Yc<- (log(CQ))*(CTM-CTO+2)
> > Xc<- ((Zc^2)*(1+(1+40/Yc)^0.5)^2)/400
> >
> > RA <- 34.992  #0.0108*3240 cal/g 02, converting weight of 02 to cal
> > RB <- 0.8   #same as 1+(-0.2) see above...
> > RQ <- 2.1
> > RTO <- 28
> > RTM <- 33
> > Za <- (log(RQ))*(RTM-RTO)
> > Ya<- (log(RQ))*(RTM-RTO+2)
> > Xa<- ((Za^2)*(1+(1+40/Ya)^0.5)^2)/400
> >
> > S <- 0.172
> >
> > FA <- 0.158
> > FB <- -0.222
> > FG <- 0.631
> >
> > UA<- 0.0253
> > UB<- 0.58
> > UG<- -0.299
> >
> > #Mass balance model parameters
> > EA <- 0.002938
> > EB <- -0.2
> > EQ <- 0.066
> > a <- 0.8
> >
> > #Specifying sex-specific parameters
> >
> > GSI<- NULL
> >
> > if (Sex==1) GSI<-0.05 else
>+ if (Sex==2) GSI<-0.17
> >
> > # Define margin of error functions
> > #merror <- function(phat,M,alpha) # (1-alpha)*100% merror for a proportion
> > #    {
> > #    z <- qnorm(1-alpha/2)
> > #    merror <- z * sqrt(phat*(1-phat)/M)  # M is (Monte Carlo) sample size
> > #    merror
> > #    }
> >
> > #Bring in temp file
> >
> > temper <- scan("temp.dat", na.strings = ".", list(Day=0, jday=0, Temp=0))
>Read 366 records
> >
> > Day<-temper$Day ; jday<-temper$jday ; Temp<-temper$Temp ;
> >
> > temp<- cbind (Day, jday, Temp)
> > #Day = number of days modelled, jday=julian day, Temp = daily avg. temp.
> > #temp [,2]
> >
> > Vc<-(CTM-(temp[,3]))/(CTM-CTO)
> > Vr<-(RTM-(temp[,3]))/(RTM-RTO)
> >
> > comp<- cbind (Day, jday, Temp, Vc, Vr)
> >
> > #comp
> >
> > bio<-matrix(NA, ncol=13, nrow=length(Day))
> > W<-NULL
> > C<-NULL
> > ASMR<-NULL
> > SMR<-NULL
> > A<-NULL
> > F<-NULL
> > U<-NULL
> > SDA<-NULL
> > Gr<-NULL
> > Hg<-NULL
> > Ed<-NULL
> > GHg<-NULL
> > K<-NULL
> > Expegk<-NULL
> > EGK<-NULL
> > p<-NULL
> > ACT<-NULL
> >
> >
> > #p <- 0.558626306252032
> > #ACT <- 1.66764519286918
> >
> > q<-c(p,ACT)
> >
> > #introduce function to solve
> > f <- function (q)
>+ {
>+
>+ M<- length(Day) #number of days iterated
>+
>+ for (i in 1:M)
>+ {
>+
>+ #Bioenergetics model
>+
>+ if (Day[i]==1) W[i] <- Wo else
>+ if (jday[i]==121 && Mat==1) W[i] <- (W[i-1]-(W[i-1]*GSI*1.2)) else
>+ W[i] <- (W[i-1]+(Gr[i-1]/Ef))
>+ #W
>+
>+ #W<-Wo
>+
>+ C[i]<- p*CA*(W[i]^CB)*((comp[i,4])^Xc)*(exp(Xc*(1-(comp[i,4]))))*Pc
>+
>+ ASMR[i]<- ACT*RA*(W[i]^RB)*((comp[i,5])^Xa)*(exp(Xa*(1-(comp[i,5]))))
>+
>+ SMR[i]<- (ASMR[i]/ACT)
>+
>+ A[i]<- (ASMR[i]-SMR[i])
>+
>+ F[i]<- (FA*((comp[i,3])^FB)*(exp(FG*p))*C[i])
>+
>+ U[i]<- (UA*((comp[i,3])^UB)*(exp(UG*p))*(C[i]-F[i]))
>+
>+ SDA[i]<- (S*(C[i]-F[i]))
>+
>+ Gr[i]<- (C[i]-(ASMR[i]+F[i]+U[i]+SDA[i]))
>+
>+ #Trudel MMBM
>+
>+ if (Day[i]==1) Hg[i] <- Hgo else Hg[i] <- 
>a*Hgp*(C[i-1]/Pc/W[i-1])/EGK[i-1]*(1-Expegk[i-1])+(Hg[i-1]*Expegk[i-1])
>+
>+ Ed[i]<- EA*(W[i]^EB)*(exp(EQ*(comp[i,3])))
>+
>+ GHg[i] <- Gr[i]/Ef/W[i]
>+
>+ if (Sex==1) 
>K[i]<-(((0.1681*(10^(1.3324+(0.000453*Hg[i])))/1000)/Hg[i])*GSI)/M else
>+ if (Sex==2) 
>K[i]<-(((0.1500*(10^(0.8840+(0.000903*Hg[i])))/1000)/Hg[i])*GSI)/M
>+ # = dw/ww conversion * gonad ~ body conc'n function(ng/g) / convert to ug/g
>+ # then express as Q times GSI gives K / M gives daily K
>+
>+ EGK[i] <- (Ed[i] + GHg[i] + (K[i]*Mat))
>+
>+ Expegk[i] <- exp(-1*EGK[i])
>+
>+ bio<- cbind(W, C, ASMR, SMR, A, F, U, SDA, Gr, Ed, GHg, EGK, Hg)
>+
>+ }
>+
>+ #warnings()
>+
>+ dimnames (bio) <-list(NULL, c("W", "C", "ASMR", "SMR", "A", "F", "U", 
>"SDA", "Gr", "Ed", "GHg", "EGK", "Hg"))
>+
>+ bioday<-cbind(jday, W, C, ASMR, SMR, A, F, U, SDA, Gr, Ed, GHg, EGK, Hg)
>+
>+ dimnames (bioday) <-list(NULL, c("jday", "W", "C", "ASMR", "SMR", "A", 
>"F", "U", "SDA", "Gr", "Ed", "GHg", "EGK", "Hg"))
>+
>+ #bioday
>+
>+ Wtmod<- bioday [length(W),2]
>+ Wtmod
>+
>+ Hgtmod<- bioday [length(Hg),14]
>+ Hgtmod
>+
>+ f <- (((Wt-Wtmod)^2 + (Hgt-Hgtmod)^2)/2)^2 ; f
>+
>+ q
>+
>+ #warnings()
>+
>+ write.table (bioday, file = "perch.csv", append = FALSE, sep=",", na = 
>NA, col.names = TRUE)
>+
>+ }
> >
> > #nlm(f,c(1,1))
> >
> > optim(c(1,1), f, method = "L-BFGS-B",
>+       lower = c(0, 0), upper=c(2, 10))
>Error in "[<-"(*tmp*, i, value = (W[i - 1] + (Gr[i - 1]/Ef))) :
>         nothing to replace with
>Execution halted
>
>
>Michael Rennie
>M.Sc. Candidate
>University of Toronto at Mississauga
>3359 Mississauga Rd. N.
>Mississauga, ON  L5L 1C6
>Ph: 905-828-5452  Fax: 905-828-3792
>	[[alternative HTML version deleted]]
>
>______________________________________________
>R-help at stat.math.ethz.ch mailing list
>https://www.stat.math.ethz.ch/mailman/listinfo/r-help
>
>
>  
>