読者です 読者をやめる 読者になる 読者になる

naoya_t@hatenablog

いわゆるチラシノウラであります

Oil FlowデータでPCA

PRML Python

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

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

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)ではどうやるんだろう、と試行錯誤したメモ。

とりあえず
f:id:n4_t:20111229201309p:plain
この図が出したいわけで

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みたいなラベルごとのマーカー指定は無理か)
f:id:n4_t:20111229205903p:plain
ちょっと違う感じなわけで

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()

みたいに3次元で描画すれば
f:id:n4_t:20111229210249p:plain
ぐりぐり探せばまあ似てなくもないアングルがあったりする感じ
(いや、このコードは後で何かを3次元プロットしたい時にコピペできるように貼ってるだけなんだからね)

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

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) を使うと今度は左右が反転したりとかね)。