2022-11-13

回帰分析の実践〜 Fast F1を使って線形回帰分析〜

以前は下記の二つの記事で単回帰・重回帰分析の概要をまとめてみました。

takarotoooooo.hatenablog.com
takarotoooooo.hatenablog.com

今回は上記で学んだことから、F1の分析用PythonパッケージであるFast F1を利用して回帰分析をしてみたいと思います。

F1の基礎知識

F1は各チーム自前のF1マシンを使ってサーキット毎に決められた周回数を一番早く走り切る人を決めるスポーツです。
F1決勝の周回中にはタイヤを交換したり、車両通しの接触などのアクシデントでセーフティーカーという先導車のもと、安全が確保されるまでペースを落として周回するといった様々な事象が発生します。

F1は基本的にマシンの事態の速さやドライバー個人のスキルが大事ですが、同じタイヤで走り続けると、1週走るのにかかる時間（ラップタイム）はだんだんと遅くなるので、いつタイヤを交換するか、どのタイヤ（タイヤにも種類がありそれぞれメリットデメリットがあります）に交換するか、交換したタイヤでどのくらい走らせるかなど、チームが考える戦略も勝つために大事な要因の一つであり、チーム力が試されるスポーツです。

今回やってみること

今回はそんなF1のデータを利用して、各状況でドライバーがフルペースで走った場合、どのくらいのラップタイムが出せるのかを回帰を用いて予測してみます。

準備

pipを使ってFast F1インストールします。

$ pip install fastf1

Pythonでモジュールのインポートをすることで利用できます。

import fastf1

今回は2022年のイタリアGPのデータを利用してみます

session = fastf1.get_session(2022, 'Monza', 'R')
session.load()

フェラーリのルクレールというドライバーの周回情報を使ってみましょう。

lec_laps = session.laps.pick_driver('LEC')

データを眺めてみる

fig = plt.figure(figsize=(16, 9))
ax = fig.add_subplot(1, 1, 1)

ax.plot(lec_laps['LapNumber'], lec_laps['LapTimeSec'], label='LEC', color='#da291c')
ax.set_xlabel('LapNumber')
ax.set_ylabel('LapTime')
ax.set_title('LEC lap times at Monza in 2022')
ax.legend()
plt.show()

上記が2022年のイタリアGPのルクレールの周回ごとのラップタイムを可視化したものです。
横軸が周回数。縦軸がその週のラップタイムです。
12, 13週辺りと33, 34週目辺り、48週以降のラップタイムが上がっていることが見えます。
12, 13週辺りと33, 34週目辺りは冒頭で説明したタイヤ履き替えるためにPITに入ったためその分ロスタイムが発生し、ラップタイムが上がっている週です。
48週目以降は他者のアクシデントにより、セーフティーカーが導入され、ペースを落として周回してレースを終了したのでラップタイムが大きくなっています。

では、今回はこのそれぞれの要因でラップタイムが上がっている周回数に、何事もなくレースペースで走っていた場合、どのくらいのラップタイムで走ることができたのかを回帰分析で予測してみます。

単回帰分析で求める

単回帰分析は目的変数と一つの説明変数の相関を導き出す分析手法でした。
今回予測したいのはラップタイムなので目的変数 $y$ をラップタイムとし、説明変数 $x$ は上で可視化した横軸の周回数としてみます。
つまり周回数 $x$ とラップタイム $y$ の回帰式 $y = ax + b$ のパラメータ $a, b$ を求めてみます。

まずは、何事もなくレースペースで走れた周回のラップタイムだけを可視化すると下記になります。

lec_clear_laps = lec_laps.query('TrackStatus == "1" and PitOutTime == "NaT" and PitInTime == "NaT"')

fig = plt.figure(figsize=(16, 9))
ax = fig.add_subplot(1, 1, 1)

ax.plot(lec_clear_laps['LapNumber'], lec_clear_laps['LapTimeSec'], label='LEC', color='#da291c')
ax.set_xlabel('LapNumber')
ax.set_ylabel('LapTime')
ax.set_title('LEC clear lap times at Monza in 2022')
ax.legend()
plt.show()

若干ですが、周回が進むとラップタイムが小さくなっているように見えます。
これを確認していきましょう。

パラメータ $a$ を求める

単回帰分析のパラメータ $a$ は下記の式で求められます。

$\begin{eqnarray} \\ a = \frac{ \sum_{i=1}^{N}{ (x_{i} - \bar{ x })(y_{i} - \bar{ y }) } }{ \sum_{i=1}^{N}{ (x_{i} - \bar{ x })^2 } } \end{eqnarray} \\$

ここで $x_i$ は $i$ 週目の周回数（つまり $i$ です）、 $y_i$ は $i$ 週目のラップタイムを表します。
$\bar{ x }$ は周回数の平均（= 24.09）、 $\bar{ y }$ はラップタイムの平均（= 85.39）です。

x_train = lec_clear_laps['LapNumber']
x_mean = x_train.mean()

y_train = lec_clear_laps['LapTimeSec']
y_mean = y_train.mean()

a = ((x_train - x_mean) * (y_train - y_mean)).sum() / ((x_train - x_mean) * (x_train - x_mean)).sum()

計算をすると $a = -0.032$ ということがわかりました。

パラメータ $b$ を求める

単回帰分析のパラメータ $b$ は $a$ を用いて下記の式で求められます。
$\begin{eqnarray} \\ b = \bar{ y } - a * \bar{ x } \end{eqnarray} \\$
計算をすると $b = 86.179$ ということがわかりました。

可視化

回帰直線引いてみると下記のように、近しい直線が引けているように見えます。

このパラメータ $a, b$ のときの評価関数の値は4.955でした。

周回数の単回帰分析によるの予測ラップタイムは結果は下記のように求めることができました。

周回数	ラップタイム
12	85.78
13	85.75
33	85.10
34	85.06
48	84.61
49	84.57
50	84.54
51	84.51
52	84.48
53	84.44

この分析で明らかになったように、周回数とラップタイムに負の相関があるかというと、答えはYesです。因果関係はありませんが、相関関係はあります。
なぜならば、周回数を重ねるにつれ、車に搭載された燃料が減り、車体重量が軽くなるためです。

以上のように、周回数とラップタイムの関係を求め、特定の周回数のラップタイムを予測することができました。

重回帰分析で求める

先ほどは周回数という一つの説明変数 $x$ を使って、ラップタイムを分析しました。
しかし、当然のように実際にはラップタイムは他の様々な要因によって変化します。
よって次は周回数以外に装着しているタイヤの種類、とタイヤを替えてからの周回数を説明変数として追加し、ラップタイムを予測してみようと思います。

カテゴリ変数の変換

装着しているタイヤの種類はソフト、ミディアム、ハードというように種類が分かれているので、数式で扱えるように今回は下記のように数値に変換します。

ソフト	1
ミディアム	2
ハード	3

スケーリング

今回はどの変数（周回数、装着しているタイヤの種類、タイヤを替えてからの周回数）がどのくらいラップタイムの良し悪しに影響するのかを知るために、各変数のスケールを揃えておきます。

x_train = lec_clear_laps[['LapNumber', 'TyreLife', 'CompoundSeq']]
scaled_x_train = ((x_train - x_train.mean()) / x_train.std()).values

そして下記の通り、バイアス項として最終列に1も追加しておきます。

scaled_x_train = np.insert(scaled_x_train, 3, 1, axis=1)

パラメータ $\boldsymbol{w}$ を求める

重回帰分析のパラメータ $\boldsymbol{w}$ は下記の式で求められます。

$\begin{eqnarray} \\ \boldsymbol{w} = ({\boldsymbol{X}^T}{\boldsymbol{X}})^{-1}{\boldsymbol{X}^T}{\boldsymbol{y}} \end{eqnarray} \\$

xt = scaled_x_train.T
xtx_inv = np.linalg.inv(np.dot(xt, scaled_x_train))
xt_y = np.dot(scaled_x_train.T, y_train)
w = np.dot(xtx_inv, xt_y)

可視化

上記で求められたパラメータによる予測を下記に表してみます。

直感的には単回帰分析よりも、近しい線が引かれるように見えます。

このパラメータ $\boldsymbol{w}$ のときの評価関数の値は2.872だったので、単回帰分析の時よりも近しい線が引けているようです。

単回帰分析の時と同様に予測ラップタイムを求めてみると下記のようになりました。

周回数	ラップタイム	単回帰分析による予測
12	86.01	85.78
13	85.43	85.75
33	85.64	85.10
34	84.61	85.06
48	84.76	84.61
49	84.18	84.57
50	84.19	84.54
51	84.20	84.51
52	84.21	84.48
53	84.22	84.44

パラメータを確認する

では、実際に求めれたパラメータを確認してみます。

$w_0$	-0.518	周回数
$w_1$	0.241	装着しているタイヤの種類
$w_2$	-0.003	タイヤを替えてからの周回数
$w_3$	85.392	バイアス項

周回数のパラメータは先ほどの通り、負の数になっているので、周回を重ねることによってラップタイムが下がることが表せている：○
装着しているタイヤの種類のパラメータは正の数になっているので、ハードタイヤよりソフトタイヤの方がラップタイムが下がることを表している：○
タイヤを替えてからの周回数は、本来は周回を重ねるごとにラップタイムは悪くなるので正の数になるはずですがそうなりませんでした：×

タイヤを替えてからの周回数を示すパラメータの値が、実際の相関とずれてしまったのは多重共線性による問題が関連しています。
タイヤを替えてからの周回数が増えると、それと同時に周回数も増加するので「タイヤを替えてからの周回数」と「周回数」には相関関係があることになります。
重回帰分析は各説明変数が独立であることを前提に行う分析なので、適切な予測ができない回帰式が出来上ってしまうので、どの説明変数を利用するのかはとても大事な要因であることがここでわかりました。

まとめ

今回は実際のデータを用いて線形の単回帰分析、重回帰分析を試してみました。
改めて回帰分析では説明変数と目的変数の相関や目的変数の予測を行うことが可能であることがわかりました。
それと同時に重回帰分析では説明変数同士が独立であることの大事さも示されました。

次は同様にロジスティック回帰に関して実際のデータを使って試してみようと思っています。

2022-11-06

回帰分析の理解まとめ〜ロジスティック回帰分析〜

これまで回帰分析のうち単回帰分析と重回帰分析についてまとめてきました。

takarotoooooo.hatenablog.com
takarotoooooo.hatenablog.com

今回はロジスティック回帰に関してまとめます。

ロジスティック回帰分析とは？

ロジスティック回帰は、事象が発生するか、発生しないかのように二値をとるような問題（コインを投げて表が出るか裏が出るかのような問題です）に関して、発生する（または発生しない）確率を予測するための分析手法です。

ロジスティック回帰ではシグモイド関数というものを利用します。

シグモイド関数

シグモイド関数は下記の式で表される関数で、下記のようなグラフを取ります。

$\begin{eqnarray} \sigma(x) = \frac{1}{1+e^{-x}} \end{eqnarray}$

ご覧の通り、シグモイド関数は、 $x$ の値が大きければ1に近づき、小さければ0に近づくような関数です。

このシグモイド関数に対して、重回帰分析の時と同様に $M$ 個の説明変数 $x_1, x_2, ... , x_M$ とパラメータ $w_1, w_2, ... , w_M$ を用いて、事象の発生確率の予測値 $\hat{y}$ を求めます。
$\begin{eqnarray} \hat{p}(x_1, x_2, ... , x_M) &=& \sigma(w_{0}x_{0} + w_{1}x_{1} + w_{2}x_{2} + ... + w_{M}x_{M}) \\ &=& \frac{1}{1+e^{-(w_{0}x_{0} + w_{1}x_{1} + w_{2}x_{2} + ... + w_{M}x_{M} )}} \end{eqnarray}$

なお、ここで $w_{0}x_{0}$ をバイアス項とし、
$\boldsymbol{x} = \begin{pmatrix} x_0 \\ x_1 \\ \vdots \\ x_M \end{pmatrix}, \boldsymbol{w} = \begin{pmatrix} w_0 \\ w_1 \\ \vdots \\ w_M \end{pmatrix}$
とすると
$\begin{eqnarray} \hat{p}(\boldsymbol{x}) &=& \sigma( \boldsymbol{x}^{T}\boldsymbol{w} ) &=& \frac{1}{1+e^{-( \boldsymbol{x}^{T}\boldsymbol{w} )}} \end{eqnarray}$
と表すことができます。

評価関数

単回帰分析、重回帰分析同様に最適なパラメータ $\boldsymbol{w}$ を求めるために評価関数を用意します。

ロジスティック回帰分析の場合、観測される目的変数 $y$ の値は事象が発生した場合=1, 発生しなかった場合=0のように二値になります。
しかしながらロジスティック回帰分析では、発生する（ $y = 1$ となる）確率を予測するため、二乗誤差がそのままでは利用することができません。

そこで $N$ を観測データの個数として、下記のような式を考えます。
$\begin{eqnarray} \mathcal{L( \boldsymbol{w} )} = \sum_{i=1}^{N}{ -y_{i}\log{ (\hat{p}(\boldsymbol{x}_i)) } - (1 - y_{i})\log{ (1 - \hat{p}(\boldsymbol{x}_i)) } } \end{eqnarray}$

この式のうち前の項である下記は、事象が発生する（ $y = 1$ である）場合の誤差を表しています。
$\begin{eqnarray} {-}y_{i}\log{ (\hat{p}(\boldsymbol{x}_i)) } \end{eqnarray}$
下記のグラフからわかるように、確率の予測値 $\hat{p}(\boldsymbol{x}_i)$ が小さい（= 発生する確率が低いと予測した）場合は誤差が大きくなり、大きい（= 発生する確率が高いと予測した）場合は誤差が小さくなります

一方後の項である下記は、事象が発生しない（ $y = 0$ である）場合の誤差を表しています。
$\begin{eqnarray} {-}(1 - y_{i}) \log{ (1 - \hat{p}(\boldsymbol{x}_i)) } \end{eqnarray}$
こちらの項は逆に、確率の予測値が小さい場合は誤差も小さくなり、大きい場合は誤差も大きくなります。

つまり、事象が発生した目的変数に対する説明変数 $\boldsymbol{x}$ による予測確率 $\hat{p}(\boldsymbol{x})$ はより大きく、逆に事象が発生しなかった目的変数に対する説明変数による予測確率はより小さくなれば、評価関数 $\mathcal{L( \boldsymbol{w} )}$ が最小になるので、それを満たすパラメータ $\boldsymbol{w}$ を求めます。

パラメータを求めるための準備

最適なパラメータ $\boldsymbol{w}$ を求めるために上記の $\mathcal{L}$ の式をまずは変形します。

$\begin{eqnarray} \mathcal{L} &=& \sum_{i=1}^{N}{ -y_{i}\log{\hat{p}(\boldsymbol{x}_i)} - (1-y_{i})\log{(1 - \hat{p}(\boldsymbol{x}_i))}} \\ &=& -\sum_{i=1}^{N}{ y_{i}\log{\hat{p}(\boldsymbol{x}_i)} + \log{(1 - \hat{p}(\boldsymbol{x}_i))} - y_{i}\log{(1 - \hat{p}(\boldsymbol{x}_i))}} \\ &=& -\sum_{i=1}^{N}{ y_{i}\{ \log{\hat{p}(\boldsymbol{x}_i)} - \log{(1 - \hat{p}(\boldsymbol{x}_i))} \} + \log{(1 - \hat{p}(\boldsymbol{x}_i))}} \\ &=& -\sum_{i=1}^{N}{ y_{i}(\log{\frac{\hat{p}(\boldsymbol{x}_i)}{(1 - \hat{p}(\boldsymbol{x}_i))}}) + \log{(1 - \hat{p}(\boldsymbol{x}_i))}} \end{eqnarray}$

ここで下記二つの式を利用することで
$\begin{eqnarray} \hat{p}(\boldsymbol{x}_i) &=& \frac{1}{1+e^{-( \boldsymbol{x}_{i}^{T}\boldsymbol{w} )}} \\ (1 - \hat{p}(\boldsymbol{x}_i)) &=& \frac{ e^{-( \boldsymbol{x}_{i}^{T}\boldsymbol{w} )} }{1+e^{-( \boldsymbol{x}_{i}^{T}\boldsymbol{w} )}} \end{eqnarray}$

下記が成り立ちます。
$\begin{eqnarray} \frac{ \hat{p}(\boldsymbol{x}_i) }{ 1 - \hat{p}(\boldsymbol{x}_i) } &=& \frac{ 1 }{ e^{-( \boldsymbol{x}_{i}^{T}\boldsymbol{w} )} } &=& e^{ \boldsymbol{x}_{i}^{T}\boldsymbol{w} } \end{eqnarray}$
この値をオッズと呼び、両辺対数をとった下記をロジットと呼びます。

$\begin{eqnarray} \log{ \frac{ \hat{p}(\boldsymbol{x}_i) }{ 1 - \hat{p}(\boldsymbol{x}_i) } } &=& \boldsymbol{x}_{i}^{T}\boldsymbol{w} \end{eqnarray}$

このロジットを利用してさらに式を変形します。
$\begin{eqnarray} \mathcal{L} &=& -\sum_{i=1}^{N}{ y_{i}(\log{\frac{\hat{p}(\boldsymbol{x}_i)}{(1 - \hat{p}(\boldsymbol{x}_i))}}) + \log{(1 - \hat{p}(\boldsymbol{x}_i))}} \\ &=& -\sum_{i=1}^{N}{ y_{i}( \boldsymbol{x}_{i}^{T} \boldsymbol{w} ) + \log{(1 - \frac{1}{ 1 + e^{(- \boldsymbol{x}_{i}^{T} \boldsymbol{w} )}})}} \\ &=& -\sum_{i=1}^{N}{ y_{i}( \boldsymbol{x}_{i}^{T} \boldsymbol{w} ) + \log{(\frac{e^{( -\boldsymbol{x}_{i}^{T} \boldsymbol{w} )}}{1 + e^{( -\boldsymbol{x}_{i}^{T} \boldsymbol{w} )}})}} \\ &=& -\sum_{i=1}^{N}{ y_{i}( \boldsymbol{x}_{i}^{T} \boldsymbol{w} ) } + \log{(e^{ ( -\boldsymbol{x}_{i}^{T} \boldsymbol{w} ) })} - \log{(1 + e^{ ( -\boldsymbol{x}_{i}^{T} \boldsymbol{w} ) })} \\ &=& -\sum_{i=1}^{N}{ y_{i}( \boldsymbol{x}_{i}^{T} \boldsymbol{w} ) } + ( -\boldsymbol{x}_{i}^{T} \boldsymbol{w} ) - \log{(1 + e^{ ( -\boldsymbol{x}_{i}^{T} \boldsymbol{w} ) })} \\ &=& -\sum_{i=1}^{N}{ (y_{i} - 1) ( \boldsymbol{x}_{i}^{T} \boldsymbol{w} ) - \log{(1 + e^{ ( -\boldsymbol{x}_{i}^{T} \boldsymbol{w} ) })}} \end{eqnarray}$

次に下記の式を考えていきます
$\begin{eqnarray} \frac{\partial}{\partial\boldsymbol{w}}{\mathcal{L(\boldsymbol{w})}} \end{eqnarray}$

勾配を求める

$j$ 番目のパラメータ $w_j$ における偏微分を考えます。

$\begin{eqnarray} \frac{\partial}{\partial\boldsymbol{w}_{j}}{\mathcal{L(\boldsymbol{w})}} &=& \frac{\partial}{\partial\boldsymbol{w}_{j}}{(-\sum_{i=1}^{N}{ (y_{i} - 1) ( \boldsymbol{x}_{i}^{T} \boldsymbol{w} ) - \log{(1 + e^{ ( -\boldsymbol{x}_{i}^{T} \boldsymbol{w} ) })}})} \\ &=& -\sum_{i=1}^{N}({ (y_i - 1)\boldsymbol{x}_{ij} - \frac{e^{(-\boldsymbol{x}_{i}^{T} \boldsymbol{w})}}{1 + e^{(-\boldsymbol{x}_{i}^{T} \boldsymbol{w})}}\boldsymbol{x}_{ij} }) \\ &=& -\sum_{i=1}^{N}({ (y_i - \frac{1}{1 + e^{(-\boldsymbol{x}_{i}^{T} \boldsymbol{w})}}){\boldsymbol{x}_{ij}} }) \\ &=& -\sum_{i=1}^{N}({ (y_i - \sigma({\boldsymbol{x}_{i}^{T} \boldsymbol{w}})){\boldsymbol{x}_{ij}} }) \\ &=& \sum_{i=1}^{N}({ (\sigma({\boldsymbol{x}_{i}^{T} \boldsymbol{w}}) - y_i){\boldsymbol{x}_{ij}} }) \\ &=& \begin{pmatrix} {\boldsymbol{x}_{1j}} & {\boldsymbol{x}_{2j}} & \dots & {\boldsymbol{x}_{Nj}} \end{pmatrix} \begin{pmatrix} \sigma({\boldsymbol{x}_{1}^{T} \boldsymbol{w}}) - y_1 \\ \sigma({\boldsymbol{x}_{2}^{T} \boldsymbol{w}}) - y_2 \\ \vdots \\ \sigma({\boldsymbol{x}_{N}^{T} \boldsymbol{w}}) - y_N \end{pmatrix} \\ &=& \begin{pmatrix} {\boldsymbol{x}_{1j}} & {\boldsymbol{x}_{2j}} & \dots & {\boldsymbol{x}_{Nj}} \end{pmatrix} \begin{pmatrix} \sigma({\boldsymbol{x}_{1}^{T} \boldsymbol{w}}) \\ \sigma({\boldsymbol{x}_{2}^{T} \boldsymbol{w}}) \\ \vdots \\ \sigma({\boldsymbol{x}_{N}^{T} \boldsymbol{w}}) \end{pmatrix} - \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_N \end{pmatrix} \end{eqnarray}$

$\begin{eqnarray} \boldsymbol{X} = \begin{pmatrix} \boldsymbol{x}_{1}^{T} \\ \boldsymbol{x}_{2}^{T} \\ \vdots \\ \boldsymbol{x}_{N}^{T} \end{pmatrix}, \quad \boldsymbol{y} = \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_N \end{pmatrix} \end{eqnarray}, \quad \boldsymbol{x}^{(j)} = \begin{pmatrix} \boldsymbol{x}_{1j} \\ \boldsymbol{x}_{2j} \\ \vdots \\ \boldsymbol{x}_{Nj} \\ \end{pmatrix}$
とすることで、下記のように表すことができます。
$\begin{eqnarray} \frac{\partial}{\partial\boldsymbol{w}_{j}}{\mathcal{L(\boldsymbol{w})}} &=& {\boldsymbol{x}^{(j)}}^T(\sigma(\boldsymbol{X}\boldsymbol{w}) - \boldsymbol{y}) \end{eqnarray}$

これが勾配の $j$ 番目の要素となるので、勾配をまとめると

$\begin{eqnarray} \frac{\partial}{\partial\boldsymbol{w}}{\mathcal{L(\boldsymbol{w})}} &=& \begin{pmatrix} {\boldsymbol{x}^{(0)}}^T(\sigma(\boldsymbol{X}\boldsymbol{w}) - \boldsymbol{y}) \\ {\boldsymbol{x}^{(1)}}^T(\sigma(\boldsymbol{X}\boldsymbol{w}) - \boldsymbol{y}) \\ {\boldsymbol{x}^{(2)}}^T(\sigma(\boldsymbol{X}\boldsymbol{w}) - \boldsymbol{y}) \\ \vdots \\ {\boldsymbol{x}^{(M)}}^T(\sigma(\boldsymbol{X}\boldsymbol{w}) - \boldsymbol{y}) \\ \end{pmatrix} \\ &=& \begin{pmatrix} {\boldsymbol{x}^{(0)}} & {\boldsymbol{x}^{(1)}} & {\boldsymbol{x}^{(2)}} & \dots & {\boldsymbol{x}^{(M)}} \end{pmatrix}^T (\sigma(\boldsymbol{X}\boldsymbol{w}) - \boldsymbol{y}) \\ &=& \boldsymbol{X}^T (\sigma(\boldsymbol{X}\boldsymbol{w}) - \boldsymbol{y}) \end{eqnarray}$

パラメータを求める

評価関数が最小値をとるパラメータ $\boldsymbol{w}$ は上記の勾配が $\boldsymbol{0}$ となるところになります。
勾配は傾きが最も大きくなる方向を示すベクトルなので、勾配が示す方向に $\boldsymbol{w}$ を変化させていけば、勾配が $\boldsymbol{0}$ の極小値にたどり着くことができます。
勾配を利用して最適なパラメータを求める方式の一つに最急降下法があるので、ここではこの方法で説明します。

最急降下法

パラメータの初期値 $\boldsymbol{w}_0$ を決定した上で、次の手順を繰り返すことで最適なパラメータを求めるのが最急降下法です。

現在のパラメータ $\boldsymbol{w}_{k}$ での勾配を求める
求められた勾配の反対方向にあらかじめ決めた学習率 $\eta$ 分、パラメータを移動させ新たなパラメータ $\boldsymbol{w}_{k+1}$ とする

つまり、下記の式を繰り返し計算することで、最適なパラメータが求められます。
$\begin{eqnarray} \boldsymbol{w}_{k+1} = \boldsymbol{w}_{k} - \eta \boldsymbol{X}^T (\sigma(\boldsymbol{X}\boldsymbol{w_k}) - \boldsymbol{y}) \end{eqnarray}$

上記の手順を繰り返してパラメータを求めたイメージが下記になります。

まとめ

ロジスティック回帰分析はある事象が発生する確率を予測する分析手法
シグモイド関数を利用する
最適なパラメータは $\boldsymbol{w}_{k+1} = \boldsymbol{w}_{k} - \eta \boldsymbol{X}^T (\sigma(\boldsymbol{X}\boldsymbol{w_k}) - \boldsymbol{y})$ を繰り返し解くことで求めることができる

今回はロジスティック回帰分析に関して概要からパラメータの求め方までをまとめてみました。
具体的なデータを使っての求め方は長くなってしまったので、別でまとめてようと思います。

2022-11-03

回帰分析の理解まとめ〜重回帰分析〜

前回こちらの記事で回帰分析のうち単回帰分析について概要とパラメータの求め方をまとめした。
takarotoooooo.hatenablog.com

今回は重回帰分析に関してまとめてみようと思います。

重回帰分析とは？

ある目的変数 $y$ を複数の説明変数 $x_1, x_2, ... , x_M$ で表す式を求めることです。
例えば、目的変数を下記の式で近似できると仮定した時に、最適なパラメータ $w_1, w_2, ... , w_M$ を求めます。

$y = w_{1}x_{1} + w_{2}x_{2} + ... + w_{M}x_{M}$
ここで、 $w_{0}x_{0}$ をバイアス項（全ての $x$ が0の時の目的変数の値）とし、
$\boldsymbol{x} = \begin{pmatrix} x_0 \\ x_1 \\ \vdots \\ x_M \end{pmatrix}, \boldsymbol{w} = \begin{pmatrix} w_0 \\ w_1 \\ \vdots \\ w_M \end{pmatrix}$
とすると
$\begin{eqnarray} y = \boldsymbol{w}^T\boldsymbol{x} \end{eqnarray}$
と表すことができます。

評価関数

重回帰分析においても、評価関数として二乗誤差を使います。
$y_i$ , $\hat{y_i}$ をそれぞれ、 $i$ 番目の実測値、予測値とすると
$\begin{eqnarray} \mathcal{L}(\boldsymbol{w}) &=& \sum_{i=1}^{N}{ ( y_i - \hat{y_i} )^2 } \\ &=& \begin{pmatrix} {{y_1} - \hat{y_1}} & {{y_2} - \hat{y_2}} & \cdots & {{y_N} - \hat{y_N}} \end{pmatrix} \begin{pmatrix} {{y_1} - \hat{y_1}} \\ {{y_2} - \hat{y_2}} \\ \vdots \\ {{y_N} - \hat{y_N}} \end{pmatrix} \end{eqnarray}$

$\boldsymbol{y} = \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_N \end{pmatrix}, \hat{ \boldsymbol{y} } = \begin{pmatrix} \hat{y_1} \\ \hat{y_2} \\ \vdots \\ \hat{y_N} \end{pmatrix}$
とすると

$\begin{eqnarray} \mathcal{L}(\boldsymbol{w}) &=& (\boldsymbol{y} - \hat{\boldsymbol{y}})^T (\boldsymbol{y} - \hat{\boldsymbol{y}}) \end{eqnarray}$

評価関数が最小値をとるパラメータを探す

まずは先ほどおいた $\hat{\boldsymbol{y}}$ を下記のように展開します
$\hat{\boldsymbol{y}} = \begin{pmatrix} { \hat{y_1} } \\ { \hat{y_2} } \\ \vdots \\ { \hat{y_N} } \end{pmatrix} = \begin{pmatrix} { \boldsymbol{x_1}^T{w_1} } \\ { \boldsymbol{x_2}^T{w_2} } \\ \vdots \\ { \boldsymbol{x_N}^T{w_N} } \end{pmatrix} = \begin{pmatrix} { \boldsymbol{x_1}^T } \\ { \boldsymbol{x_2}^T } \\ \vdots \\ { \boldsymbol{x_N}^T } \end{pmatrix} \boldsymbol{w} = \begin{pmatrix} { x_{10} } & { x_{11} } & { x_{12} } & \cdots & { x_{1M} } \\ { x_{20} } & { x_{21} } & { x_{22} } & \cdots & { x_{2M} } \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ { x_{N0} } & { x_{N1} } & { x _{N2} } & \cdots & { x_{NM} } \end{pmatrix} \boldsymbol{w}$

ここで

$\boldsymbol{X} = \begin{pmatrix} { x_{10} } & { x_{11} } & { x_{12} } & \cdots & { x_{1M} } \\ { x_{20} } & { x_{21} } & { x_{22} } & \cdots & { x_{2M} } \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ { x_{N0} } & { x_{N1} } & { x _{N2} } & \cdots & { x_{NM} } \end{pmatrix}$

とすると、下記のように表せます。

$\hat{\boldsymbol{y}} = \boldsymbol{X}\boldsymbol{w}$

これを評価関数に代入し、式を展開すると下記のように表せます。

$\begin{eqnarray} \mathcal{L}(\boldsymbol{w}) &=& (\boldsymbol{y} - \hat{\boldsymbol{y}})^T (\boldsymbol{y} - \hat{\boldsymbol{y}}) \\ &=& (\boldsymbol{y} - \boldsymbol{X}\boldsymbol{w})^T (\boldsymbol{y} - \boldsymbol{X}\boldsymbol{w}) \\ &=& (\boldsymbol{y}^T - \boldsymbol{w}^T\boldsymbol{X}^T) (\boldsymbol{y} - \boldsymbol{X}\boldsymbol{w}) \\ &=& \boldsymbol{y}^T \boldsymbol{y} - \boldsymbol{w}^T\boldsymbol{X}^T \boldsymbol{y} - \boldsymbol{y}^T \boldsymbol{X}\boldsymbol{w} + \boldsymbol{w}^T\boldsymbol{X}^T \boldsymbol{X}\boldsymbol{w} \\ &=& \boldsymbol{y}^T \boldsymbol{y} - 2\boldsymbol{y}^T \boldsymbol{X}\boldsymbol{w} + \boldsymbol{w}^T\boldsymbol{X}^T \boldsymbol{X}\boldsymbol{w} \end{eqnarray}$

ここで評価関数を最小にする $\boldsymbol{w}$ を求めます。

$\begin{eqnarray} \frac{\partial}{\partial{\boldsymbol{w}}}\mathcal{L}(\boldsymbol{w}) &=& \frac{\partial}{\partial{\boldsymbol{w}}}{({\boldsymbol{y}^T}{\boldsymbol{y}} - 2{\boldsymbol{y}^T}{\boldsymbol{X}}{\boldsymbol{w}} + {\boldsymbol{w}^T\boldsymbol{X}^T}{\boldsymbol{X}\boldsymbol{w}})} \\ &=& \frac{\partial}{\partial{\boldsymbol{w}}}{({\boldsymbol{y}^T}{\boldsymbol{y}})} - \frac{\partial}{\partial{\boldsymbol{w}}}{(2{\boldsymbol{y}^T}{\boldsymbol{X}}{\boldsymbol{w}})} + \frac{\partial}{\partial{\boldsymbol{w}}}{({\boldsymbol{w}^T\boldsymbol{X}^T}{\boldsymbol{X}\boldsymbol{w}})} \\ &=& \boldsymbol{0} - 2{\boldsymbol{X}^T}{\boldsymbol{y}} + ({\boldsymbol{X}^T}{\boldsymbol{X}}+({\boldsymbol{X}^T}{\boldsymbol{X}})^T)\boldsymbol{w} \\ &=& - 2{\boldsymbol{X}^T}{\boldsymbol{y}} + 2{\boldsymbol{X}^T}{\boldsymbol{X}}\boldsymbol{w} \end{eqnarray}$

つまり

$\begin{eqnarray} \frac{\partial}{\partial{\boldsymbol{w}}}\mathcal{L}(\boldsymbol{w}) &=& 0 \\ {- 2{ \boldsymbol{X}^T }{ \boldsymbol{y} } } + 2{ \boldsymbol{X}^T }{ \boldsymbol{X} }\boldsymbol{w} &=& 0 \\ { \boldsymbol{X}^T }{ \boldsymbol{X} }\boldsymbol{w} &=& { \boldsymbol{X}^T }{ \boldsymbol{y} } \\ \end{eqnarray}$

つまり、下記が成り立ちます

$\begin{eqnarray} \boldsymbol{w} &=& ({ \boldsymbol{X}^T }{ \boldsymbol{X} })^{-1}{ \boldsymbol{X}^T }{ \boldsymbol{y} } \end{eqnarray}$

実際に求めてみる

パラメータの求め方がわかったので、実際に以下のサンプルデータを使って求めてみます。

$x_1$	$x_2$	$x_3$	$x_4$	$y$
0	0	0	0	-5.092
1	1	9	1	55.14354986
・・・	・・・	・・・	・・・	・・・
48	3	2	0	89.69
49	4	1	1	101.47423445

このデータを横軸を $x_1$ 、縦軸を $x_2$ としてプロットしたものが下記です。

このデータをもとに $\boldsymbol{w}$ を下記の式から求めると
$\begin{eqnarray} ({ \boldsymbol{X}^T }{ \boldsymbol{X} })^{-1}{ \boldsymbol{X}^T }{ \boldsymbol{y} } \end{eqnarray}$

$w_0$	$w_1$	$w_2$	$w_3$	$w_4$
-0.78	2.069	-4.28	5.41	11.69

と、なり
$\hat{\boldsymbol{y}} = \boldsymbol{X}\boldsymbol{w}$
を計算して予測データとしてプロットしたものが下記となり、おおよそ観測データを近似できているように見えます。

$\boldsymbol{w}$ が示す意味

$\boldsymbol{w}$ のそれぞれの値 $w_1, ... , w_M$ は単回帰分析の時と同様に $x_1, x_2, ... , x_M$ が $y$ に与える影響度の大きさを表しています。
ただし、 $x_1, x_2, ... , x_M$ それぞれの値のスケールが違うことで、 $w_1 < w_2$ の時に $x_1$ よりも $x_2$ の方が $y$ により大きな影響を与えるということが言えません。

パラメータ $\boldsymbol{w}$ の値を見るだけで、どの変数 $x_i$ が $y$ に対して、どの程度の影響があるのかを知れるようにするために、各 $x_i$ の値を標準化した上でパラメータ $\boldsymbol{w}$ を求めます。

標準化

各変数の標準化の方法は変数が従う確率分布によって異なりますが、ここでは各変数が正規分布に従うと仮定して標準化してみます。
各 $x_i$ に関して平均を引き、標準偏差で割ることで平均0、分散1の値として標準化することができます。

$x_1$	$x_2$	$x_3$	$x_4$
-1.69	-1.41	-1.56	-1
-1.62	-0.70	1.56	1
・・・	・・・	・・・	・・・
1.62	0.70	-0.87	-1
1.69	1.41	-1.21	1

この標準化された $\boldsymbol{X}_{scaling}$ をもとに算出されたパラメータ $\boldsymbol{w}_{scaling}$ が下記です。

$w_0$	$w_1$	$w_2$	$w_3$	$w_4$
71.54	29.85	-6.05	15.54	5.84

これを見ると $w_1$ の絶対値が大きいので、 $x_1$ が $y$ への影響度が一番大きいことが判断できます。

まとめ

重回帰分析は、ある目的変数 $y$ を複数の説明変数 $x_1, x_2, ... , x_M$ で表す式を求めること
最適なパラメータ $w$ は行列を用いて $({ \boldsymbol{X}^T }{ \boldsymbol{X} })^{-1}{ \boldsymbol{X}^T }{ \boldsymbol{y} }$ で求められる
パラメータ $w$ は説明変数を標準化することで、説明変数毎の目的変数への影響度を表すことができる。

今回は重回帰分析に関して概要からパラメータの求め方までをまとめてみました。
次回はロジスティック回帰分析に関してまとめてみようと思います。

2022-10-27

回帰分析の理解まとめ〜単回帰分析〜

回帰分析とは？

ある目的変数 $y$ と説明変数 $x$ の関係性を明らかにする統計手法です。
単回帰、重回帰、ロジスティック回帰、主成分回帰、部分的最小二乗回帰などさまざまな種類があります。
この記事では単回帰分析に関してまとめてみます。

単回帰分析とは？

ある目的変数 $y$ を一つの説明変数 $x$ で表す式を求めることです。
例えば、身長と体重の関係。年齢と年収の関係などが、下記の式で近似できると仮定した時に、最適なパラメータ $a$ , $b$ を求めます。
$\begin{eqnarray} \\ y = ax + b \end{eqnarray} \\$

最適なパラメータとは？

「最適」とはどんなものになるでしょうか？

下のようなデータを考えてみます。

このデータに対して
$y = ax + b$ における $a$ , $b$ をそれぞれ変えると下記の通りさまざまな直線を書くことができます。

直感的にどの $a$ , $b$ による直線が一番観測されたデータと近いでしょうか？

それぞれ直線がどのくらい観測データに近しいかを表す関数を評価関数と呼びます。

評価関数

評価関数として二乗誤差というものがよく使われます。
$i$ 番目の観測値 $x_i$ に対する $y_i$ と、あるパラメータ $a$ , $b$ によって表される推定値 $\hat{y}_i ( = ax_i + b)$ の差分 $y_i - \hat{y}_i$ の二乗を全ての $x$ において計算し総和を求めます。

$N$ を観測されたデータの個数とすると下記の式で表されます

$\begin{eqnarray} \\ \mathcal{L}(a, b) &=& \sum_{i=1}^{N}{ (y_{i} - \hat{y_{i}})^2 } \\ &=& \sum_{i=1}^{N}{ (y_{i} - (ax_{i} + b))^2 } \\ \end{eqnarray}$

この二乗誤差の総和を評価関数とすると、評価関数の値がより小さければ観測データをよく表しているように見えます。
よって評価関数の値をより小さくするようなパラメータ $a$ , $b$ を見つけてみます。

評価関数が最小値をとるパラメータを探す

評価関数の値が最小になる点は微分を用いることで、見つけることができます。

今回は評価関数を二乗誤差の総和としているので、下記を満たすパラメータ $a$ , $b$ を求めれば良いことになります。
$\begin{eqnarray} \\ \frac{\partial}{\partial a}{ \sum_{i=1}^{N}{ (y_{i} - (ax_{i} + b))^2 } } = 0 \\ \frac{\partial}{\partial b}{ \sum_{i=1}^{N}{ (y_{i} - (ax_{i} + b))^2 } } = 0 \\ \end{eqnarray}$

パラメータ $b$ の導出

$\begin{eqnarray} \\ \frac{\partial}{\partial a}{ \sum_{i=1}^{N}{ (y_{i} - (ax_{i} + b))^2 } } &=& -2\sum_{i=1}^{N}{ (y_{i} - (ax_{i} + b)) } \\ &=& 2\sum_{i=1}^{N}{ (b) } -2\sum_{i=1}^{N}{ (y_{i} - ax_{i}) } \\ &=& 2Nb -2\sum_{i=1}^{N}{ (y_{i} - ax_{i}) } \\ \end{eqnarray}$
つまり
$\begin{eqnarray} \\ 2Nb -2\sum_{i=1}^{N}{ (y_{i} - ax_{i}) } = 0 \\ => b &=& \frac{ \sum_{i=1}^{N}{ (y_{i} - ax_{i}) } }{ N } \end{eqnarray}$

ここで、 $\frac{ \sum_{i=1}^{N}{ y_i } }{ N }$ と $\frac{ \sum_{i=1}^{N}{ x_i } }{ N }$ はそれぞれ観測データの $y$ と $x$ の平均となるので、それぞれを $\bar{ y }$ , $\bar{ x }$ と表すと
$\begin{eqnarray} \\ b &=& \bar{ y } - a\bar{ x } \end{eqnarray} \\$

パラメータ $a$ の導出

$\begin{eqnarray} \\ \frac{\partial}{\partial a}{ \sum_{i=1}^{N}{ (y_{i} - (ax_{i} + b))^2 } } &=& \frac{\partial}{\partial a}{ \sum_{i=1}^{N}{ (y_{i} - (ax_{i} + \bar{ y } - a\bar{ x }))^2 } } \\ &=& \frac{\partial}{\partial a}{ \sum_{i=1}^{N}{ (y_{i} - \bar{ y } - a(x_{i} - \bar{ x }))^2 } } \\ &=& \sum_{i=1}^{N}{ (x_{i} - \bar{ x })(y_{i} - \bar{ y } - a(x_{i} - \bar{ x })) } \\ &=& \sum_{i=1}^{N}{ (x_{i} - \bar{ x })(y_{i} - \bar{ y }) } - a \sum_{i=1}^{N}{ (x_{i} - \bar{ x }))^2 } \end{eqnarray}$
つまり
$\begin{eqnarray} \\ \sum_{i=1}^{N}{ (x_{i} - \bar{ x })(y_{i} - \bar{ y }) } - a \sum_{i=1}^{N}{ (x_{i} - \bar{ x }))^2 } = 0 \\ => a = \frac{ \sum_{i=1}^{N}{ (x_{i} - \bar{ x })(y_{i} - \bar{ y }) } }{ \sum_{i=1}^{N}{ (x_{i} - \bar{ x })^2 } } \end{eqnarray}$
と $a$ を観測データのみで求めることができることがわかります。
なお、 $y_{i} - \bar{ y }$ , $x_{i} - \bar{ x }$ は、観測データを中心化したデータであり、 $\tilde{ y_i } = y_{i} - \bar{ y }$ , $\tilde{ x_i } = x_{i} - \bar{ x }$ と表すと下記の関係になります。

実際に求めてみる

パラメータの算出方法がわかったところで、実際のデータを使って求めてみます。

以下がサンプルデータの一部です。

$x$	$y$	$\tilde{x}$	$\tilde{y}$
0	45.010787	-24.5	-102.428666
1	62.315182	-23.5	-85.124271
2	66.829897	-22.5	-80.609556
・・・	・・・	・・・	・・・
47	245.844440	22.5	98.404987
48	238.553507	23.5	91.114054
49	246.184165	24.5	98.744713

$\begin{eqnarray} \\ a &=& \frac{ (-24.5 * -102.428666) + ... + (24.5 * 98.744713) }{ (-24.5)^2 + ... + (24.5)^2 } \\ &=& 3.8460887159267987 \end{eqnarray}$
なお、
$\bar{x} = 24.5, \bar{y} = 147.4394527007682$
なので、
$\begin{eqnarray} \\ b &=& \bar{ y } - a\bar{ x } \\ &=& 53.21027916056161 \end{eqnarray} \\$

これらの $a$ , $b$ を用いた直線を表してみると下記になります。

実際に観測されたデータに近しい直線で引けており、 $x$ , $y$ の表すことができました。

まとめ

単回帰分析は目的変数 $y$ を一つの説明変数 $x$ で表す式を求めること
最適な直線は観測データと直線の差を表す評価関数を最小にするようなパラメータを見つけることで求められる

今回は回帰分析のうち、単回帰分析に関して概要からパラメータの求め方までをまとめてみました。
次回は重回帰分析に関してまとめてみようと思います。

F1の基礎知識

今回やってみること

準備

データを眺めてみる

単回帰分析で求める

パラメータ を求める

パラメータ を求める

可視化

重回帰分析で求める

カテゴリ変数の変換

スケーリング

パラメータ を求める

可視化

パラメータを確認する

まとめ

ロジスティック回帰分析とは？

シグモイド関数

評価関数

パラメータを求めるための準備

勾配を求める

パラメータを求める

最急降下法

まとめ

重回帰分析とは？

評価関数

評価関数が最小値をとるパラメータを探す

実際に求めてみる

が示す意味

標準化

まとめ

回帰分析とは？

単回帰分析とは？

最適なパラメータとは？

評価関数

評価関数が最小値をとるパラメータを探す

パラメータ の導出

パラメータ の導出

実際に求めてみる

まとめ

パラメータ $a$ を求める

パラメータ $b$ を求める

パラメータ $\boldsymbol{w}$ を求める

$\boldsymbol{w}$ が示す意味

パラメータ $b$ の導出

パラメータ $a$ の導出