今回は「Pythonで散布図」の続きで、散布図に使った二系統のデータから、最小二乗法で求めた単回帰直線を描いてみたいと思います。
単回帰直線を求めるためのライブラリ SciPy は、前回の「Mac に SciPy をインストール」をみてインストールしてください。
単回帰直線を求める
単回帰直線は、scipy.stats.linregress(x, y) で計算します。
戻り値として、直線の傾き、切片、相関係数などが返ってきます。
from scipy import stats slope, intercept, r_value, _, _ = stats.linregress(x, y) x: データ1 y: データ2 slope: 単回帰直線の傾き intercept: 単回帰直線の切片（y = ax + b の b に相当） r_value: 相関係数
直線を描く
(x1, y1), (x2, y2), ... という点を通る直線は、matplotlib.lines.Line2D を使って、
from matplotlib import lines from matplotlib import pyplot as plt fig = plt.figure() ax = fig.add_subplot(111) line = lines.Line2D([x1, x2, ...], [y1, y2, ...]) ax.add_line(line)
と書くことができます。
なので、例えば一番シンプルに (0, 150), (50, 0) の2点間の線を書くとすると、
line = lines.Line2D([0, 50], [150, 0])
となります。
（2013.8.25追記）
また、直線の傾き（slope）とY切片（intercept）を使い、Y切片(0, intercept)、X切片(50, 50 * slope + intercept) の2点間の直線を引くことで、単回帰直線を描くこともできます。
line = lines.Line2D([0, 50], [intercept, 50 * slope + intercept])
となります。でも、次のように、lambda式を使って関数を作った方が分かりやすくていいですね。
func = lambda x: x * slope + intercept line = lines.Line2D([0, 50], [func(0), func(50)])
サンプル
scipy.stats.linregress(x, y) で求めた傾き（slope）とY切片（intercept）を使い、
Y切片(0, intercept)、X切片(50, 50 * slope + intercept) の2点間の直線を引くことで、単回帰直線を描いています。
test_linear_regression.py
#! /usr/bin/env python # -*- coding: utf-8 -*- import numpy import xlrd from matplotlib import lines from matplotlib import mlab from matplotlib import pyplot as plt from scipy import stats def get_data(sheet, rowx, colx): data = [] for row in range(rowx, sheet.nrows): value = sheet.cell(row, colx).value if value != '': data.append(value) data = numpy.array(data) return data if __name__ == '__main__': book = xlrd.open_workbook('/Users/akiyoko/Documents/temp/2nd_test.xls') sheet = book.sheet_by_name('Statistics (total score)') data1 = get_data(sheet, 9, 5) # データの起点はF10 data2 = get_data(sheet, 9, 6) # データの起点はG10 print 'data1=%s' % data1 print 'size of data 1=%d' % len(data1) print 'mean value of data 1=%.1f' % numpy.mean(data1) print 'data2=%s' % data2 print 'size of data 2=%d' % len(data2) print 'mean value of data 2=%.0f' % numpy.mean(data2) # 線形回帰 slope, intercept, r, _, _ = stats.linregress(data1, data2) print 'slope=%.2f' % slope print 'intercept=%.2f' % intercept print 'coefficient of correlation=%.2f' % r # 共通初期設定 plt.rc('font', **{'family': 'serif'}) # キャンバス fig = plt.figure() # プロット領域（1x1分割の1番目に領域を配置せよという意味） ax = fig.add_subplot(111) # 散布図 sc = ax.scatter(data1, data2, s=25, marker='x', color='b') # 単回帰直線 func = lambda x: x * slope + intercept line = lines.Line2D([0, 50], [func(0), func(50)], color='r') ax.add_line(line) # X, Y方向の表示範囲 ax.set_xlim(0, 50) ax.set_ylim(0, 6000) # タイトル ax.set_title('Scatter Graph: r=%.2f' % r, size=16) ax.set_xlabel('Score', size=14) ax.set_ylabel('Error Severity', size=14) # 凡例 ax.legend((sc,), ('2nd Test',), scatterpoints=1, loc='upper left', fontsize=10) # グリッド表示 ax.grid(True) plt.show()
実行結果
$ python test_linear_regression.py data1=[ 42. 42. 40. 38. 38. 38. 38. 35. 35. 35. 35. 34. 34. 34. 34. 33. 33. 33. 33. 32. 32. 32. 30. 30. 30. 30. 30. 30. 30. 30. 29. 29. 29. 28. 28. 28. 28. 28. 27. 27. 27. 27. 27. 27. 26. 26. 26. 26. 26. 26. 25. 24. 24. 24. 24. 23. 23. 23. 23. 23. 23. 23. 23. 23. 22. 22. 22. 22. 22. 22. 22. 22. 21. 21. 20. 20. 20. 20. 20. 19. 19. 19. 19. 19. 18. 18. 17. 17. 16. 16. 16. 16. 15. 15. 15. 15. 14. 13. 13. 13. 13. 12. 12. 11. 9.] size of data 1=105 mean value of data 1=24.9 data2=[ 576. 668. 682. 795. 814. 1046. 1414. 1084. 1107. 1196. 1533. 1205. 1327. 1427. 1434. 1575. 1578. 1683. 1772. 1522. 1585. 1665. 1542. 1642. 1677. 1795. 1813. 1824. 1865. 1980. 1483. 1830. 2257. 1763. 2120. 2207. 2712. 2828. 1793. 1846. 1959. 2049. 2124. 2677. 2223. 2236. 2511. 2910. 3028. 3371. 2528. 2196. 2274. 2469. 2539. 2155. 2158. 2474. 2601. 2733. 2793. 2863. 3398. 3585. 1984. 2413. 2574. 2752. 2892. 3212. 3318. 3794. 2839. 3406. 2435. 2571. 3011. 3133. 3265. 2506. 2931. 2983. 3455. 3609. 2895. 3586. 3722. 3846. 2883. 2937. 3213. 3456. 3161. 3623. 3888. 4035. 4400. 4051. 4364. 4533. 4825. 4063. 4289. 3762. 4453.] size of data 2=105 mean value of data 2=2510 slope=-120.83 intercept=5513.78 coefficient of correlation=-0.92
ちなみに、直線を描くコードは、
x_range = numpy.array([0, 50]) line = lines.Line2D(x_range, x_range * gradient + intercept, color='r')
としてもよいと思います。
ただし、
x_range = [0, 50] line = lines.Line2D(x_range, x_range * gradient + intercept, color='r')
とするとエラーになってしまいます。
[0, 50] * 3 と numpy.array([0, 50]) * 3 の動きが違うからですね。
>>> [0, 50] * 3 [0, 50, 0, 50, 0, 50] >>> import numpy >>> numpy.array([0, 50]) * 3 array([ 0, 150])
