R+Pythonでマルコフ連鎖モンテカルロ

バイト先で広告の効果測定を頼まれたので、前々から興味があったMCMCを使って測定を行った。

PythonにはPyMCという専用のパッケージがあるけど、そっちはどうも小難しい感じがしたので

Rpyを使って使えるようにしてみた。

使い方はRで回帰分析をするときとほとんど同じでデータフレームと数式を入れると結果を出すという形式にした。

用いたデータは２００４年1月〜２００９年9月までの円ドルレートとアメリカの失業率。（http://www.mediafire.com/?mfo5mmezow3）

MCMCで解析した結果とその際に行われた推定の過程をPDFファイルで吐き出してくれるようにした。

解析した結果はと言うと、失業率と円/ドルはそこまで関係が強くないらしい（回帰分析でR2が0.63くらい）

単変量であたるくらい単純なものだったら、FXで損する人はいないか。

以下、今回用いたソースコード。

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

#rpy2が使える環境であることを確認
import numpy as NP
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のMCMCパッケージであるMCMCpackをダウンロード
try:
	robjects.r["library"]("MCMCpack")
except:
	print "CRANのミラーサイトを使ってMCMCパッケージをダウンロードしてください。"
	robjects.r["install.packages"]("MCMCpack")

#Numpyのarrayをdata.frameに変換する関数
def MyAsDataframe(array):

	DataArray = NP.array(Dataset)
	Dataframe = {}
	for data in NP.transpose(DataArray):
		Dataframe[data[0]] = robjects.FloatVector(data[1:])
	return robjects.r["data.frame"](**Dataframe)
#Regress(回帰),poisson(ピアソン回帰)
#デフォルトの値が用いられています。もしも、パラメータをいじくりたい場合は変更が必要です。
def MCMC(fomula,Dataframe,MCtype):
	if MCtype == "regress":
		MCMCResult =  robjects.r["MCMCregress"](fomula,RDataframe)
	elif MCtype == "poisson":
		MCMCResult =  robjects.r["MCMCpoisson"](fomula,RDataframe)
	else:
		"Type Error!"
	return MCMCResult
#単なる回帰分析
def LinerRegress(fomula,Dataframe):
	return robjects.r["lm"](fomula,RDataframe)
#データのハンドリング
Dataset = []
for dataLine in open("sample.txt"):
	Dataset.append(dataLine.strip().split("\t"))
RDataframe = MyAsDataframe(Dataset)

#MCMC開始
#数式を作る
fomula = robjects.RFormula('YD ~ jobless')
#MCMCの結果を表示する。
MCMCResult = MCMC(fomula,RDataframe,"regress")
robjects.r["pdf"]("Report.pdf")
robjects.r["plot"](MCMCResult)
robjects.r["dev.off"]()
print "MCMC-Result"
print robjects.r["summary"](MCMCResult)

#回帰分析の結果を表示する。
LMResult = LinerRegress(fomula,RDataframe)
print "LM-Result"
print robjects.r["summary"](LMResult)