今回から何回かに渡って機械学習の基本的なモデルをご紹介したいと思います。
これらの内容ついては本やUdemyの講座などで素晴らしい解説がありますので、そちらを参照していただいても良いかと思います。
ただ、本を買ったりUdemyの講座を受ける前に、これらがどんなものかという参考になればと思っています。
ですので、数式を交えながらも、企業分析などを通じて感覚的に理解できるように解説していきたいと思います。
ということで、今回は線形回帰(Linear Regression)モデルを解説したいと思います。
勾配ブースティングなどの機械学習手法やディープラーニングにスポットが当たることが多いですが、実際にはビッグデータではなくほんの数サンプルだけの超スモールデータで予測をするといったことも多くあります。
特に企業分析などはビッグデータと呼べるほどデータがないこともたくさんあります。
そういった場合に、データから答えを導くというよりは、人の感覚もたくさん入れて結果を導く線形回帰モデルは非常に有効です。
非常にシンプルで結果の解釈がしやすく、どの変数の影響がどの程度といったことがわかるので、うまくいっているかどうかの判断も可能です。
最後にはPythonで実装し、Scikit-Learnを使った回帰方法について解説したいと思います。
線形回帰モデル
回帰モデルとは
線形回帰モデルの前に、まず一般的な回帰モデルとはというところを説明します。
例として、トヨタ自動車の自動車販売台数とトヨタ自動車の売上高の関係を考えてみましょう。
実際に、販売台数と売上は以下のような関係になっています。
これではまだわかりにくいかもしれませんが、何となく自動車販売台数が増えると売上高も増えている感じがします。
これは、トヨタ自動車の1年間の自動車販売台数を\(x\)、トヨタ自動車の年次の売上高を\(y\)とすると、\(x\)が増えると\(y\)も増えるということになります。
統計学なのか機械学習の分野なのかによって多少呼び方が異なりますが、\(x\)が原因、\(y\)が結果になっているので、一般的にこの\(x\)を説明変数や原因変数などとと呼びます。
\(y\)は目的変数、被説明変数、結果変数などと呼ばれます。
また、自動車販売台数\(x\)がある値であった場合のトヨタ自動車の売上高\(y\)の平均値を条件付き平均値と呼びます。
この条件付き平均値を\(\hat{y}\)とすると、
$$\hat{y}=f(x)$$
と書くことができます。
\(f\)は任意の関数で、必ずしも線形関数でなくても構いません。
このように説明変数\(x\)の値と目的変数の条件付き平均値\(\hat{y}\)の関係を回帰関係と言います。
そしてこのような回帰関係を求めることを、\(y\)を\(x\)の上に回帰する、\(y\)の\(x\)に対して回帰すると言います。
線形回帰モデル
線形回帰モデルは\(f\)が説明変数\(x\)の線形関数になっている場合のことを言います。
ですので、重要な点として線形回帰では以下の仮定をおいています。
目的変数が説明変数の線形和で表される。
線形和というのは、\(x_1,x_2,\cdots\)を説明変数、\(a_0,a_1,a_2,\cdots\)を係数とすると、
$$a_0 + a_1x_1+a_2x_2+\cdots$$
というような足し算になっているという意味です。
これが当てはまらないデータだと線形回帰モデルの予測は適切な予測になりません。
例えば1変数の場合に以下のように非線形な(直線でない)関係があった場合、線形モデルでモデリングするのは妥当でありません。
こういった場合は、説明変数や目的変数を適切に処理することによって、うまく線形の関係に持っていくことも可能な場合もあります。
単回帰
説明変数が1つの場合には単回帰モデルと呼ばれます。
単回帰は、目的変数を\(y\)、説明変数を\(x_1\)として、以下の式で表されます。
$$y=\beta_0 + \beta_1 x_1$$
\(\beta_0\)、\(\beta_1\)はパラメータですが、一般的に\(\beta_0\)は切片(項)、\(\beta_1\)を係数と呼びます(本サイトでは切片項も含めて係数と呼ぶ場合もありますのでご留意ください)。
先ほどの例だと、トヨタ自動車の販売台数が1台増加すれば、ざっくりその価格分の売り上げが増加するはずなので、係数\(\beta_1\)は平均的なトヨタ車の価格200万円~300万といった水準に近くなると考えられます。
実際にどうなるかは後程見ていきたいと思います。
係数の求め方
では、モデルの考え方、形はわかったとして、どうやって係数\(\beta_0\)、\(\beta_1\)を求めるかを考えていきます。
係数を求めるためには“どのような結果になるような係数が欲しいか”を決める必要があります。
それを目的関数と言い、機械学習では損失関数としても使われます。
目的関数 – 最小二乗法
どのように係数を決めるかという問題でまず考えられるのは、実際の値\(y\)と推定された値\(\hat{y}\)の誤差がもっとも小さくなるように決めるというのです。
ただし、誤差は正だったり負だったりするので、単純に誤差を足すというのは良くありません。
そこで、誤差の絶対値をもっとも小さくすることが考えられますが、絶対値は数学的に扱いにくいので、数学的に取り扱いやすい誤差を二乗したものをもっとも小さくしようという発想になります。
これが線形回帰モデルで恐らくもっとも一般的な目的関数である“二乗誤差の最小化”で、この方法により係数を求めることを最小二乗法と呼びます。
\(i\)番目のインプットデータ\(x^{(i)}\)に対するモデルの予測結果を\(\hat{y}^{(i)}\)とすると、二乗誤差は
$$\left(y^{(i)}-\hat{y}^{(i)}\right)^2$$
と表されます。
この平均を平均二乗誤差(Mean Squared Error; MSE)と呼びます(偏差平方和など色々な呼び方があります)。
$$\text{MSE}=\frac{1}{N}\sum_{i=1}^{N}\left(y^{(i)}-\hat{y}^{(i)}\right)^2$$
最小二乗法は、この平均を最小化するように係数を求める方法をです。
係数の導出
では、上記のMSEを最小化して係数を求めたいと思います。
ここで話を簡単にするため、切片項のない単回帰を考えます。
$$y_i = \beta_1 x_1$$
この場合、MSEは
$$\text{MSE}=\frac{1}{N}\sum_{i=1}^{N}\left(y^{(i)}-\beta_1 x_1^{(i)}\right)^2$$
となります。
このMSEを係数\(\beta_1\)で微分しましょう。
2乗は微分が簡単なので、絶対値を使うより取り扱いやすいということです。
$$\frac{\partial \text{MSE}}{\partial \beta_1} =-\frac{1}{N}\sum_{i=1}^{N}2\left(y^{(i)}-\beta_1 x_1^{(i)}\right)x_1^{(i)}$$
最小値を求めるために、この微分をゼロと置いて、\(\beta_1\)について解きます。
$$\begin{align}
-\frac{1}{N}\sum_{i=1}^{N}2\left(y^{(i)}-\beta_1 x_1^{(i)}\right)x^{(i)}&=-\beta_1\frac{1}{N}\sum_{i=1}^{N} {x_1^{(i)}}^2 +\frac{1}{N}\sum_{i=1}^{N}x_1^{(i)} y^{(i)}\\
&=0
\end{align}$$
なので、\(\beta_1\)について解いてそれを\(\beta_1\)の推定値\(\hat{\beta_1}\)とすると、
$$\hat{\beta}_1=\frac{\sum_{i=1}^{N}x^{(i)}y^{(i)}}{\sum_{i=1}^{N}{x^{(i)}}^2}$$
となります。
\(x_1^{(i)}\)や\(y^{(i)}\)はデータなので、それに上の式にデータを代入することで係数が求まります。
切片項がある場合も同じで、\(\beta_0\)、\(\beta_1\)で微分をして、方程式を解くだけです。
行列で表現
機械学習では行列で表現すると表記上も計算上も便利になることが多いです。
そこで、単回帰モデルを行列で表現してみましょう。
$$y=\beta_0+\beta_1x_1$$
というモデルを行列で表します。
\(\beta_0\)と\(\beta_1\)を連結した行列(ベクトル)を\(\beta\)とします。
$$\beta=\left(\begin{array}{c}\beta_0\\ \beta_1\end{array}\right)$$
そして、\(1\)と\(x_1\)を連結した行列を\({\bf{x}}\)とします。
$${\bf{x}}=\left(\begin{array}{c}1\\ x_1\end{array}\right)$$
すると、
$$y=\beta_0+\beta_1x_1$$
は、
$$y=\left(\begin{array}{c}\beta_0& \beta_1\end{array}\right) \left(\begin{array}{c}1\\ x_1\end{array}\right) =\beta ^T {\bf{x}}$$
と表され、すっきりとします。
これは説明変数が複数の重回帰になるとさらに力を発揮し、いくつ説明変数があっても上の式で記述できます。
このときMSEは
$$\text{MSE}=\frac{1}{N}\sum_{i=1}^{N}(y_i-\hat{\beta}^T {\bf{x}}^{(i)})^2$$
と表されます。
さらにデータ\(x_1^{(i)}\)をすべて縦に並べた行列を\(X\)とし、\(y^{(i)}\)を縦に並べたベクトルを\({\bf{y}}\)とします。
$$\begin{align}
X=\left(\begin{array}{cc}1&x_1^{(1)}\\
\vdots&\vdots\\
1&x_1^{(N)}\\
\end{array}\right)\\
y=\left(\begin{array}{c}y^{(1)}\\ \ldots \\ y^{(N)}\end{array}\right)\\
\end{align}$$
するとMSEは
$$\text{MSE}=\frac{1}{N}({\bf{y}}- X\beta )^T({\bf{y}}-X \beta)$$
と変形できます。
ちなみに、\(X\)は\(N\times 2\)行列、\(\beta\)は\(2\times 1\)行列(ベクトル)なので、\(X\beta\)は\(N\times 1\)になります。
そのため\(\left({\bf{y}}- X\beta\right)^T\)は\(1\times N\)で、\({\bf{y}}- X\beta\)は\(1\times N\)なのでMSEはスカラー(\(1\times 1\))になります。
さきほどと同様にMSEを係数ベクトル\(\beta\)で微分してゼロとおき、\(\beta\)について解くと、
$$\hat{\beta}=\left(X^TX\right)^{-1}X^T{\bf{y}}$$
とすっきりした形で書くことができます。
\(\left(X^TX\right)^{-1}\)は\(X^TX\)の逆行列を表します(ようは逆数みたいなものです)。
重回帰
説明変数が複数ある場合の線形回帰を重回帰モデルと呼びます。
トヨタ自動車の売上の例だと、説明変数が自動車販売台数だけでなく、例えば為替レートも入れる場合です。
この場合、自動車販売台数を\(x_1\)、為替レートを\(x_2\)とすると、
$$y=\beta_0 + \beta_1 x_1 + \beta_2 x_2$$
というモデルで表されます。
この場合も、係数を並べた\(\beta\)と、説明変数に1という定数を追加して\(x_1\)、\(x_2\)を並べた\({\bf{x}}\)を使うと、単回帰のときと同様に
$$y=\left(\begin{array}{c}\beta_0& \beta_1 &\beta_2\end{array}\right) \left(\begin{array}{c}1\\ x_1 \\ x_2\end{array}\right) =\beta ^T {\bf{x}}$$
と表されます。
単回帰のときと同じなので、係数も同じように求めることができます。
説明変数を\(p\)個と以下のように一般化しても
$$y=\beta_0 + \beta_1 x_1 + \beta_2 x_2+\cdots+\beta_p x_p$$
さきほどと同様に行列で表すと同じ結果になります。
Pythonで企業分析を実装
では、トヨタ自動車の自動車販売台数と売上高の関係について見てみましょう。
まずは、データをプロットしてみます。
確かに自動車販売台数が増加すると売上高も増加しているようです。
ちなみに、いきなりモデルを作るのではなく、こういった形でグラフを見て、データを眺める、関係性を見ることは非常に重要です。
Pythonで実装
ではPythonを使って線形回帰を実装してみましょう。
まずnumpyと逆関数を求めるためにlinalgをインポートしましょう。
import numpy as np from numpy import linalg
データは以下のように作成します。
car_sales = np.array([8139584, 8423984, 8334485, 9692948, 10133095, 10167489, 10093781, 10250567, 10441126, 10602559, 10456617, 9919759]) sales = np.array([189509, 189936, 185836, 220641, 256919, 272345, 284031, 275971, 293795, 302256, 299299, 272145]) years = np.arange(2009, 2021)
続いて、係数\(\beta\)を求める関数を作成します。
def linear_regression(x, y): n_samples = len(x) x = x.reshape(n_samples, -1) X = np.concatenate([np.ones((n_samples, 1)), x], axis=1) beta = np.matmul(np.matmul(linalg.inv(np.matmul(X.T, X)), X.T), y) return beta
3行目はベクトルを行列の形に変形しており、3行目は1というベクトルと説明変数\({\bf{x}}\)を連結しています。
そして、5行目で以下の式に従って、行列計算で係数ベクトル\(\beta\)を求めています。
$$\hat{\beta}=\left(X^TX\right)^{-1}X^T{\bf{y}}$$
では、この関数の\(x\)に自動車販売台数を、\(y\)にトヨタの売上高を渡して、係数を求めてみます。
beta = linear_regression(x=car_sales, y=sales)
この結果、係数は\(\hat{\beta}_0=-212907\), \(\hat{\beta}_1=0.048\)となりました。
販売台数の係数が0.048で、売上高が億円単位なので、販売台数が1台増加すると480万円売上が増加するという結果なりました。
1台480万円もすることになってしまいますので、係数を見ると「本当にこの結果があっているか？」と疑問になります。
また、切片項が-212,907なので1台も売れなければ売上高がものすごい額のマイナスになっていまいます。
気になるので、実際にプロットして確認してみましょう。こういった確認は必須です。
予測結果\(\hat{y}\)を求める関数を作成します。
def linear_predict(x, beta): beta = beta.reshape(-1, 1) n_samples = len(x) x = x.reshape(n_samples, -1) X = np.concatenate([np.ones((n_samples, 1)), x], axis=1) return np.matmul(X, beta)
作成した関数に自動車販売台数と係数を渡して、実際に予測します。
prediction = linear_predict(car_sales, beta)
では、これをプロットしてみましょう。
それっぽくなっているのでプログラムは間違っていないようです。
ではどこが良くないのでしょうか？
その辺りを機械学習をするにあたって非常に便利なライブラリであるScikit-learnを使いながら確認しましょう。
Scikit-learnを使う
Scikit-learnのインストールについてはたくさんの解説サイトがありますので、そちらをご参照ください。
まず、sklearn.linear_modelモジュールからLinearRegressionクラスをインポートします。
from sklearn.linear_model import LinearRegression
そして、LenearRegressionをインスタンス化し、あとはfitで説明変数、目的変数を渡すだけです。
linear_model = LinearRegression() x = car_sales.reshape(len(car_sales), -1) linear_model.fit(x, sales)
2行目は説明変数はベクトルではなく行列でないといけないので、\(N\)ベクトルから\(N\times 1\)行列に変換しています。
係数は学習したモデルのintercept_プロパティ、coef_プロパティで確認できます。
以下のコードで見てみると、
print(linear_model.intercept_) print(linear_model.coef_)
さきほどのPythonで自作した結果と一致していると思います。
先ほどと同じ切片項がマイナスになるという結果ですが、ここで販売台数がゼロであれば売上高はほぼゼロと考えられるので、切片項をゼロにするようにしてみましょう。
LinearRegressionモデルはいくつか設定が可能で、切片項を含めない場合はfit_intercept=Falseとすれば良いです。
linear_model = LinearRegression(fit_intercept=False)
この場合、販売台数がゼロの場合、売上もゼロになるような予測になります。
その結果、係数は\(\beta_1=0.026\)なりました。
1台販売台数が増えると250万円売上が増加することになります。
こちらの方が感覚に合っていますね。
このように実務では係数を見て結果の妥当性を判断するというのは非常に重要です。
特に、係数の水準や係数の符号は必ず見ないといけません。
次回解説するp値などの統計量も重要ですが、何より係数の妥当性の方が重要です。
では回帰した結果をプロットして見てみます。
一部大き目に外れているところがありますね。
この要因としては、他にも売上に影響を与える変数があることなどが挙げられます。
自動車を国内だけで販売しているわけではないので、同じ台数を販売しても為替レートにより売上高が増減します。
円高になれば円ベースの売上高は減りますし、円安になれば円ベースの売上高は増えます。
長くなってきたので、このあたりは次回の記事で確認したいと思います。
(補足) Sklearn.linear_model.LinearRegressionの引数について
簡単にLinearRegressionモデルの重要な引数を載せておきます。
fit_intercept
先ほど見たように切片項を使用するかどうかを指定します。
Trueなら使用し、Falseなら使用しません。
デフォルトはTrueです。
normalize
Trueにすると、説明変数から平均を引いて標準偏差で割ることで、平均がゼロ、標準偏差が1にしたもの説明変数とします(標準化と言います)。
重回帰などで説明変数ごとに水準が大きく異なる場合は、使用した方が良いです。
新しいバージョンではなくなるようなので、StandardScalerを使った方が良さそうです。
positive
Trueを指定すると係数を正に限定します。
重回帰だと複数の説明変数が影響しあって、本来は正の影響がある変数の係数が負になってしまうことがあります。
例えば、複数の説明変数を含めたせいで、自動車販売台数の係数が負になってしまうような場合です。
この場合、自動車販売台数が増えたときに売上が減るような予測になるので良くありません。
そこでpositiveをTrueに設定して、そのような場合は係数がゼロになるように(その変数を無視するように)します。
Scikit-LearnのLinearRegressionの公式ドキュメントはこちらです。
『sklearn.linear_model.LinearRegression』
まとめ
今回は機械学習の基礎として線形回帰モデルを見ていきました。
単純なモデルといっても非常に奥が深く、まだまだ深堀することは可能です。
次回は重回帰を行って、係数の有意性を検証したり、モデルの評価をしたりといったことをしたいと思います。
皆さまが機械学習に興味を持ったり、楽しみながら学んでいただければ幸いです。
最後に、この記事を読んでいただいて、数学的なところがわからない、統計の話がわからない、数学・統計からもう少しきちんと理解したいという方向けに参考となる本を載せておきます。
では！
参考本
数学についての人気の本としては以下のような本があります。
(無責任ですみませんが、私自身は読んでいないのであくまでご参考です)
統計の基礎に関しては以下がわかりやすいのでオススメです。
今回の記事もこの本を参考にしています。
あと、この記事で作成しているグラフはPlotlyというPythonのライブラリを使っています。
こちらもオススメですので以下の記事をチェックしてみてください！
コメント