[R] Excel can do what R can't?????
Michael Rennie
mrennie at utm.utoronto.ca
Tue Jul 15 22:24:47 CEST 2003
At 11:47 AM 7/15/03 -0700, Jerome Asselin wrote:
>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.
Hi, Jerome
I don;t think I can vectorize it, since in the iteration loop, the value
for each [i] is dependent on the value of [i-1], so I require the loop to
go through each [i] before I can get my values for any particular vector
(variable). I actually had my program operating this way in the first
place, but I get all sorts of warnings and the 'optim' function especially
doesn't seem to appreciate it.
>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!
The calls to Wtmod, q, and Hgtmod are all just residual from the
development of the loop inside function f. I would like to get the last
line of 'bioday' reported from within the loop, had I figured out the
optimization, but that point is rather moot unless I can get the
optimization functioning.
>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.
Thanks for the advice- every bit helps if I eventually get this thing to
work.....
Mike
>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
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