No.14682　svm()の出力から手計算で予測値を得るためには？　　【波音】　2011/05/25(Wed) 08:19

Rのe1071パッケージに含まれるsvm()を使ってサポートベクタマシンを実行しています。

library(e1071)
data(iris)
result <- svm(Species ~ ., data=iris)

例えばこれ（http://www.orsj.or.jp/~archive/pdf/bul/Vol.46_05_225.pdf）における式(24)でいうαyの部分に該当するのが

> result$coefs[1:5, ]
[,1] [,2]
[1,] 0.0890967 0.0000000
[2,] 0.0000000 0.1454778
[3,] 0.8651998 0.9486997
[4,] 0.0000000 0.1315259
[5,] 0.0000000 0.2761224

として得られます。カーネルのサポートベクタ部分は

> result$SV[1:5,]
Sepal.Length Sepal.Width Petal.Length Petal.Width
9 -1.7430170 -0.3609670 -1.335752 -1.311052
14 -1.8637803 -0.1315388 -1.505695 -1.442245
16 -0.1730941 3.0804554 -1.279104 -1.048667
21 -0.5353840 0.7861738 -1.165809 -1.311052
23 -1.5014904 1.2450302 -1.562342 -1.311052

として得られます。

応答変数の水準が2つの場合は1つの（1列の）判別予測値（result$decision.valuesで得られる）が得られますが，3つ以上の場合はαyが応答変数の水準数-1だけでてきます。

こうした応答変数の水準数が3つ以上の場合に，どのようにしたらresult$decision.valuesで得られる判別予測値を計算できるでしょうか？

練習のために手作業でresult$decision.valuesを導きたいためです。
No.14702　Re: svm()の出力から手計算で予測値を得るためには？　　【波音】　2011/05/28(Sat) 11:01

自己解決しましたが，情報共有のために手計算のコードを載せておきます（計算過程を把握するために，あえて汎用性がないようにしてあります）。
setwd("c://users/ユーザー名/documents/R/workdir")

data(iris)              # アヤメのデータ
x <- iris[, 1:4]	        # 説明変数部分のデータを抽出
library(e1071)　　　　　# パッケージの読み込み

# svmを使って判別分析
result <- svm(Species ~ ., data=iris, probability=TRUE, kernel="polynomial", cross=1)

sv.num <- c(1, cumsum(result$nSV))		# サポートベクタの区切り番号：[1]  1  6 32 54
sv1 <- result$SV[1:sv.num[2], ]			# 1行〜6行までのサポートベクタ
sv2 <- result$SV[(sv.num[2]+1):sv.num[3], ]	# 7行〜32行までのサポートベクタ
sv3 <- result$SV[(sv.num[3]+1):sv.num[4], ]	# 33行〜54行までのサポートベクタ


	# combn(3, 2)	# この組み合わせ数だけiとjを変更する

	i <- 1		# ここは手入力
	j <- 2		# ここも手入力

	start <- c(1, cumsum(result$nSV)+1)
	start <- start[-length(start)]

	ri <- start[i] : (start[i] + result$nSV[i] - 1)
	rj <- start[j] : (start[j] + result$nSV[j] - 1)

	# coefs for (i,j):
	coef1 <- result$coefs[ri, j-1]
	coef2 <- result$coefs[rj, i]

	sv1.size <- nrow(sv1)
	sv2.size <- nrow(sv2)
	sv3.size <- nrow(sv3)
	all.size <- c(sv1.size, sv2.size, sv3.size)
	which(all.size == length(coef1))	            # 返される番号が該当するcoef1に対応するsvの番号
	which(all.size == length(coef2))	            # coef2に対応するsvの番号

g <- result$gamma
b <- result$rho
dgr <- result$degree

scale.x <- scale(x)

#### 1個目の計算 ####
zero <- numeric(length(scale.x[,1])*length(sv1[,1]))                            # sv.**の部分を変更する
kanel <- matrix(zero, nrow=length(scale.x[,1]), ncol=length(sv1[,1]), byrow=F)  # 同じく
omomi <- matrix(zero, nrow=length(scale.x[,1]), ncol=length(sv1[,1]), byrow=F)  # 同じく
# iは元データの説明変数部分の行数（オブジェクトscale.xの行数）
# sはサポートベクタの数（オブジェクトsvの行数）
for(i in 1:length(scale.x[,1])){
for(s in 1:length(sv1[,1])){					# sv.**を変更
kanel[i, s] <- (g*t(scale.x[i, ])%*%(sv1[s, ]))^dgr		# 同じく
omomi[i, s] <- coef1[s]*kanel[i, s]
}
}
my.sum <-apply(omomi, 1, sum)		# 係数で重みづけされたカーネル部分の行和を計算

#### 2個目の計算 ####
zero2 <- numeric(length(scale.x[,1])*length(sv2[,1]))                           # sv.**の部分を変更する
kanel2 <- matrix(zero2, nrow=length(scale.x[,1]), ncol=length(sv2[,1]), byrow=F)	# 同じく
omomi2 <- matrix(zero2, nrow=length(scale.x[,1]), ncol=length(sv2[,1]), byrow=F)	# 同じく
# iは元データの説明変数部分の行数（オブジェクトscale.xの行数）
# sはサポートベクタの数（オブジェクトsvの行数）
for(i in 1:length(scale.x[,1])){
for(s in 1:length(sv2[,1])){                            # sv.**を変更
kanel2[i, s] <- (g*t(scale.x[i, ])%*%(sv2[s, ]))^dgr    # 同じく
omomi2[i, s] <- coef2[s]*kanel2[i, s]
}
}
my.sum2 <-apply(omomi2, 1, sum)		# 係数で重みづけされたカーネル部分の行和を計算

ans <- (my.sum + my.sum2) - b[1]	# b[]の添え字を変えること
No.14703　Re: svm()の出力から手計算で予測値を得るためには？　　【さくらとお城】　2011/05/28(Sat) 15:38

便乗質問です。
　
　水準2つの場合もよく理解できていません。例えば，result$coefsとirisの原データなどから，どのようにresult$decision.valuesや分類結果を計算されますか？

　どこかに分かりやすい説明はありませんか？どうぞよろしくお願いします。
No.14704　Re: svm()の出力から手計算で予測値を得るためには？　　【波音】　2011/05/28(Sat) 16:35

kernlabパッケージのkesvm()を使った例であれば

　線形svmの識別関数を求める-ksvmパッケージ

というページがありますよ。ここではsetosaに該当する50行を削除して，応答変数（Species）が2水準である場合を例示しています。

応答変数が2水準の場合は1つの分割予測式(combn(3, 2))が得られるだけなので，カーネル部分と係数の掛け算の部分も単純です。

data(iris)
iris2 <- iris[iris$Species != "setosa",]
library(e1071)
result <- svm(Species ~ ., data=iris2, kernel="polynomial")

x <- iris2[, 1:4]	# 説明変数部分のデータ
x <- scale(x)		# 標準化する

# 以下結果オブジェクトから必要なものを取得
cf <- result$coefs
sv <- result$SV
g <- result$gamma
b <- result$rho
dgr <- result$degree

comb.num <- ncol(combn(result$nclasses, 2))
ans <- matrix(numeric(nrow(x)*comb.num), nrow=nrow(x), ncol=comb.num)

x.rows <- nrow(x)	# 元データの行数
sv.rows <- nrow(sv)	# サポートベクタの数（result$SVの行数）

kanel <- matrix(numeric(x.rows * sv.rows), nrow=x.rows, ncol=sv.rows, byrow=FALSE)
omomi <- matrix(numeric(x.rows * sv.rows), nrow=x.rows, ncol=sv.rows, byrow=FALSE)

# カーネル関数（多項式カーネル）
my.kernel <- function(x, sv){
	res <- (g * x %*% sv)^dgr
	res
}

# kは元データの説明変数部分の行数（オブジェクトxの行数）
# sはサポートベクタの数（オブジェクトsvの行数）
for(k in 1:x.rows){
	for(s in 1:sv.rows){
		kanel[k, s] <- my.kernel(x[k, ], sv[s, ])	# カーネル部分の計算
		omomi[k, s] <- cf[s] * kanel[k, s]			# カーネル部分を係数で重み付け
	}
}

my.sum <-apply(omomi, 1, sum)	# 係数で重みづけされたカーネル部分の行和を計算
								# 　つまりxの行数分だけの要素を持つ
my.sum - b[1]					# 各サンプルの判別予測値からbを引く
								# これがresult$decision.valuesと一致する

そして私も分かりやすい解説ページを散々見つけまわりましたがなかなか見つけだすことができませんでした(^_^;)
No.14705　Re: svm()の出力から手計算で予測値を得るためには？　　【さくらとお城】　2011/05/28(Sat) 17:21

　波音さん，ありがとうございます。

　どこまで理解できるか，ちょっと不安ですが，このコードをじっくり解析してみます。
　もし，その原理について分りやすい解説が見つかりましたら，教えてください。ただ，あまり高度な数式ですと，余計わからなくなります。
No.14728　Re: svm()の出力から手計算で予測値を得るためには？　　【さくらとお城】　2011/06/06(Mon) 09:34

波音さん

　お陰さまで，”線形svmの識別関数を求める-ksvmパッケージ”の内容をほぼ理解できました。そこで，質問です。

＊No. 14704の内容も，”pred3 <- iris3 %*% coeff - b(model) # 検算の式”のようにできますか？というのは，このように表現できれば，変数選択に使えるかもと考えました。
　もしかしてカーネル関数が線形の場合だけ，上の単純な式になるかな。

＊No. 14702の内容は，まだ理解できていません。2水準の場合，判別直線が1本あればできるとのイメージですが，3水準の場合，どのようなイメージですか？

　どうぞよろしくお願いします。

● 「統計学関連なんでもあり」の過去ログ--- 044 の目次へジャンプ
● 「統計学関連なんでもあり」の目次へジャンプ
● 直前のページへ戻る