目的

使用法

all.quad.disc(data, group)
print.all.quad.disc(obj)

引数

data        データ行列（行がケース，列が変数）
group       各ケースがどの群であるかを表す変数
            データフレーム中の変数の場合には，iris[, 5] ではなく iris[5] のように
           （1 列だけのデータフレームになるように）指定すること（変数名を関数に渡せるように）

ソース

以下のほかに，quad.disc が必要。

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

# 総当たり法による二次の判別分析
all.quad.disc <- function(   data,                                           # 説明変数のデータフレーム（データ行列）
                        group)                                          # 群を表す変数（データフレームから指定するときには iris[,5] ではなく，iris[5] のように）
{
        BinConv <- function(nv)
        {
               n <- 2^nv                                             # 独立変数を取り出す取り出し方
                bincomb <- matrix(FALSE, nrow=n, ncol=nv)            # e1071 パッケージの bincombinations より
                for (j in 1:nv) {
                        bincomb[, j] <- rep(c(rep(FALSE, n/2^j), rep(TRUE, n/2^j)), length = n)
                }
                bincomb <- bincomb[-1,]
                return(bincomb)
        }
        nv <- ncol(data)
        vname <- colnames(data)                                              # 変数名（なければ作る）
        if (is.null(vname)) {
                vname <- colnames(data) <- paste("x", 1:nv, sep="")
        }
        gname <- names(group)
        if (is.null(gname)) {                                           # group を，データフレームから iris[, 5] のようにすると，名前がなくなる
                gname <- ""                                          # iris[5] のように 1 列のデータフレームとして指定すること
        }
        group <- factor(as.matrix(group))
        ok <- complete.cases(data, group)                            # 欠損値を含まないケースだけを対象にする
        data <- data[ok,]
        group <- group[ok]
        n <- nrow(data)
        bincomb <- BinConv(nv)
        nr <- nrow(bincomb)
        correct <- numeric(nr)                                               # 正準判別関数の正判別率を基準とする
        for (i in 1:nr) {
                dat <- data[, bincomb[i,], drop=FALSE]
                a <- quad.disc(dat, group)                                   # 別途用意してある quad.disc.R で定義
                correct[i] <- sum(a$correct)/n               # 正判別率を記録しておく
        }
        ans <- data.frame(correct, bincomb)
        colnames(ans) <- c("correct rate", vname)
        return(structure(list(ans=ans, name=vname, gname=gname), class="all.quad.disc"))
}
# print メソッド
print.all.quad.disc <- function(obj)
{
        ans <- obj$ans
        name <- obj$name
        gname <- obj$gname
        o <- order(ans[, 1], decreasing=TRUE)
        ans <- ans[o,]
        nc <- ncol(ans)
        cat("\ncorrect rate   formula\n")
        for (i in 1:nrow(ans)) {
                cat(sprintf("%10.5f     %s ~ %s\n",ans[i, 1], gname, paste(name[as.matrix(ans[i, 2:nc])], collapse=" + ")))
        }
        invisible(ans)                                                  # 結果をソートしただけのものを返す
}


使用例

> all.quad.disc(as.matrix(iris[1:4]), iris[5])

correct rate   formula
   0.98667     Species ~ Sepal.Length + Petal.Length + Petal.Width
   0.98000     Species ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width
   0.97333     Species ~ Petal.Length + Petal.Width
   0.96667     Species ~ Sepal.Width + Petal.Length + Petal.Width
   0.96667     Species ~ Sepal.Length + Petal.Width
   0.96000     Species ~ Petal.Width
   0.96000     Species ~ Sepal.Width + Petal.Length
   0.96000     Species ~ Sepal.Length + Petal.Length
   0.96000     Species ~ Sepal.Length + Sepal.Width + Petal.Length
   0.95333     Species ~ Petal.Length
   0.95333     Species ~ Sepal.Length + Sepal.Width + Petal.Width
   0.94667     Species ~ Sepal.Width + Petal.Width
   0.80667     Species ~ Sepal.Length + Sepal.Width
   0.72000     Species ~ Sepal.Length
   0.55333     Species ~ Sepal.Width