Hi all,
I am not used to using the computer to do calculus and have up to
now done my differentiation "by hand" , calling on skills I learned
many years ago and some standard cheat sheets.
My interest at present is in getting the second derivative of a
gaussian, which I did by hand and results in a somewhat messy
result involving terms in sigma^5 .. I have done some spot checks
against R on the first derivative and my algebra seems OK, for the
second derivative I get inflection points where expected via my
algebra and R but I have not yet done a detailed numerical
comparison (so I cannot be 100% sure of my algebra).
What I would like is for the output of D to look at least somewhat
similar to my hand algebra. It occurred to me that if expressions
could be constructed from subexpressions, I could get a more
readable result (from D), as opposed to
> fx2
-(invroot2pi/sigma) * (exp(-0.5 * (((x - mu)/sigma)^2)) * (0.5 *
(2 * (1/sigma * (1/sigma)))) - exp(-0.5 * (((x - mu)/sigma)^2)) *
(0.5 * (2 * (1/sigma * ((x - mu)/sigma)))) * (0.5 * (2 *
(1/sigma * ((x - mu)/sigma)))))
So I tried, in the interests of simplification..
> expr1<-expression( ((x-mu)/sigma)^2)
> expr1
expression(((x - mu)/sigma)^2)> x
[1] 5 15 20> expr2<-expression(exp(-0.5*expr1))
> expr2
expression(exp(-0.5 * expr1))> eval(expr2)
Error in -0.5 * expr1 : non-numeric argument to binary operator
I thought this might mean that I needed to define expr2 as calling
expr1(x) , as follows> expr2<-expression(exp(-0.5*expr1(x)))
> expr2
expression(exp(-0.5 * expr1(x)))> eval(expr2)
Error in eval(expr, envir, enclos) : couldn't find function
"expr1"
which seems to indicate that I am lacking in understanding of
some issues of scope and calling enviroments.
Is what I am trying to do ( get a "readable" machine generated
expression/formula for a higher derivative) possible? and if so, could
someone give me some clues as to how to go about it?
Thank you
----------------------------------------- the original code
> invroot2pi
[1] 0.3989423> fx
expression((invroot2pi/sigma) * exp(-0.5 * (((x -
mu)/sigma)^2)))> fx1<-D(fx,"x") # first derivative
> fx1
-(invroot2pi/sigma) * (exp(-0.5 * (((x - mu)/sigma)^2)) * (0.5 *
(2 * (1/sigma * ((x - mu)/sigma)))))> fx2<-D(fx1,"x") # second derivative
> fx2
-(invroot2pi/sigma) * (exp(-0.5 * (((x - mu)/sigma)^2)) * (0.5 *
(2 * (1/sigma * (1/sigma)))) - exp(-0.5 * (((x - mu)/sigma)^2)) *
(0.5 * (2 * (1/sigma * ((x - mu)/sigma)))) * (0.5 * (2 *
(1/sigma * ((x - mu)/sigma)))))> mu<-0
> sigma<-1
> x<-c(-1,0,1)
> eval(fx1)
[1] 0.2419707 0.0000000 -0.2419707> eval(fx2)
[1] 0.0000000 -0.3989423 0.0000000> # seems OK at first blush .. inflection points should be at mu+-
sigma> x<-seq(-3,3,by=0.1)
> x
[1] -3.0 -2.9 -2.8 -2.7 -2.6 -2.5 -2.4 -2.3 -2.2 -2.1 -2.0 -1.9 -1.8 -1.7 -
1.6
[16] -1.5 -1.4 -1.3 -1.2 -1.1 -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2
-0.1
[31] 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3
1.4
[46] 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8
2.9
[61] 3.0> plot(eval(fx1))
> plot(x,eval(fx1))
> plot(x,eval(fx2))
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at
stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
