Oil FlowデータでPCA

では、データを使って PCA を行ってみよう。

oilflow <- read.table("DataTrn.txt");
result <- prcomp(oilflow)   # Rなら主成分分析が１行で書ける！！！

oilflow.labels <- read.table("DataTrnLbls.txt");
col <- colSums(t(oilflow.labels) * c(4,3,2));  # ラベルごとに色を指定
pch <- colSums(t(oilflow.labels) * c(3,1,4));  # ラベルごとにマーカーを指定
plot(result$x[,1:2], col=col, pch=pch, xlim=c(-3,3),ylim=c(-3,3))

PRML 12章 主成分分析を試す(棒読み - Mi manca qualche giovedi`?

こういうのってPython（というかnumpy/scipy）ではどうやるんだろう、と試行錯誤したメモ。

とりあえず

この図が出したいわけで

matplotlib.mlab.PCA() というそのまんまっぽい関数があるんだけど、これで

from pylab import *
import matplotlib.mlab 
import numpy

X = loadtxt('./DataTrn.txt')
pca = matplotlib.mlab.PCA(X)
proj = pca.Y

labels = loadtxt('./DataTrnLbls.txt', dtype=int)
colors = ['','black','red','green','blue','cyan','magenta','yellow','gray']
labelc = map(lambda _:colors[_],map(sum,labels*[4,3,2]))

scatter(proj[:,0],proj[:,1], color=labelc, marker='x')
show()

みたいにすれば出るんじゃないかと期待したけど（Rみたいなラベルごとのマーカー指定は無理か）

ちょっと違う感じなわけで

from pylab import *
import matplotlib.mlab
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

X = loadtxt('./DataTrn.txt')
pca = matplotlib.mlab.PCA(X)
proj = pca.Y

labels = loadtxt('./DataTrnLbls.txt', dtype=int)
colors = ['','black','red','green','blue','cyan','magenta','yellow','gray']
labelc = map(lambda _:colors[_],map(sum,labels*[4,3,2]))

fig = plt.figure()
ax3d = fig.add_subplot(111, projection='3d')
ax3d.scatter(proj[:,0], proj[:,1], proj[:,2], color=labelc, marker='x')

show()

みたいに３次元で描画すれば

ぐりぐり探せばまあ似てなくもないアングルがあったりする感じ
（いや、このコードは後で何かを３次元プロットしたい時にコピペできるように貼ってるだけなんだからね）

で、最初の出力はどう出したかというと

from pylab import *
import numpy

def pca(X):
  M = X - X.mean(axis=0)
  [latent, coeff] = numpy.linalg.eig(cov(M.T))
  score = dot(M, coeff)
  return score

X = loadtxt('./DataTrn.txt')
proj = pca(X)

labels = loadtxt('./DataTrnLbls.txt', dtype=int)
colors = ['','black','red','green','blue','cyan','magenta','yellow','gray']
labelc = map(lambda _:colors[_],map(sum,labels*[4,3,2]))

axis([-3, 3, -3, 3])
scatter(proj[:,0], -proj[:,1], color=labelc, marker='x')
show()

こんな感じ（、とさらっと言ってみたけど、eig() でやる方法とか svd() でやる方法とか色々試した結果これ）。

上下反転して出てきたのでscatter()に渡す時に補正してます。（cov(M.T) の代わりに dot(M.T, M) を使うと今度は左右が反転したりとかね）。