Re: [R] Excel can do what R can't?????
On Sun, 3 Aug 2003, Fan wrote:
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.
Thanks in advance
--
Fan
Damon Wischik wrote:
Michael Rennie wrote:
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 solution, 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.
1. Either your function or the Excel solver is wrong. I executed your
source code (which defines f), then ran it over a grid of points, and
plotted the answer, using this code:
xvals - seq(.2,2,by=.2)
yvals - seq(1,3,by=.2)
z - matrix(NA,nrow=length(xvals),ncol=length(yvals))
for (i in 1:length(xvals)) for (j in 1:length(yvals)) {
x - xvals[i]
y - yvals[j]
z[i,j] - f(c(x,y))
}
filled.contour(x=xvals,y=yvals,z=log(z))
Your solution from Excel evaluates to
f(c(.558626306252032,1.66764519286918)) == 0.3866079
while I easily found a point which was much better,
f(c(.4,1)) = 7.83029e-05
You should have tried executing your function over a grid of points, and
plotting the results in a contour plot, to see if optim was working
sensibly. You could do the same grid in Excel and R to verify that the
function you've defined does the same thing in each.
Since your optimization is only over a 2D parameter space, it is easy for
you to plot the results, to see at a glance what the optimum is, and to
work out what is going wrong.
2. Your code executes very slowly because it is programmed inefficiently.
You need to iterate a function to get your final solution, but you don't
need to keep track of all the states you visit on the way. The way R
works, whenever you assign a value to a certain index in a vector, as in
A[i] - 10,
the system actually copies the entire vector. So, in every iteration, you
are copying very many vectors, and this is needlessly slowing down the
program. Also, at the end of each iteration, you define
bio - cbind(W, C, ASMR, SMR, A, F, U, SDA, Gr, Ed, GHg, EGK, Hg)
which creates a matrix. But you only ever use this matrix right at the
end, and so there is no need to create this 365*14 matrix at every single
iteration.
It looks to me as if you took some Excel code and translated it directly
into R. This will not produce efficient R code. Your iterative loop would
be more naturally expressed in R as
f - function(q) {
p - q[[1]]
ACT - q[[2]]
# cat(paste(Trying p=,p, ACT=,ACT,\n,sep=))
state - c(W=Wo,Hg=Hgo)
numdays - length(temps)
for (i in 1:numdays)
state - updateState(state,
jday=temps$jday[i],temp=temps$Temp[i],M=numdays,
p=p,ACT=ACT)
Wtmod - state[[W]]
Hgtmod - state[[Hg]]
(Wt-Wtmod)^2/Wt + (Hgt-Hgtmod)^2/Hgt
}
updateState - function(state,jday,temp,M,p,ACT) {
# Given W[i-1] and Hg[i-1], want to compute W[i] and Hg[i]
W - state[[W]]
Hg - state[[Hg]]
# First compute certain parameters: Vc[i-1] ... Expegk[i-1]
Vc - (CTM-temp)/(CTM-CTO)
Vr - (RTM-temp)/(RTM-RTO)
C - p * CA * W^CB * Vc^Xc * exp(Xc*(1-Vc)) * Pc
ASMR - ACT * RA * W^RB * Vr^Xa * exp(Xa*(1-Vr))
...
# Now find W[i] and Hg[i]
Wnew - if (!(jday==121 Mat==1)) W+Gr/Ef
elseW * (1-GSI*1.2)
Hgnew - a*Hgp*C*(1-Expegk)/(Pc*W*EGK) + Hg*Expegk
c(W=Wnew,Hg=Hgnew)
}
In this code, I do not attempt to keep the entire array in memory. All I
need to know at each iteration is the value of state=(W,Hg) at time i-1,
and from this I compute the new value at time i.
3. You use some thoroughly weird code to read in a table. You should add a
row to the top of your table with variable names, then just use
temps - read.table(TEMP.DAT, header=TRUE)
temps$Vc - (CTM-temps$temp)/(CTM-CTO)
This would also avoid leaving global variables (like Day) hanging around
the
Re: [R] Excel can do what R can't?????
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. __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] Excel can do what R can't?????
On Mon, 4 Aug 2003, Jean Fan wrote: 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 ? Not at all. I meant what I actually said (and I said nothing about which was better). To spell it out: First optim() has five different methods, so Excel's could be at most essentially the same as *one* of them. Second, none of them are called GRG2 by their authors. Third, GRG2 by that name does not appear in the index of my books. 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. So now you will be able to work out how it differs, it at all, from the methods of R (which are fully documented). But looking at just the first reference suggests that GRG is not intended for unconstrained or box-constrained problems, those covered by optim(), and that it is a class of methods rather than an actual algorithm. -- 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. -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: Re: [R] Excel can do what R can't?????
On Mon, 4 Aug 2003, [iso-8859-1] Jean Fan wrote: 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. The code in Excel is actually called GRG2 (the 2 does matter). Unlike any of the methods for optim(), it can handle nonlinear inequality constraints and does not need a feasible initial solution. There's a blurb about it in the NEOS optimisation guide: http://www-fp.mcs.anl.gov/otc/Guide/SoftwareGuide/Blurbs/grg2.html Judging from this blurb, it will be similar to L-BFGS-B for problems with no constraints or box constraints. -thomas Thomas Lumley Assoc. Professor, Biostatistics [EMAIL PROTECTED] University of Washington, Seattle __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] Excel can do what R can't?????
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 ?
Thanks in advance
--
Fan
Damon Wischik wrote:
Michael Rennie wrote:
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 solution, 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.
1. Either your function or the Excel solver is wrong. I executed your
source code (which defines f), then ran it over a grid of points, and
plotted the answer, using this code:
xvals - seq(.2,2,by=.2)
yvals - seq(1,3,by=.2)
z - matrix(NA,nrow=length(xvals),ncol=length(yvals))
for (i in 1:length(xvals)) for (j in 1:length(yvals)) {
x - xvals[i]
y - yvals[j]
z[i,j] - f(c(x,y))
}
filled.contour(x=xvals,y=yvals,z=log(z))
Your solution from Excel evaluates to
f(c(.558626306252032,1.66764519286918)) == 0.3866079
while I easily found a point which was much better,
f(c(.4,1)) = 7.83029e-05
You should have tried executing your function over a grid of points, and
plotting the results in a contour plot, to see if optim was working
sensibly. You could do the same grid in Excel and R to verify that the
function you've defined does the same thing in each.
Since your optimization is only over a 2D parameter space, it is easy for
you to plot the results, to see at a glance what the optimum is, and to
work out what is going wrong.
2. Your code executes very slowly because it is programmed inefficiently.
You need to iterate a function to get your final solution, but you don't
need to keep track of all the states you visit on the way. The way R
works, whenever you assign a value to a certain index in a vector, as in
A[i] - 10,
the system actually copies the entire vector. So, in every iteration, you
are copying very many vectors, and this is needlessly slowing down the
program. Also, at the end of each iteration, you define
bio - cbind(W, C, ASMR, SMR, A, F, U, SDA, Gr, Ed, GHg, EGK, Hg)
which creates a matrix. But you only ever use this matrix right at the
end, and so there is no need to create this 365*14 matrix at every single
iteration.
It looks to me as if you took some Excel code and translated it directly
into R. This will not produce efficient R code. Your iterative loop would
be more naturally expressed in R as
f - function(q) {
p - q[[1]]
ACT - q[[2]]
# cat(paste(Trying p=,p, ACT=,ACT,\n,sep=))
state - c(W=Wo,Hg=Hgo)
numdays - length(temps)
for (i in 1:numdays)
state - updateState(state,
jday=temps$jday[i],temp=temps$Temp[i],M=numdays,
p=p,ACT=ACT)
Wtmod - state[[W]]
Hgtmod - state[[Hg]]
(Wt-Wtmod)^2/Wt + (Hgt-Hgtmod)^2/Hgt
}
updateState - function(state,jday,temp,M,p,ACT) {
# Given W[i-1] and Hg[i-1], want to compute W[i] and Hg[i]
W - state[[W]]
Hg - state[[Hg]]
# First compute certain parameters: Vc[i-1] ... Expegk[i-1]
Vc - (CTM-temp)/(CTM-CTO)
Vr - (RTM-temp)/(RTM-RTO)
C - p * CA * W^CB * Vc^Xc * exp(Xc*(1-Vc)) * Pc
ASMR - ACT * RA * W^RB * Vr^Xa * exp(Xa*(1-Vr))
...
# Now find W[i] and Hg[i]
Wnew - if (!(jday==121 Mat==1)) W+Gr/Ef
elseW * (1-GSI*1.2)
Hgnew - a*Hgp*C*(1-Expegk)/(Pc*W*EGK) + Hg*Expegk
c(W=Wnew,Hg=Hgnew)
}
In this code, I do not attempt to keep the entire array in memory. All I
need to know at each iteration is the value of state=(W,Hg) at time i-1,
and from this I compute the new value at time i.
3. You use some thoroughly weird code to read in a table. You should add a
row to the top of your table with variable names, then just use
temps - read.table(TEMP.DAT, header=TRUE)
temps$Vc - (CTM-temps$temp)/(CTM-CTO)
This would also avoid leaving global variables (like Day) hanging around
the place. Global variables cause confusion: see the next point.
4. Here are some lines taken from your code.
p - NULL
ACT - NULL
#starting values for p, ACT
p - 1
ACT - 2
f - function (q)
{
F[i]- (FA*((comp[i,3])^FB)*(exp(FG*p))*C[i])
# (and ACT is never referred to)
}
Why did you define p-NULL and ACT-NULL at the top? Those definitions are
irrelevant, because they are overridden by p-1 and ACT-2.
In the
Re: [R] Excel can do what R can't?????
Quoting Roger D. Peng [EMAIL PROTECTED]: I'm having a little difficulty understanding this thread. If Excel can do the job correctly and suits your needs, why not just use Excel? Primarily because I don't know how to automate this in excel. The reason for me doing this is I eventually need it to go through 1000 rows of input variables so that I can get output with error associated with it. Plus, people think your cool if you can do it in R. As far as I know, 'optim' cannot optimize a function subject to arbitrary equality constraints. The 'constrOptim' function allows for linear inequality constraints but in general, you will either have to reformulate the problem or add your own penalty into the objective function. Also, just a small note, but using lexical scoping in your problem would eliminate the need to have all those variables defined in the global environment (but otherwise it won't change anything). What's lexical scoping? Mike -- 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 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] Excel can do what R can't?????
For lexical scoping, read some of the better documentation on R. At minimum, check www.r-project.org - search - R site search. spencer graves Michael Rennie wrote: Quoting Roger D. Peng [EMAIL PROTECTED]: I'm having a little difficulty understanding this thread. If Excel can do the job correctly and suits your needs, why not just use Excel? Primarily because I don't know how to automate this in excel. The reason for me doing this is I eventually need it to go through 1000 rows of input variables so that I can get output with error associated with it. Plus, people think your cool if you can do it in R. As far as I know, 'optim' cannot optimize a function subject to arbitrary equality constraints. The 'constrOptim' function allows for linear inequality constraints but in general, you will either have to reformulate the problem or add your own penalty into the objective function. Also, just a small note, but using lexical scoping in your problem would eliminate the need to have all those variables defined in the global environment (but otherwise it won't change anything). What's lexical scoping? Mike __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] Excel can do what R can't?????
Michael Rennie wrote:
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 solution, 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.
1. Either your function or the Excel solver is wrong. I executed your
source code (which defines f), then ran it over a grid of points, and
plotted the answer, using this code:
xvals - seq(.2,2,by=.2)
yvals - seq(1,3,by=.2)
z - matrix(NA,nrow=length(xvals),ncol=length(yvals))
for (i in 1:length(xvals)) for (j in 1:length(yvals)) {
x - xvals[i]
y - yvals[j]
z[i,j] - f(c(x,y))
}
filled.contour(x=xvals,y=yvals,z=log(z))
Your solution from Excel evaluates to
f(c(.558626306252032,1.66764519286918)) == 0.3866079
while I easily found a point which was much better,
f(c(.4,1)) = 7.83029e-05
You should have tried executing your function over a grid of points, and
plotting the results in a contour plot, to see if optim was working
sensibly. You could do the same grid in Excel and R to verify that the
function you've defined does the same thing in each.
Since your optimization is only over a 2D parameter space, it is easy for
you to plot the results, to see at a glance what the optimum is, and to
work out what is going wrong.
2. Your code executes very slowly because it is programmed inefficiently.
You need to iterate a function to get your final solution, but you don't
need to keep track of all the states you visit on the way. The way R
works, whenever you assign a value to a certain index in a vector, as in
A[i] - 10,
the system actually copies the entire vector. So, in every iteration, you
are copying very many vectors, and this is needlessly slowing down the
program. Also, at the end of each iteration, you define
bio - cbind(W, C, ASMR, SMR, A, F, U, SDA, Gr, Ed, GHg, EGK, Hg)
which creates a matrix. But you only ever use this matrix right at the
end, and so there is no need to create this 365*14 matrix at every single
iteration.
It looks to me as if you took some Excel code and translated it directly
into R. This will not produce efficient R code. Your iterative loop would
be more naturally expressed in R as
f - function(q) {
p - q[[1]]
ACT - q[[2]]
# cat(paste(Trying p=,p, ACT=,ACT,\n,sep=))
state - c(W=Wo,Hg=Hgo)
numdays - length(temps)
for (i in 1:numdays)
state - updateState(state,
jday=temps$jday[i],temp=temps$Temp[i],M=numdays,
p=p,ACT=ACT)
Wtmod - state[[W]]
Hgtmod - state[[Hg]]
(Wt-Wtmod)^2/Wt + (Hgt-Hgtmod)^2/Hgt
}
updateState - function(state,jday,temp,M,p,ACT) {
# Given W[i-1] and Hg[i-1], want to compute W[i] and Hg[i]
W - state[[W]]
Hg - state[[Hg]]
# First compute certain parameters: Vc[i-1] ... Expegk[i-1]
Vc - (CTM-temp)/(CTM-CTO)
Vr - (RTM-temp)/(RTM-RTO)
C - p * CA * W^CB * Vc^Xc * exp(Xc*(1-Vc)) * Pc
ASMR - ACT * RA * W^RB * Vr^Xa * exp(Xa*(1-Vr))
...
# Now find W[i] and Hg[i]
Wnew - if (!(jday==121 Mat==1)) W+Gr/Ef
elseW * (1-GSI*1.2)
Hgnew - a*Hgp*C*(1-Expegk)/(Pc*W*EGK) + Hg*Expegk
c(W=Wnew,Hg=Hgnew)
}
In this code, I do not attempt to keep the entire array in memory. All I
need to know at each iteration is the value of state=(W,Hg) at time i-1,
and from this I compute the new value at time i.
3. You use some thoroughly weird code to read in a table. You should add a
row to the top of your table with variable names, then just use
temps - read.table(TEMP.DAT, header=TRUE)
temps$Vc - (CTM-temps$temp)/(CTM-CTO)
This would also avoid leaving global variables (like Day) hanging around
the place. Global variables cause confusion: see the next point.
4. Here are some lines taken from your code.
p - NULL
ACT - NULL
#starting values for p, ACT
p - 1
ACT - 2
f - function (q)
{
F[i]- (FA*((comp[i,3])^FB)*(exp(FG*p))*C[i])
# (and ACT is never referred to)
}
Why did you define p-NULL and ACT-NULL at the top? Those definitions are
irrelevant, because they are overridden by p-1 and ACT-2.
In the body of your function f, in defining F[i], you refer to the
variable p. The only assignment to p is in the line p-1. I strongly
suspect this is an error. Probably you want to refer to q[1]. The best way
to do this (as you can see from my code above) is to define p and ACT at
the beginning of f.
5. Some minor comments on code. It's unwise to use T or F as variable
names in R,
Re: [R] Excel can do what R can't?????
?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 - 10*(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. __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] Excel can do what R can't?????
The phrase:
f - 10*(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]
10*(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]
10*(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 - 10*(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.
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
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 [EMAIL PROTECTED]:
The phrase:
f - 10*(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]
10*(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]
10*(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 - 10*(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.
__
[EMAIL PROTECTED] 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
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
RE: [R] Excel can do what R can't?????
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:[EMAIL PROTECTED] 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:[EMAIL PROTECTED] 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
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:[EMAIL PROTECTED] 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:[EMAIL PROTECTED] 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
Re: [R] Excel can do what R can't?????
Clearly there is much I don't understand about your problem. I suspect you are doing something that isn't correct in R, and the complexity of the problem makes it difficult for you or anyone else to isolate the likely error. Have you tried discarding Hg (or W) and simplifying the problem to estimating a single parameter Hgmod with only 3 observations, say? Then you can program those three numbers in Excel and R and have a better chance of finding the problem. If you still can't find it, give that to the world and maybe someone else can see it quickly. You can't expect others to search a huge haystack for your needle, but if you cut it to a small handful of straw, it is much easier to find the needle just by squeezing a little. If you get the toy problem to work and still can't see why the full problem doesn't, add complexity back into your function a little at a time until you figure it out. In this process, I suggest you modularize your code: Don't write one huge complicated function. Write a series of relatively simple functions, with the more complicated tasks being done by simpler functions. Example: Define the following functions: X1 - function(X, Xmod)(((X-Xmod)^2)/X) Wt.Hg - function(Xmod, Wt, Hgt)(X1(Wt, Xmod[1])+X1(Hgt, Xmod[2])) Then (1e9*Wt.Hg(Xmod, Wt, Hgt)) computes the following line: 10*(Wt-Wtmod)^2)/Wt) + (((Hgt-Hgtmod)^2)/Hgt))2) This may not be a good example, but you need to invent concepts for portions of your code. hope this helps. spencer graves Huntsinger, Reid wrote: 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:[EMAIL PROTECTED] 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:[EMAIL PROTECTED] 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
RE: [R] Excel can do what R can't?????
Hi, Reid and Spencer- I think I've figured something out pretty critical to the problem. Loking at my 'solver' options, I have a condition added that 'Hgtmod = Hgt'. Without this conditional statement, I have to run solver 3-4 times before I get a final solution. MEANING- solver and R, when left to their own devices, suck equally at finding a solution with similar starting points. BUT, given a conditional statement that demands Hgt = Hgtmod, it gives it somewhere to look withing the given parameter space. So, the millon dollar quesiton: Is there any way of setting up a contitional statement like this in 'optim', to specify a solution such that Hgtmod = Hgt? Or, write it into function f? The control statements 'fnscale' and 'parscale' look like candidates, but will only help me if I build Hgtmod into the optimization statement- can I do that? How do you specify multiple parameters? Or should I specify two functions to optimize- one, fucntion f, and a second, something like g- (Hgt = Hgtmod)^2 Can I do that? If this is all I need to do, I am off to the races, and owe you both a beer! I'm going to try some stuff. All is not lost! If you have any ideas, I'd love to hear them. THanks again for everything. Quoting Huntsinger, Reid [EMAIL PROTECTED]: 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. Yes- this was what I was trying to communicate (albeit unsucessfully) in my last e-mail. values of 1,2 or 1.5,1.5 or whatever I try, I get the same answers for f, Wtmod, Hgtmod in Excel as I do in R. Understand, though, that in order to determine the agreement for these values, I'm not using an optimization function of any sort- just plugging in numbers, and seeing what I get back. THat's what I meant about taking the 365-day loop out of the optimization function; even though I am still calculating a value for f, I am no longer defining it at the top of the loop with 'f-function(q)'; simply performing the calculation f - 10*Wt-Wtmod)^2)/Wt) + (((Hgt-Hgtmod)^2)/Hgt)) ; f. In order to find out if the value of f is the same in excel and in R given the same starting parameters of p, ACT. Sorry- I didn't mean to be unclear. 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. This is next on the list- I suspect that this is where I may be running into problems. But, if this is were I am running into problems, then does that mean that Excel is using a better optimization process than R? Since they are recognizing the same parameter space over p, ACT, as evidenced that they offer identical solutions for the same values of p, ACT, then shouldn't they be reaching the same solution? AH! Ureka moment- see top of e-mail Thanks again for your suggestions. Mike Reid Huntsinger -Original Message- From: Spencer Graves [mailto:[EMAIL PROTECTED] 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
Re: [R] Excel can do what R can't?????
I don't expect you to have a complete solution from the simplifications.
I expect you to learn something from the toy problems that can help
you solve the real problems.
Michael Rennie wrote:
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 - 10*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
Re: [R] Excel can do what R can't?????
I'm having a little difficulty understanding this thread. If Excel can do the job correctly and suits your needs, why not just use Excel? As far as I know, 'optim' cannot optimize a function subject to arbitrary equality constraints. The 'constrOptim' function allows for linear inequality constraints but in general, you will either have to reformulate the problem or add your own penalty into the objective function. Also, just a small note, but using lexical scoping in your problem would eliminate the need to have all those variables defined in the global environment (but otherwise it won't change anything). -roger Michael Rennie wrote: Hi, Reid and Spencer- I think I've figured something out pretty critical to the problem. Loking at my 'solver' options, I have a condition added that 'Hgtmod = Hgt'. Without this conditional statement, I have to run solver 3-4 times before I get a final solution. MEANING- solver and R, when left to their own devices, suck equally at finding a solution with similar starting points. BUT, given a conditional statement that demands Hgt = Hgtmod, it gives it somewhere to look withing the given parameter space. So, the millon dollar quesiton: Is there any way of setting up a contitional statement like this in 'optim', to specify a solution such that Hgtmod = Hgt? Or, write it into function f? The control statements 'fnscale' and 'parscale' I don't think 'fnscale' or 'parscale' will help you here at all. If you want a constraint you need to specify it in the objective function directly (something like an additive penalty). look like candidates, but will only help me if I build Hgtmod into the optimization statement- can I do that? How do you specify multiple parameters? Or should I specify two functions to optimize- one, fucntion f, and a second, something like g- (Hgt = Hgtmod)^2 Can I do that? If this is all I need to do, I am off to the races, and owe you both a beer! I'm going to try some stuff. All is not lost! If you have any ideas, I'd love to hear them. THanks again for everything. __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] Excel can do what R can't?????
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
Re: [R] Excel can do what R can't?????
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 -
Re: [R] Excel can do what R can't?????
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 -
Re: [R] Excel can do what R can't?????
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
Re: [R] Excel can do what R can't?????
I thought that you could simplify your code by using something like c(0,W[-length(W)]) as opposed to W[i-1] in a loop, but now I understand it's not that easy. Unless you can analytically simplify the calculation of W in order to vectorize it, it's going to be slow. However, many of the lines don't depend on [i] and not on [i-1]. Therefore you could simplify those as they don't need to be calculated within the loop. HTH, Jerome On July 15, 2003 01:24 pm, Michael Rennie wrote: 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 # #
