Recursive algorithms for computing variance and mean

These algorithms are used in kernel_adapt() to simplify variance-covariance recalculation at every step of the algorithm.

cov_recursive(
  X_t,
  Cov_t,
  Mean_t_prev,
  t.,
  Mean_t = NULL,
  eps = 0,
  Sd = 1,
  Ik = diag(ncol(Cov_t))
)

mean_recursive(X_t, Mean_t_prev, t.)

Arguments

X_t

Last value of the sample

Cov_t

Covariance in t

t.

Sample size up to t-1.

Mean_t, Mean_t_prev

Vectors of averages in time t and t-1 respectively.

Sd, eps, Ik

See kernel_adapt().

Details

The variance covariance algorithm was described in Haario, Saksman and Tamminen (2002).

References

Haario, H., Saksman, E., & Tamminen, J. (2001). An adaptive Metropolis algorithm. Bernoulli, 7(2), 223–242. https://projecteuclid.org/euclid.bj/1080222083

Examples

# Generating random data (only four points to see the difference)
set.seed(1231)
n <- 3
X <- matrix(rnorm(n*4), ncol = 4)

# These two should be equal
mean_recursive(
  X_t         = X[1,],
  Mean_t_prev = colMeans(X[-1,]),
  t.          = n - 1
)
#> [1] -0.3815965  0.6974209  0.2725030  0.3904914
colMeans(X)
#> [1] -0.3815965  0.6974209  0.2725030  0.3904914

# These two should be equal
cov_recursive(
  X_t         = X[1, ], 
  Cov_t       = cov(X[-1,]), 
  Mean_t      = colMeans(X),
  Mean_t_prev = colMeans(X[-1, ]),
  t           = n-1
)
#>           [,1]      [,2]      [,3]      [,4]
#> [1,] 1.5745595 1.4188822 0.7010067 0.9173480
#> [2,] 1.4188822 1.5675687 1.0409572 0.8027636
#> [3,] 0.7010067 1.0409572 0.8917114 0.3745821
#> [4,] 0.9173480 0.8027636 0.3745821 0.5364269
cov(X)
#>           [,1]      [,2]      [,3]      [,4]
#> [1,] 1.5745595 1.4188822 0.7010067 0.9173480
#> [2,] 1.4188822 1.5675687 1.0409572 0.8027636
#> [3,] 0.7010067 1.0409572 0.8917114 0.3745821
#> [4,] 0.9173480 0.8027636 0.3745821 0.5364269

# Speed example -------------------------------------------------------------
set.seed(13155511)
X <- matrix(rnorm(1e3*100), ncol = 100)

ans0 <- cov(X[-1,])
t0 <- system.time({
  ans1 <- cov(X)
})

t1 <- system.time(ans2 <- cov_recursive(
  X[1, ], ans0,
  Mean_t      = colMeans(X),
  Mean_t_prev = colMeans(X[-1,]),
  t. = 1e3 - 1
))

# Comparing accuracy and speed
range(ans1 - ans2)
#> [1] -2.220446e-16  2.220446e-16
t0/t1
#>    user  system elapsed 
#>       8     NaN       8