akiyoko blog

akiyoko の IT技術系ブログです

Pythonで単回帰直線

今回は「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

f:id:akiyoko:20130616003219p:plain


ちなみに、直線を描くコードは、

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])