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]]
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
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
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-helpMichael 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
R is good at automating specific kinds of complex loops, namely those that can be vectorized, or that can be written to draw on otherwise built-in facilities. It's usually reasonable for other kinds of loops but not spectacularly fast. You can write this part in C, though, quite easily, and R provides very convenient utilities for this. As for your code: You seem to have a system of equations that relates W and Hg to their one-period-ago values. It might clarify things if you coded this as a function: input time t values and q, output time t + 1 values. (You wouldn't need any arrays.) Then f would just iterate this function and calculate the criterion. Does the trajectory of (W, Hg) for given q in R seem correct? Does it agree with Excel? What does the criterion function look like? You could plot it in R and perhaps see if the surface is complicated, in which case a simple grid search might work for you. Reid Huntsinger -----Original Message----- From: Michael Rennie [mailto:mrennie at utm.utoronto.ca] Sent: Wednesday, July 16, 2003 11:18 AM To: Spencer Graves Cc: R-Help; M.Kondrin Subject: Re: [R] Excel can do what R can't????? Hi, Spencer I know I submitted a beastly ammount of code, but I'm not sure how to simplify it much further, and still sucessfully address the problem that i am having. The reason being is that the funciton begins f<- function (q) At the top of the iterative loop. This is what takes q and generates Wtmod, Hgtmod at the end of the iterative loop. the assignment to f occurs at the bottom of the iterative loop. So, yes, the call to f is performing an immediate computation, but based on arguments that are coming out of the iterative loop above it, arguments which depend on q<-(p, ACT). Maybe this is the problem; I've got too much going on between my function defenition and it's assignment, but I don't know how to get around it. So, I'm not sure if your example will work- the output from the iterative process is Wtmod, Hgtmod, and I want to minimize the difference between them and my observed endpoints (Wt, Hgt). The numbers I am varying to reach this optimization are in the iterative loop (p, ACT), so re-defining these outputs as x's and getting it to vary these doesn't do me much good unless they are directly linked to the output of the iterative loop above it. Last, it's not even that I'm getting error messages anymore- I just can't get the solution that I get from Excel. If I try to let R find the solution, and give it starting values of c(1,2), it gives me an optimization sulution, but an extremely poor one. However, if I give it the answer I got from excel, it comes right back with the same answer and solutions I get from excel. Using the 'trace' function, I can see that R gets stuck in a specific region of parameter space in looking for the optimization and just appears to give up. Even when it re-set itself, it keeps going back to this region, and thus doesn't even try a full range of the parameter space I've defined before it stops and gives me the wrong answer. I can try cleaning up the code and see if I can re-submit it, but what I am trying to program is so parameter heavy that 90% of it is just defining these at the top of the file. Thank you for the suggestions, Mike Quoting Spencer Graves <spencer.graves at PDF.COM>:> The phrase: > > f <- 1000000000*(((((Wt-Wtmod)^2)/Wt) + (((Hgt-Hgtmod)^2)/Hgt))2) ; f > > is an immediate computation, not a function. If you want a function, > try something like the following: > > f <- function(x){ > Wt <- x[1] > Wtmod <- x[2] > Hgt <- x[3] > Hgtmod <- x[4] > 1000000000*(((((Wt-Wtmod)^2)/Wt) + (((Hgt-Hgtmod)^2)/Hgt))2) > } > > OR > > f <- function(x, X){ > Wt <- X[,1] > Hgt <- X[,2] > Wtmod <- x[1] > Hgtmod <- x[2] > 1000000000*(((((Wt-Wtmod)^2)/Wt) + (((Hgt-Hgtmod)^2)/Hgt))2) > } > > "par" in "optim" is the starting values for "x". Pass "X" to "f" via > "..." in the call to "optim". > > If you can't make this work, please submit a toy example with the > code and error messages. Please limit your example to 3 observations, > preferably whole numbers so someone else can read your question in > seconds. If it is any longer than that, it should be ignored. > > hope this helps. > Spencer Graves > > M.Kondrin wrote: > > >?optim > > > > optim(par, fn, gr = NULL, > > method = c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN"), > > lower = -Inf, upper = Inf, > > control = list(), hessian = FALSE, ...) > > > > ..... > > fn: A function to be minimized (or maximized), with first > > argument the vector of parameters over which minimization is > > to take place. It should return a scalar result. > > > > Your fn defined as: > > f <- 1000000000*(((((Wt-Wtmod)^2)/Wt) + (((Hgt-Hgtmod)^2)/Hgt))2) ; f > > What is its first argument I wonder? > > I think you have just an ill-defined R function (although for Excel it > > may be OK - do not know) and optim just chokes on it. > > > > ______________________________________________ > > 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 ______________________________________________ R-help at stat.math.ethz.ch mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help ------------------------------------------------------------------------------ Notice: This e-mail message, together with any attachments, ...{{dropped}}
Hi, Reid>Do the values of W and Hg over time for a given q agree between R and Excel? >Not the optimal value of q, just the trajectories for fixed q (trying >several values for q).If I take the iterative loop out of the function, and ask for values of Hgmod, Hgtmod, and f, then I get EXACTLY what I get out of Excel. It's the optimization that seems to be the problem. If I trace the solutions, R isn't even exploring the full parameter space I tell it to look in. SO, the iterative loop is correct, and doing what it's supposed to, since values of p, ACT match exactly what they do in excel- it's something about how R is examining the possibilities in the optimization process that is giving me different answers between the two. I dunno- I'm going to tinker with it some more tonight. Mike>Reid > >-----Original Message----- >From: Michael Rennie [mailto:mrennie at utm.utoronto.ca] >Sent: Wednesday, July 16, 2003 2:47 PM >To: Huntsinger, Reid >Subject: RE: [R] Excel can do what R can't????? > > > >Hi, Reid > >At 02:09 PM 7/16/03 -0400, you wrote: > >R is good at automating specific kinds of complex loops, namely those that > >can be vectorized, or that can be written to draw on otherwise built-in > >facilities. It's usually reasonable for other kinds of loops but not > >spectacularly fast. You can write this part in C, though, quite easily, and > >R provides very convenient utilities for this. > > > >As for your code: You seem to have a system of equations that relates W and > >Hg to their one-period-ago values. It might clarify things if you coded >this > >as a function: input time t values and q, output time t + 1 values. (You > >wouldn't need any arrays.) Then f would just iterate this function and > >calculate the criterion. > >Wouldn't I still need to loop this function to get it through 365 days? Is >there a big difference, then, between this and what I've got? > > >Does the trajectory of (W, Hg) for given q in R seem correct? Does it agree > >with Excel? What does the criterion function look like? You could plot it >in > >R and perhaps see if the surface is complicated, in which case a simple >grid > >search might work for you. > >When I give R the values that I get in excel for p, ACT, the function >solution is actually more precise than what I get in Excel; the values for >p, ACT come back identical (then again, they are exactly what I >assigned..) But, if I leave R on it's own to find the solution, it keeps >getting jammed in a particular region. I've never done any function >plotting in R, but it would help if I could see what kind of surface I get >for f as a function of p, ACT- this would at least force R to examine the >full range of values specified by the upper and lower limits I've set >(which it isn't doing under the 'optim' command). > >Mike > > > >Reid Huntsinger > > > > > > > > > > > >-----Original Message----- > >From: Michael Rennie [mailto:mrennie at utm.utoronto.ca] > >Sent: Wednesday, July 16, 2003 11:18 AM > >To: Spencer Graves > >Cc: R-Help; M.Kondrin > >Subject: Re: [R] Excel can do what R can't????? > > > > > > > >Hi, Spencer > > > >I know I submitted a beastly ammount of code, but I'm not sure how to > >simplify > >it much further, and still sucessfully address the problem that i am >having. > > > >The reason being is that the funciton begins > > > >f<- function (q) > > > >At the top of the iterative loop. This is what takes q and generates >Wtmod, > > > >Hgtmod at the end of the iterative loop. the assignment to f occurs at the > >bottom of the iterative loop. So, yes, the call to f is performing an > >immediate > >computation, but based on arguments that are coming out of the iterative > >loop > >above it, arguments which depend on q<-(p, ACT). Maybe this is the >problem; > > > >I've got too much going on between my function defenition and it's > >assignment, > >but I don't know how to get around it. > > > >So, I'm not sure if your example will work- the output from the iterative > >process is Wtmod, Hgtmod, and I want to minimize the difference between >them > > > >and my observed endpoints (Wt, Hgt). The numbers I am varying to reach >this > > > >optimization are in the iterative loop (p, ACT), so re-defining these > >outputs > >as x's and getting it to vary these doesn't do me much good unless they are > >directly linked to the output of the iterative loop above it. > > > >Last, it's not even that I'm getting error messages anymore- I just can't > >get > >the solution that I get from Excel. If I try to let R find the solution, > >and > >give it starting values of c(1,2), it gives me an optimization sulution, >but > >an > >extremely poor one. However, if I give it the answer I got from excel, it > >comes right back with the same answer and solutions I get from excel. > > > >Using the 'trace' function, I can see that R gets stuck in a specific >region > >of > >parameter space in looking for the optimization and just appears to give >up. > > > >Even when it re-set itself, it keeps going back to this region, and thus > >doesn't even try a full range of the parameter space I've defined before it > >stops and gives me the wrong answer. > > > >I can try cleaning up the code and see if I can re-submit it, but what I am > >trying to program is so parameter heavy that 90% of it is just defining > >these > >at the top of the file. > > > >Thank you for the suggestions, > > > >Mike > > > > > >Quoting Spencer Graves <spencer.graves at PDF.COM>: > > > > > The phrase: > > > > > > f <- 1000000000*(((((Wt-Wtmod)^2)/Wt) + (((Hgt-Hgtmod)^2)/Hgt))2) ; >f > > > > > > is an immediate computation, not a function. If you want a function, > > > try something like the following: > > > > > > f <- function(x){ > > > Wt <- x[1] > > > Wtmod <- x[2] > > > Hgt <- x[3] > > > Hgtmod <- x[4] > > > 1000000000*(((((Wt-Wtmod)^2)/Wt) + (((Hgt-Hgtmod)^2)/Hgt))2) > > > } > > > > > > OR > > > > > > f <- function(x, X){ > > > Wt <- X[,1] > > > Hgt <- X[,2] > > > Wtmod <- x[1] > > > Hgtmod <- x[2] > > > 1000000000*(((((Wt-Wtmod)^2)/Wt) + (((Hgt-Hgtmod)^2)/Hgt))2) > > > } > > > > > > "par" in "optim" is the starting values for "x". Pass "X" to "f" via > > > "..." in the call to "optim". > > > > > > If you can't make this work, please submit a toy example with >the > > > code and error messages. Please limit your example to 3 observations, > > > preferably whole numbers so someone else can read your question in > > > seconds. If it is any longer than that, it should be ignored. > > > > > > hope this helps. > > > Spencer Graves > > > > > > M.Kondrin wrote: > > > > >?optim > > > > > > > > optim(par, fn, gr = NULL, > > > > method = c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", >"SANN"), > > > > lower = -Inf, upper = Inf, > > > > control = list(), hessian = FALSE, ...) > > > > > > > > ..... > > > > fn: A function to be minimized (or maximized), with first > > > > argument the vector of parameters over which minimization is > > > > to take place. It should return a scalar result. > > > > > > > > Your fn defined as: > > > > f <- 1000000000*(((((Wt-Wtmod)^2)/Wt) + (((Hgt-Hgtmod)^2)/Hgt))2) ; f > > > > What is its first argument I wonder? > > > > I think you have just an ill-defined R function (although for Excel it > > > > may be OK - do not know) and optim just chokes on it. > > > > > > > > ______________________________________________ > > > > 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 > > > >______________________________________________ > >R-help at stat.math.ethz.ch mailing list > >https://www.stat.math.ethz.ch/mailman/listinfo/r-help > > > >--------------------------------------------------------------------------- >--- > >Notice: This e-mail message, together with any attachments, contains > >information of Merck & Co., Inc. (Whitehouse Station, New Jersey, > >USA) that may be confidential, proprietary copyrighted and/or legally > >privileged, and is intended solely for the use of the individual or entity > >named on this message. If you are not the intended recipient, and > >have received this message in error, please immediately return this by > >e-mail and then delete it. > >--------------------------------------------------------------------------- >--- > >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 > > >------------------------------------------------------------------------------ >Notice: This e-mail message, together with any attachments, contains >information of Merck & Co., Inc. (Whitehouse Station, New Jersey, >USA) that may be confidential, proprietary copyrighted and/or legally >privileged, and is intended solely for the use of the individual or entity >named on this message. If you are not the intended recipient, and >have received this message in error, please immediately return this by >e-mail and then delete it. >------------------------------------------------------------------------------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
I don't understand how you can take the loop out of the function and still get values for the final timepoint. And whether the optimal parameter values agree wasn't my question. First I'd like to determine whether Hg and W (in your code) have the same value in R as they do in Excel, for a few possible values of the parameter q. Then, if so, whether the surface over which you're minimizing is complicated or not. As I said, if it is you can't expect local optimizers to work very well without good starting values, and they really don't explore the whole parameter space--that would be too slow for their usual applications. Optimization approaches like grid search or simulated annealing do try to cover the parameter space and may be better suited to your use. I would certainly try plotting and grid search just to see what's happening since that's clearly possible with two parameters. Reid Huntsinger -----Original Message----- From: Spencer Graves [mailto:spencer.graves at PDF.COM] Sent: Wednesday, July 16, 2003 5:20 PM To: Michael Rennie Cc: Huntsinger, Reid; R-Help Subject: Re: [R] Excel can do what R can't????? I'm confused: I've done this type of thing by programming the same objective function in R (or S-Plus) and Excel. After the answers from my objective function in R match the answers in Excel, then I pass that objective function to something like "optim", which then finds the same answers as "solver" in Excel. Your latest description makes me wonder if the function you pass to "optim" tries to do some of the optimization that "optim" is supposed to do. ??? hope this helps. spencer graves Michael Rennie wrote:> Hi, Reid > >> Do the values of W and Hg over time for a given q agree between R and >> Excel? >> Not the optimal value of q, just the trajectories for fixed q (trying >> several values for q). > > > If I take the iterative loop out of the function, and ask for values of > Hgmod, Hgtmod, and f, then I get EXACTLY what I get out of Excel. It's > the optimization that seems to be the problem. If I trace the > solutions, R isn't even exploring the full parameter space I tell it to > look in. SO, the iterative loop is correct, and doing what it's > supposed to, since values of p, ACT match exactly what they do in excel- > it's something about how R is examining the possibilities in the > optimization process that is giving me different answers between the two. > > I dunno- I'm going to tinker with it some more tonight. > > Mike > >> Reid >> >> -----Original Message----- >> From: Michael Rennie [mailto:mrennie at utm.utoronto.ca] >> Sent: Wednesday, July 16, 2003 2:47 PM >> To: Huntsinger, Reid >> Subject: RE: [R] Excel can do what R can't????? >> >> >> >> Hi, Reid >> >> At 02:09 PM 7/16/03 -0400, you wrote: >> >R is good at automating specific kinds of complex loops, namely those >> that >> >can be vectorized, or that can be written to draw on otherwise built-in >> >facilities. It's usually reasonable for other kinds of loops but not >> >spectacularly fast. You can write this part in C, though, quite >> easily, and >> >R provides very convenient utilities for this. >> > >> >As for your code: You seem to have a system of equations that relates >> W and >> >Hg to their one-period-ago values. It might clarify things if you coded >> this >> >as a function: input time t values and q, output time t + 1 values. (You >> >wouldn't need any arrays.) Then f would just iterate this function and >> >calculate the criterion. >> >> Wouldn't I still need to loop this function to get it through 365 >> days? Is >> there a big difference, then, between this and what I've got? >> >> >Does the trajectory of (W, Hg) for given q in R seem correct? Does it >> agree >> >with Excel? What does the criterion function look like? You could >> plot it >> in >> >R and perhaps see if the surface is complicated, in which case a simple >> grid >> >search might work for you. >> >> When I give R the values that I get in excel for p, ACT, the function >> solution is actually more precise than what I get in Excel; the values >> for >> p, ACT come back identical (then again, they are exactly what I >> assigned..) But, if I leave R on it's own to find the solution, it keeps >> getting jammed in a particular region. I've never done any function >> plotting in R, but it would help if I could see what kind of surface I >> get >> for f as a function of p, ACT- this would at least force R to examine the >> full range of values specified by the upper and lower limits I've set >> (which it isn't doing under the 'optim' command). >> >> Mike >> >> >> >Reid Huntsinger >> > >> > >> > >> > >> > >> >-----Original Message----- >> >From: Michael Rennie [mailto:mrennie at utm.utoronto.ca] >> >Sent: Wednesday, July 16, 2003 11:18 AM >> >To: Spencer Graves >> >Cc: R-Help; M.Kondrin >> >Subject: Re: [R] Excel can do what R can't????? >> > >> > >> > >> >Hi, Spencer >> > >> >I know I submitted a beastly ammount of code, but I'm not sure how to >> >simplify >> >it much further, and still sucessfully address the problem that i am >> having. >> > >> >The reason being is that the funciton begins >> > >> >f<- function (q) >> > >> >At the top of the iterative loop. This is what takes q and generates >> Wtmod, >> > >> >Hgtmod at the end of the iterative loop. the assignment to f occurs >> at the >> >bottom of the iterative loop. So, yes, the call to f is performing an >> >immediate >> >computation, but based on arguments that are coming out of the iterative >> >loop >> >above it, arguments which depend on q<-(p, ACT). Maybe this is the >> problem; >> > >> >I've got too much going on between my function defenition and it's >> >assignment, >> >but I don't know how to get around it. >> > >> >So, I'm not sure if your example will work- the output from the >> iterative >> >process is Wtmod, Hgtmod, and I want to minimize the difference between >> them >> > >> >and my observed endpoints (Wt, Hgt). The numbers I am varying to reach >> this >> > >> >optimization are in the iterative loop (p, ACT), so re-defining these >> >outputs >> >as x's and getting it to vary these doesn't do me much good unless >> they are >> >directly linked to the output of the iterative loop above it. >> > >> >Last, it's not even that I'm getting error messages anymore- I just >> can't >> >get >> >the solution that I get from Excel. If I try to let R find the >> solution, >> >and >> >give it starting values of c(1,2), it gives me an optimization sulution, >> but >> >an >> >extremely poor one. However, if I give it the answer I got from >> excel, it >> >comes right back with the same answer and solutions I get from excel. >> > >> >Using the 'trace' function, I can see that R gets stuck in a specific >> region >> >of >> >parameter space in looking for the optimization and just appears to give >> up. >> > >> >Even when it re-set itself, it keeps going back to this region, and thus >> >doesn't even try a full range of the parameter space I've defined >> before it >> >stops and gives me the wrong answer. >> > >> >I can try cleaning up the code and see if I can re-submit it, but >> what I am >> >trying to program is so parameter heavy that 90% of it is just defining >> >these >> >at the top of the file. >> > >> >Thank you for the suggestions, >> > >> >Mike >> > >> > >> >Quoting Spencer Graves <spencer.graves at PDF.COM>: >> > >> > > The phrase: >> > > >> > > f <- 1000000000*(((((Wt-Wtmod)^2)/Wt) + >> (((Hgt-Hgtmod)^2)/Hgt))2) ; >> f >> > > >> > > is an immediate computation, not a function. If you want a function, >> > > try something like the following: >> > > >> > > f <- function(x){ >> > > Wt <- x[1] >> > > Wtmod <- x[2] >> > > Hgt <- x[3] >> > > Hgtmod <- x[4] >> > > 1000000000*(((((Wt-Wtmod)^2)/Wt) + (((Hgt-Hgtmod)^2)/Hgt))2) >> > > } >> > > >> > > OR >> > > >> > > f <- function(x, X){ >> > > Wt <- X[,1] >> > > Hgt <- X[,2] >> > > Wtmod <- x[1] >> > > Hgtmod <- x[2] >> > > 1000000000*(((((Wt-Wtmod)^2)/Wt) + (((Hgt-Hgtmod)^2)/Hgt))2) >> > > } >> > > >> > > "par" in "optim" is the starting values for "x". Pass "X" to "f" via >> > > "..." in the call to "optim". >> > > >> > > If you can't make this work, please submit a toy example with >> the >> > > code and error messages. Please limit your example to 3 >> observations, >> > > preferably whole numbers so someone else can read your question in >> > > seconds. If it is any longer than that, it should be ignored. >> > > >> > > hope this helps. >> > > Spencer Graves >> > > >> > > M.Kondrin wrote: >> > > > >?optim >> > > > >> > > > optim(par, fn, gr = NULL, >> > > > method = c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", >> "SANN"), >> > > > lower = -Inf, upper = Inf, >> > > > control = list(), hessian = FALSE, ...) >> > > > >> > > > ..... >> > > > fn: A function to be minimized (or maximized), with first >> > > > argument the vector of parameters over which >> minimization is >> > > > to take place. It should return a scalar result. >> > > > >> > > > Your fn defined as: >> > > > f <- 1000000000*(((((Wt-Wtmod)^2)/Wt) + >> (((Hgt-Hgtmod)^2)/Hgt))2) ; f >> > > > What is its first argument I wonder? >> > > > I think you have just an ill-defined R function (although for >> Excel it >> > > > may be OK - do not know) and optim just chokes on it. >> > > > >> > > > ______________________________________________ >> > > > 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 >> > >> >______________________________________________ >> >R-help at stat.math.ethz.ch mailing list >> >https://www.stat.math.ethz.ch/mailman/listinfo/r-help >> > >> >--------------------------------------------------------------------------->> >> --- >> >Notice: This e-mail message, together with any attachments, contains >> >information of Merck & Co., Inc. (Whitehouse Station, New Jersey, >> >USA) that may be confidential, proprietary copyrighted and/or legally >> >privileged, and is intended solely for the use of the individual or >> entity >> >named on this message. If you are not the intended recipient, and >> >have received this message in error, please immediately return this by >> >e-mail and then delete it. >> >--------------------------------------------------------------------------->> >> --- >> >> 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 >> >> >>---------------------------------------------------------------------------- -->> >> Notice: This e-mail message, together with any attachments, contains >> information of Merck & Co., Inc. (Whitehouse Station, New Jersey, >> USA) that may be confidential, proprietary copyrighted and/or legally >> privileged, and is intended solely for the use of the individual or >> entity >> named on this message. If you are not the intended recipient, and >> have received this message in error, please immediately return this by >> e-mail and then delete it. >>---------------------------------------------------------------------------- -->> > > > 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 > > ______________________________________________ > R-help at stat.math.ethz.ch mailing list > https://www.stat.math.ethz.ch/mailman/listinfo/r-help------------------------------------------------------------------------------ Notice: This e-mail message, together with any attachments, ...{{dropped}}
Dear Professor Ripley, I'm little confused of your reply: do you mean that GRG would not be a standard optimization algorithm, so it couldn't be better than what exist in R ? I'm not a specialist of numerical optimization algorithms, but it seems that GRG is actually implemented in several specialized optimisation toolbox (sure generally commercial), not only the limited one in Excel. And with google, search "GRG generalized reduced gradient" is giving 424 links. -- Fan> I've found that the discussions are interesting, generally > speaking, peoples seem equally confident on R's optim/nlm and > Excel's solver. > > The authors of the algorithm GRG2 (Generalized Reduced Gradient) > are not cited in the documentation of optim(), so I'm wondering if > the optimization algorithm implemented in Excel is "fondamentally" > the same than that in R ?I don't suppose Excel cites the method*s* used in optim() either, but GRG2 is not in the index of my copies of any of the standard texts on numerical optimization. ********** L'ADSL A 20 EUR/MOIS********** Tiscali propose l'ADSL le moins cher du march? : 20 EUR/mois et le modem ADSL offert ! Pour profiter de cette offre exceptionnelle, cliquez ici : http://register.tiscali.fr/adsl/ Offre soumise ? conditions.