sage/PRML - SageでRのPCAを試す

FrontPage

2011/07/06からのアクセス回数 7215

ここで紹介したSageワークシートは、以下のURLからダウンロードできます。

http://sage.math.canterbury.ac.nz/home/pub/109/

また、Sageのサーバを公開しているサイト（http://sage.math.canterbury.ac.nz/ 、http://www.sagenb.org/）にユーザIDを作成することで、ダウンロードしたワークシートを アップロードし、実行したり、変更していろいろ動きを試すことができます。

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関数でイメージ読み込みします。

sageへの入力:

# Rでデータを読み込みPCAを計算
fileName = DATA + 'DataTrn.txt'
oilflow = r("oilflow <- read.table('%s')" %fileName)
result = r("result <- prcomp(oilflow)")

sageへの入力:

#　ラベルの読み込み
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行目を出力しています。

sageへの入力:

# Rの主成分分析の結果からλ1, λ2成分, λ3成分をsageの変数rsに取り込む
rs = sageobj(r("result$x[,1:3]"))
print rs[0]

sageからの出力:

(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への入力:

#ラベル情報をsageの形式に変換すると'DATA'ディクショナリに列単位でV1, V2のようにセットされる
lb = sageobj(labels); lb

sageからの出力:

{'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,
以下省略

sageへの入力:

# これをzipでまとめて使う
lbs = zip(lb['DATA']['V1'],lb['DATA']['V2'],lb['DATA']['V3'])

sageへの入力:

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

コメント †

皆様のご意見、ご希望をお待ちしております。

• ノートブックで fileName = DATA + 'DataTrn.txt' を実行する場合、データの有り場所を指定する必要はないのでしょうか？ターミナルで入力する場合、pwd, ls があるので可能と思います。 -- ysato? 2012-09-06 (木) 18:33:15
• ysatoさま、DATA変数にdataディレクトリのパスがセットされていますので、ワークシートではDATA + ファイル名の形式で指定することでワークシート関連するファイルにアクセスすることができるようになります。 -- 竹本　浩? 2012-09-06 (木) 18:49:08

お名前:

(Input image string)