判別分析(二次の判別関数)     Last modified: Aug 27, 2009

目的

二次の判別関数を使った判別分析を行う。
(R の MASS パッケージに qda, predict.qda がある)

使用法

quad.disc(data, group, func.name=c("solve", "ginv"))
print.quad.disc(obj, digits=5)
summary.quad.disc(obj, digits=5)
plot.quad.disc(obj, which=c("boxplot", "barplot", "scatterplot"), nclass=20, pch=1:obj$ngroup, col=1:obj$ngroup, xpos="topright", ypos=NULL, ...)

引数

data        データ行列(行がケース,列が変数)
group       各ケースがどの群であるかを表す変数
func.name   逆行列を計算する関数名(solve か ginv)。
            普通に逆行列が求まる場合にはどちらを使っても結果は同じ
            省略時は solve
            ムーア・ペンローズ型一般化逆行列を使うときは ginv

obj         quad.disc 関数が返すオブジェクト
digits      結果の表示桁数

which       2 群判別の場合に,描画するグラフの種類
nclass      barplot のときのおよその階級数
pch         scatterplot を描く記号
col         scatterplot の記号の色
xpos, ypos  scatterplot の凡例の位置

...         boxplot, barplot, scatterplot に引き渡す引数

ソース

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

# 二次の判別関数
quad.disc <- function(       data,                                           # 説明変数データ行列
                        group,                                          # グループを表すベクトル
                        func.name=c("solve", "ginv"))                   # 逆行列を求める関数
{
        inverse <- if (match.arg(func.name) == "solve") solve else { library(MASS); ginv}
        data <- as.data.frame(data)
        if (is.null(colnames(data))) {
                colnames(data) <- paste("Var", 1:p, sep="")
        }
        vname <- colnames(data)
        group <- as.factor(as.matrix(group))                         # 群を表す変数は factor にする
        OK <- complete.cases(data, group)                            # 欠損値を持つケースを除く
        data <- as.matrix(data[OK,])
        group <- group[OK]
        p <- ncol(data)                                                      # 説明変数の個数
        n <- nrow(data)                                                      # データの個数
        n.i <- table(group)                                          # 各群の例数
        g.name <- names(n.i)                                         # 群の名前
        k <- length(n.i)                                             # 群の個数
        group.means <- matrix(unlist(by(data, group, colMeans)), p)  # 各群・各変数の平均
        vars <- array(unlist(by(data, group, var)), c(p, p, k))              # 各群の分散・共分散行列
        inv.vars <- array(apply(vars, 3, inverse), c(p, p, k))               # 各群の分散・共分散行列の逆行列

        scores <- sapply(1:k, function(i) {                          # 各ケースの各群からの距離
                temp <- t(data)-group.means[,i];
                sapply(1:n, function (j) temp[,j]%*%inv.vars[,,i]%*%temp[,j])
                }
        )
        p.values <- pchisq(scores, p, lower.tail=FALSE)                      # 各ケースが各群に属する確率
        prediction <- as.factor(g.name[apply(p.values, 1, order)[k,]])       # どの群に属するか判別
        correct <- ifelse(prediction == group, TRUE, FALSE)          # 判別が正しいか?
        correct.table <- table(group, prediction)                    # 判別結果概括表
        correct.rate <- sum(diag(correct.table))/n*100                       # 正判別率
# add names
        colnames(group.means) <- colnames(scores) <- colnames(p.values) <- dimnames(vars)[[3]] <- dimnames(inv.vars)[[3]] <- g.name
        colnames(vars) <- rownames(vars) <- colnames(inv.vars) <- rownames(inv.vars) <- rownames(group.means) <- vname
        rownames(scores) <- rownames(p.values) <- paste("case", 1:n)
        return(structure(list(group.means=group.means, vars=vars, inv.vars=inv.vars, scores=scores, p.values=p.values, prediction=prediction, correct=correct, correct.table=correct.table, correct.rate=correct.rate, group=group, ngroup=k), class="quad.disc"))
}
# print メソッド
print.quad.disc <- function( obj,                                    # quad.disc 関数が返すオブジェクト
                                digits=5)                               # 結果の表示桁数
{
        cat("\n判別結果\n\n")
        print(obj$correct.table)
        cat(sprintf("\n正判別率 = %.1f %%\n", obj$correct.rate))
}
# summary メソッド                                                      # すべての結果を表示する
summary.quad.disc <- function(       obj,                                    # quad.disc が返すオブジェクト
                                digits=5)                               # 結果の表示桁数
{
        print.default(obj, digits=digits)
}
# plot メソッド
plot.quad.disc <- function(  obj,                                    # quad.disc 関数が返すオブジェクトの
                                which=c("boxplot", "barplot", "scatterplot"),   # 箱髭図か棒グラフか散布図かの選択
                                nclass=20,                              # barplot の場合のおよその階級数
                                pch=1:obj$ngroup,                       # scatterplot を描く記号
                                col=1:obj$ngroup,                       # scatterplot の記号の色
                                xpos="topright", ypos=NULL,             # scatterplot の凡例の位置
                                ...)                                    # boxplot, barplot, plot に引き渡す引数
{
        if (obj$ngroup == 2) {
                group <- obj$group
                score <- obj$score[,1]-obj$score[,2]
                which <- match.arg(which)
                if (which == "boxplot") {                               # boxplot
                        plot(score ~ group, xlab="群", ylab="二乗距離の差", ...)
                }
                else if (which == "barplot") {                          # barplot
                        tbl <- table(group, cut(score,
                                        breaks=pretty(score, n=nclass)))
                        barplot(tbl, beside=TRUE, legend=TRUE, xlab="二乗距離の差", ...)
                }
                else {                                                  # scatterplot 各群の重心までの二乗距離
                        group <- obj$group
                        group.levels <- levels(group)
                        score1 <- obj$score[,1]
                        score2 <- obj$score[,2]
                        max1 <- max(score1)
                        max2 <- max(score2)
                        max0 <- max(max1, max2)
                        plot(score1, score2, col=col[as.integer(group)], pch=col[as.integer(group)],
                                xlim=c(0, max0), xlab=paste(group.levels[1], "の重心への二乗距離"),
                                ylim=c(0, max0), ylab=paste(group.levels[2], "の重心への二乗距離"), asp=1, ...)
                        abline(0, 1, lty=2)
                        text(max1, max2/2, paste(group.levels[2], "に判別される領域"), pos=2)
                        text(0, max2+strheight("H")*1.5, paste(group.levels[1], "に判別される領域"), pos=4)
                        legend(x=xpos, y=ypos, legend=group.levels, col=col, pch=pch)
                }
        }
        else {
                warning("3群以上の場合にはグラフ表示は用意されていません")
        }
}


使用例

data(iris)	# Fisher のアヤメのデータ
result <- quad.disc(iris[1:4], iris[5])	# quad.disc 関数からの出力は自動的には出力されないので注意
result	# 結果を出力する

出力結果例

> (result <- quad.disc(iris[1:4], iris[5]) ) # print による結果出力

判別結果

            prediction
group        setosa versicolor virginica
  setosa         50          0         0
  versicolor      0         47         3
  virginica       0          0        50

正判別率 = 98.0 %

> summary(result) # summary による結果出力

$group.means	# 各群ごとの平均値

             setosa versicolor virginica
Sepal.Length  5.006      5.936     6.588
Sepal.Width   3.428      2.770     2.974
Petal.Length  1.462      4.260     5.552
Petal.Width   0.246      1.326     2.026


$vars	# 各群の分散・共分散行列

, , setosa

             Sepal.Length Sepal.Width Petal.Length Petal.Width
Sepal.Length     0.124249    0.099216    0.0163551   0.0103306
Sepal.Width      0.099216    0.143690    0.0116980   0.0092980
Petal.Length     0.016355    0.011698    0.0301592   0.0060694
Petal.Width      0.010331    0.009298    0.0060694   0.0111061

, , versicolor

             Sepal.Length Sepal.Width Petal.Length Petal.Width
Sepal.Length     0.266433    0.085184     0.182898    0.055780
Sepal.Width      0.085184    0.098469     0.082653    0.041204
Petal.Length     0.182898    0.082653     0.220816    0.073102
Petal.Width      0.055780    0.041204     0.073102    0.039106

, , virginica

             Sepal.Length Sepal.Width Petal.Length Petal.Width
Sepal.Length     0.404343    0.093763     0.303290    0.049094
Sepal.Width      0.093763    0.104004     0.071380    0.047629
Petal.Length     0.303290    0.071380     0.304588    0.048824
Petal.Width      0.049094    0.047629     0.048824    0.075433

$inv.vars	# 各群の分散・共分散行列の逆行列

, , setosa

             Sepal.Length Sepal.Width Petal.Length Petal.Width
Sepal.Length      18.9434    -12.4048      -4.5002     -4.7761
Sepal.Width      -12.4048     15.5705       1.1111     -2.1041
Petal.Length      -4.5002      1.1111      38.7762    -17.9350
Petal.Width       -4.7761     -2.1041     -17.9350    106.0459

, , versicolor

             Sepal.Length Sepal.Width Petal.Length Petal.Width
Sepal.Length       9.5028     -3.6762      -8.6317      6.4545
Sepal.Width       -3.6762     19.7110       2.1160    -19.4803
Petal.Length      -8.6317      2.1160      19.8038    -26.9372
Petal.Width        6.4545    -19.4803     -26.9372     87.2448

, , virginica

             Sepal.Length Sepal.Width Petal.Length Petal.Width
Sepal.Length      10.5339     -3.4797      -9.9601      1.7882
Sepal.Width       -3.4797     15.8754       1.1027     -8.4729
Petal.Length      -9.9601      1.1027      13.4058     -2.8909
Petal.Width        1.7882     -8.4729      -2.8909     19.3141


$scores	# ユークリッドの二乗距離

             setosa versicolor virginica
case 1      0.44911  114.80449 182.93591
case 2      2.08109   83.31536 153.97495
case 3      1.28434   94.92043 160.49415
case 4      1.70621   82.77880 140.64139
case 5      0.76169  120.48102 184.03695
case 6      3.71265  120.48007 183.29809
	【中略】
case 147  578.01838   23.36413   4.00770
case 148  618.18551   16.74095   1.11181
case 149  720.69230   33.20289   3.94189
case 150  550.78799   10.11250   2.69108

$p.values	# 各群への所属の確率

              setosa versicolor  virginica
case 1    9.7826e-01 6.8699e-24 1.7457e-38
case 2    7.2085e-01 3.4538e-17 2.8628e-32
case 3    8.6403e-01 1.1849e-19 1.1454e-33
case 4    7.8959e-01 4.4882e-17 2.0574e-29
case 5    9.4351e-01 4.2162e-25 1.0126e-38
case 6    4.4629e-01 4.2181e-25 1.4594e-38
	【中略】
case 147 8.8576e-124 1.0709e-04 4.0497e-01
case 148 1.7956e-132 2.1703e-03 8.9239e-01
case 149 1.1523e-154 1.0855e-06 4.1393e-01
case 150 6.9093e-118 3.8575e-02 6.1078e-01

$prediction	# 判別結果

  [1] setosa     setosa     setosa     setosa     setosa    
  [6] setosa     setosa     setosa     setosa     setosa    
 [11] setosa     setosa     setosa     setosa     setosa    
 [16] setosa     setosa     setosa     setosa     setosa    
 [21] setosa     setosa     setosa     setosa     setosa    
 [26] setosa     setosa     setosa     setosa     setosa    
 [31] setosa     setosa     setosa     setosa     setosa    
 [36] setosa     setosa     setosa     setosa     setosa    
 [41] setosa     setosa     setosa     setosa     setosa    
 [46] setosa     setosa     setosa     setosa     setosa    
 [51] versicolor versicolor versicolor versicolor versicolor
 [56] versicolor versicolor versicolor versicolor versicolor
 [61] versicolor versicolor versicolor versicolor versicolor
 [66] versicolor versicolor versicolor versicolor versicolor
 [71] virginica  versicolor virginica  versicolor versicolor
 [76] versicolor versicolor versicolor versicolor versicolor
 [81] versicolor versicolor versicolor virginica  versicolor
 [86] versicolor versicolor versicolor versicolor versicolor
 [91] versicolor versicolor versicolor versicolor versicolor
 [96] versicolor versicolor versicolor versicolor versicolor
[101] virginica  virginica  virginica  virginica  virginica 
[106] virginica  virginica  virginica  virginica  virginica 
[111] virginica  virginica  virginica  virginica  virginica 
[116] virginica  virginica  virginica  virginica  virginica 
[121] virginica  virginica  virginica  virginica  virginica 
[126] virginica  virginica  virginica  virginica  virginica 
[131] virginica  virginica  virginica  virginica  virginica 
[136] virginica  virginica  virginica  virginica  virginica 
[141] virginica  virginica  virginica  virginica  virginica 
[146] virginica  virginica  virginica  virginica  virginica 
Levels: setosa versicolor virginica

$correct	# 判別の正誤

  [1]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
 [11]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
 [21]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
 [31]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
 [41]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
 [51]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
 [61]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
 [71] FALSE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
 [81]  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
 [91]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
[101]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
[111]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
[121]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
[131]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
[141]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE

$correct.table	# 判別結果表

            prediction
group        setosa versicolor virginica
  setosa     50      0          0       
  versicolor  0     47          3       
  virginica   0      0         50       

$correct.ratio	# 正判別率

[1] 98

$group

  [1] setosa     setosa     setosa     setosa     setosa    
  [6] setosa     setosa     setosa     setosa     setosa    
 [11] setosa     setosa     setosa     setosa     setosa    
 [16] setosa     setosa     setosa     setosa     setosa    
 [21] setosa     setosa     setosa     setosa     setosa    
 [26] setosa     setosa     setosa     setosa     setosa    
 [31] setosa     setosa     setosa     setosa     setosa    
 [36] setosa     setosa     setosa     setosa     setosa    
 [41] setosa     setosa     setosa     setosa     setosa    
 [46] setosa     setosa     setosa     setosa     setosa    
 [51] versicolor versicolor versicolor versicolor versicolor
 [56] versicolor versicolor versicolor versicolor versicolor
 [61] versicolor versicolor versicolor versicolor versicolor
 [66] versicolor versicolor versicolor versicolor versicolor
 [71] versicolor versicolor versicolor versicolor versicolor
 [76] versicolor versicolor versicolor versicolor versicolor
 [81] versicolor versicolor versicolor versicolor versicolor
 [86] versicolor versicolor versicolor versicolor versicolor
 [91] versicolor versicolor versicolor versicolor versicolor
 [96] versicolor versicolor versicolor versicolor versicolor
[101] virginica  virginica  virginica  virginica  virginica 
[106] virginica  virginica  virginica  virginica  virginica 
[111] virginica  virginica  virginica  virginica  virginica 
[116] virginica  virginica  virginica  virginica  virginica 
[121] virginica  virginica  virginica  virginica  virginica 
[126] virginica  virginica  virginica  virginica  virginica 
[131] virginica  virginica  virginica  virginica  virginica 
[136] virginica  virginica  virginica  virginica  virginica 
[141] virginica  virginica  virginica  virginica  virginica 
[146] virginica  virginica  virginica  virginica  virginica 
Levels: setosa versicolor virginica

$ngroup
[1] 3


二群のときには,二群の重心への距離の差を求め,グラフを描く > (result <- quad.disc(iris[51:150, 1:4], iris[51:150, 5])) 判別結果 prediction group versicolor virginica versicolor 47 3 virginica 0 50 正判別率 = 97.0 % > plot(result) fig > plot(result, which="barplot") fig > plot(result, which="s", xpos="right") fig ・ 解説ページ


・ 直前のページへ戻る  ・ E-mail to Shigenobu AOKI

Made with Macintosh