2.2.1. 線形回帰の基礎

2020/3/25公開

導入

ここではまず簡単な例として1次元の線形回帰について解説する。扱う問題自体は中学校で習うような簡単なものだが、モデリングや数式の表記は手加減せずに論文で使われているレベルのものを用いるので、機械学習の難しい数式に慣れることに集中できるだろう。

流れとしては直線$y=ax$による回帰を生成モデルによって捉え直し、最尤推定を行うことで最小二乗法を導出する。今後の議論の基礎になるので、ぜひ時間をかけて取り組んでみてほしい。

問題

中学校でオームの法則を習ったことを覚えているだろうか。物体の電気抵抗$R$と物体にかける電圧$V$、流れる電流$I$の間には

$$V=RI \tag{2.2.1.1}$$

の関係が成り立つ。この法則を利用すれば物体の電気抵抗を測定でき、物体に所望の電圧をかけたり、電流を流したりできるようになる。

より一般化して$(2.2.1.1)$式を直線の方程式で書き直せば

$$y=ax \tag{2.2.1.2}$$

となる。電圧、電気抵抗、電流をそれぞれ$y,a,x$で置いただけである。機械学習の文化では$x$を入力(input)、$y$を出力(output)またはラベル(label)という。また$a$は重み(weight)と呼ばれ、本来なら$w$で表記することが多い。

$(2.2.1.2)$式のような関係性を持つことがわかっている、またはそのような関係性を持ちそうだと思われる現象は、重み$a$さえ求めれば入力から出力を予測できる。いまの状況ではどれくらいの電流を流したときにどれくらいの電圧がかかっているかが推測できる。

中学校や高校の範囲では、適当に1点$(x,y)$を測定してそれを代表とし、

$$a = \frac{y}{x} \tag{2.2.1.3}$$

とすれば求まるということになっていた。しかし電気抵抗を実際に何回か測定してみると、どの点を代表とするかによって求まる$a$の値は若干変わり、毎回同じ値になるとは限らない。これは物体の電気抵抗が温度に依存することや電気回路の接触などに起因するが、何が原因かはよくわからないことがほとんどだろう。

こういった「よくわからない観測値の揺らぎ」の扱いを得意とするのが確率論である。中学高校の範囲でならいくつかの観測点について求めた$a$を平均する、大学の学部教養レベルであれば最小二乗法を用いて$a$を求める、といった方法がある。

どちらも十分に実用的な方法であり、実際それで片付いてしまうことも多いが、いまは機械学習への導入に適した形で問題を解き直そう。

モデリング

出力$y$と入力$x$の間には以下の関係があると仮定する。

$$y = a x + \varepsilon \tag{2.2.1.4}$$

この式は「$y$は基本的に$ax$で求まるが、よくわからない誤差$\varepsilon$を含む」という意味であり、理論値からのズレを$\varepsilon$という補正項に吸収させたことになる。ここで$\varepsilon$は平均$0$、分散$\sigma ^ 2$のガウス分布に従うものとする。

$$\varepsilon \sim \mathcal{N} (\varepsilon | 0, \sigma ^ 2) \tag{2.2.1.5}$$

この段階はモデリング(modeling)といって、現象に対して人間が勝手に仮定を置くことであり、誤差が本当にガウス分布に従うかどうかはわからない。しかし人間は現象に何かしら仮定を置かなければ対象を数学的に扱うことはできないし、また仮定の置き方に自由度があることで、人間やコンピュータが計算しやすい形の数式が導出されることもある。仮定がどれだけ正しそうかは、得られた結果を見てから考えればよく、もし満足いかない結果なら満足いくまで仮定を変更すればよい。

ガウス分布(Gaussian distribution)または正規分布(normal distribution)は以下の式で定義される連続確率分布である。

$$\mathcal{N} (x | \mu, \sigma ^ 2) = \frac{1}{\sqrt{2 \pi \sigma ^ 2}} \exp \left(- \frac{(x - \mu) ^ 2}{2 \sigma ^ 2} \right) \tag{2.2.1.6}$$

出力$y$は$ax$に$\varepsilon$を足したものなので、平均$ax$、分散$\sigma ^ 2$のガウス分布にしたがうものとみなせるから、その確率分布は

$$p(y | x, a) = \mathcal{N} (y | ax, \sigma ^ 2) \tag{2.2.1.7}$$

と表せる。最初なので一応丁寧に書いておくと、

$$p(y | x , a) = \frac{1}{\sqrt{2 \pi \sigma ^ 2}} \exp \left(- \frac{(y - ax) ^ 2}{2 \sigma ^ 2} \right) \tag{2.2.1.8}$$

である。パラメータ$\sigma ^ 2$は今後の計算にはあまり関わってこないので、議論の簡略化のため既知の定数とした。ちなみに、今回のように出力$y$がある確率分布から生成されるモデルを生成モデル（generative model）という。

観測した$N$個のデータ点の集合を

$$\mathcal{D} = \left\{ (x ^ {(1)}, y ^ {(1)}), (x ^ {(2)}, y ^ {(2)}), \ldots, (x ^ {(N)}, y ^ {(N)}) \right \} \tag{2.2.1.9}$$

とおく。$\mathcal{D} ^ {(i)} = (x ^ {(i)}, y ^ {(i)})$は$i$番目に観測されたデータ点を表す。また出力と入力を分けるときは

$$\begin{aligned} \mathcal{D}_y = \left\{ y ^ {(1)}, y ^ {(2)}, \ldots, y ^ {(N)} \right \} \\ \mathcal{D}_x = \left\{ x ^ {(1)}, x ^ {(2)}, \ldots, x ^ {(N)}\right \} \end{aligned} \tag{2.2.1.10}$$

と書くことにする。

いま、簡単のために1番目のデータ$\mathcal{D} ^ {(1)} = (x ^ {(1)}, y ^ {(1)})$のみが観測されたとして$a$を求めてみよう。$(2.2.1.8)$式では$y$が確率変数、$x$は確率変数だが観測されて（今回は人間が勝手に定めたことで）実現値となっており、$a$はとりあえず与えられた固定値とみなしていた。しかしいまは$a$を求めたいので$a$が変数であり、また$y$は観測されたことで実現値（固定値）となっている。そこで以下の尤度関数（likelihood function）を定義する。

$$\begin{aligned} \mathcal{L}(a|\mathcal{D} ^ {(1)}) &= \mathcal{L}(a | (x ^ {(1)}, y ^ {(1)})) \\ &= \frac{1}{\sqrt{2 \pi \sigma ^ 2}} \exp \left(- \frac{(y ^ {(1)} - ax ^ {(1)}) ^ 2}{2 \sigma ^ 2} \right) \end{aligned} \tag{2.2.1.11}$$

尤度関数とは確率分布のパラメータを変数とみなし、確率変数を定数とみなした関数である。「尤度」は「ゆうど」と読み「尤もらしい」と書いて「もっともらしい」と読むので、尤度は「もっともらしさ」という意味である。

尤度関数は確率分布から生成したデータのもっともらしさを表す指標のひとつである。尤度関数を最大化する$a$を求める、つまり最も尤もらしい$a$の値を求める手法を最尤推定(さいゆうすいてい：maximum likelihood estimation)という。

$(2.2.1.10)$式が最大となるのは$\exp$の中が 0 になるときであるから、$y ^ {(1)} - a x ^ {(1)} = 0$、すなわち

$$a = \frac{y ^ {(1)}}{x ^ {(1)}} \tag{2.2.1.11}$$

のときである。これは原点と観測した1点を正確に通る直線を引くときの$(2.2.1.3)$式と一致している。

注意点は尤度関数は確率分布ではないということである。つまり一般には

$$\mathcal{L}(a|\mathcal{D} ^ {(1)}) \neq p(a | \mathcal{D} ^ {(1)}) \tag{2.2.1.12}$$

である。なぜならばいま、そもそも$a$は確率変数ではないとみなしているし、仮に確率変数とみなしても、ベイズの定理により

$$p(a | \mathcal{D} ^ {(1)}) = p(a| y , x) = \frac{p(y|x,a)p(a)}{p(y)} \tag{2.2.1.13}$$

となるので、$(2.2.1.12)$式の両辺が形式的にでも等号で結ばれるのは$p(y) = p(a) = {\rm const.}$という特殊な状況だけだからである。

では$N$個のデータ点を観測した場合に話を戻そう。各データ点の観測が独立で、同じ確率分布からのサンプリングであると仮定する。このとき$N$個のデータの確率分布は

$$\begin{aligned} p(\mathcal{D} _ y | \mathcal{D} _ x, a) &= p(y ^ {(1)}| x ^ {(1)}, a) \times \cdots \times p(y^{(N)}|x^{(N)}, a) \\ &= \prod _ {i=1} ^ N p(y ^ {(i)}| x ^ {(i)}, a) \\ &= \prod _ {i=1} ^ N \frac{1}{\sqrt{2 \pi \sigma ^ 2}} \exp \left(- \frac{(y ^ {(i)} - ax ^ {(i)}) ^ 2}{2 \sigma ^ 2} \right) \end{aligned} \tag{2.2.1.14}$$

と表せる。具体的な状況としては、入力$x ^ {(1)}, x^{(2)}, \ldots, x ^ {(N)}$のときの出力を観測する、と決めたがまだ実際には出力$y ^ {(1)}, y ^ {(2)}, \ldots, y ^ {(N)}$の観測を行っていない状態（まだ確率変数）である。

出力$y ^ {(1)}, y ^ {(2)}, \ldots, y ^ {(N)}$を観測する（実現値を代入して定数とみなす）ことで再び尤度関数が定義できる。すなわち

$$\mathcal{L}(a|\mathcal{D}) = \prod _ {i=1} ^ N \frac{1}{\sqrt{2 \pi \sigma ^ 2}} \exp \left(- \frac{(y ^ {(i)} - ax ^ {(i)}) ^ 2}{2 \sigma ^ 2} \right) \tag{2.2.1.15}$$

である。この尤度関数を最大化することで$a$を求めるから、この最適化問題は

$$\underset{a \in \mathbb{R}}{\operatorname{arg}\operatorname{max}} \,\mathcal{L}(a | \mathcal{D}) \tag{2.2.1.16}$$

と書ける。

$(2.2.1.16)$式のような最適化問題の表記は論文や書籍でよく用いられるので覚えておくとよい。最大化問題（maximization problem）に使われる$\operatorname{arg} \operatorname{max}$は argument of the maximum の略で日本語に訳すと「関数の最大値を与える引数」の意味だ。$\operatorname{arg} \operatorname{max}$の下についている$a \in \mathbb{R}$は「実数全体の範囲で」の意味であり、$(2.2.1.16)$式は「関数$\mathcal{L}(a|\mathcal{D})$の値を最大化する$a$を実数全体の範囲で求めよ」という意味になる。この問題の最適解$a ^ \ast$は

$$a ^ \ast = \underset{a \in \mathbb{R}}{\operatorname{arg} \operatorname{max}} \,\, \mathcal{L}(a|\mathcal{D})\tag{2.2.1.17}$$

のように表記される（$a ^ \ast$は$a _ {\rm opt.}$などの表記も好まれる）。厳密には最大値を与える$a$はひとつとは限らず右辺は集合を表しているものと解釈されるため、左辺が実数であるという意図と矛盾するので、上式のような表記は本来であれば好ましくないのだが、論文や書籍ではよく用いられる。

同様に最小化問題（minimization problem）では$\operatorname{arg} \operatorname{min}$を用い、これは argument of the minimum の略である。また、最大化問題や最小化問題の対象となる関数、ここでは$\mathcal{L}(a | \mathcal{D})$のことを目的関数（objective function）という。

最適化問題を解く

最適化問題（optimization problem）が解けるかどうかは異なるし、またその解き方にもバリエーションがあるものだが、$(2.2.1.16)$式は対数尤度関数の最大化というよく使われる方法で解ける。

対数尤度関数とは尤度関数の対数をとったもの、すなわち

$$\ln \mathcal{L}(a|\mathcal{D}) \tag{2.2.1.18}$$

のことである。対数変換は単調変換であるから、対数変換した関数の最大化問題ともとの関数の最大化問題の最適解は一致する。したがって

$$\underset{a \in \mathbb{R}}{\operatorname{arg} \operatorname{max}} \,\, \ln \mathcal{L}(a|\mathcal{D}) = \underset{a \in \mathbb{R}}{\operatorname{arg} \operatorname{max}} \,\, \mathcal{L} (a|\mathcal{D}) \tag{2.2.1.19}$$

であり、解きやすい左辺のほうを右辺の代わりに解いてもよいことになる。

$(2.2.1.19)$式の左辺の$\ln \mathcal{L}(a | \mathcal{D})$を変形していこう。対数を取ると掛け算が足し算になり、さらに$\exp$も外れるという部分がミソだ。式変形に使った式番号も書いておくので、少し大変だが追ってみてほしい。

$$\begin{aligned} \ln \mathcal{L}(a|\mathcal{D}) &= \ln \prod _ {i=1} ^ N \frac{1}{\sqrt{2 \pi \sigma ^ 2}} \exp \left(- \frac{(y ^ {(i)} - ax ^ {(i)}) ^ 2}{2 \sigma ^ 2} \right) \quad \leftarrow (2.2.1.15) \\ &= \sum _ {i = 1} ^ N \ln \left\{ \frac{1}{\sqrt{2 \pi \sigma ^ 2}} \exp \left(- \frac{(y ^ {(i)}- ax^{(i)}) ^ 2}{2 \sigma ^ 2} \right) \right\} \\ &= \sum _ {i = 1} ^ N \left[ \ln \left\{ \frac{1}{\sqrt{2 \pi \sigma ^ 2}} \right\} +\ln \left\{ \exp \left(- \frac{(y ^ {(i)}- ax^{(i)}) ^ 2}{2 \sigma ^ 2} \right) \right\} \right] \\ &= \sum _ {i=1} ^ N \left(- \frac{(y ^ {(i)}- ax^{(i)}) ^ 2}{2 \sigma ^ 2} \right) + \sum _ {i = 1} ^ N \ln \left\{ \frac{1}{\sqrt{2 \pi \sigma ^ 2}} \right\} \\ &= - \frac{1}{2 \sigma ^ 2} \sum _ {i=1} ^ N (y ^ {(i)}- ax^{(i)}) ^ 2 + N \cdot \ln \left\{ \frac{1}{\sqrt{2 \pi \sigma ^ 2}} \right\} \end{aligned}$$

したがって

$$\ln \mathcal{L}(a|\mathcal{D}) = - \frac{1}{2 \sigma ^ 2} \sum _ {i=1} ^ N (y ^ {(i)}- ax^{(i)}) ^ 2 + N \cdot \ln \left\{ \frac{1}{\sqrt{2 \pi \sigma ^ 2}} \right\} \tag{2.2.1.20}$$

が求まった。$a$に関する最大化問題なので、第1項に関する正数倍、および$a$に関係ない第2項は問題の解に影響しないから、

$$\begin{aligned} \underset{a \in \mathbb{R}}{\operatorname{arg} \operatorname{max}} \,\, \ln \mathcal{L}(a|\mathcal{D}) &= \underset{a \in \mathbb{R}}{\operatorname{arg} \operatorname{max}} \,\, \left\{ -\sum _ {i=1} ^ N (y ^ {(i)}- ax^{(i)}) ^ 2 \right\} \\ &=\underset{a \in \mathbb{R}}{\operatorname{arg} \operatorname{min}} \,\, \left\{ \sum _ {i=1} ^ N (y ^ {(i)}- ax^{(i)}) ^ 2 \right\} \end{aligned} \tag{2.2.1.21}$$

である。$(2.2.1.21)$式は最小二乗法（least squares method）の式そのものである。

あとは$a$に関して整理して

$$\sum _ {i=1} ^ N (y ^ {(i)}- ax^{(i)}) ^ 2 = \left( \sum _ {i=1} ^ N \left(x ^ {(i)}\right) ^ 2 \right) a ^ 2 - 2 \left( \sum _ {i=1} ^ N x ^ {(i)} y ^ {(i)}\right) a + \left(y ^ {(i)}\right) ^ 2$$

であり、また二次関数の最小値を取る点は微分して$0$になる点だから、

$$\frac{d}{da} \sum _ {i=1} ^ N (y ^ {(i)}- ax^{(i)}) ^ 2 = 2\left( \sum _ {i=1} ^ N \left(x ^ {(i)}\right) ^ 2 \right) a - 2 \left( \sum _ {i=1} ^ N x ^ {(i)} y ^ {(i)}\right) = 0$$

より、

$$a = \frac{ \sum _ {i=1} ^ N x ^ {(i)} y ^ {(i)}}{\sum _ {i=1} ^ N \left(x ^ {(i)}\right) ^ 2} \tag{2.2.1.22}$$

が解となる。データ点をひとつだけ用いるとき、すなわち$N=1$のとき

$$a = \frac{x ^ {(1)} y ^ {(1)}}{\left(x ^ {(1)}\right) ^ 2} = \frac{y ^ {(1)}}{x ^ {(1)}} \tag{2.2.1.23}$$

であり、この結果は$(2.2.1.22)$式が、代表点をひとつ決めて$a$を求める方法の一般化になっていることを示している。

これで私たちは最小二乗法の裏にあるモデルとその数理を理解したことになる。これから登場する回帰モデルはほとんど今回の議論の拡張に過ぎないので、ここまでの話を理解できたのであればあとはぐっと楽になるであろう。