単回帰

1. 単回帰#

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches
import matplotlib.colors
import matplotlib.animation
from IPython.display import HTML
import japanize_matplotlib

1.1. 回帰とは#

\(\def\bm{\boldsymbol}\)「夏になるとかき氷が食べたくなる」「寒くなると温かい鍋を食べたくなる」「雨が降るとお店の客足が鈍る」など、我々の生活や行動は気象と密接に連動している。気象庁は気候リスク管理技術に関する調査（スーパーマーケット及びコンビニエンスストア分野）：販売促進策に2週間先までの気温予測を活用するという資料を公開し、気象情報の活用による効率的な在庫管理や商品陳列を紹介している。

ここでは、新規オープンする売店の店長になったつもりで、アイスクリームの在庫管理について考えてみたい。売店の店長としては、商品の在庫をあまり抱えたくない。在庫は売れれば利益となるが、その売店の資金がモノに変化したことになるため、在庫を抱えている間はキャッシュフローが減少する。アイスクリームには賞味期限がないため、売れるまで辛抱強く待てばよいと思われるかもしれないが、冷凍庫で保管するにも保管スペースとコストがかかる。また、最中アイスなどは製造から時間が経過しすぎると、サクサク感が損なわれてしまう。一方で、暑い日にアイスクリームが売れすぎて、在庫を切らしてしまうと、せっかくの商機を失ってしまうことになる。

そこで、気象とアイスクリームの売れ行きの関係を調べ、アイスクリームの在庫の量を調整することを検討したい。例えば、天気予報で最高気温が摂氏30度と予想された日のアイスクリームの売れ行きを推測できれば、それに合わせてアイスクリームの発注を増やしたり、商品を陳列するスペースを増やすことが可能である。ところが、その売店は新しく商売を始めたため、アイスクリームの売れ行きに関する過去の販売データを持ち合わせていない。

そこで、2011年から2020年の東京における

\(x\): 最高気温（東京都の日最高気温の月平均）
\(y\): アイスクリーム・シャーベットへの支出額（一つの世帯が支出する金額の月平均）

の関係を表すデータを用いて、最高気温からアイスクリーム・シャーベットの支出額（以降「支出額」と略す）を近似的に計算する式を求め、売れ行きの傾向を把握したい。以下は、気象庁と総務省統計局が公開している統計データを加工し、グラフ化したものである。

Show code cell source Hide code cell source

X = np.array([
     9.1, 11.2, 12.3, 18.9, 22.2, 26. , 30.9, 31.2, 28.8, 23. , 18.3,
    11.1,  8.3,  9.1, 12.5, 18.5, 23.6, 24.8, 30.1, 33.1, 29.8, 23. ,
    16.3, 11.2,  9.6, 10.3, 16.4, 19.2, 24.1, 26.5, 31.4, 33.2, 28.8,
    23. , 17.4, 12.1, 10.6,  9.8, 14.5, 19.6, 24.7, 26.9, 30.5, 31.2,
    26.9, 23. , 17.4, 11. , 10.4, 10.4, 15.5, 19.3, 26.4, 26.4, 30.1,
    30.5, 26.4, 22.7, 17.8, 13.4, 10.6, 12.2, 14.9, 20.3, 25.2, 26.3,
    29.7, 31.6, 27.7, 22.6, 15.5, 13.8, 10.8, 12.1, 13.4, 19.9, 25.1,
    26.4, 31.8, 30.4, 26.8, 20.1, 16.6, 11.1,  9.4, 10.1, 16.9, 22.1,
    24.6, 26.6, 32.7, 32.5, 26.6, 23. , 17.7, 12.1, 10.3, 11.6, 15.4,
    19. , 25.3, 25.8, 27.5, 32.8, 29.4, 23.3, 17.7, 12.6, 11.1, 13.3,
    16. , 18.2, 24. , 27.5, 27.7, 34.1, 28.1, 21.4, 18.6, 12.3])

Y = np.array([
    463.,  360.,  380.,  584.,  763.,  886., 1168., 1325.,  847.,
    542.,  441.,  499.,  363.,  327.,  414.,  545.,  726.,  847.,
   1122., 1355.,  916.,  571.,  377.,  465.,  377.,  362.,  518.,
    683.,  838., 1012., 1267., 1464., 1000.,  629.,  448.,  466.,
    404.,  343.,  493.,  575.,  921., 1019., 1149., 1303.,  805.,
    739.,  587.,  561.,  486.,  470.,  564.,  609.,  899.,  946.,
   1295., 1325.,  760.,  667.,  564.,  633.,  478.,  450.,  567.,
    611.,  947.,  962., 1309., 1307.,  930.,  668.,  496.,  650.,
    506.,  423.,  531.,  672.,  871.,  986., 1368., 1319.,  924.,
    716.,  651.,  708.,  609.,  535.,  717.,  890., 1054., 1077.,
   1425., 1378.,  900.,  725.,  554.,  542.,  561.,  459.,  604.,
    745., 1105.,  973., 1263., 1533., 1044.,  821.,  621.,  601.,
    549.,  572.,  711.,  819., 1141., 1350., 1285., 1643., 1133.,
    784.,  682.,  587.])

icecream_data = (X, Y)

fig, ax = plt.subplots(dpi=100)
ax.scatter(X, Y, marker='.')
ax.set_title('最高気温とアイスクリーム・シャーベットの支出額')
ax.set_xlabel('最高気温の月平均（℃）')
ax.set_ylabel('支出額（円）')
ax.set_xlim(0, 35)
ax.set_ylim(-250, 2000)
ax.grid()
plt.show()

../_images/4534b55e8a5e576e9f24efc62be937f0da764b1e04c611a8f0c994a5cef8d978.png

グラフから２次関数や指数関数の形が思い浮かぶが、まずは最も単純な場合として、最高気温（\(x \in \mathbb{R}\)）と支出額（\(y \in \mathbb{R}\)）の関係を１次関数に当てはめる（フィッティングさせる）ことを考える。すなわち、１次関数の一般形、

(1.1)#\[ \begin{align} y = f_{a,b}(x) = ax + b \end{align} \]

を仮定し、\(x\)と\(y\)の関係を最も適切に表現するように、傾き\(a \in \mathbb{R}\)と切片\(b \in \mathbb{R}\)を決定したい。１次関数の一般形は\(a\)と\(b\)の値を変えることで無数の関数を表現できる（下図参照）。１次関数の一般形のように、説明したい対象（この場合はアイスクリーム・シャーベットの支出額）の形や構造を抽象化したものをモデル（model）と呼ぶ。モデルの中で説明したい対象を定量的に表現する変数（１次関数の一般形の例では\(a\)と\(b\)）をパラメータ（parameter）と呼ぶ。データに合うようにモデルのパラメータ（この例では\(a\)と\(b\)）を求めることをパラメータ推定（parameter estimation）と呼ぶ。

../_images/d878e0ebd8b12804f40ee3c8eaaaff875c8d6c95d9e02e0488f161bcf76d5bed.png

ところで、人間は外気温が高くなるとアイスクリームやシャーベットなどの冷たいものを食べたくなる。逆に、人間がアイスクリームやシャーベットを食べるから、外気温が上昇すると考えるのは不自然である。

icecream

ゆえに、最高気温からアイスクリーム・シャーベットの支出額の方向に因果関係があると考える。そこで、最高気温で支出額を説明することを考える。支出額\(y\)のように、説明したい変数は目的変数（target variable）や被説明変数（explained variable）、従属変数（dependent variable）などと呼ばれる。また、最高気温\(x\)のように、目的変数を説明するための変数は説明変数（explanatory variable）や独立変数（independent variable）などと呼ばれる。目的変数が量的変数（quantitative variable）であるとき、すなわち説明変数から目的変数の数や量の大小を表現することを回帰（regression）と呼ぶ。１次関数へのフィッティングのように、説明変数が１つの回帰を単回帰（simple regression）と呼び、単回帰によって求められた直線を回帰直線（regression line）と呼ぶ。説明変数が２つ以上の回帰を重回帰（multiple regression）と呼び、次章以降で取り扱う。

ところで、回帰は目的変数が数値で表される問題であれば、様々な局面で応用できる。

採用予定の従業員の学歴、スキル、経験などから妥当な年俸額を算定する
地価、最寄り駅からの距離、築年数、床面積、間取りなどから妥当な家賃を算定する
最寄り駅、最寄り駅からの距離、周辺の人口などから、出店時の売上高を予想する
職業、年齢、年収、勤続年数、貯蓄額などから与信限度額を算定する
テレビ、ラジオ、インターネットなどに配信した広告に対する支出が売上高に与える影響を検証する
あるエリアにおいて使われていた携帯電話の台数から、そのエリアの滞在人口を推定する

なお、最後に挙げた例では、人間が携帯電話の位置を動かしているので、滞在人口から携帯電話の台数の方向に因果関係があると考えるべきである。ところが、この回帰分析では説明変数が滞在人口（因果関係における原因側）となっている。滞在人口は直接観測する（数える）のが難しいため、基地局やGPSなどで携帯電話の台数を観測することにして、携帯電話の台数を説明変数、滞在人口を目的変数としている。

1.2. データの表現#

先ほどの例は観測点が多いため、まずは簡単な例として以下のグラフで示された関係を考える。これは\(x\)を説明変数、\(y\)を目的変数として、\(xy\)平面上の4点\((1, 3), (3, 6), (6, 5), (8, 7)\)で構成されている。

D = np.array([[1, 3], [3, 6], [6, 5], [8, 7]])

../_images/5b74f2e8d0f842b18fc1c46ae7e15854c487d55428b05473d30278c3570b50ac.png

ここでは、説明変数と目的変数の一つの組を事例（instance）として表現する。１個以上の事例をまとめたものをデータ（data）と呼ぶ。 \(1\)番目の事例を\((x_1, y_1)\)、\(2\)番目の事例を\((x_2, y_2)\)、\(i\)番目の事例を\((x_i, y_i)\)と表すことにすると、\(N\)個の事例からなるデータ\(\mathcal{D}\)は次のように表される。

(1.2)#\[\begin{align} \mathcal{D} = \left\{(x_1, y_1), (x_2, y_2), \dots, (x_N, y_N)\right\} = \left\{(x_i, y_i)\right\}_{i=1}^{N} \end{align}\]

先ほどのグラフの例（\(N=4\)）のデータ\(\mathcal{D}_{s}\)は次のように表現される。

(1.3)#\[\begin{align} \mathcal{D}_s = \left\{(x_i, y_i)\right\}_{i=1}^{4} = \left\{(1, 3), (3, 6), (6, 5), (8, 7)\right\} \end{align}\]

データの説明変数\(x\)と目的変数\(y\)を対応付ける関数\(y=f(x)\)を回帰で求めることは、教師データ（supervision data）（もしくは学習データや訓練データ（training data）とも呼ばれる）を用いて、関数\(y = f(x)\)を教師あり学習（supervised learning）で求めると見なすことができる。関数を学習することで、説明変数と目的変数を対応付けるメカニズムの理解が深まるとともに、訓練データには無かった新しい事例（説明変数の値）に対し、その出力（目的変数の値）を予測することに役立つ。

1.3. 残差#

1次関数のモデル\(f_{a,b}(x) = ax + b\)において、傾き\(a\)と切片\(b\)に何らかの値を設定した状況を考える。与えられた\(x\)に1次関数を適用して計算された値を\(\hat{y} \in \mathbb{R}\)と書くことにする。説明変数\(x_i\)に対して、目的変数の実測値（観測値）\(y_i\)と、1次関数\(f_{a,b}(x_i)\)（モデル）から計算された予測値\(\hat{y}_i\)をハット（^）の有無で区別していることに注意されたい。

(1.4)#\[\begin{align} \hat{y} = f_{a,b}(x) = ax + b \end{align}\]

例えば、\(a = 1, b = 1\)に設定して\(y=x+1\)の回帰直線を考える。データ\(\mathcal{D}_s\)の各点\((x_i, y_i)\)と、回帰で計算された点\((x_i, \hat{y}_i)\)をプロットすると、以下のようになる。

../_images/afd7967732272d9f48be45b7e54ea4089ec6b50c03cf73203e8de79e756defbb.png

さて、設定した\(a\)と\(b\)がどのくらいデータ\(\mathcal{D}\)をうまく反映しているのか測定したい。そこで、データ\(\mathcal{D}\)中に含まれる各点\((x_i, y_i)\)に関して、1次関数\(\hat{y}_i = f_{a,b}(x_i) = ax_i + b\)で近似した場合の実測値からのずれ、すなわち残差（residual）\(\epsilon_i \in \mathbb{R}\)を次のように定義する。

(1.5)#\[\begin{align} \epsilon_i = y_i - \hat{y}_i = y_i - (ax_i + b) \end{align}\]

この残差\(\epsilon_i\)は、プロットされた点から回帰直線まで、\(y\)軸方向に差を取ったものである。例えば、データ\(\mathcal{D}_s\)に対して\(a = 1, b = 1\)としたとき、

(1.6)#\[\begin{align} \epsilon_1 &= 3 - (1 \times 1 + 1) = 1 \\ \epsilon_2 &= 6 - (1 \times 3 + 1) = 2 \\ \epsilon_3 &= 5 - (1 \times 6 + 1) = -2 \\ \epsilon_4 &= 7 - (1 \times 8 + 1) = -2 \end{align}\]

以下に、各事例の残差\(\epsilon_i\)を点線で示した。

../_images/eca50a657eee13f276b2e74400f1817d6ae09b6423a23df62fa29b251776bbdf.png

各事例における残差\(\epsilon_i\)を定義したので、次はデータ\(\mathcal{D}\)全体の残差を定義し、設定した\(a\)と\(b\)がどのくらいデータ\(\mathcal{D}\)をうまく反映しているのか、定量的に表現したい。ここで、各事例の残差の平均\(\frac{1}{N}\sum_{i=1}^N \epsilon_i\)を計算してしまうと、正と負の残差が相殺してしまう。ところが、各事例の残差の絶対値の平均\(\frac{1}{N}\sum_{i=1}^N |\epsilon_i|\)を取ると、数学的な取り扱いが煩雑になる。ここでは、平均二乗残差 (MSR: Mean Square Residual) を用いる。

(1.7)#\[ \begin{align} \hat{L} = \frac{1}{N} \sum_{i=1}^{N} \epsilon_i^2 = \frac{1}{N} \sum_{i=1}^{N} (y_i - \hat{y}_i)^2 = \frac{1}{N} \sum_{i=1}^{N} \left\{y_i - (ax_i + b)\right\}^2 \end{align} \]

ちなみに、先ほどのデータ\(\mathcal{D}_s\)に対する１次関数\(y = x + 1\)の平均二乗残差は、

(1.8)#\[\begin{align} \hat{L} = \frac{1}{4} \sum_{i=1}^{4} \epsilon_i^2 = \frac{1}{4} \left\{1^2 + 2^2 + (-2)^2 + (-2)^2\right\} = \frac{13}{4} \end{align}\]

なお、平均二乗残差は平均二乗誤差（MSE: Mean Square Error）と呼ばれることも多い。

1.4. 平均二乗残差は\(a,b\)の関数#

これまでの議論では、傾き\(a\)と切片\(b\)に具体的な値を与えていたが、\(a\)と\(b\)の値を変えることで異なる回帰直線を描画できる。以下は、\(a\)と\(b\)の値をランダムに試行しながら回帰直線を描画し、その平均二乗残差を計算するアニメーションである。\(a\)と\(b\)の値により回帰直線が決定され、その回帰直線がデータ点に近い所を通るとき、平均二乗残差（MSR）の値が小さくなる。過去に試行した\(a\)と\(b\)のうち、最も平均二乗残差が小さい回帰直線を赤色で示している。

単回帰

目次

1. 単回帰#

1.1. 回帰とは#

1.2. データの表現#

1.3. 残差#

1.4. 平均二乗残差は\(a,b\)の関数#