[R] ?subexpressions, D, deriv

[John Aitchison]

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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._