[R] Excel can do what R can't?????
Michael Rennie
mrennie at utm.utoronto.ca
Thu Jul 17 02:51:23 CEST 2003
Hmmm.
I tried entering 'Hgtmod = Hgt' at the end of my 'optim' function, but that
didn't help me any- still getting poor optimizations. Perhaps this isn't
working to set a condition, as I was hoping it to. I think that if I can set
the condition Hgmod = Hgt, then it should be able to find reasonable solutions
to this problem set, since that seems to be the trick in my Excel 'solver'
function.
Also, I'm a little hesitant to simplify this too much in terms of reducing the
model, becuase I need this thing to work over 365 iterations, as it does in
excel. I've at least cleaned up some of the commenting so it should be a bit
more straightforward. I tried making the arrays more clear with tabs, but lost
them upon pasting them into this file.
I've included below the code I am currently using with my temp.dat file
attached (it's also below so someone can copy and paste it into a text file
named 'temp.dat')- that's all anyone should need to play around with this if
they are feeling so inclined by cutting and pasting into R.
Thanks again,
Mike
#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
#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
p <- 1 # 0.558626306252032
ACT <- 2 # 1.66764519286918
q<-c(p,ACT)
#introduce function to solve
f <- function (q, Hgtmod)
{
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))
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)
}
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"))
Wtmod<- bioday [length(W),2]
Wtmod
Hgtmod<- bioday [length(Hg),14]
Hgtmod
q
f <- 1000000000*((((Wt-Wtmod)^2)/Wt) + (((Hgt-Hgtmod)^2)/Hgt)) ; f
}
optim(q, f, method = "L-BFGS-B",
lower = c(0.2, 2), upper=c(2, 3),
Hgtmod = Hgt)
#-----------------------------
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
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: temp.dat
Url: https://stat.ethz.ch/pipermail/r-help/attachments/20030716/1285b704/temp.pl
More information about the R-help
mailing list