#freeze
[[FrontPage]]

#contents

2011/07/06からのアクセス回数 &counter;

ここで紹介した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を解く [#i8fe1197]

はやり主成分分析(PCA)は、Rの得意とする分野で、
[[PRML 12章 主成分分析を試す(棒読み>http://d.hatena.ne.jp/n_shuyo/20100524/pca]]
の実装をSageで表示する方法とRの計算結果をSageに渡す方法について紹介します。	


*** Rなら主成分分析が1行で書ける [#sb5539ba]

Oil Flowのデータを使って主成分分析をします。r関数を使ってRのコマンドを記述します。

各行のデータが「一様」、「環状」、「層状」のいずれかを示すラベルファイル
DataTrnLbls.txtを読み込み、「青」、「緑」、「赤」でプロットしています。

Rのプロット結果をSageで表示するには、一度ファイルに出力する必要があります。
r.pdf, r.dev_offで囲んでプロット処理を記述します。
結果の出力は、html関数でイメージ読み込みします。

sageへの入力:
#pre{{
# Rでデータを読み込みPCAを計算
fileName = DATA + 'DataTrn.txt'
oilflow = r("oilflow <- read.table('%s')" %fileName)
result = r("result <- prcomp(oilflow)")
}}

sageへの入力:
#pre{{
# ラベルの読み込み
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">')
}}

&ref(fig1.png);

** Rの結果をsageに渡す [#d7a29e22]

Rの計算結果をsageに渡してためにsageobj関数または、_save_メソッドが提供されています。

以下の例では、主成分分析の結果のλ1, λ2成分, λ3成分をsageの変数rsに取り込んでいます。
この場合の結果は、sageのマトリックスとして返され、rs[0]で第1行目を出力しています。


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


*** sageの3次元プロットで表示 [#x4535360]

Rで読み込んだlabelsは、辞書形式に変換されます。'DATA'のV1, V2, V3に列単位でラベル情報が
セットされます。これをまとめてlbsにセットします。

3次元プロットは、λ1, λ2成分, λ3成分を[x, y, z]にセットし、point3dでプロットするだけです。

sageでは他のライブラリとのインタフェースを用意し、データの変換ができるのでこれらを組み合わせ
ることが簡単にできます。

sageへの入力:
#pre{{
#ラベル情報をsageの形式に変換すると'DATA'ディクショナリに列単位でV1, V2のようにセットされる
lb = sageobj(labels); lb
}}
sageからの出力:
#pre{{
{'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への入力:
#pre{{
# これをzipでまとめて使う
lbs = zip(lb['DATA']['V1'],lb['DATA']['V2'],lb['DATA']['V3'])
}}


sageへの入力:
#pre{{
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()
}}

&ref(fig2.png);

** コメント [#mf349858]
#vote(おもしろかった[4],そうでもない[0],わかりずらい[0])
#vote(おもしろかった[5],そうでもない[0],わかりずらい[0])

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

#comment_kcaptcha

トップ   編集 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
SmartDoc