[R] The gradient of a multivariate normal density with respect to its parameters

Ravi Varadhan RVaradhan at jhmi.edu
Mon Jun 22 16:48:29 CEST 2009


Hi Karl,

I am not aware of any.  May I ask what your purpose is?  You don't really
need this if you are going to use it in optimization, since most optimizers
use a simple finite-difference approximation if you don't provide the
gradient.  Using the numerical approximation from "numDeriv" will be quite
time-consuming in an optimization routine, since numDeriv uses a high-order
Richardosn extrapolation to compute an accurate approximation of the
gradient.

Regardless of your purpose, there is a small bug in your function.  You
should change   `dmvnorm(cbind(x,y),mu,sig)'  to
`dmvnorm(cbind(xx,yy),mu,sig)'.  Also, taking the derivative of log-density
might be better.

Ravi.


----------------------------------------------------------------------------
-------

Ravi Varadhan, Ph.D.

Assistant Professor, The Center on Aging and Health

Division of Geriatric Medicine and Gerontology 

Johns Hopkins University

Ph: (410) 502-2619

Fax: (410) 614-9625

Email: rvaradhan at jhmi.edu

Webpage:
http://www.jhsph.edu/agingandhealth/People/Faculty_personal_pages/Varadhan.h
tml

 

----------------------------------------------------------------------------
--------


-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On
Behalf Of Karl Ove Hufthammer
Sent: Monday, June 22, 2009 10:10 AM
To: r-help at stat.math.ethz.ch
Subject: [R] The gradient of a multivariate normal density with respect to
its parameters

Does anybody know of a function that implements the derivative (gradient) of
the multivariate normal density with respect to the *parameters*?

It's easy enough to implement myself, but I'd like to avoid reinventing the
wheel (with some bugs) if possible. Here's a simple example of the result
I'd like, using numerical differentiation:

library(mvtnorm)
library(numDeriv)
f=function(pars, xx, yy)
{
  mu=pars[1:2]
  sig1=pars[3]
  sig2=pars[4]
  rho=pars[5]
  sig=matrix(c(sig1^2,rho*sig1*sig2,rho*sig1*sig2,sig2^2),2)
  dmvnorm(cbind(x,y),mu,sig)
}

mu1=1
mu2=2
sig1=3
sig2=4
rho=.5
x=2 # or a x vector
y=3 # or a y vector
jacobian(f,c(mu1,mu2,sig1,sig2,rho),xx=x,yy=y)
# (Can replace 'jacobian' with 'grad' if x and y have length 1.)

--
Karl Ove Hufthammer

______________________________________________
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