目的

主座標分析を行う
R には cmdscale という関数名で用意されている（距離行列を与える）。

使用法

princo(s)
print.princo(obj, ax=res$ax, digits=5)
plot.princo(obj, labels=FALSE, text.cex=0.7, ...)

引数

s         類似度行列（正方行列。対称行列でなくてはならない）
ax        出力する次元数
digits    結果の出力桁数
labels    ラベルを描くかどうか
text.cex  ラベルのフォントの大きさ
...        plot への引数

ソース

インストールは，以下の 1 行をコピーし，R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/princo.R", encoding="euc-jp")

# 主座標分析を行う
princo <- function(  s)                              # 類似度行列
{
        if (is.data.frame(s)) {
                s <- as.matrix(s)
        }
        n <- nrow(s)                                 # 行列のサイズ
        object.names <- colnames(s)
        if (is.null(object.names)) {
                object.names <- paste("対象", 1:n, sep="")
        }
        m <- colSums(s)/n
        h <- s+sum(s)/n^2-outer(m, m, "+")
        res <- eigen(h)                                      # 固有値・固有ベクトルを求める
        values <- res$values[res$values > 1e-5]              # 解の個数を決める
        ax <- length(values)
        vectors <- t(t(res$vectors[,1:ax])*sqrt(values))# 対象に数値を与える
        colnames(vectors) <- names(values) <- paste("解", 1:ax, sep="")
        rownames(vectors) <- object.names
        return(structure(list(ax=ax, n=n, values=values, vectors=vectors), class="princo"))
}
# print メソッド
print.princo <- function(    res,                    # princo が返すオブジェクト
                                ax=res$ax,              # 何次元までの解を出力するか
                                digits=5)               # 表示桁数
{
        ax <- min(ax, res$ax)
        val <- res$values
        val2 <- val/sum(val)
        val <- rbind(val, val2, cumsum(val2))
        rownames(val) <- c("固有値", "寄与率", "累積寄与率")
        print(val[, 1:ax], digits=digits)
        cat("\nベクトル\n\n")
        print(res$vectors[, 1:ax], digits=digits)
}
# plot メソッド
plot.princo <- function(res,                         # princo が返すオブジェクト
                        labels=FALSE,                   # ラベルを描くかどうか
                        text.cex=0.7,                   # ラベルのフォントの大きさ
                        ...)                            # plot への引数
{
        if (res$ax >= 2) {                              # 二次元以上の解が得られたら，
                plot(res$vectors[,1:2], ...)            # 二次元の図を描く
                abline(v=0, h=0)
                old <- par(xpd=TRUE)
                if (labels) {
                        text(res$vectors[,1], res$vectors[,2], rownames(res$vectors),
                             pos=4, offset=.2, cex=text.cex)
                }
                par(old)
        }
        else {
                warning("解が一次元なので，二次元配置図は描けません。")
        }
}


使用例

> s <- matrix(c(	# 4行4列の類似度行列例（ファイルから読んでも良い）
	0, -1, -2, -3,
	-1, 0, -3, -4,
	-2, -3, 0, -1,
	-3, -4, -1, 0
), byrow=TRUE, ncol=4)

> (ans <- princo(s)) # print メソッドで結果を出力する

               解1     解2     解3
固有値     5.23607 1.00000 0.76393
寄与率     0.74801 0.14286 0.10913
累積寄与率 0.74801 0.89087 1.00000

ベクトル

           解1  解2      解3
対象1 -0.85065  0.5  0.52573
対象2 -1.37638 -0.5 -0.32492
対象3  0.85065  0.5 -0.52573
対象4  1.37638 -0.5  0.32492

インストールは，以下の 1 行をコピーし，R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/similarity.matrix.R", encoding="euc-jp")

# データ行列から類似度行列を作る
similarity.matrix <- function(       dat,                            # データ行列
                                power=1,                        # 距離のべき乗数（ユークリッド二乗距離なら 2 を指定）
                                ...)                            # dist 関数への引数（距離の種類などの指定）
{
        d <- as.matrix(dist(dat, ...))
        if (power != 1) {
                d <- d^power
        }
        d <- -d
        diag(d) <- 0
        return(d)
}


# plot メソッドでグラフ表示する

> a <- similarity.matrix(iris[1:4]) # 類似度行列を作る
> b <- princo(a)
> plot(b, labels=TRUE, col=(1:3)[as.integer(iris[,5])])



> x <- princo(similarity.matrix(iris[1:4], power=2)) # ユークリッドの二乗距離から類似度行列を作る
> plot(x, labels=TRUE, col=(1:3)[as.integer(iris[,5])], xlim=c(-5, 6), cex=0.6)



 解説ページ