はじめに

今回は回帰分析に関するお勉強と簡単な実装の回です。

回帰分析とは

　回帰分析とはある結果とその結果をもたらす要因との関係を明らかにしようとする分析です。回帰分析には単回帰分析と重回帰分析というものがあります。
　まず、単回帰分析について、例えば年齢と年収の関係を思い浮かべてみましょう。無作為に労働人口からn人の標本を選び出して、横軸に年齡、縦軸に年収を設定した領域に各標本をプロットしていきます。このプロットしたものを散布図と呼びます。さて、年齢と年収の間に綺麗な線形関係が見えるかもしれません。この時、年収という結果を生み出す要因は年齢であると言えます。さて、この関係を回帰分析に当てはめて考えた場合、年収を目的変数（従属変数）、年齢を説明変数（独立変数）と呼びます。では年齢と年収の関係はどのようなものになるのでしょうか？さっきちらっと出ましたね、線形関係です。ではこの直線、回帰直線を式で表しましょう。
$$y = \alpha + \beta{x}-(1)$$
$y$が年収で$x$が年齢です。$\alpha$はいわゆる切片でバイアスとよび、$\beta$は回帰係数とよびます。回帰分析のゴールは、実際の目的変数と説明変数の関係を精度よく表現するような$\alpha$と$\beta$の値を求めることです。$\alpha$と$\beta$さえわかれば未知の説明変数に対する目的変数を予測することが可能です。
重回帰分析とは単回帰分析ではひとつだった説明変数に、さらに説明変数を付け加えて目的変数を表現して、その関係におけるパラメータを求める分析です。追加する説明変数には年齢の他に役職や資格などのスキルなどが挙げられるでしょう。回帰分析とは目的変数と説明変数の間の関係式を求めるものだということがわかりました。では、どのような基準で線形関係式のパラメータを決定するのでしょうか？

最小二乗法

回帰問題では直線の式を求めます。しかし実際の値はその直線上にあるとは限りません。したがって、できるだけ直線上の値と実際の値の誤差が小さくなるようにバイアスと回帰係数を決定する必要があります。ここで登場するのが最小二乗法です。最小二乗法ではこの直線上の値$y_i$と実際の値$Y_i$の誤差を残差と呼び,$Y_i-y_i$で表します。$Y_i$は$y_i$より大きかったり小さかったりする場合があり、単に残差を足しあ合わせて評価してしまうと、正確な評価ができません。そこで残差の平方和（残差を２乗して足し合わせたもの）について評価を行います。これが最小二乗法の考え方です。残差平方和$Q$は以下の式で表します。
$$Q=\sum_{i=1}^{n}(Y_i-y_i)^2=\sum_{i=1}^{n}(Y_i-(\alpha+\beta{x_i}))^2-(2)$$
ではQを最小にする$\alpha$と$\beta$はどのようにして求めるのでしょうか？少し紛らわしいのですが、式（２）においては、$\alpha$と$\beta$を変数として考える必要があります。$\alpha$と$\beta$を変数としてみた場合、$Q$は２変数関数です。$\alpha$と$\beta$に関する残差平方和の偏微分の値が0のとき$Q$は最小値であると言えるでしょう。似たような話をこの記事の勾配法の導入、という部分で展開しているのでよかったらのぞいてみてください。式（２）を少しいじりましょう。
$$Q=\sum_{i=1}^{n}(Y_i-(\alpha+\beta{x_i}))^2=\sum_{i=1}^{n}(Y_i^2-2Y_i\alpha-2Y_i\beta{x_i}+\alpha^2+2\alpha\beta{x_i}+\beta^2x_i^2)-(3)$$
これで偏微分を求めやすくなりました。では実際に計算してみましょう。
$$\frac{\partial Q}{\partial \alpha}=-2\sum_{i=1}^nY_i+2n\alpha+2\beta\sum_{i=1}^nx_i=0-(4)$$
$$\frac{\partial Q}{\partial \beta}=-2\sum_{i=1}^nY_ix_i+2\alpha\sum_{i=1}^nx_i+2\beta\sum_{i=1}^nx_i^2=0-(5)$$
式(4),(5)について連立方程式をとけば残差$Q$を最小にするような$\alpha$,$\beta$を求めることができます。またこの連立方程式を正規方程式と呼びます。この計算はむずかしいことを考えずに、連立方程式の基本に立ち返って考えてください。まずは$\beta$から求めます。式(4),(5)より、
$$\beta=\frac{\sum_{i=1}^nx_i\sum_{i=1}^nY_i-n\sum_{i=1}^nx_iY_i}{(\sum_{i=1}^nx_i)^2-n\sum_{i=1}^nx_i^2}=
\frac{\sum_{i=1}^nx_iY_i-\frac{1}{n}\sum_{i=1}^nx_i\sum_{i=1}^nY_i}{\sum_{i=1}^nx_i^2-\frac{1}{n}(\sum_{i=1}^nx_i)^2}-(6)$$
一旦ここで式変形を止めておきます。ここから先に進むには分散について説明しなくてはなりません。

分散について

分散とは、無作為に選んだ標本の各値$x_i$からの標本の平均値$\bar{x}$を引き、それらの値を２乗して全て足し合わせたものの平均と定義されています。式で表すと以下の通りです。
$$s^2=\frac{1}{n}\sum_{i=1}^n(x_i-\bar{x})^2-(7)$$
ここで$s^2$が分散です。この分散は標本集団に関する分散です。ではこの値を元に母集団の分散についても考えてみましょう。標本分散に$\frac{n}{n-1}$をかけたものが母集団の分散です。詳しくは次回以降に紹介します。
母分散は以下の式で表すことができます。
$$\sigma^2=\frac{1}{n-1}\sum_{i=1}^{n}(x_i-\bar{x})^2-(8)$$
ここで$\sigma^2$は不偏分散と呼びます。標本分散はあくまでその標本における分散なので母集団の推定に関しては不偏分散を用います。あとで使いやすいように式(8)を以下のように変形しておきます。
$$\sigma^2=\frac{1}{n-1}\sum_{i=1}^{n}(x_i^2-2x_i\bar{x}+\bar{x}^2)=\frac{1}{n-1}(\sum_{i=1}^{n}x_i^2-2\bar{x}\sum_{i=1}^{n}x_i+n\bar{x}^2)=\frac{1}{n-1}(\sum_{i=1}^{n}x_i^2-\frac{1}{n}(\sum_{i=1}^{n}x_i)^2)-(9)$$
次に共分散というものについて考えます。年齢と年収の例を思い出しましょう。年齢と年収の組、$(x,y)$について、n組のデータが観測できました。ここで年齢$x$と年収$y$の平均をそれぞれ$\bar{x}$,$\bar{y}$と表します。共分散$s_{xy}$は以下の式で表せます。
$$s_{xy}=\frac{1}{n-1}\sum_{i=1}^{n}(x_i-\bar{x})(y_i-\bar{y})=\frac{1}{n-1}(\sum_{i=1}^{n}x_iy_i-n\bar{x}\bar{y})=\frac{1}{n-1}(\sum_{i=1}^{n}x_iy_i-\frac{1}{n}(\sum_{i=1}^{n}x_i)(\sum_{i=1}^{n}y_i))-(10)$$
本題である式(6)に戻りましょう。
式(9)(10)を式（６）に代入します。すると母集団における回帰係数$\beta'$は以下の式で表すことができます。
$$\beta'=\frac{s_{xy}}{s_x^2}-(11)$$また、式(4)より、母集団の回帰直線におけるバイアス$\alpha'$は、$$\alpha'=\bar{Y_i}-\beta'\bar{x_i}-(12)$$となります。
さて、母集団中の標本${Y'}$について標本の説明変数、回帰係数、そしてバイアスを用いて表現してみます。式(1)に$Y=Y'$,$\alpha=\alpha'$,$\beta=\beta'$を代入してみると、$$Y'=\alpha'+\beta'x-(13)$$となります。ここに式(12)を代入します。
すると未知の説明変数$x$に関する目的関数は$$Y'=\beta'(x-\bar{x})-(14)$$として表現できます。ここで(11)を合わせてみてみると、母集団中の標本${Y'}$は標本中の値のみで表現できることがわかります。つまり母集団全てに対して調査を行わずとも任意の標本さえあれば未知の値を予測できてしまうのです。
それではこの考えに基づいて、一定期間の株の値動きから、今後の値動きを予想してみることにします。

株価予想プログラムの実装

それでは早速、回帰分析を用いて株価を予想してみます。今回説明変数にはある日の株価の終値、目的変数に30日後の終値を設定します。回帰直線を求めるために今までたくさんの数式が登場しました。しかしPythonにはsklearnという非常に便利なモジュールが備わっています。標本の値を代入するだけで回帰直線を導出してくれるのです。では導入してみましょう。

import pandas_datareader.data as web
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import sklearn
from sklearn import preprocessing, model_selection, linear_model
import datetime

class ComtureKabu:
    #株価データの取得、データ整形
    def __init__(self):
        self.df_comture = web.DataReader('3844.T', 'yahoo', '2018-10-01')
        self.df_comture['SMA'] = self.df_comture['Close'].rolling(window=14).mean()
        self.df_comture['label'] = self.df_comture['Close'].shift(-30)
        x = np.array(self.df_comture.drop(['label', 'SMA'], axis=1))
        x = sklearn.preprocessing.scale(x)
        self.predict_data = x[-30:]
        x = x[:-30]
        y = np.array(self.df_comture['label'])
        y = y[:-30]
        self.x_train, self.x_test, self.y_train, self.y_test = sklearn.model_selection.train_test_split(x, y, test_size = 0.2,)
    #学習モデル適用
    def learn(self, x_train, y_train):
        lr = sklearn.linear_model.LinearRegression()
        models = lr.fit(x_train, y_train)
        return lr
    #テストデータに対するモデル適合度出力
    def accuracy(self):
        x_test, y_test = self.x_test, self.y_test
        lr = self.learn(x_test, y_test)
        accuracy = lr.score(x_test, y_test)
        return accuracy
    #得られたモデルを元に30日後までの株価予想
    def predict(self):
        df_comture = self.df_comture
        x_test, y_test = self.x_test, self.y_test
        lr = self.learn(x_test, y_test)
        predict_data = self.predict_data
        predicted_data = lr.predict(predict_data)
        df_comture['Predicted'] = np.nan

        last_data = df_comture.iloc[-1].name
        oneday = 86400
        next_unix = last_data.timestamp() + oneday
        for data in predicted_data:
            next_data = datetime.datetime.fromtimestamp(next_unix)
            next_unix += oneday
            df_comture.loc[next_data] = np.append([np.nan] * (len(df_comture.columns)-1), data)
        df_comture = df_comture[-45:]
        df_comture['Close'].plot(figsize=(15, 6), color='green')
        df_comture['Predicted'].plot(figsize=(15, 6), color='orange')
        plt.show()

df_comture = ComtureKabu()
accuracy = df_comture.accuracy()
df_comture.predict()
print(accuracy)

0.7521657632651639

黄色で示した部分が予想した株価のグラフです。
テストデータに対する予想精度はおよそ0.75とあまり良い数字ではありませんね。これに関しては株の勉強をすることで説明変数にどのような値を用いるべきかもっと考える必要がありそうです。あとは学習モデルに線形回帰以外のものを用いたものと比較してより良いモデルを選択すべきですね。

おわりに

他にも様々な学習モデルがあるのでそれについてもいつか紹介します。