Генерация матрицы корреляции населения с ошибкой модели на Такера, Купмана, Линна (1969) и Хонга (1991)

Question

Я пытаюсь воссоздать процесс создания реалистичной корреляционной матрицы для исследования Монте-Карло, включающего факторный анализ.Я попытался выполнить шаги, изложенные Хонгом (1999), который можно посмотреть здесь: https://link.springer.com/content/pdf/10.3758/BF03200754.pdf.

Используя R, я до сих пор придумал это, но я продолжаю получать популяционную корреляциюматрица с корреляциями> 1, и это довольно далеко от тех, о которых сообщил Хонг.Ниже приведен код R, который я использовал:

#factor loadings
V <- matrix(c(.673, .020, -.039, .606, .136, .008, .711, .012, .068, .638, -.085, -.007, -.058, .529, -.019, .008, .727, .061, .035, .609, -.053, .197, .618, .174, -.002, -.002, .601, -.003, .001, .626), ncol = 3, nrow = 10, byrow = TRUE)

#phi matrix (major factor intercorrelations = 0.30)
diag <- rep(1.0, 3)
phi <- matrix(rep(0.3, 9), ncol = 3, nrow = 3)
diag(phi) <- diag

#Y matrix (correlation between major and minor factors = 0.30)
Y <- matrix(rep(0.3, 150), nrow = 3, ncol = 50)

#gamma matrix (minor factor intercorrelations = 0.30)
gamma <- matrix(rep(0.3, 2500), nrow = 50, ncol = 50)
diag <- rep(1.0, 50)
diag(gamma) <- diag

#minor factor loadings (std. dev decreases with each loading)
W <- matrix(NA, nrow = 10, ncol = 50)
for(i in 1:50) {
  mu = 0
  n = i - 1
  sigma = 0.80^n
  W[, i] <- rnorm(10, mu, sigma)
}

#square minor factor loadings
Wsq <- matrix(sapply(W, W2Wsq), nrow = 10, ncol = 50)

#ensure minor factors account for 10% of variance
sums <- matrix(rowSums(Wsq), ncol = 1)
sums <- sapply(sums, normalize)
Wsq <- matrix(sums*Wsq, nrow = 10, ncol = 50)
W <- matrix(sapply(Wsq, Wsq2W), nrow = 10, ncol = 50)

#J matrix
J <- cbind(V, W)

#B matrix
Bt <- cbind(phi, Y)
Bb <- cbind(t(Y), gamma)
B <- rbind(Bt, Bb)

#Get D matrix
D <- diag(10)
Dp <- J%*%B%*%t(J)
Dm <- matrix(rep(0, 100), ncol = 10, nrow = 10)
diag(Dm) <- diag(Dp)
D <- D - Dm

#Get population correlation matrix
P <- Dp + D

А функции, используемые в приведенном выше коде, следующие:

W2Wsq <- function(x) {
x^2
}

Wsq2W <- function(x) {
sqrt(x)
}

normalize <- function(x) {
1/(10*x)
}