[R] how to get higher derivatives with "deriv"
Ravi Varadhan
rvaradhan at jhmi.edu
Fri Jul 30 16:28:34 CEST 2010
I am not as lazy as Peter (but neither am I as good as him)! So, please
allow me to take a stab at completing his exercise:
Deriv <- function(fn, order=1) {
ddf <- fn
body(ddf) <- D(body(f),"x")
if (order > 1) {
for (i in 2:order) {
body(ddf) <- D(body(ddf),"x")
}
}
ddf
}
f <- function (x,alpha) x^alpha
x <- seq(0, 3, length=100)
plot(x, f(x, alpha=1.5), type="l", ylab="")
lines(x, Deriv(f, 1)(x, alpha=1.5), col=2)
lines(x, Deriv(f, 2)(x, alpha=1.5), col=3)
Thanks for the elegant solution, Peter.
Best,
Ravi.
-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On
Behalf Of peter dalgaard
Sent: Friday, July 30, 2010 9:52 AM
To: Wu Gong
Cc: r-help at r-project.org
Subject: Re: [R] how to get higher derivatives with "deriv"
On Jul 29, 2010, at 11:27 PM, Wu Gong wrote:
>
> ## specify the function string
> f.str <- "x^alpha"
>
> ## higher derivatives
> DD <- function(f.str, x = 2, alpha=3,order = 1){
> expr.s <- parse(text=f.str)
> while(order>=1) {
> expr.s <- D(expr.s,"x")
> order <- order-1}
> eval(expr.s)
> }
>
> compute
> DD(f.str,x=1,alpha=0.5,order=1)
>
Somewhat more elegant things can be done, consider for instance
> f
function (x,alpha) x^alpha
> ddf <- f
> body(ddf) <- D(D(body(f),"x"),"x")
> ddf
function (x, alpha)
x^((alpha - 1) - 1) * (alpha - 1) * alpha
> ddf(1,.5)
[1] -0.25
(wrapping this as a function which accepts an order argument is left as an
exercise, because I'm too lazy....)
--
Peter Dalgaard
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: pd.mes at cbs.dk Priv: PDalgd at gmail.com
______________________________________________
R-help at 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.
More information about the R-help
mailing list