R＋PythonでカーネルPCA（おまけ付き）

昨日の勉強会で「Rは便利だけど、図がかっこよくない」という話になったので

Rpyで出力した結果をmatplotlibで出力する場所まで書いてみた。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

#rpy2が使える環境であることを確認
import numpy as NP
import random
try:
	import rpy2.robjects as robjects
except:
	print "Rpy2をインストールしてください。"
	print "UBUNTUの場合、apt-get install python-rpy2　でインストールできます。"
	print "easy_install が入っている場合は easy_install rpy2　でインストールできます。"
	print "その他の場合は　http://sourceforge.net/projects/rpy/files/rpy2/　からダウンロードして対応ください。"


#Rのカーネル法のパッケージであるkernlabをダウンロード
try:
	robjects.r["library"]("kernlab")
except:
	print "CRANのミラーサイトを使ってkernlabパッケージをダウンロードしてください。"
	robjects.r["install.packages"]("kernlab")

#Numpyのarrayをdata.frameに変換する関数
def MyAsDataframe(Dataset):
	DataArray = NP.array(Dataset)
	Dataframe = {}
	for data in NP.transpose(DataArray):
		if data[0] == "Species":
			Dataframe[data[0]] = robjects.StrVector(data[1:])
		else:
			Dataframe[data[0]] = robjects.FloatVector(data[1:])
	return robjects.r["data.frame"](**Dataframe)
#Hebbian PCA
#デフォルトの値が用いられています。もしも、パラメータをいじくりたい場合は変更が必要です。
def kernelPCA(formula, dataset,dot,f):
	RFResult =  robjects.r["kpca"](formula, dataset,kernel=dot,features = f)
	return RFResult
Dataset = []
label = []
for dataLine in open("iris.txt"):
	SplitedLine = dataLine.strip().split("\t")
	Dataset.append(SplitedLine[:-1])
	if SplitedLine[-1] != "Species":
		label.append(SplitedLine[-1])
RDataframe = MyAsDataframe(Dataset)
#用いる式を作成。この場合は全変数を入れるので以下の様に記載する。
formula = robjects.RFormula('~. ')

#カーネルを指定する。今回はANOVAカーネルを使用。
Anovadot = robjects.r["anovadot"](sigma = 1, degree = 1)
KPCAResult = kernelPCA(formula, RDataframe,Anovadot,3)

#主成分得点を記録する。
pcv = NP.array(robjects.r["pcv"](KPCAResult))
TransValue = NP.transpose(pcv)
#ここからグラフ作成に移る。
import numpy as NP
try:
	from mpl_toolkits.mplot3d import Axes3D
	import matplotlib.pyplot as plt
except:
	print "matplotlibをインストールしてください。"
	print "UBUNTUの場合、apt-get install python-matplotlib　でインストールできます。"
	print "easy_install が入っている場合は easy_install matplotlib　でインストールできます。"
	print "その他の場合は　http://matplotlib.sourceforge.net/　からダウンロードして対応ください。"

#キャンパスの設定。(3Dプロット編)
"""
fig = plt.figure()
ax = Axes3D(fig)

x = TransValue[0]
y = TransValue[1]
z = TransValue[2]
ax.scatter(x, y, z)

ax.set_xlabel('KPC1')
ax.set_ylabel('KPC2')
ax.set_zlabel('KPC3')

plt.show()
"""
#キャンパスの設定。(2Dプロット編)
fig = plt.figure()
ax = fig.add_subplot(111)

x = TransValue[0]
y = TransValue[1]
Rand = random.random()
color = []
for l in label:
	color.append(len(l)*100*Rand)
ax.scatter(x, y,c=color,alpha=0.50)
ax.set_title('Kernel PCA')

ax.set_xlabel('KPC1')
ax.set_ylabel('KPC2')

plt.show()

これで出力されるのが以下の図二つ。
ま、そこまでぱっとしないか。（色彩センスがないので。）

Exampleはそこそこかっこいいのになぁ。