2018.11.13
数理統計学演習

判別分析

東京大学大学院農学生命科学研究科　大森宏

http://lbm.ab.a.u-tokyo.ac.jp/~omori/kensyu/discriminant.htm

１．判別分析とは

　データを分類する手法である．分類のわかっているデータ（トレイニングデータ，教師データ）があり，データの属性値に基づいて，この分類を判別する手法である．機械学習の分野では，教師あり学習（supervised learning）と呼ばれている．データ分類手法として知られるクラスター分析は，分類がわかっていないので，教師無し学習（unsupervised learning）と呼ばれる．
　応用分野は多岐にわたり，病気診断，スニップデータからの発現分類，スパムメールフィルター，などがあげられる．

２．２群の判別

２−１．空間の分割と判別規則

　母集団が２つの部分母集団 A，B の混合からなっている場合を考える． 部分母集団 A，B の確率密度をそれぞれ fA(x )，fB(x ) とし，その混合比率をπA，πB とすると，母集団全体の確率密度は，
f　(x ) ＝ πAfA(x ) ＋ πBfB(x )                (1)
となる．いま，帰属未知の対象 X の属性値 x に基づき部分母集団への帰属を決めるとする．このとき， x の取り得る値の全体（標本空間） R　は，RA と RB に分割され，
x ∈ RA なら X は A に帰属，
x ∈ RB なら X は B に帰属． となる．これを判別規則という．

２−２．最適な判別規則

　x の空間 R の分割が判別規則に対応するが，どのような空間分割を行うのがよいであろうか． まず，ベイズの定理による事後確率の最大化という考え方がある．
ベイズ（Bayes）の定理
P(A) を事象 A が起こる確率，P(A|B) を事象 B が生起したもとでの事象 A が起こる条件付き確率， P(A∩B) を事象 A と事象 B が同時に生起する確率とすると，
P(A∩B) ＝ P(A|B)P(B) ＝ P(B|A)P(A)
が成り立つので
P(A|B) ＝ P(B|A)P(A)/P(B)               
である．これをベイズの定理と言う． 　いま，標本空間 Ω が互いに排反な事象，
A1，…，An （P(Ai∩Aj，i≠j， P(A1∪…∪An) ＝ P(A1)＋…＋P(An)＝1），
に分割されたとすると，Ω の任意の部分集合 B の生起する確率 P(B) は，
P(B) ＝ P(B|A1)P(A1)＋…＋P(B|An)P(An) ＝ ΣiP(B|Ai)P(Ai)
と表現される．これより，事象 B が生起したときの事象 Ai が生起する事後確率は，
              
と表せる．
ベイズ判別規則
　対象 X が A に帰属したとき，対象 X の値が x のごく近傍を取る確率は， P(x |A) ＝ fA(x )Δx であり，対象 X が B に帰属したとき， x のごく近傍を取る確率 は，P(x |B) ＝ fB(x )Δx とあらわせる． 　対象 X の値が x の近傍であるとき，X が部分母集団 A 由来である事後確率は， P(A) ＝ πA に注意すると，
であり，X が部分母集団 B 由来である事後確率は，
となる．事後確率の高い部分集団に対象 X を帰属させるとすると，

対象 X は A 由来と判定 ⇔ P(A|x ) ＞ P(B|x ) ⇔ πAfA(x ) ＞ πBfB(x ) ⇔ logπA ＋ logfA(x ) ＞ logπB ＋ logfB(x )，

すなわち，
RA ＝ { x   |   logπA ＋ logfA(x ) ＞ logπB ＋ logfB(x ) }               
となる．
誤判別確率最小判別
　判別規則により，対象 X を判別したとする．このときは次のような誤判別がある． 　すると，誤判別の確率は，
となる．これより，RA として，
の部分のみを含み，
の部分をふくまなければ誤判別確率は最小となる．よって，式(2)と同じ判別規則が得られる．
１次元分布での判別
　母集団が１次元の分布の場合の判別による空間（数直線）分割の例を以下に示す．
missclass1 missclass2
図１．空間の分割と誤判別確率

図２．分散のみが異なるときの空間の分割
missclass3

２−３．正規母集団のもとでの判別

２次判別
　部分母集団が正規分布で記述されているときは，簡単な判別規則をつくることができる． 平均 μ，分散共分散行列 Σ の p 次元正規分布の密度関数は，
              
と表せ，その対数は，
とかける．部分母集団 A，B の平均を μA，μB，分散共分散行列を ΣA，ΣB とすると，対象 X は A 由来と判定する領域は，
              
と表せる．これを２次判別（Quadratic Discriminant Analysis）という．
線形判別
　ところで，部分母集団 A，B が共通の共分散行列 Σ をもち，その存在比率が等しい （πA ＝ πB ＝ 1/2）とすると，A 由来と判定する領域は，
と簡略化される．
マハラノビス距離
　ところで DA，
              
を，部分母集団 A と対象の値 x とのマハラノビス（Mahalanobis）距離という．この距離は，通常 のユークリッド（Euclidean）距離と異なり，分布を考慮にいれた距離で，マハラノビス距離が小さい方の部分 母集団の方が x の近傍のような値が得られる確率が高い（起こりやすい）ことを意味 している．つまり，マハラノビス距離の近い部分母集団に対象を割り付ける，という判別規則となる．
線形判別関数の導出
　部分母集団 A，B と対象の値 x とのマハラノビス距離， DA，DB，を用いると，
           
という x の線形関数ができる．これを線形判別関数（linear discriminant function）という． 判別規則は，

w ＞ 0 ⇔ 対象 X は A 由来と判定，w ＜ 0 ⇔ 対象 X は B 由来と判定，

となる．つまり，超平面 w ＝ 0 により，標本空間が２つに分割される．
線形判別関数の例
　２変量正規母集団を考える．部分母集団 A，B の平均 μA，μB と共通の 分散共分散行列 Σ が，
であったとする．このとき，部分母集団 A，B が定義されている平面を分割する線形判別関数は，
と計算できる．
図３．２つの２次元正規母集団に対する線形判別関数
母集団比率を考えた線形判別関数
　母集団構成比率 πA，πB が既知であるときや，未知でもこの 比率に大きな偏りが認められる時は構成比率の値を考慮に入れた判別関数をつくるべきである． 式(4)の判別領域で，log|ΣA|，log|ΣB| の項を消去すると，対象 X を A 由来 と判別する領域は，
           
となる．
母集団母数が未知の場合
　部分母集団 A，B の母数が未知であるのが普通である．このときは部分母集団への帰属がわかっている トレイニングデータ（training data），教師データともいう，から母集団母数を推定する． 部分母集団 A，B それぞれから大きさ nA，nB のトレイニングデータが得られた とすると，それぞれの標本平均，
で母集団平均 μA，μB が推定される．
　部分母集団 A，B の分散共分散はそれぞれの群内標本分散共分散行列，
で推定され，２次判別ではこれを用いる。
　両者が等しいと考えられるときは，共通の分散共分散行列が重みづけ平均，
で推定される．これを用いて，分類未知の対象 x に対する線形判別関数が，
と推定できる．また，部分母集団の構成比率は，トレイニングデータが完全にランダムに抽出された とすれば
で推定できる．

2-4．ロジスティック判別

ロジスティック回帰
　ある事象の生起確率 p が，変数により制御されている場合を考える．よく知られた例では，薬品の投与量 x と死亡率 p の関係が考えられる．死亡オッズの対数（p のロジット関数）が投与量に比例すると考え，
としたモデルがロジスティック回帰である．この式を変形すると死亡率 p は，
投与量 x のシグモイド関数で表現できる．
　一般には，m 個の説明変数を考えると，
となる．
判別分析への適用
　いま，２つの群 A，B を記述する変数が m 個あり，B が生起する確率を p，A が生起する確率を 1 - p とすると，ロジスティック回帰モデルで生起確率が記述できる．ここで，p > 0.5 なら B，p < 0.5 なら A と判別 すると決めれば判別分析が行える．判別の境界線は p = 0.5 となる直線である．R では，glm() 関数でロジスティック判別が行える．
　ロジスティック判別では，各群の確率分布を記述する必要がないので，正規母集団モデルに適合しないようなデータにも適用できる．

2-5．２群の判別例

　R の MASS ライブラリィーにある Cushings データを用いて解説する。これは，コルチゾールの過剰分泌による高血圧疾患であるクッシング症候群（Cushing's syndrome）のデータである。クッシング症候群には，腺腫，両側過形成，腫瘍の３つの型があり，組織病理学的診断結果がそれぞれ，a，b，c，で表示されている。なお，u は未知のタイプである。この型診断をステロイド代謝産物であるテトラヒドロコルチゾン（Tetrahydrocortisone）とプレグナントリオール（Pregnanetriol）の腎排泄率（mg/24h）のデータで判別することを試みる。人数の多い腺腫（a）と両側過形成（b）の判別の様子を，線形判別，２次（quadratic）判別，ロジスティック判別で比較してみる。 
library(MASS)			# MASS ライブラリィーの呼び出し
Cushings			# Cushings データ
cush2 <- log(Cushings[1:16,1:2])
tp2 <- factor(Cushings$Type[1:16])		# カテゴリー変数化
cush2.lda <- lda(cush2, tp2)			# 線形判別
cush2.qda <- qda(cush2, tp2)			# ２次判別
Tetrahydrocortisone <- cush2[,1]
Pregnanetriol <- cush2[,2]
# ロジスティック回帰
cush2.glm <- glm(tp2 ~ Tetrahydrocortisone + Pregnanetriol, family=binomial)
#
# 散布図上の判別線表示
len <- 100
plot(cush2, type="n", xlim=c(0,3), ylim=c(-4,3),
	xlab="log(Tetrahydrocortisone)", ylab="log(Pregnanetriol)")
text(cush2, as.character(tp2), col=as.numeric(tp2)+1)
title(main="２群の判別")
xp <- seq(0, 3, length = len) 
yp <- seq(-4, 3, length = len)
# 100 × 100 のグリッドの定義
cushT <- expand.grid(Tetrahydrocortisone=xp, Pregnanetriol=yp)
Z <- predict(cush2.lda, cushT)		# 線形判別でのグリッド上の事後確率
zp <- Z$posterior[,1] - Z$posterior[,2]		# a の事後確率 − b の事後確率
contour(xp, yp, matrix(zp, len), add=T, levels=0)  # 等事後確率線
#
Z <- predict(cush2.qda, cushT)		# ２次判別でのグリッド上の事後確率
zp <- Z$posterior[,1] - Z$posterior[,2]
contour(xp, yp, matrix(zp, len), add=T, levels=0, col="blue")
#
Z <- predict(cush2.glm, cushT, type="response")			# a の確率
contour(xp, yp, matrix(Z, len), add=T, levels=0.5, col="purple")	# a が 0.5 になる線
legend("bottomright", legend=c("線形判別","２次判別","ロジスティック判別"), 
	lty=1, col=c("black","blue","purple"), cex=0.8)
discriexample

３．多群の判別

3-1．一般規則

母集団が g 個の部分母集団にわかれていて，各部分母集団の構成比率が πi で，その 密度関数が fi(x ) であるとすると，母集団全体の確率密度は，

f (x ) ＝ π1f1(x ) ＋…＋ πgfg(x )               

と表せる．２群の判別のときと同様に考えると，帰属未知の対象 X の値 x が得られたとき

対象 X は群 i 由来と判定 ⇔ πifi(x ) が最大

という判別規則となり，x の取り得る値の全体空間が g 個に分割される．

3-2．多項ロジスティック判別

　ロジスティック判別は２群の判別であるが，g 群の判別では g(g - 1)/2 のロジステッィク判別を行えば判別できる．R では，nnet ライブラリィーの multinom() 関数で多項ロジスティック判別が行える．

3-3．決定木，分類木

　変数ごとに領域を分割して２分木に枝分かれさせて分類する手法．若干荒っぽい手法であるが，解釈が容易なので使いやすい．介護保険の要介護度判定に用いられている． R では，rpart ライブラリィーの rpart() 関数で行える．

3-4．ナイーブベイズ

　3-1 の一般規則は，πifi(x ) が最大である群 i に分類未知の対象 x を割り付けるものであった．このためには，群 i での x の分布 fi(x ) が推定できる必要がある．x が多変量になるにつれて，共分散を正しく推定することが困難になってくる．そこで，変量間の相関構造を無視し，互いに独立な変量と仮定して fi(x ) を近似的に推定するのがナイーブベイズである．このとき，
naive
と仮定する．R では，e1071 ライブラリィーの naiveBayes() 関数でナイーブベイズ判別が行える．
　ナイーブベイズは，文章の分類など変数の数が非常に多い（判別に重要な単語の種類数）ときに高速化のため用いられる．

3-5．サポートベクトルマシン（Support Vector Machine SVM）

　データを判別に重要なデータ（サポートベクトル）とそうでないデータに分け，サポートベクトルで判別規則を構築する手法である．また，データを表現する関数を内積で定義するので，非常に柔軟な判別規則を構成（カーネルトリック）することができる．
分離超平面とマージン最大化
　いま，p 次元空間において，(a・b) をベクトル a，b 間の内積としたとき， 任意の p 次元変数 x に対して分離超平面を

g(x) = (w・x) ＋ b ＝ 0

とする．ここで，関数 sgn(u) を

sgn(u) = 1，if u ≧ 0，sign(u) = -1，if u ＜ 0

となる符号関数として，判別関数を

f(x) ＝ sgn((w・x) ＋ b)

と定義する．すなわち，f(x) ＝ 1，のときは x を群１に割り付け， f(x) ＝ -1， のときは群２に割り付けるものとする．
　ここで，m 個のトレイニングデータ xi，i ＝ 1，…，m， すべてがこの超平面で分離可能であったとする．すなわち，i が群１由来のときは， yi ＝ f(xi) ＝ 1，であり， 群２由来のときは，yi ＝ f(xi) ＝ -1，であったとする． このとき，群間のマージンを最大にする判別関数求める．これは，トレイニングデータを群に分ける超平面 のうち，なるべく間があいた余裕のある超平面を求めようとするものである． この最適な超平面は下の図に示すように，
　
の最適化問題を解くことによって得られる．
　非負のラグランジュ（Lagrange）乗数 α＝ (α1，…，αm) ≧ 0 を用いて 目的関数を
　
として，双対問題
　
を解くことに帰着する．すなわち，この式を w と b で最小化し， ラグランジュ乗数 αi で最大化する．
　この式を w と b で偏微分して 0 とおく，すなわち，
とすると，その解は，
となるので，これを目的関数に代入すると，これは以下の最大化問題
　
を満たすラグランジュ乗数 αi を 求めることに帰着する．
　この解のうち，αi ＞ 0 となるトレイニングデータは，図のマージン， (w・x) ＋ b ＝ 1 か (w・x) ＋ b ＝ -1 上に乗っている．このようなデータを サポートベクター（Support Vector）と呼んでいる．残りの多くの αi ＝ 0 となるトレイニング データは式(17)の分離超平面や式(18)の判別関数の構成に関与しない．
カーネルトリック
　分離超平面の識別関数，

g(x) = (w・x) ＋ b = x'w ＋ b

において，非線形関数 φ(x) に対し， カーネル関数 K(xi, xj) を用いて，

g(x) =φ(x)'w ＋ b ＝ ΣαiyiK(x, xi) + b

と書ける．こうすると，カーネル関数を定義すれば，非線形関数 φ(x) を具体的に定義する必要はない．これをカーネルトリックという．
　R では，e1071 ライブラリィーの svm() 関数で SVM が行える．

3-6．ランダムフォレスト

　分類木が与える分類は，直線や平面の組み合わせで行うのでいかにも荒っぽい．そこで，多くの分類木を生成させ，それらを組み合わせて分類規則を生成する方法がある．これを，バギング（bagging : bootstrap aggregating）と言う．
　データから重複を許したブートストラップサンプルデータを生成し，それぞれで分類木を生成する．最終分類は，それぞれの分類木の分類結果の多数決で決める．ランダムフォレストでは，分類に使用する変数の組み合わせもランダムに決めている．このような弱学習を統合したバギングは予測誤差を減少させることが知られている．
　R では，randomForest ライブラリィーの randomForest() 関数でランダムフォレスト判別が行える．

3-7．３群の判別例

　先ほど使用したデータで，クッシング症候群の３つの型，腺腫（a），両側過形成（b），腫瘍（c）の判別を色々な手法で行って見る。 
cush <- log(Cushings[1:21,1:2])
tp <- factor(Cushings$Type[1:21])
len <- 100
xp <- seq(0, 4.5, length = len) 
yp <- seq(-4, 3, length = len)
cushT <- expand.grid(Tetrahydrocortisone=xp, Pregnanetriol=yp)
#
# 判別線表示関数定義
predplot <- function(pred, main=""){
  plot(cush, type="n", xlim=c(0,4.5), ylim=c(-4,3),
	xlab="log(Tetrahydrocortisone)", ylab="log(Pregnanetriol)", main=main)
  text(cush, as.character(tp), col=as.numeric(tp)+1)
  zp <- pred[,3] - pmax(pred[,2], pred[,1]) 
  contour(xp, yp, matrix(zp, len), add=T, levels=0) 
  zp <- pred[,1] - pmax(pred[,2], pred[,3]) 
  contour(xp, yp, matrix(zp, len), add=T, levels=0)
}
#
#
# 線形判別
cush.lda <- lda(cush, tp)			# 線形判別関数
Z <- predict(cush.lda, cushT)
main <- "線形判別"
predplot(Z$posterior, main=main)
#
# ２次判別
cush.qda <- qda(cush, tp)
Z <- predict(cush.qda, cushT)
main <- "２次判別"
predplot(Z$posterior, main=main)
#
# 多項ロジスティック判別
Tetrahydrocortisone <- cush[,1]
Pregnanetriol <- cush[,2]
library(nnet) 
cush.multinom <- multinom(tp ~ Tetrahydrocortisone + Pregnanetriol, maxit=250)
zz <- predict(cush.multinom, cushT, type="probs")
main <- "多項ロジスティック判別"
predplot(zz, main=main)
#
# 決定木，分類木
library(rpart)
cush.tree <- rpart(tp ~ Tetrahydrocortisone + Pregnanetriol)
plot(cush.tree); text(cush.tree)
zz <- predict(cush.tree, cushT, type="prob")
main <- "分類木"
predplot(zz, main=main)
#
#
library(e1071)
# ナイーブベイズ（naive bayes）
cush.bayes <- naiveBayes(cush, tp)
zz <- predict(cush.bayes, cushT, type="raw")
main <- "ナイーブベイズ"
predplot(zz, main=main)
#
# サポートベクトルマシン（SVM）
cush.svm <- svm(tp ~ Tetrahydrocortisone + Pregnanetriol, probability = TRUE)
zz <- predict(cush.svm, cushT, decision.values = TRUE, probability = TRUE)
zz.dec <- attr(zz, "decision.values")
zz.prob <- attr(zz, "probabilities")
main <- "サポートベクトルマシン"
predplot(zz.prob, main=main)
#
# ランダムフォレスト
library(randomForest)
set.seed(71)
cush.rf <- randomForest(cush, tp)
zz <- predict(cush.rf, cushT, type="prob")
main <- "ランダムフォレスト"
predplot(zz, main=main)
linear quadro
logit tree
naive svm
rf

４．正準判別

　母集団分布が多変量のとき，判別のようすを視覚化することは難しい．そこで，主成分分析と似たような考え方 で次元を縮約して，判別を比較的少数の次元で行う手法がある．これが正準判別である．
　いま，p 次元空間に g 個の部分母集団があり，i 番目の部分母集団から大きさ ni の トレイニングデータ，
が得られたとする．トレイニングデータ全体での偏差積和行列 T は，データの各群への割付により， 群間偏差積和行列 B と群内偏差積和行列 W に，

T ＝ B ＋ W                (13)

のように分解される．なお，
である．ここで，分類未知の対象 X の値 x が与えられたとする． このとき，x の線形結合，
              
を考える．
　x の分散が Var[x] ＝ Σ のとき，u の分散は，
Var[u } ＝ Var[a'x] ＝ a'Σa
となることから，トレイニングデータ全体での u の総平方和は a'Ta と あらわされ，これが

a'Ta ＝ a'Ba ＋ aWa               

のように分解される．a'Ba は u の群間平方和で， a'Wa は u の群内平方和である．
　判別に有効と考えられる線形結合 u は，群の差をなるべく大きくするのがよいので，群間と群内の比
がなるべく大きくなる係数 a を求めればよい．
　一意的な a を求めるため，a'Wa ＝ 1 と基準化すれば， a はラグランジュの未定係数法を用いて，
Q ＝ a'Ba − λ (a'Wa − 1 )
を最大化すればよい．Q を a で偏微分して 0 とおくと，
となるので，
              
という方程式が得られる．これより，a は W-1B の固有 ベクトルであることがわかる．固有値 λ は，λ ＝ a'Ba となり，線形結合 u の群間平方和となっている．
　W-1B の階数 m は，m ≦ min(p, g−1) であり，全部で m 個の正の固有値 が得られ，残りは 0 固有値となる．正の固有値を大きい順に並べて，
λ1 ≧ λ2 ≧ … ≧ λm
とおく．個々の固有値 λi に対応する固有ベクトル ai により， 線形結合 ui ＝ ai'x が得られる． u1 が最もよく群間の変動を説明し，u2 がその次に群間の変動を説明する， というようになる．この ui を正準変量と呼ぶこともある．
　通常は，u1，u2 で群間の変動のかなりの部分が説明されるので，データを この２つの正順得点で (u1，u2 ) 平面上にプロットし，判別関数で平面 を分割すれば，判別の様子が視覚化される．
　また，i≠j のとき，ai'Waj ＝ 0，が成り立つ ので，u1 と u2 は群内で無相関となる．しかし， ai'aj ≠ 0 なので，ui 軸と uj 軸は直交しない．この点が主成分分析による次元の減少とは異なっている．

５．多変量・多群の判別例

　MASS ライブラリィーのガラス破片データ（fgl）を用いて多変量，多群の判別例を解説する。このデータは，出所により6つに分類された鑑識データの214個のガラス破片に関して，反射特性と成分に関する9つの項目に対するデータである。 
library(MASS)
head(fgl)
nf <- which(fgl$type == "WinF")
nn <- which(fgl$type == "WinNF")
nv <- which(fgl$type == "Ven")
nc <- which(fgl$type == "Con")
nt <- which(fgl$type == "Tabl")
nh <- which(fgl$type == "Head")
type2 <- rep("", nrow(fgl))
type2[nf] <- "F"
type2[nn] <- "N"
type2[nv] <- "V"
type2[nc] <- "C"
type2[nt] <- "T"
type2[nh] <- "H"
typecol <- rep(0, nrow(fgl))
typecol[nf] <- 1
typecol[nn] <- 2
typecol[nv] <- 3
typecol[nc] <- 4
typecol[nt] <- 5
typecol[nh] <- 6

fgl.pca <- prcomp(scale(fgl[,-10]))
summary(fgl.pca)
plot(fgl.pca$x[,1:2], type="n")
text(fgl.pca$x[,1:2], type2, col=typecol+1, cex=0.7)
title(main="主成分スコア")

fgl.lda <- lda(type ~ ., data=fgl)
fgl.lda
fgl.pred <- predict(fgl.lda)
#
#
xp <- seq(-3.6,7.4, length=len)
yp <- seq(-8,3.8, length=len)
fglT <- expand.grid(LD1=xp, LD2=yp)
mld1 <- tapply(fgl.pred$x[,1], fgl$type, mean)
mld2 <- tapply(fgl.pred$x[,2], fgl$type, mean)

pred <- NULL
for(i in 1:nrow(fglT)){
  d1 <- dnorm(fglT[i,1],mean=mld1, sd=1)  
  d2 <- dnorm(fglT[i,2],mean=mld2, sd=1) 
  pred <- rbind(pred, fgl.lda$prior*d1*d2/sum(fgl.lda$prior*d1*d2))
}
#
# 
plot(fgl.pred$x[,1:2], type="n")
text(fgl.pred$x[,1:2], type2, col=typecol+1, cex=0.7)
title(main="正準判別スコアと判別境界")
#
# 判別境界表示
zp <- pred[,1] - pmax(pred[,2], pred[,3], pred[,4], pred[,5], pred[,6]) 
  contour(xp, yp, matrix(zp, len), add=T, levels=0) 
zp <- pred[,2] - pmax(pred[,1], pred[,3], pred[,4], pred[,5], pred[,6]) 
  contour(xp, yp, matrix(zp, len), add=T, levels=0)
zp <- pred[,3] - pmax(pred[,1], pred[,2], pred[,4], pred[,5], pred[,6]) 
  contour(xp, yp, matrix(zp, len), add=T, levels=0)
zp <- pred[,4] - pmax(pred[,1], pred[,2], pred[,3], pred[,5], pred[,6]) 
  contour(xp, yp, matrix(zp, len), add=T, levels=0)
zp <- pred[,5] - pmax(pred[,1], pred[,2], pred[,3], pred[,4], pred[,6]) 
  contour(xp, yp, matrix(zp, len), add=T, levels=0)
#
# データ，判別表
fgltab.lda <- table(fgl$type, fgl.pred$class)
fgltab.lda
1 - sum(diag(fgltab.lda))/sum(fgltab.lda)
#
#
# force predict to return class labels only
mypredict.lda <- function(object, newdata)
  predict(object, newdata = newdata)$class
# 10-fold CV 10 回
library(ipred)
set.seed(41)
error.lda <- numeric(10)
for(i in 1:10) error.lda[i] <- 
  errorest(type ~ ., data = fgl,
  model = lda, predict= mypredict.lda)$error
summary(error.lda)
#
#
# SVM
library(e1071)
fgl.svm <- svm(type ~ ., data = fgl)
# データ，判別表
fgltab.svm <- table(fgl$type, fitted(fgl.svm))
fgltab.svm
1 - sum(diag(fgltab.svm))/sum(fgltab.svm)
# 10-fold CV 10 回
set.seed(563)
error.SVM <- numeric(10)
for (i in 1:10) error.SVM[i] <-
  errorest(type ~ ., data = fgl,
  model = svm)$error
summary(error.SVM)
#
#
library(randomForest)
set.seed(17)
fgl.rf <- randomForest(type ~ ., data = fgl,
	mtry = 2, importance = TRUE,
	do.trace = 100)
print(fgl.rf)
fgltab.rf <- table(fgl$type, fgl.rf$predicted)
fgltab.rf
1 - sum(diag(fgltab.rf))/sum(fgltab.rf)
#
# 10-fold CV 10 回
library(ipred)
set.seed(131)
error.RF <- numeric(10)
for(i in 1:10) error.RF[i] <-
  errorest(type ~ ., data = fgl,
  model = randomForest, mtry = 2)$error
summary(error.RF)
#

主成分分析と正準判別の比較

　主成分はデータのちらばりを表現するが，正準判別はクラスがなるべく分離するようにデータを配置する．
pca lda

正準判別，SVM，ランダムフォレストの予測誤差

　10 fold クロスバリデーション（CV）の10回の平均で予測誤差を比較．ランダムフォレストが良く，しかも，判別誤差と予測誤差の差が小さかった． 
　正準判別
        WinF WinNF Veh Con Tabl Head
  WinF    52    15   3   0    0    0
  WinNF   17    54   0   3    2    0
  Veh     11     6   0   0    0    0
  Con      0     5   0   7    0    1
  Tabl     1     2   0   0    6    0
  Head     1     2   0   1    0   25
判別誤差：0.3271，予測誤差：0.3818

　SVM
        WinF WinNF Veh Con Tabl Head
  WinF    59    11   0   0    0    0
  WinNF   15    61   0   0    0    0
  Veh     10     7   0   0    0    0
  Con      0     0   0  13    0    0
  Tabl     0     1   0   0    8    0
  Head     1     0   0   0    0   28
判別誤差：0.2103，予測誤差：0.2935

　ランダムフォレスト
       WinF WinNF Veh Con Tabl Head
  WinF    63     6   1   0    0    0
  WinNF   11    61   1   1    2    0
  Veh      7     4   6   0    0    0
  Con      0     3   0   9    0    1
  Tabl     0     1   0   0    8    0
  Head     1     3   0   0    0   25
判別誤差：0.1916，予測誤差：0.2065

参考文献

  1. Anderson T. W. (1958). An Introduction to Multivariate Statistical Analysis, John Wiley & Sons, Inc.
  2. Morrison, D. F. (1976). Multivariate Statistical Methods, McGraw-Hill.
  3. 奥野忠一，久米均，芳賀敏郎，吉沢正(1971)．多変量解析法，日科技連．
  4. 竹内啓，柳井晴夫(1972)．多変量解析の基礎，東洋経済新報社．
  5. Venables, W. N. and Ripley, B. D. (1999). Modern Applied Statistics with S-Plus, Springer-Verlag N.Y.
  6. 甘利俊一，金明哲，村上征勝，永田昌明，大津起夫，山西健司（2003）． 言語と心理の統計−ことばと行動の確率モデルによる分析　統計科学のフロンティア10，岩波書店
Copyright (C) 2007, Hiroshi Omori. 最終更新：2017年10月 3日