[R] Excel can do what R can't?????
Jerome Asselin
jerome at hivnet.ubc.ca
Tue Jul 15 20:47:47 CEST 2003
Mike,
The definition of your function f() seems quite inefficient. You could
vectorize it, which would shorten and speed up your code, especially if M
is large. See the R introduction file available online to learn how to do
it if you don't already know how. Also, you have to return only one
argument. Unless I'm wrong, your function wants to return Wtmod, Hgtmod, q
and f. I'm don't think this would change anything in this case, but you
should definitely clean this up!
Another advice... If you can simplify your example into a few lines of
"ready-to-execute" code with a toy dataset, then it's easy for everyone to
try it out and you can get more feedback. The code you've included is
quite large and cumbersome. For one thing, you could easily have removed
the lines of code that were commented out.
Meanwhile, I would suggest that you go back to the basics of R to clean up
your code.
Sorry I can't be more helpful.
Jerome
On July 15, 2003 10:46 am, Michael Rennie wrote:
> Hi there
>
> I thought this would be of particular interest to people using 'optim'
> functions and perhaps people involved with R development.
>
> I've been beaten down by R trying to get it to perform an optimization
> on a mass-balance model. I've written the same program in excel, and
> using the 'solver' function, it comes up with an answer for my variables
> (p, ACT, which I've assigned to q in R) that gives a solution to the
> function "f" in about 3 seconds, with a value of the function around
> 0.0004. R, on the other hand, appears to get stuck in local minima, and
> spits back an approximation that is close the the p, ACT values excel
> does, but not nearly precise enough for my needs, and not nearly as
> precise as excel, and it takes about 3 minutes. Also, the solution for
> the value it returns for the function is about 8 orders of magnitude
> greater than the excel version, so I can't really say the function is
> approximating zero. I was able to determine this using a "trace"
> command on function f, which is listed below.
>
> This is very likely due to the fact that I've made some coding error
> along the way, or have done something else wrong, but I don't know.
> Either way, I am shocked and surprised that a program like excel is
> outperforming R. I've attached my command file and the dataset
> "temp.dat" at the bottom of this e-mail for anyone who would like to
> fiddle around with it, and if you come up with something, PLEASE let me
> know- In the meantime, I've got to start fiddling with excel and
> figuring out how to automate the solver calculation.....
>
> Briefly, the point of the program is to approximate the model output
> from an iterative calculation, Wtmod and Hgtmod, to user-specified
> endpoints Wt and Hgt, by seeking the optimal values of p, ACT involved
> in the iterative process.
>
> Also, if your interested in recent correspondence that explains the
> point of the program a bit, and how the function ties in to the
> iterative process, search the R help forum for e-mails entitled "[R]
> problem with coding for 'optim' in R". Thanks also to Roger Peng and
> numerous others for helping me get this far.
>
> The whole point of me doing this in R was because it's supposed to be
> spectacularly fast at automating complex loops, but seems to be falling
> short for this application. Hopefully it's something wrong with my
> coding and not with R itself.
>
> Mike
>
> R COMMAND FILE:
>
> ####################################
> # 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))
>
> 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
>
> #starting values for p, ACT
> p <- 1 # 0.558626306252032 #solution set for p, ACT from excel 'solver'
> f'n ACT <- 2 # 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 <- 1000000000*(((((Wt-Wtmod)^2)/Wt) + (((Hgt-Hgtmod)^2)/Hgt))^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.2, 1), upper=c(2, 3),
> control = list(fnscale = 0.001))
>
>
> TRACE FUNCTION USED TO DETERMINE WHERE R IS GETTING STUCK;
>
> trace("f", quote(print(q)), at = 1, print = FALSE)
>
> o <- optim(c(1,2), f, method = "L-BFGS-B", lower = c(0.2,1), upper =
> c(2,3))
>
> DATA FOR TEMP.DAT:
>
> 1 153 9.4
> 2 154 9.6
> 3 155 9.8
> 4 156 10
> 5 157 10.2
> 6 158 10.4
> 7 159 10.6
> 8 160 10.8
> 9 161 11
> 10 162 11.2
> 11 163 11.4
> 12 164 11.6
> 13 165 11.8
> 14 166 12
> 15 167 12.3
> 16 168 12.5
> 17 169 12.7
> 18 170 12.9
> 19 171 13.1
> 20 172 13.4
> 21 173 13.6
> 22 174 13.8
> 23 175 14
> 24 176 14.2
> 25 177 14.5
> 26 178 14.7
> 27 179 14.9
> 28 180 15.1
> 29 181 15.4
> 30 182 15.6
> 31 183 15.8
> 32 184 16
> 33 185 16.2
> 34 186 16.5
> 35 187 16.7
> 36 188 16.9
> 37 189 17.1
> 38 190 17.3
> 39 191 17.5
> 40 192 17.7
> 41 193 17.9
> 42 194 18.1
> 43 195 18.3
> 44 196 18.5
> 45 197 18.7
> 46 198 18.9
> 47 199 19
> 48 200 19.2
> 49 201 19.4
> 50 202 19.5
> 51 203 19.7
> 52 204 19.9
> 53 205 20
> 54 206 20.2
> 55 207 20.3
> 56 208 20.4
> 57 209 20.5
> 58 210 20.7
> 59 211 20.8
> 60 212 20.9
> 61 213 21
> 62 214 21.1
> 63 215 21.2
> 64 216 21.3
> 65 217 21.3
> 66 218 21.4
> 67 219 21.5
> 68 220 21.5
> 69 221 21.6
> 70 222 21.6
> 71 223 21.6
> 72 224 21.7
> 73 225 21.7
> 74 226 21.7
> 75 227 21.7
> 76 228 21.7
> 77 229 21.7
> 78 230 21.7
> 79 231 21.6
> 80 232 21.6
> 81 233 21.6
> 82 234 21.5
> 83 235 21.5
> 84 236 21.4
> 85 237 21.3
> 86 238 21.3
> 87 239 21.2
> 88 240 21.1
> 89 241 21
> 90 242 20.9
> 91 243 20.8
> 92 244 20.7
> 93 245 20.5
> 94 246 20.4
> 95 247 20.3
> 96 248 20.2
> 97 249 20
> 98 250 19.9
> 99 251 19.7
> 100 252 19.5
> 101 253 19.4
> 102 254 19.2
> 103 255 19
> 104 256 18.9
> 105 257 18.7
> 106 258 18.5
> 107 259 18.3
> 108 260 18.1
> 109 261 17.9
> 110 262 17.7
> 111 263 17.5
> 112 264 17.3
> 113 265 17.1
> 114 266 16.9
> 115 267 16.7
> 116 268 16.5
> 117 269 16.2
> 118 270 16
> 119 271 15.8
> 120 272 15.6
> 121 273 15.4
> 122 274 15.1
> 123 275 14.9
> 124 276 14.7
> 125 277 14.5
> 126 278 14.2
> 127 279 14
> 128 280 13.8
> 129 281 13.6
> 130 282 13.4
> 131 283 13.1
> 132 284 12.9
> 133 285 12.7
> 134 286 12.5
> 135 287 12.3
> 136 288 12
> 137 289 11.8
> 138 290 11.6
> 139 291 11.4
> 140 292 11.2
> 141 293 11
> 142 294 10.8
> 143 295 10.6
> 144 296 10.4
> 145 297 10.2
> 146 298 10
> 147 299 9.8
> 148 300 9.6
> 149 301 9.4
> 150 302 9.3
> 151 303 9.1
> 152 304 8.9
> 153 305 8.7
> 154 306 8.6
> 155 307 8.4
> 156 308 8.2
> 157 309 8.1
> 158 310 7.9
> 159 311 7.8
> 160 312 7.6
> 161 313 7.5
> 162 314 7.3
> 163 315 7.2
> 164 316 7
> 165 317 6.9
> 166 318 6.8
> 167 319 6.7
> 168 320 6.5
> 169 321 6.4
> 170 322 6.3
> 171 323 6.2
> 172 324 6.1
> 173 325 6
> 174 326 5.8
> 175 327 5.7
> 176 328 5.6
> 177 329 5.5
> 178 330 5.5
> 179 331 5.4
> 180 332 5.3
> 181 333 5.2
> 182 334 5.1
> 183 335 5
> 184 336 5
> 185 337 4.9
> 186 338 4.8
> 187 339 4.7
> 188 340 4.7
> 189 341 4.6
> 190 342 4.5
> 191 343 4.5
> 192 344 4.4
> 193 345 4.4
> 194 346 4.3
> 195 347 4.3
> 196 348 4.2
> 197 349 4.2
> 198 350 4.1
> 199 351 4.1
> 200 352 4
> 201 353 4
> 202 354 4
> 203 355 3.9
> 204 356 3.9
> 205 357 3.8
> 206 358 3.8
> 207 359 3.8
> 208 360 3.8
> 209 361 3.7
> 210 362 3.7
> 211 363 3.7
> 212 364 3.6
> 213 365 3.6
> 214 366 3.6
> 215 1 3.2
> 216 2 3.2
> 217 3 3.2
> 218 4 3.2
> 219 5 3.2
> 220 6 3.2
> 221 7 3.2
> 222 8 3.2
> 223 9 3.2
> 224 10 3.2
> 225 11 3.2
> 226 12 3.2
> 227 13 3.2
> 228 14 3.2
> 229 15 3.2
> 230 16 3.2
> 231 17 3.2
> 232 18 3.2
> 233 19 3.2
> 234 20 3.2
> 235 21 3.2
> 236 22 3.2
> 237 23 3.2
> 238 24 3.2
> 239 25 3.2
> 240 26 3.2
> 241 27 3.2
> 242 28 3.2
> 243 29 3.2
> 244 30 3.2
> 245 31 3.2
> 246 32 3.2
> 247 33 3.2
> 248 34 3.2
> 249 35 3.2
> 250 36 3.2
> 251 37 3.2
> 252 38 3.2
> 253 39 3.2
> 254 40 3.2
> 255 41 3.2
> 256 42 3.2
> 257 43 3.2
> 258 44 3.2
> 259 45 3.2
> 260 46 3.2
> 261 47 3.2
> 262 48 3.2
> 263 49 3.2
> 264 50 3.2
> 265 51 3.2
> 266 52 3.2
> 267 53 3.2
> 268 54 3.3
> 269 55 3.3
> 270 56 3.3
> 271 57 3.3
> 272 58 3.3
> 273 59 3.3
> 274 60 3.3
> 275 61 3.3
> 276 62 3.3
> 277 63 3.3
> 278 64 3.3
> 279 65 3.3
> 280 66 3.3
> 281 67 3.3
> 282 68 3.3
> 283 69 3.3
> 284 70 3.3
> 285 71 3.4
> 286 72 3.4
> 287 73 3.4
> 288 74 3.4
> 289 75 3.4
> 290 76 3.4
> 291 77 3.4
> 292 78 3.4
> 293 79 3.5
> 294 80 3.5
> 295 81 3.5
> 296 82 3.5
> 297 83 3.5
> 298 84 3.5
> 299 85 3.6
> 300 86 3.6
> 301 87 3.6
> 302 88 3.6
> 303 89 3.6
> 304 90 3.7
> 305 91 3.7
> 306 92 3.7
> 307 93 3.8
> 308 94 3.8
> 309 95 3.8
> 310 96 3.8
> 311 97 3.9
> 312 98 3.9
> 313 99 4
> 314 100 4
> 315 101 4
> 316 102 4.1
> 317 103 4.1
> 318 104 4.2
> 319 105 4.2
> 320 106 4.3
> 321 107 4.3
> 322 108 4.4
> 323 109 4.4
> 324 110 4.5
> 325 111 4.5
> 326 112 4.6
> 327 113 4.7
> 328 114 4.7
> 329 115 4.8
> 330 116 4.9
> 331 117 5
> 332 118 5
> 333 119 5.1
> 334 120 5.2
> 335 121 5.3
> 336 122 5.4
> 337 123 5.5
> 338 124 5.5
> 339 125 5.6
> 340 126 5.7
> 341 127 5.8
> 342 128 6
> 343 129 6.1
> 344 130 6.2
> 345 131 6.3
> 346 132 6.4
> 347 133 6.5
> 348 134 6.7
> 349 135 6.8
> 350 136 6.9
> 351 137 7
> 352 138 7.2
> 353 139 7.3
> 354 140 7.5
> 355 141 7.6
> 356 142 7.8
> 357 143 7.9
> 358 144 8.1
> 359 145 8.2
> 360 146 8.4
> 361 147 8.6
> 362 148 8.7
> 363 149 8.9
> 364 150 9.1
> 365 151 9.3
> 366 152 9.3
>
>
>
> 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
More information about the R-help
mailing list