1 files changed, 50 insertions, 0 deletions
diff --git a/R/matrix_metrics.R b/R/matrix_metrics.R
new file mode 100644
index 0000000..08521b0
--- /dev/null
+++ b/R/matrix_metrics.R
@@ -0,0 +1,50 @@
+#' Matrix Metrics
+#'
+#' Measure the performance of a covariance or correlation
+#' matrix estimate with various metrics.
+#'
+#' @param S_true The true covariance matrix.
+#' @param R_hat The estimate correlation matrix.
+#' @param strong Eigenvalues strength threshold (default 0.75).
+#'
+#' @return A data.frame with 5 metrics.
+#'
+#' @examples
+#' X_obj <- simulate_data(n = 10, p = 10)
+#' X <- X_obj$X
+#' S_true <- X_obj$S
+#' R_hat <- cor(X)
+#' results <- matrix_metrics(S_true, R_hat)
+#'
+#' @export
+matrix_metrics <- function(S_true, R_hat, strong = 0.75){
+	p <- nrow(S_true)
+	D <- S_true - R_hat
+
+	rmse <- sqrt(mean(D^2))
+	# NOTE: Serra's text is ambiguous on squared-vs-not. This uses ||.||_F / sqrt(p),
+	# the reading under which the identity matrix normalises to exactly 1. If you
+	# want the squared version, use sum(D^2)/p instead. Verify against RSC source.
+	L_F <- sqrt(sum(D^2)) / sqrt(p)
+
+	ev_t <- sort(eigen(S_true, symmetric = TRUE, only.values = TRUE)$values)
+	ev_h <- sort(eigen(R_hat,  symmetric = TRUE, only.values = TRUE)$values)
+	L_S  <- sum(abs(ev_t - ev_h))
+
+	lt <- lower.tri(S_true)          # unique off-diagonal pairs (v < l)
+	t0 <- S_true[lt]
+	th <- R_hat[lt]
+
+	E1 <- sum( (t0 == 0) & (abs(th) > strong) )
+	E2 <- sum( (th == 0) & (abs(t0) > strong) )
+	SS <- sum( (t0 > 0 & th < 0) | (t0 < 0 & th > 0) )
+
+	data.frame(
+		rmse = rmse,
+		L_F = L_F,
+		L_S = L_S,
+		E1 = E1,
+		E2 = E2,
+		SS = SS
+	)
+}