You can use the numDeriv
package to numerically estimate the gradient,
or Ryacas
if you want an exact value.
library(numDeriv)
f <- function(x) c( x, x[1] + x[2], ( x[1] + x[2] ) / x[3] )
x0 <- c(1,1,1)
V <- diag(1:3)
J <- jacobian(f, x0)
J
# [,1] [,2] [,3]
# [1,] 1 0 0
# [2,] 0 1 0
# [3,] 0 0 1
# [4,] 1 1 0
# [5,] 1 1 -2
f(x0) # (Biased) estimator of the mean of f(X)
J %*% V %*% t(J) # Estimator of the variance of f(X)
Including the hessian will give a better approximation.