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]]
______________________________________________ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
______________________________________________ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help