Karl,
I have written a function to compute the gradient of a multivariate normal
density (derivatives at a given point, with respect to the parameters of the
density). It is in the atttached file. I would appreciate any feedback.
There are couple of interesting points:
1. Note is that there is a minor adjustment needed to the formulas given by
Dwyer (Table 2, Eqs. (11.7 and 11.8) to account for the symmetry in the
derivatives of the covariance parameters. Without this adjustment the
derivatives of covariances are off by a factor of 2.
2. I couldn't figure out how to avoid inverting the covariance matrix for
Eq. (11.7), using `solve(sigma)'. I would appreciate knowing if there is a
trick for doing this.
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: rvarad...@jhmi.edu
Webpage:
http://www.jhsph.edu/agingandhealth/People/Faculty_personal_pages/Varadhan.h
tml
----------------------------------------------------------------------------
--------
-----Original Message-----
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
Behalf Of Karl Ove Hufthammer
Sent: Wednesday, June 24, 2009 3:21 AM
To: r-h...@stat.math.ethz.ch
Subject: Re: [R] The gradient of a multivariate normal density with
respecttoits parameters
Ravi Varadhan:
> You may want to look at the paper by Dwyer on "Some applications of
> matrix derivatives in multivariate analysis" (JASA 1967), especially
> the Table 2 on p. 617.
Thanks! That's a *very* useful paper. Note there also is a small corrigenda
available:
Corrigenda: Some Applications of Matrix Derivatives in Multivariate Analysis
Paul S. Dwyer Journal of the American Statistical Association, Vol. 62, No.
320 (Dec., 1967), p. 1518 Stable URL: http://www.jstor.org/stable/2283811
--
Karl Ove Hufthammer
______________________________________________
R-help@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.
library(mvtnorm)
library(numDeriv)
####################################################
mvnorm.grad <- function(x, mu, sigma, log=FALSE) {
# x = a vector denoting location at which gradient is computed
# mu = a vector denoting mean of multivariate normal
# sigma = covariance matrix
# log = logical variable indicating whether density or log-density should be
used
# Note: gradient is computed with respect to location parameters, and with
respect to variances and covariances
#
# Author: Ravi Varadhan, Johns Hopkins University, June 24, 2009
#
siginvmu <- c(solve(sigma, x-mu))
mu.grad <- siginvmu
sigma.grad <- tcrossprod(siginvmu) - solve(sigma) # is there a way to avoid
inverting `sigma'?
diag(sigma.grad) <- 0.5 * diag(sigma.grad) # this is required to account for
symmetry of covariance matrix
if (log) list(mu.grad=mu.grad, sigma.grad=sigma.grad) else
{
f <- dmvnorm(x, mean=mu, sigma=sigma, log=FALSE)
list(mu.grad=mu.grad*f, sigma.grad=sigma.grad*f)
}
}
####################################################
# application to a bivariate normal
#
f=function(pars, xx, yy, log=FALSE)
# Numerical gradient using "numDeriv" for bivariate normal
#
{
mu = pars[1:2]
sig1 = pars[3]
sig12 = pars[4]
sig2 = pars[5]
sig = matrix(c(sig1,sig12,sig12,sig2),2)
dmvnorm(c(xx,yy),mean=mu,sigma=sig, log=log)
}
mu1=1
mu2=2
sig1=3^2
sig2=4^2
sig12=.5 * sqrt(sig1 * sig2)
sig <- matrix(c(sig1,sig12,sig12,sig2), 2, 2)
xx=rnorm(1) # or a x vector
yy=rnorm(1) # or a y vector
ans1 <- jacobian(func=f, x=c(mu1,mu2,c(sig)[-3]), xx=xx, yy=yy)
ans2 <- mvnorm.grad(x=c(xx, yy), mu=c(mu1, mu2), sigma=sig)
ans1
ans2
all.equal(c(ans1), c(ans2$mu, ans2$sigma.grad[lower.tri(ans2$sigma.grad,
diag=TRUE)]))
#########
# application to a 4-dim multivariate normal
mu <- rnorm(4, mean=c(1,2,4,8))
sigma <- matrix(rexp(16), 4, 4)
sigma <- sigma %*% t(sigma)
x0 <- rnorm(4, mean=c(1,2,4,8))
ans <- mvnorm.grad(x=x0, mu=mu, sigma=sigma)
ans
______________________________________________
R-help@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.
