covariance-example.R

n <- 50   # number of observations
d <- 4    # number of dimensions

X <- matrix(rnorm(n * d), nrow = n, ncol = d)

covariance_helper <- function(x, y) {
  x_centered <- x - mean(x)
  y_centered <- y - mean(y)
  sum(x_centered * y_centered) / (length(x) - 1)
}

covariance <- function(X) {

  result <- matrix(ncol(X) * ncol(X), nrow = ncol(X), ncol = ncol(X))

  for (col1 in 1:ncol(X)) {
    for (col2 in 1:ncol(X)) {

      x <- X[, col1]
      y <- X[, col2]

      result[col1, col2] <- covariance_helper(x, y)

    }
  }

  result
}

covariance(X)

### check for correctness
all.equal(cov(X), covariance(X))