[R] problem with coding for 'optim' in R
Michael Rennie
mrennie at utm.utoronto.ca
Tue Jul 15 01:04:36 CEST 2003
Hi,
My argument q is related to "Wtmod" and "Hgtmod" in function f (see below).
What I attempted to do was make function f the squared difference between my
modelled inputs and my observed endpoints of the 365-day outcomes of the
iteration process. The modelled outcomes 'Wtmod' and 'Hgtmod' are dependent
upon two particular variables, p and ACT. I've tried to combine the two of
those parameters into q by q<-(c(p, ACT)). Thus, the function f does rely on
the argument q, but only after performing the iterative process; I am trying to
minimize the difference between the outcome of the iterative process and the
observed endpoints I have specified in function f. Because there are 365
calculations that use p, ACT that need to happen before I get the iterative
endpoints, I need to get to the end of this iterative process before I can
compare my endpoints, in order to figure out how to adjust p and ACT. Upon
doing so, it needs to return to the beginning of the iterative process and
attempt different starting points for p, ACT (components of q).
It seems as though this should work, since when I specify the known solution
set for p, ACT, I appear to get a solution back for function f (It takes a hell
of a long time to figure out though, which is strange since I have literally
given it the answer- see OUTPUT #1). However, it does give me a wierd message
at the end, so something is still wrong.
But, when I try a different set of starting points (1,1), it doesn't work.
(OUTPUT # 2), and it gives me some strange error, stating that "L-BFGS-B needs
finite values of fn", which it didn't seem to need when I gave it the values
for p, ACT that I know yield the correct answer.
So, if the only time I can get 'optim' to work is when I give it the actual
solution, then it's pretty much useless, and not terribly sucessful at
optimizing.
Any ideas?
Mike
Quoting "Roger D. Peng" <rpeng at stat.ucla.edu>:
> Heuristically, 'optim' works by changing the imputs to your function 'f'
> and tries to find a minimum of the function. But your 'f' function does
> not actually reference the imputs that are passed in (i.e. your 'q'
> parameter). Therefore, changing the imput values will do nothing and I
> would be surprised if 'optim' worked at all. You appear to have another
> variable 'q' defined in the global workspace but this is not related to
> the 'q' which is passed as an argument to 'f'. 'optim' is modifying the
> value of 'q' that is passed to 'f', not the one stored in the global
> workspace.
>
> It seems your interpretation of how 'optim' works is correct, but you
> have to rewrite 'f' so that it actually uses the arguments passed to it.
>
> -roger
>
>
OUTPUT 1- using known solution set for p, ACT, the 'optim' function appears to
work (but not for any other starting point, see output 2 below), and I get a
message back, which I'm not sure what it means ("CONVERGENCE:
REL_REDUCTION_OF_F <= FACTR*EPSMCH").
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)
>
> #specify sttarting values
> #q0<-c(p = 1, ACT = 1)
>
> #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]<- q[1]*CA*(W[i]^CB)*((comp[i,4])^Xc)*(exp(Xc*(1-(comp[i,4]))))*Pc
+
+ ASMR[i]<- q[2]*RA*(W[i]^RB)*((comp[i,5])^Xa)*(exp(Xa*(1-(comp[i,5]))))
+
+ SMR[i]<- (ASMR[i]/q[2])
+
+ 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
+
+ q
+
+ f <- (((Wt-Wtmod)^2 + (Hgt-Hgtmod)^2)/2)^2 ; f
+
+ #warnings()
+
+ #write.table (bioday, file = "perch.csv", append = FALSE, sep=",", na = NA,
col.names = TRUE)
+
+
+
+ #nlm(f,c(1,1))
+ }
>
> optim(q, f, method = "L-BFGS-B",
+ lower = c(0, 0), upper=c(2, 10))
$par
[1] 0.5586205 1.6676453
$value
[1] 8.703706e-10
$counts
function gradient
14 14
$convergence
[1] 0
$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
>
> write.table (bioday, file = "perch.csv", append = FALSE, sep=",", na = NA,
col.names = TRUE)
Error in inherits(x, "data.frame") : Object "bioday" not found
Execution halted
OUTPUT 2- trying starting points of p, ACT = 1, the program no longer works.
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 <- 1 #0.558626306252032
> ACT <- 1 #1.66764519286918
>
> q<-c(p,ACT)
>
> #specify sttarting values
> #q0<-c(p = 1, ACT = 1)
>
> #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]<- q[1]*CA*(W[i]^CB)*((comp[i,4])^Xc)*(exp(Xc*(1-(comp[i,4]))))*Pc
+
+ ASMR[i]<- q[2]*RA*(W[i]^RB)*((comp[i,5])^Xa)*(exp(Xa*(1-(comp[i,5]))))
+
+ SMR[i]<- (ASMR[i]/q[2])
+
+ 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
+
+ q
+
+ f <- (((Wt-Wtmod)^2 + (Hgt-Hgtmod)^2)/2)^2 ; f
+
+ #warnings()
+
+ #write.table (bioday, file = "perch.csv", append = FALSE, sep=",", na = NA,
col.names = TRUE)
+
+
+
+ #nlm(f,c(1,1))
+ }
>
> optim(q, f, method = "L-BFGS-B",
+ lower = c(0, 0), upper=c(2, 10))
Error in optim(q, f, method = "L-BFGS-B", lower = c(0, 0), upper = c(2, :
L-BFGS-B needs finite values of fn
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
More information about the R-help
mailing list