princo <- function(s)
{
n <- nrow(s)
m <- apply(s, 2, sum)/n
h <- s+sum(s)/n^2-outer(m, m, "+")
res <- eigen(h)
res$vectors <- res$vectors*matrix(sqrt(res$values), n, n, byrow=TRUE)
colnames(res$vectors) <- names(res$values) <- paste("Axis", 1:n)
rownames(res$vectors) <- paste("Object", 1:n)
res$values <- res$values[res$values > 1e-5]
ax <- length(res$values)
res$vectors <- res$vectors[,1:ax]
if (ax > 2) {
plot(res$vectors[,1:2])
abline(v=0, h=0)
text(res$vectors[,1], res$vectors[,2], 1:n, pos=4, offset=.2)
}
res
}
使用例
s <- matrix(c(
0, -1, -2, -3,
-1, 0, -3, -4,
-2, -3, 0, -1,
-3, -4, -1, 0
), byrow=TRUE, ncol=4)
princo(s)
結果
$values
Axis 1 Axis 2 Axis 3
5.236068 1.000000 0.763932
$vectors
Axis 1 Axis 2 Axis 3
Object 1 -0.8506508 -0.5 0.5257311
Object 2 -1.3763819 0.5 -0.3249197
Object 3 0.8506508 -0.5 -0.5257311
Object 4 1.3763819 0.5 0.3249197
このほか，対象の二次元プロット