#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