R - functions of an estimated vector and its covariance matrix

Question

Given a estimates for a vector of variables x[1:3] and its covariance matrix V[1:3,1:3], I seek a utility that will apply linear or nonlinear functions to two or more variables, such as

x[4] = x[1] + x[2]

x[5] = x[4] / x[3].

The estimated values of x[4] and x[5] require simple algebra.

The covariance matrix including the linear transformation x[4] is simply (H * V * H'), where

H =

| 1 0 0 | 
| 0 1 0 |
| 0 0 1 |
| 1 1 0 |

The covariance matrix adding x[5] could be estimated with the first-order Taylor series approximation for terms a and b:

H =

| 1 0 0 0 | 
| 0 1 0 0 |
| 0 0 1 0 |
| 1 1 0 0 |
| 0 0 a b |

In concept, I know how the algorithm should work. But it will take a lot of coding, especially if I attempt some kind of general-purpose equation parser in the user-interface.

Are there any existing R libraries that work for this problem?

Answer 1

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.