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