10章 Rとの連携でPCAを解く
はやり主成分分析（PCA）は、Rの得意とする分野で、 PRML 12章 主成分分析を試す(棒読み の実装をSageで表示する方法とRの計算結果をSageに渡す方法について紹介します。
Rなら主成分分析が１行で書ける
Oil Flowのデータを使って主成分分析をします。r関数を使ってRのコマンドを記述します。
各行のデータが「層状」、「環状」、「一様」（順序は要確認です）のいずれかを示すラベルファイル DataTrnLbls.txtを読み込み、「青」、「緑」、「赤」でプロットしています。
Rのプロット結果をSageで表示するには、一度ファイルに出力する必要があります。 r.pdf, r.dev_offで囲んでプロット処理を記述します。 結果の出力は、html関数でイメージ読み込みします。
# Rでデータを読み込みPCAを計算
fileName = DATA + 'DataTrn.txt'
oilflow = r("oilflow <- read.table('%s')" %fileName)
result = r("result <- prcomp(oilflow)")
|
|
# ラベルの読み込み
fileName = DATA + 'DataTrnLbls.txt'
labels = r("oilflow.labels <- read.table('%s')" %fileName)
# プロットファイル名の設定
filename = DATA+'pca.pdf'
r.pdf(file='"%s"' %filename)
# 結果のプロット
r("col <- colSums(t(oilflow.labels) * c(4,3,2))")
r("pch <- colSums(t(oilflow.labels) * c(3,1,4))")
r("plot(result$x[,1:2], col=col, pch=pch, xlim=c(-3,3), ylim=c(-3,3))")
r.dev_off()
# 式を変えたときにはブラウザーで再読込必要
html('<img src="pca.pdf">')
|
|
Rの結果をsageに渡す
Rの計算結果をsageに渡してためにsageobj関数または、_save_メソッドが提供されています。
以下の例では、主成分分析の結果のλ1, λ2成分, λ3成分をsageの変数rsに取り込んでいます。 この場合の結果は、sageのマトリックスとして返され、rs[0]で第1行目を出力しています。
# Rの主成分分析の結果からλ1, λ2成分, λ3成分をsageの変数rsに取り込む
rs = sageobj(r("result$x[,1:3]"))
print rs[0]
(0.853313762157, 0.409097298967, 0.499457870747) (0.853313762157, 0.409097298967, 0.499457870747) |
sageの3次元プロットで表示
Rで読み込んだlabelsは、辞書形式に変換されます。'DATA'のV1, V2, V3に列単位でラベル情報が セットされます。これをまとめてlbsにセットします。
3次元プロットは、λ1, λ2成分, λ3成分を[x, y, z]にセットし、point3dでプロットするだけです。
sageでは他のライブラリとのインタフェースを用意し、データの変換ができるのでこれらを組み合わせ ることが簡単にできます。
#ラベル情報をsageの形式に変換すると'DATA'ディクショナリに列単位でV1, V2のようにセットされる
lb = sageobj(labels); lb
{'row_names': [None, -1000], '_r_class': 'data.frame', '_Names': ['V1',
'V2', 'V3'],
'DATA': {'V1': [1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 1,
0, 0, 1, 0,
以下省略
{'row_names': [None, -1000], '_r_class': 'data.frame', '_Names': ['V1', 'V2', 'V3'],
'DATA': {'V1': [1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0,
以下省略
|
# これをzipでまとめて使う
lbs = zip(lb['DATA']['V1'],lb['DATA']['V2'],lb['DATA']['V3'])
|
|
N = len(lbs)
plt = Graphics()
for n in range(N):
[x, y, z] = rs[n]
if lbs[n][0] == 1:
plt += point3d([x, y, z], rgbcolor='blue')
elif lbs[n][1] == 1:
plt += point3d([x, y, z], rgbcolor='green')
else:
plt += point3d([x, y, z], rgbcolor='red')
plt.show()
|
Sleeping...
 |
|
|