1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
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
)
}
|
