回帰分析の実践〜 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	バイアス項

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

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

まとめ

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

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

F1の基礎知識

今回やってみること

準備

データを眺めてみる

単回帰分析で求める

パラメータ を求める

パラメータ を求める

可視化

重回帰分析で求める

カテゴリ変数の変換

スケーリング

パラメータ を求める

可視化

パラメータを確認する

まとめ

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

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

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