Create a gist now

Instantly share code, notes, and snippets.

Python + PypeRでPythonからRをつかってみる
# -*- coding:utf-8 -*-
import numpy
import pandas
import pylab
import pyper
n = 200
# データの生成
score_x = numpy.random.normal(171.77, 5.54, n)
score_y = numpy.random.normal(62.49, 7.89, n)
score_x.sort()
score_x = numpy.around(score_x + numpy.random.normal(scale=3.0, size=n), 2)
score_y.sort()
score_y = numpy.around(score_y + numpy.random.normal(size=n), 2)
# 散布図を描く
pylab.scatter(score_x, score_y, marker='.', linewidths=0)
pylab.grid(True)
pylab.xlabel('X')
pylab.ylabel('Y')
# Rで回帰分析
df = {'X': score_x, 'Y': score_y}
df = pandas.DataFrame(df)
r = pyper.R(use_pandas='True')
# Rへデータ渡す
r.assign('df', df)
# Rのコマンド実行
print(r("summary(df)"))
r("result <- lm(Y~X, data=df)")
print(r("summary(result)"))
#予測区間と信頼区間を算出するため
new_x = numpy.arange(155, 190, 0.1)
new_df = pandas.DataFrame({'X': new_x})
r.assign('new', new_df)
# 予測区間(R)
r("prediction <- predict(result, new, interval='prediction')")
# 信頼区間(R)
r("confidence <- predict(result, new, interval='confidence')")
# Python側にとってくる
lm_result = r.get('result$fitted.values')
prediction = pandas.DataFrame(r.get('prediction'))
confidence = pandas.DataFrame(r.get('confidence'))
# 回帰直線, 予測区間, 信頼区間を描く
pylab.plot(score_x, lm_result, 'r', linewidth=2)
pylab.plot(new_x, prediction[1], 'g', linewidth=1)
pylab.plot(new_x, prediction[2], 'g', linewidth=1)
pylab.plot(new_x, confidence[1], 'c', linewidth=1)
pylab.plot(new_x, confidence[2], 'c', linewidth=1)
pylab.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment