Hi, I posted this message earlier in "Rmetrics" and I don't know whether I posted in the wrong place, so I'm posting it again in Rhelp. I have a function in x and y and let's call it f(x,y). I need to get the Hessian matrix. i.e I need (d^2f/dx^2), (d^2f/dxdy), (d^2f/dydx), (d^2f/dy^2).I can get these using the D function. now I need to evaluste the hessian matrix for -5<x<5 and -2<y<2 and I also need to print the output. I tried using a for loop,but it is giving me an error. Thanks -- View this message in context: http://www.nabble.com/derivatives-in-R-tp16265419p16265419.html Sent from the R help mailing list archive at Nabble.com.
Lavan <rsumithran at yahoo.com> wrote in news:16265419.post at talk.nabble.com:> > Hi, I posted this message earlier in "Rmetrics" and I don't know > whether I posted in the wrong place, so I'm posting it again in > Rhelp. > > I have a function in x and y and let's call it f(x,y). I need to get > the Hessian matrix. i.e I need (d^2f/dx^2), (d^2f/dxdy), > (d^2f/dydx), (d^2f/dy^2).I can get these using the D function. now I > need to evaluste the hessian matrix for -5<x<5 and -2<y<2 and I also > need to print the output. > > I tried using a for loop,but it is giving me an error.I was wondering if the fact that you provide neither an example nor even the text of the error that could be the explanation for no response. -- David Winsemius
It's pretty hard to debug something when all we are told is "it is giving me an error". Take a specific example, (e.g. which we know will be dead flat at (0,0), just as a check). Here's one way to do it:> fxy <- quote(x^2*log(x^2+y^2+1) + sin(x+y)) > fxyx^2 * log(x^2 + y^2 + 1) + sin(x + y)> Hxy <- function(x, y) NULL > body(Hxy) <- call('cbind',fxx = D(D(fxy, "x"), "x"), fxy = D(D(fxy, "x"), "y"), fyy = D(D(fxy, "y"), "y"))> Hxyfunction (x, y) cbind(fxx = 2 * log(x^2 + y^2 + 1) + 2 * x * (2 * x/(x^2 + y^2 + 1)) + (2 * x * (2 * x/(x^2 + y^2 + 1)) + x^2 * (2/(x^2 + y^2 + 1) - 2 * x * (2 * x)/(x^2 + y^2 + 1)^2)) - sin(x + y), fxy = 2 * x * (2 * y/(x^2 + y^2 + 1)) - x^2 * (2 * x * (2 * y)/(x^2 + y^2 + 1)^2) - sin(x + y), fyy = x^2 * (2/(x^2 + y^2 + 1) - 2 * y * (2 * y)/(x^2 + y^2 + 1)^2) - sin(x + y))> xy <- expand.grid(x = (-5):5, y = (-2):2) > cbind(xy, do.call(Hxy, xy))x y fxx fxy fyy 1 -5 -2 13.0149369 0.87920882 1.87920882 2 -4 -2 11.1066815 0.08339629 0.66389516 3 -3 -2 9.0947006 -0.34667938 -0.40790387 4 -2 -2 7.2919676 0.23085183 -0.65803706 5 -1 -2 5.2801945 1.25223112 0.03000890 6 0 -2 4.1281733 0.90929743 0.90929743 7 1 -2 5.9805455 -0.26964013 0.73035987 8 2 -2 8.0487701 -0.98765432 0.09876543 9 3 -2 9.2121539 -1.45371588 -0.29045058 10 4 -2 10.4767996 -1.27210922 0.03401323 11 5 -2 12.2168303 -0.36334223 1.08110221 12 -5 -1 12.1421622 -0.22454581 1.43526214 13 -4 -1 10.5502143 -0.86015884 0.62132264 14 -3 -1 9.5431203 -0.55845539 0.58204048 15 -2 -1 8.6135278 0.58556445 1.03000890 16 -1 -1 5.9954109 1.79818632 1.13151965 17 0 -1 2.2277653 0.84147098 0.84147098 18 1 -1 5.0861135 -0.88888889 0.22222222 19 2 -1 7.6309368 -1.28591543 0.04741790 20 3 -1 9.3906254 -1.10764453 0.42954555 21 4 -1 11.3680186 -0.23988544 1.43912691 22 5 -1 13.1783802 0.70193281 2.47148014 23 -5 0 11.4744286 -0.95892427 0.96415265 24 -4 0 10.7781363 -0.75680250 1.12555045 25 -3 0 10.5062902 0.14112001 1.94112001 26 -2 0 9.5681733 0.90929743 2.50929743 27 -1 0 6.2277653 0.84147098 1.84147098 28 0 0 0.0000000 0.00000000 0.00000000 ### check 29 1 0 4.5448234 -0.84147098 0.15852902 30 2 0 7.7495784 -0.90929743 0.69070257 31 3 0 10.2240502 -0.14112001 1.65887999 32 4 0 12.2917413 0.75680250 2.63915544 33 5 0 13.3922771 0.95892427 2.88200120 34 -5 1 11.6647752 -0.81167218 0.95787515 35 -4 1 11.6502586 0.04235458 1.72136692 36 -3 1 11.2092202 0.71095032 2.24814040 37 -2 1 9.3138788 0.39702654 1.73035987 38 -1 1 5.0861135 -0.88888889 0.22222222 39 0 1 0.5448234 -0.84147098 -0.84147098 40 1 1 4.1768160 -0.02040854 -0.68707520 41 2 1 8.3312878 0.30332444 0.74776888 42 3 1 11.0567253 0.95514960 2.09564547 43 4 1 12.4680629 1.05768971 2.53917119 44 5 1 12.7009932 0.33428518 1.99409314 45 -5 2 12.4990703 -0.08110221 1.36334223 46 -4 2 12.2953945 0.54648564 1.85260808 47 -3 2 10.8950958 0.22922609 1.39249139 48 -2 2 8.0487701 -0.98765432 0.09876543 49 -1 2 4.2976035 -1.95258210 -0.95258210 50 0 2 2.3095784 -0.90929743 -0.90929743 51 1 2 4.9979545 0.96999110 -0.25223112 52 2 2 8.8055726 1.74445682 0.85556793 53 3 2 11.0125491 1.57116917 1.50994468 54 4 2 11.6655125 0.64222729 1.22272616 55 5 2 11.7009637 -0.43476438 0.56523562>Bill Venables CSIRO Laboratories PO Box 120, Cleveland, 4163 AUSTRALIA Office Phone (email preferred): +61 7 3826 7251 Fax (if absolutely necessary): +61 7 3826 7304 Mobile: +61 4 8819 4402 Home Phone: +61 7 3286 7700 mailto:Bill.Venables at csiro.au http://www.cmis.csiro.au/bill.venables/ -----Original Message----- From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On Behalf Of Lavan Sent: Tuesday, 25 March 2008 10:00 AM To: r-help at r-project.org Subject: [R] derivatives in R Hi, I posted this message earlier in "Rmetrics" and I don't know whether I posted in the wrong place, so I'm posting it again in Rhelp. I have a function in x and y and let's call it f(x,y). I need to get the Hessian matrix. i.e I need (d^2f/dx^2), (d^2f/dxdy), (d^2f/dydx), (d^2f/dy^2). I can get these using the D function. now I need to evaluste the hessian matrix for -5<x<5 and -2<y<2 and I also need to print the output. I tried using a for loop, but it is giving me an error. Thanks -- View this message in context: http://www.nabble.com/derivatives-in-R-tp16265419p16265419.html Sent from the R help mailing list archive at Nabble.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.
hi thanks,
pls see the code below. the elements of matrix Il should be the derivatives
in matrix h (in Boldface) plus some additional terms, here I did calculate
derivatoves in the matrix, but I want R to calculate it for me.
pid<- function (tau1,tau2,sigma=10,n=100,x_bar=20) {
library(akima)
library(lattice)
library(rgl)
a1 <- 1.5
a2 <- 4
t1 <-seq(-1,1);
t2 <-seq(3,4);
l <-length(t1);
m <-length(t2);
N1 <- matrix(NA,l,m);
N2 <- matrix(NA,l,m);
wp <- matrix(NA,l,m);
rr <- matrix(NA,l,m);
N <- matrix(NA,l,m);
mu <- matrix(NA,l,m);
Il <- matrix(NA,l,m);
iIl <- matrix(NA,l,m);
theta <-expression(x+y);
dx <-D(theta, "x");
dy <-D(theta, "y");
d2x <-D(dx, "x");
dxy <-D(dx, "y");
dyx <-D(dy, "x");
d2y <-D(dy, "y");
h <- matrix(c(d2x,dxy,dyx,d2y),nrow=2,ncol=2);
for (i in 1:l)
{
for (j in 1:m)
{
wp[i,j]<-((t1[i]-a1)/tau1)^2+((t2[j]-a2)/tau2)^2
mu[i,j] <-(t1[i] + t2[j]);
Il[i,j]
<-matrix(c(n/(sigma^2)+1/(tau1^2),n/(sigma^2),n/(sigma^2),n/(sigma^2)+1/
tau2^2)),2,2);
iIl <-solve(Il);
ul[i,j] <-
(c(n*(x_bar-mu[i,j])/(sigma^2)-(t1[i]-a1)/(tau1^2),n*(x_bar-mu[i,j])/
(sigma^2)-(t2[j]-a2)/(tau2^2)));
ul <-matrix(c(n*(x_bar-mu)/(sigma^2)-(theta1-a1)/(tau1^2),n*(x_bar-mu)/
(sigma^2)-(theta2-a2)/(tau2^2)),2,1);
ult <-t(ul);
wl[i,j] <-t(ul)%*%solve(Il,ul);
rr[i,j] <-(wp[i,j]/wl[i,j]);
N[i,j] <-(wl[i,j]-wp[i,j])/wl[i,j]
}
}
list("Test Statistics Wp"= wp,"Stat Ratio"=rr,
"N"=N,"N1"=N1,
"N2"=N2, "mean"=mu, "Il"=Il,
"inverse of Il"=iIl)
}
--
View this message in context:
http://www.nabble.com/derivatives-in-R-tp16265419p16267397.html
Sent from the R help mailing list archive at Nabble.com.