[R] Excel can do what R can't?????
Spencer Graves
spencer.graves at pdf.com
Tue Jul 15 20:25:55 CEST 2003
I've programmed many things like this in both Excel and R. When I
did not get the same answer from both, it was because I had an error in
one (or both). I do this routinely as part of debugging: I catch many
mistakes this way, and I often feel I can not trust my answers without
this level of checking.
I use Excel with Solver if I only need one solution or if I'm working
with someone who doesn't have R or S-Plus. Otherwise, I prefer S-Plus
of R.
First forget about "optim": Do you get the same numbers from your
function "f" and from Excel? Have you plotted the function to be sure
you even have local minima? Naively, I would expect Excel to be more
likely to get stuck in local minima than "optim".
I'm sorry you've had such a frustrating experience with R. The S
language is very powerful but does have a steep learning curve. I had
S-Plus for 3-5 years before I was finally able to do anything useful
with it. Now, it is an integral part of how I do much of what I do.
hope this helps.
spencer graves
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