主成分分析（R）

PRML(12章)のまとめです．

観測値 $\{x_n\}_{n=1}^N\subset \mathbf{R}^D$ があるとき，$u_1\in\mathbf{R}^D$への観測値の射影が最大になるような $u_1\in\mathbf{R}^D$ を探すことを考えます．

$u_1$ は観測値 $\{x_n\}_{n=1}^N$ の一番の【特徴】を表す方向であると考えられます．極端な例でいうと，$u_1$ 方向に射影したとき観測値がすべて一点に集約されてしまっては，観測値の持っている情報は無くなってしまうと考えられます．そこで射影したときに射影した観測値の分散が最大になる方向は，射影したときの観測値の情報を最大限保持する方向であると考えられます．

$u_1$だけでなく，$u_2,\ldots, u_M$まで同様に考えて，$D$次元の観測値$\{x_n\}_{n=1}^N\subset \mathbf{R}^D$を$M$次元の$u_1,\ldots,u_M$の張る空間への射影を考えたのが主成分分析になります．

高次元データを$2$次元に射影して可視化するときや，データの圧縮，データの白色化に使われます．

Rによる例

主成分分析はRの関数prcomp()を使うことでできます．以下ではFisherのあやめのデータ(iris)で主成分分析した結果を載せてあります．主部分空間(2次元)で見たときに種によって結構分かれている気がします．

#iris data
#pca(prcomp)

data(iris)
iris.data <- iris[, 1:4]
iris.lab <- iris[, 5]
result.pca <- prcomp(data.matrix(iris.data), scale = TRUE)

#label
lab <- as.integer(iris.lab)
lev <- levels(iris.lab)

#plot
plot(result.pca$x[, 1:2], col=lab, pch=lab)
par(new = TRUE)
legend("topright", legend = lev, col = 1:3, pch = 1:3)

ちなみに第3, 4主成分（寄与率の最も低い二つの主成分）に関して描画すると以下のようになります:

定性的な理解通り，寄与率の低い方向に関して描画しても，あまり元のデータの特徴は捉えられていないかと思います．

分散最大化による定式化

具体的に $u_1$ を求める方法は以下のようになります:

（射影されたデータの分散）$=\frac{1}{N}\sum_{n=1}^N(u_1^{\top}(x_n-\bar{x}))$

これを $u_1^{\top}u_1=1$ の制約条件（方向だけに興味があるので規格化の制約をつけます）のもとで解くと（ラグランジュの未定乗数法）

\[
u_1=共分散行列の最大固有値に対応する固有ベクトル
\]

となります．上位$M$個の固有ベクトルを考えれば $u_1,u_2,\ldots,u_M$ に相当します．射影誤差の最小化としても $u_1,\ldots,u_M$ が得られます．射影誤差の最小化はデータの近似（データ圧縮）の考え方に基づきます．

確率的主成分分析（最尤推定量としての解釈）

上で述べた主成分分析の結果は，ある確率モデルの最尤推定量として表現することができます．最尤推定量として表現できること自体も面白い性質ですが，この表現により主成分分析のベイズ的扱いが可能になるという実際的な利点もあります．

主成分分析は次のような確率モデルの最尤推定量として得ることができます：

$z$を潜在変数とします．$z$ は$M$次元の標準正規分布 $\mathcal{N}(z|0, I)$ に従うとします:
\[
p(z)=\mathcal{N}(z|0, I)
\]
$D\times M$ 行列 $W$ と $\mu, \sigma^2$ をパラメータとして，$z$ が与えられた下での $x$ の確率分布を
\[
p(x|z)=\mathcal{N}(x|Wz+\mu, \sigma^2I)
\]
とします．この確率モデルのパラメータ $W$ の最尤推定量 $\hat{W}$ の$M$本の列ベクトルが主成分分析で得られる$u_1,\ldots, u_M$に一致します．

この確率モデルの解釈は次のようになると思います：

$\mu$は観測値$\{x_n\}_{n=1}^N$の平均$\bar{x}$で，$z$は$\mu$からの "$M$" 次元空間内でのデータの分布を表します．$W$はデータの分布している潜在的な低次元の線形空間（多様体）から$D$次元への写像を表します（$M$次元空間の標準基底の行先（$D$次元空間の直交系）を選べば写像が決まります）．$z$は重要で，ただの潜在変数というだけでなく，観測データの中には低次元（$M$次元）の自由度しかないという気持ちを表したものであると考えられます．

潜在変数を導入した確率モデルを考えることで，EMアルゴリズムを使って主成分分析を行うことができます．これはデータの次元$D$が非常に大きく，固有ベクトルを求めることが計算的に厳しい時に役に立ちます（確率的解釈をしたことで得られる利点の一つです）．

ベイズ的主成分分析

上にも述べたように主成分分析の確率的解釈のおかげでベイズ推定を行うことができます．ベイズ的に考えることの利点の一つに関連自動決定ということがあげられます．上で述べた主成分分析では，主部分空間 $W$ の次元 $M$ を固定していましたが，データの可視化（$M=2$）などの場合以外は$M$をどう決めればよいかは不明です．ベイズでは事前分布を決めることで自動的に"枝刈り"して次元 $M$ を調節します（モデルの複雑度とデータへの当てはまりの良さのトレードオフ）．とPRMLには説明がありますが，これも事前分布の入れ方によって変わるのかなと思います．つまり自分で主部分空間の次元を決める代わりに，事前分布を決める必要が出てくると思います．

カーネル主成分分析

観測データの主成分を見るのではなく，一回特徴ベクトルでデータを変換してから主成分分析を行う方法がカーネル主成分分析になります．
こうすることで，非線形の部分空間（多様体）に射影することができます．この特徴ベクトルを陽に指定しなくても，内積の形だけがわかれば特徴ベクトルで変換した方で主成分分析ができます（カーネルトリック）．

実際にRでカーネル主成分分析をやってみます．
ガウシアンカーネルでカーネルのパラメータsigmaを0.0003にした場合です．同様にsigmaを0.1, 1, 10, 100に変えて主成分分析した結果です．

#package 'kernlab' = kernel pca
library(kernlab)

data(iris)
iris.data <- iris[, 1:4]
iris.lab <- iris[, 5]

#label
lab <- as.integer(iris.lab)
lev <- levels(iris.lab)

#kernel pca
#kernel = gaussian, sigma=0.0003
#feature dim =2
result <- kpca(~., data=iris.data, features=2, kpar=list(sigma=0.0003))
plot(rotated(result), col=lab, pch=lab)

#legend
legend("topright", legend=lev, pch=1:3, col=1:3)

sigmaによって結果が大きく違うことがわかります．