[Work/Class/Python3の基礎とデータ処理/2_PythonLibraries]

CSVデータからの単純線形回帰分析

気象庁オープンデータ

気象庁|過去の気象データ・ダウンロード」ページから,必要なデータを絞ってCSV形式でダウンロードすることができる.これを分析対象にする.

線形単純回帰分析

xとyの関係が例えば「y = ax + b」という「線形」の「モデル」すなわち一次関数の直線で表現できると仮定し,実際に観測されたデータを用いて,その直線を表現するaとbを求めることを,線形回帰分析という.(xを説明変数(Explanatory Variable)もしくは独立変数,yを目的変数(Response Variable)もしくは従属変数,応答変数などと呼ぶ)

実際には「y = ax1 + bx2 + cx3 + dx4 + e」というモデルも考えられる.これは説明変数が5つ,目的変数が1つの場合である.

説明変数が1つの時「単純回帰分析」と呼び,説明変数が2以上の場合は「重回帰分析」と呼ぶ.とりあえず本項では「単純回帰分析」,特に「y = ax + b」のモデルで捉える単純線形回帰分析を扱う.

(なぜ「特に」かというと,「y = ax^2 + b」というモデルも,求めるべき係数aが説明変数「x二乗」に対して線形であることから,こういうモデルも「単純線形回帰分析」と呼べる,という理由による.この場合は「x二乗」をそのまま説明変数として使うため,計算上そのまま扱えるように説明変数の配列をを二乗になっている状態で作成しておく.)

最小二乗法

現実に観測されたデータには「ノイズ」が乗っている.すなわちどんなデータも綺麗なy = ax + bの直線上にあることは殆どなく,そこからズレていることが前提である.そこで観測された全てのデータの,直線で表されるyの量との「ズレ」の総量を最小にする「最小二乗法」が一般的に用いられる.

numpyライブラリには,この「最小二乗法を用いた単純線形回帰分析」を行う関数が用意されている.polyfit()関数で計算を行い,その結果をpoly1d()に関数に入れてやると「引数に説明変数の配列を取り,戻り値が目的変数の配列」という関数オブジェクトが生成される.この関数オブジェクトを使いプロットしてやる.

ちなみに英語名は「Ordinary Least Squares」であり略称がOLSである.

対象データ

上記の気象庁のオープンデータサイトより,東京都八王子市の1977年〜2016年の8月の月平均気温,日の最高気温の月平均,日の最低気温の月平均,3つが入っているデータをダウンロードした.八王子市でも夏の温暖化が進んでいるとすれば,説明変数xを年(データがある1977年〜2016年を0〜にする),応答変数yを日中最高気温の月平均として,単純線形回帰分析をかけると「y = ax + b」が成立するかもしれない.

ヘッダを削って文字コードUTF-8に直して保存しておく.加工したCSVデータ

コード例

# -*- coding: utf-8 -*-
# OLSFromCSV.py
import numpy, pandas
import codecs
from PyQt5.QtWidgets import QApplication
import matplotlib.pyplot

class OLSFromCSV():
    def __init__(self, csv_file_name, column_index):
        self.__csv_file_name = csv_file_name
        # CSVファイルの何列目を目的変数にするか
        self.__column_index = column_index
        return None

    def read_csv_then_push_stack(self):
        # pandasのCSV読み込み機能を使って読み込み,
        # ndarrayに2次元配列として入れていく
        # Pandsライブラリを使ってCSVファイルを読み込む
        df = pandas.read_csv(self.__csv_file_name, header=-1)
        print(df)
        
        
        # Pandsの読み込みデータはNumPyの2次元arrayでくるので,必要な列のみを取り出す
        self.__response_variable_array = df[self.__column_index]
        self.__num_of_data = self.__response_variable_array.shape[0]
        # 説明変数の配列を作成
        self.__explanatory_variable_array = numpy.arange(self.__num_of_data)
        return self
            
    def make_plot_csv_data_2d(self):
        # Matplotlib.plotで図を描く
        # plotに使うデータを設定: x軸とy軸にそれぞれ配列を指定する.
        matplotlib.pyplot.plot(self.__explanatory_variable_array, self.__response_variable_array, color='gray')

    def plot_show(self, x_label, y_label):
        # plotの軸の名前を設定
        matplotlib.pyplot.xlabel(x_label, fontsize=10, fontname='serif')
        matplotlib.pyplot.ylabel(y_label, fontsize=10, fontname='serif')
        # plotする
        matplotlib.pyplot.show()

    def doOLS(self):
        # x=説明変数の配列
        # y=目的変数の配列
        # intercept=直線が原点(0,0)を通るか切片を持つかどうかのパラメータ
        # Trueだと原点を通らない,つまり「y = ax + b」の「b」を計算できる
        # numpy.polyfit(x, y, n): n次式で説明変数1,目的変数1の回帰分析
        z = numpy.polyfit(self.__explanatory_variable_array, self.__response_variable_array, 1)
        # 計算が終わったpolyfitの結果をpoly1dに入れると1変数1次の式を出してくれる
        p = numpy.poly1d(z)
        # y = ax + bの時,p[0]がa, p[1]がbになる
        # が,同時にpという関数オブジェクトも生成してくれるので,
        # p(x軸の配列)を実行すると y軸用の配列も出力してくれる
	# それをそのままplotする
        matplotlib.pyplot.plot(self.__explanatory_variable_array, p(self.__explanatory_variable_array), color='blue')
	# 式の出力もできる
        print(p)

if __name__ == "__main__":
    ols_from_csv = OLSFromCSV('TempDataAug.csv', 4)
    ols_from_csv.read_csv_then_push_stack()
    ols_from_csv.make_plot_csv_data_2d()
    ols_from_csv.doOLS()
    ols_from_csv.plot_show('Year(1977~2016)', 'Average Max Temp. of Aug.')

結果

式の結果は「0.05397 x + 29.62」であった.これをプロットすると図のようになる.

numpyライブラリの単純線形回帰による八王子市の8月の平均気温と線形回帰の結果

なんとなく合っているような気がする.

Numpy以外のライブラリ

pandasライブラリや, その分析機能だけ取り出して開発を継続しているstatsmodelsなど,簡単に回帰分析ができるライブラリが,Pythonでは多く提供されている(全てnumpyの配列を取る).これらのライブラリは今回やった単純な回帰分析の他に,その計算結果の統計的な正当性なども計算してくれる.

今回は「なんとなく合っているような気がする」だったが,高機能なライブラリは「統計的な正当性」(例えば検定など)も同時に行ってくれるため,論文のデータなどに使う場合は,そちらを用いる,もしくはNumpyで別に検定を行う,などを行う必要がある.

参考: Python で回帰分析 を行う 6 つ の方法 ~ 回帰分析した結果の t値・p値・(自由度修正済み)決定係数・D.W.値などを出力できる方法 / 出力できない方法 - Qiita