On Thu, 17 Sep 2009, [ISO-8859-1] Noela Sánchez wrote:
Hi all,
I need to automate a process in order to prepare a a big loop in the future
but I have a problem with the *command function*
First I fit a model with lm
model1<-lm(data2[,2]~data2[,1]+I(data2[,1]^2)+I(data2[,1]^3)+I(data2[,1]^4))
I extract the coefficients to build the polynomial.
coef<-as.matrix(model1$coefficients)
In the next step I need to define the polynomial to derive it. If I write
the coefficients manually (writing the numbers by hand) the deriv command
works fine!
bb<-deriv(~2847.22015 -463.06063*x+ 25.43829*x^2 -0.17896*x^3, namevec="x",
function.arg=40)
I would like to automate this step by being able to extract the coefficients
from the linear model and adding them into the polynomial (and not write
them by hand)!
But if I build the polynomial with the function(x) command calling the *
coef* values, the numeric values are not interpreted, the command function
does not read properly the coefficients from the linear model.
fun<-function(x) coef[1]+coef[2]*x+coef[3]*x^2+coef[4]*x^3
fun
function(x) coef[1]+coef[2]*x+coef[3]*x^2+coef[4]*x^3
How can i avoid to write the values of the coefficients by hand??
As is often the case there are two parts to this answer: how to do what you
asked, and how to do what you want.
You can use bquote() to get the value of coef[1] into the expression, so
deriv(bquote(.(coef[1])+.(coef[2])*x+.(coef[3])*x^2+.(coef[4])*x^3),"x",
function.arg="x")
will give the derivative. Bill Venables would tell you to use Horner's Rule
and write the polynomial as
coef[1]+x*(coef[2]+x*(coef[3]+x*coef[4]))
to get better speed and numerical stability, but the same trick still works.
However, you really don't need deriv() to differentiate a polynomial, and so
you can use ordinary lexical scope rather than substitution, for a much tidier
answer
derivative <- function(coef){
function(x) coef[2]+x*(2*coef[3]+x*3*coef[4])
}
-thomas
Thomas Lumley Assoc. Professor, Biostatistics
tlum...@u.washington.edu University of Washington, Seattle
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
