PRML 12章 確率的主成分分析を試す

PCAを試す に続いて確率的主成分分析(Probability Principal Component Analysis)。
解析的に解けてしまって、閉形式の解がわかっているので実装としてはたいしておもしろくない(いや、いいことなんですけどね)。

M <- 2;
directory <- ".";

argv <- commandArgs(T);
if (length(argv)>0) directory <- commandArgs(T)[1];
if (length(argv)>1) M <- as.integer(commandArgs(T)[2]);

oilflow <- as.matrix(read.table(sprintf("%s/DataTrn.txt", directory)));
oilflow.labels <- read.table(sprintf("%s/DataTrnLbls.txt", directory));
D <- ncol(oilflow);

# mu = mean x_bar (PRML 12.1)
mu <- colMeans(oilflow);

# eigenvalues and eigenvectors of covariance S (PRML 12.3)
e <- eigen(var(oilflow))

# sigma^2 = sum(rest of eigenvalues) / (D - M) (PRML 12.46)	
sigma2 <- mean(e$values[-(1:M)]);

# W_ML = U_M(L_M - sigma^2 I)R (PRML 12.45)
W_ML <- e$vectors[,1:M] %*% diag(e$values[1:M] - sigma2) %*% diag(c(1,-1))

# M = W^T W + sigma^2 I (PRML 12.41)
M_inv <- solve(t(W_ML) %*% W_ML + sigma2 * diag(M));

# projection into principal subspace (PRML 12.48)
z <- t(M_inv %*% t(W_ML) %*% (t(oilflow) - mu))

# draw chart
col <- colSums(t(oilflow.labels) * c(4,3,2));
pch <- colSums(t(oilflow.labels) * c(3,1,4));
plot(z, col=col, pch=pch, xlim=c(-2,4),ylim=c(-4,2))

一番めんどくさいのはデータ共分散行列を求める \bf{S}=\frac{1}{N}\sum_{n=1}^N (\bf{x}_n-\bar{\bf{x}})(\bf{x}_n-\bar{\bf{x}})^T だが、それすら var(oilflow) で済んでしまうのだから、R にはかなわない。
座標系の回転を司る R は任意の直交行列であるが、単位行列にしたところ PCA の図を上下反転したものが出てきたので、見た目を合わせるために第2成分を反転する行列を指定しておいた。


ちょっとずれてはいるものの、数式を見ていて予想できてしまった通り、PCA とほぼ同じ図ができた。
これなら PCA でいいや……じゃあなくて、ここから EM アルゴリズムベイズアプローチへの発展があるのが 確率的 PCA のいいところ、なので、次回はEM。