最尤法とは、ベイズとは （１．最小二乗法と最尤法）

# はじめに

<!---1. ラフ　　2. 文章化　　3．修正　　4. 最終#-->
　この記事は、統計学を勉強し始めた人に向けて書いています。最尤法やベイズという用語は聞いたけれども、具体的にイメージがしにくいという人に向けて、それらの手法の一面を解説する内容です。

　最尤法やベイズといった手法の根幹はとてもシンプルです。そして多くの教科書では最初にその数式を見せ、考え方の基礎を解説します。しかし統計の初学者にとって、それらシンプルな式が意味することを見出すのはとても難しいです。そこでこの記事では、統計学を勉強し始めた頃の自分を思い出しながら、まず確率論を使わない統計手法の解説を行い、それとの対比として最尤法を紹介し、最尤法を使うことによる利点を明確にしたいと思います。

　なおこの記事内での数式や記号は、基本的に機械学習の教科書「パターン認識と機械学習」（C. M. Bishop）を踏襲しています。

# 統計モデルの基本

<!---1. ラフ　　2. 文章化　　3．修正#　　4. 最終-->
　統計モデルについて考えるために、まずは簡単な例を通じて必要最低限の道具を揃えていきます。今、あなたは植物の大きさと水やりの量の関係を調べる実験をしているとします。鉢植えに種を蒔いて、毎日決まった量の水をやり、30日後の植物の大きさを測るという簡単な実験です。

　あなたが興味を持っている**変数**は２種類あります。毎日やる水の量と、播種後30日目の植物の大きさです。これらをそれぞれ、$x$ と $t$ で表すことにします。あなたの興味があるのは、$x$ が $t$ にもたらす影響です。このとき、説明「する」方の変数 $x$ を**独立変数**、説明「される」方の変数 $t$ を**従属変数**と呼びます。

　$x$ と $t$ は変数なので、どんな値も取ることができます。もし無限に実験を繰り返せば、$x$ と $t$ の関係性について限りなく正確な知識を得られるでしょう。しかし実際にできるのは、有限の回数だけ実験を繰り返し、$x$ と $t$ の**データ**を集めることです。データを表す記号は $x_i$ や $t_i$ のように書きます。$i$ は何番目のデータかを表す記号で、例えば実験を $10$ 回繰り返せば $i=1,...,10$ の値を取ります。データを集めて図示するとFig 1のようになりました（このデータはサンプルです）。

<img width="400" alt="Plant_fig1.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/607418/66e16aaa-70c1-a218-961f-2d5c2b201933.png">

<!---1. ラフ　　2. 文章化　　3．修正#　　4. 最終-->
　Fig 1 を見ると、$x$ の増加に伴って $t$ も増加する直線的な関係があるように見えます。そこで次にやることは、$x$ と $t$ の関係性についての仮説を数式化すること、すなわち**モデル**を考えることです。今回あなたは、水の量と植物体の重さの間に以下のような関係があるという仮説を立てました。

```math
y = w_0 + w_1 x \tag{1}
```

　新しく $y$ と $w_0, w_1$ という値が登場しました。$y$ はモデルによる**推定値**です。植物体の大きさに対応する点では $t$ と同じですが、$t$ が実際に観測された値を表すのに対し、$y$ はモデルによって説明される値を表します。もしモデルが完璧に $x$ と $t$ の関係性を説明できるのであれば $t=y$ となりますが、現実にそのようなケースはほぼありません。また、$w_0, w_1$ はモデルの特徴を表す値で、**パラメータ**と呼ばれます。この場合、$w_0, w_1$ はそれぞれ直線の切片と傾きに対応します。

<!---1. ラフ　　2. 文章化　　3．修正#　　4. 最終-->
　先ほども触れた通り、どんなにうまいモデルを作っても現象を完璧に表現することはできません。水やりの量と植物体の関係性は厳密には直線的ではないでしょうし、計測の誤差なども影響するでしょう。そこで、モデルによる推定値 $y$ と計測された値 $t_i$ の差を $e_i$ で表します。

```math
t_i = y_i + e_i \tag{2}
```

　ここで、$y_i$ は式 $(1)$ の変数 $x$ にデータ $x_i$ を代入した値、すなわち $y_i = w_0 + w_1 x_i$ です。$e_i$ はモデルによって説明できない量を表し、**残差**と呼びます。


### ここまでのまとめ

<!---1. ラフ　　2. 文章化　　3．修正#　　4. 最終-->

　ここまでに登場した人物を整理します。

- $x, t$ は**変数**で、いかなる値も取り得ます。$x$ が独立変数、$t$ が従属変数です。
- $x_i, t_i$ は**データ**を表し、それらの値は（70 mL, 53.2 g などの）計測された値です。
- これらの関係性を表すために $y = w_0 + w_1 x$ という**モデル**を用意しました。
- $w_0, w_1$ は**パラメータ**です。
- $y$ はモデルによる**推定値**です。実際にデータ $x_i$ を代入された値は $y_i$ と表します。
- $e_i$ は**残差**で、推定値とデータの差を表します。


　さて、ようやく準備が完了しました。次のセクションでは、パラメータが実際にどういった値を取るのかという点について説明します。

　なお、ここまでの話では確率の話は一切出てきていないことに注意してください！ 確率が登場するのはもう少し先の話です。



# 最小二乗法

<!---1. ラフ　　2. 文章化　　3．修正#　　4. 最終-->

　前のセクションでは、$x$ と $t$ の関係性を表すモデルを導入しました。しかし、モデルの特徴を表すパラメータ $w_0, w_1$ の具体的な決め方については論じませんでした。このままでは、モデルが表す関係性を一つに決定することができません。

<img width="400" alt="Plant_fig2.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/607418/873a16b3-4051-94bd-58ee-b0650771dec2.png">

　Fig 2には、５組のパラメータの値に対応する５本の直線が引いてあります。どの線が最もふさわしいモデルと言えるでしょうか？

　ここで鍵となるのが残差 $e_i$ です。私たちが欲しいのはデータをうまく説明できるようなモデルです。つまり、モデルに説明できない部分（残差）ができるだけ小さくなるような基準でパラメータの値を決めるのは妥当と言えるでしょう。そうした手法の一つに**最小二乗法**があります。あるパラメータの値が与えられると、残差 $e_i$ の具体的な値が計算できます。次に残差を二乗し全て足し合わせた値を計算します。

```math
\sum_{i=1}^{10} e_i ^ 2 = \sum_{i=1}^{10} \{ t_i - (w_0 + w_1 x_i) \} ^ 2 \tag{3}
```

　そして、この値が最小となる $w_0, w_1$ を見つけ出し、パラメータの推定値とします。Fig 2に登場した５本の残差二乗和（Sum of squared errors, SSE）を計算すると、赤い線が最もデータを良く表していることが分かります。実際は残差二乗和を最小化する $w_0, w_1$ は計算で導き出すことができ、その結果は黒の点線で示してあります。

<img width="400" alt="Plant_fig3.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/607418/d34e10e1-d8b7-03db-31de-a82ac3f5d7b4.png">

　こうして最小二乗法を導入することで、無事にパラメータの値を決定させることができました。しかしここで疑問が残ります。この記事では天啓のように二乗和を最小化するという発想をしましたが、例えば残差の絶対値の和や、四乗の和を使ったり、あるいは全く別の方法を使っても良いはずです。実は最小二乗法を正当化する理由は、それらしい直線が引けるということ以外にはここでは与えられていません。

　そこで次のセクションでは、モデルに**確率**の考え方を導入します。そして更に次のセクションにて、確率の考え方によって裏付けされた方法である**最尤法**でパラメータを推定します。



# モデルに確率分布を導入する

<!---1. ラフ　　2. 文章化　　3．修正#　　4. 最終-->

　ここで再び登場人物を整理します。

- $x, t$ は**変数**で、いかなる値も取り得ます。$x$ が独立変数、$t$ が従属変数です。
- $x_i, t_i$ は**データ**を表し、それらの値は（250mL, 50gなどの）計測された値です。
- これらの関係性を表すために $y = w_0 + w_1 x$ という**モデル**を用意しました。
- $w_0, w_1$ は**パラメータ**です。
- $y$ はモデルによる**推定値**です。実際にデータ $x_i$ を代入された値は $y_i$ と表します。
- $e_i$ は**残差**で、推定値とデータの差を表します。

　さて、ここからは確率の考え方を導入するわけですが、具体的に登場人物の誰に確率を当てはめればよいでしょうか？ 着目するのは**残差** $e_i$ です。もしモデル $y = w_0 + w_1x$ が水やりの量と植物体の重さの関係性を十分に表現できているならば、そのモデルをもってしてなお表現できない残差は、重さの測定誤差や植物体の成長の偶発的なばらつきといった**ランダムな**要素であると考えられます。そしてランダムな数を扱うツールこそが、確率の考え方です。

　ここで新しい表記を加えます。今まで残差は変数ではなく、各データ点に対応した値と考えていたため、必ず添え字 $i$ をつけて $e_i$ と表記してきました。これからは残差も $x, t$ と同じ変数 $e$ と表記し、ランダムに変動する $e$ を実験で観測できた値として $e_i (i = 1, ..., 10)$ があると仮定します。

　次にすべきことは、確率変数 $e$ の従う**確率分布**の設定です。$e$ はどのような確率でどのような値を取るのでしょうか。例として、Fig 4には正規分布、ラプラス分布、一様分布の３つの分布を描きました。

<img width="400" alt="Plant_fig4.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/607418/7517c940-93c2-8b07-50b6-ce247fa18b1b.png">

　まず一様分布から見てみます。一様分布は、残差はある一定以上・以下の値は取らない取らないと仮定することになります。私たちは残差の取る範囲があらかじめ分かるわけではないので、これは仮定が強すぎるように思えます。

　正規分布とラプラス分布は、０に近い値を取る確率が高く設定され、０から離れるごとにその値を取る確率が下がる分布です。これは自然な仮定に見えます。これら二つの中からどちらを選ぶかとなりますが、今回は**正規分布**を採用します。理由は色々とありますが、分かりやすい理由として計算の簡単さがあります。ラプラス分布は尖っている分布なので微分不可能な点がありますが、正規分布は滑らかで微分の計算も簡単です。この記事では計算式は登場しませんが、パラメータの推定などにおいてこの性質は非常によい働きをします。

　また、正規分布それ自体にも、平均と分散の二つのパラメータがあります。平均は分布の中心の値、分散は裾野の広がり方を決定します。今回は平均は０に固定し、分散は未知の値 $\sigma_e^2$ として扱います。これを式で表すとこうなります。

```math
p(e) = N(e | 0, \sigma_e^2) \tag{4}
```

　確率変数を導入したところで、登場人物を整理しましょう。

- $x$ は**独立変数**で、いかなる値も取り得ます。
- $x_i, t_i$ は**データ**を表し、それらの値は（250mL, 50gなどの）計測された値です。
- これらの関係性を表すために $y = w_0 + w_1 x$ という**モデル**を用意しました。
- $y$ はモデルによる**推定値**です。実際にデータ $x_i$ を代入された値は $y_i$ と表します。
- $e$ は**残差**で、推定値とデータの差を表します。分布 $p(e) = N(e | 0, \sigma_e^2)$ に従う**確率変数**です。
- $t$ は**従属変数**であると同時に、分布 $p(t|x) = N(t | w_0 + w_1 x, \sigma_e^2)$ に従う**確率変数**です。
- $w_0, w_1, \sigma_e^2$ は**パラメータ**です。

　気づきましたか？ $e$ だけでなく、特に言及しなかった $t$ も**確率変数**になっています。今の仮定において、$t$ は以下のように表されます。

```math
t = y + e = w_0 + w_1 x + e \tag{5}
```

　$w_0, w_1$ はパラメータ、 $x$ は変数で、これらは確率変数ではありません。$e$ を確率変数に設定したので、それらに数を足した $t$ も確率変数となります。$t$ の従う確率分布は、$e$ の確率分布を $w_0 + w_1 x$ だけずらしたものです。すなわち、平均 $w_0 + w_1 x$ 、分散 $\sigma_e^2$ の正規分布 $N(t | w_0 + w_1 x, \sigma_e^2)$ です。

　もう一つの変化にも気づいたでしょうか。パラメータが $w_0, w_1, \sigma_e^2$ の３つに増えています。$\sigma_e^2$ がどれほどの大きさになるのか、事前には分からないので、データから推定することになります。

　長くなりましたが、モデルに確率変数の導入することができました。遂に次のセクションで、確率の考え方を元にパラメータの推定を行います。



# 最尤法

<!---1. ラフ　　2. 文章化　　3．修正#　　4. 最終-->

　さっそく、導入された確率を使ってパラメータを推定してみましょう。先の議論によって、$t$ は確率変数で正規分布 $N(t | w_0 + w_1 x, \sigma_e^2)$ に従うことが分かりました。すなわち、３つのパラメータ $(w_0, w_1, \sigma_e)$ の値を決めれば、データ $t_i (i = 1, \ldots, 10)$ が生成される確率を計算できることになります。

　例えば $(w_0, w_1, \sigma_e^2) = (91, 0, 25)$ と決めたときの $t$ の確率分布を図にしてみます。データが生成される確率の高い部分は青色、低い部分は白色で表しています。

<img width="400" alt="Plant_fig5.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/607418/7bb74604-5ef8-e5e7-026a-27df2e1d4bdf.png">

　この場合、$70 \leq x \leq 90$ の５つの点は、確率が高い場所に位置しています。一方で $50 \leq x \leq 65$, $x = 95$ の５つの点は、確率の低い場所に位置しています。では実際どれくらいの確率でこれら10個の点が生成されるのか、計算してみましょう。

　$t$ は連続的に変化する変数なので、本来 $t$ がある一定の値（例えば50ちょうど）を取る確率を計算することはできません。０になってしまいます。ここでは代わりに、**確率密度**を計算することにしましょう。これは Fig 4の縦軸の値で、例えばFig 4の正規分布の場合、値が０となるときの確率密度は $0.399$ のように計算することができます。

　数式で表すと、$t_i$ の生成される確率密度は $p(t_i | x_i, w_0, w_1, \sigma_e^2)$ と書きます。これは**条件付き確率**と呼ばれるもので、$p(A|B)$ と書くと「$B$が与えられた時の$A$の確率密度」を表します。この場合、$t_i$ の確率を計算するときは、対応するデータ $x_i$ とパラメータ $w_0, w_1, \sigma_e^2$ を固定して考えるので、条件付確率の中で「与えられる」側の数として扱われます。

　具体的に計算すると、Fig 5の左の点から順に  $0.00162$, $0.00083$, $0.02348$, $0.00014$, $0.07740$, $0.07081$, $0.07524$, $0.06923$, $0.07685$, $0.00200$ となります。青い部分に乗っている５点は、確率密度も高い値を示すことが分かります。

　データ点それぞれが生成される確率密度は計算できました。次に欲しいのは10点のデータ全体が生成される確率密度です。ここで、各データ点が**独立に**生成されていると仮定します。例えば $t_1$ が推定値 $y_1$ より低い値だったとき、$t_2$ の値もそれに引っ張られて推定値 $y_2$ より低くなる、というような現象があれば、$t_1$ と $t_2$ は独立ではないと言えます。逆に $t_1$ と $t_2$ の値の出方に関係がないとき、これらを独立であると言います。

　データ点同士が独立であると仮定すると、10点のデータ全体が生成される確率密度を計算するには、各データ点の確率密度を全てかけ合わせれば良いと**積の法則**から分かります。数式で表すとこうなります。

```math
p(t_1, ..., t_{10} | x_1, ..., x_{10}, w_0, w_1, \sigma_e^2) = \prod_{i=1}^{10} p(t_i | x_i, w_0, w_1, \sigma_e^2) \tag{6}
```

　これに従ってFig 5のデータ10点が生成される確率密度を計算すると、$2.04946 \times 10 ^ {-20}$ という値になります。

　もう一つ例を見てみましょう。今度は $(w_0, w_1, \sigma_e^2) = (45, 0.6, 9)$ です。

<img width="400" alt="Plant_fig6.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/607418/8cda211d-e69b-805b-b185-0ce8a45467ac.png">

　Fig 5よりも多くの点をカバーできているように見えますが、３つほど青い帯からはみ出した点が見えます。各点の確率密度は、左から $0.10538$, $0.10372$, $0.10214$, $0.00022$, $0.08687$, $0.06884$, $0.13237$, $0.00509$, $0.01156$, $0.09187$ で、これらを全てかけ合わせた値は $1.05882 \times 10 ^ {-15}$ です。

　それでは上に挙げた２つのパラメータの値、$(w_0, w_1, \sigma_e^2) = (91, 0, 5)$ と $(45, 0.6, 3)$ はどちらを採用するのがよさそうでしょうか。ここまでの議論してきた通り、確率の考え方を導入したことで、データ点が生成される確率密度を計算することができました。今回は確率密度の高い後者の方が、データ点を表す上で**もっともらしい**値と言えます。

　今回は２通りのパラメータ値しか検討しませんでしたが、実際はパラメータは連続的にどのような値も取ることができます。どうやって最もふさわしいパラメータの値を選べばよいか、もうお分かりでしょう。確率密度を最大にするパラメータの値こそが、データを説明する上で最ももっともらしい値と言えるでしょう。この方法こそが、**最尤法**と呼ばれる、今日の統計学の中でも頻繁に採用される由緒正しき手法です。

　１つだけ用語の説明を加えましょう。最尤法では「データが生成される確率密度」を最大化すると言ってきましたが、この最大化される値を**尤度**と呼びます。数式ではこう書きます。

```math
p(\boldsymbol{t}|\boldsymbol{x}, \boldsymbol{\theta}) \tag{7}
```

　ベクトルの表記を導入したので見慣れない形に見えますが、式(6)と同じことを言っています。

- 従属変数　  $\boldsymbol{t} = (t_1, ..., t_{10})$
- 独立変数　  $\boldsymbol{x} = (x_1, ..., x_{10})$
- パラメータ  $\boldsymbol{\theta} = (w_0, w_1, \sigma_e)$

　最尤法の説明の最後に、最尤法でサンプルデータを説明するモデルのパラメータを推定した結果を見てみましょう。$(w_0, w_1, \sigma_e^2) = (47.727, 0.54665, 28.858)$ で、尤度は $9.344612 \times 10 ^ {-14}$ でした。

<img width="400" alt="Plant_fig6.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/607418/092aef49-f2d1-b81c-4ed6-38da4958a704.png">



# 最尤法と最小二乗法

　ここまでで、パラメータを推定する手法として**最小二乗法**と**最尤法**を紹介してきました。サンプルデータに対して最小二乗法を用いてモデルパラメータを推定した結果がFig 3の黒点線、最尤法を用いてモデルパラメータを推定した結果がFig 7です。見比べると、非常に近い位置に直線が引かれていることが分かります。実はこういったケースにおいて、最小二乗法と最尤法で推定されるパラメータの値は一致することが知られています。

　最小二乗法の説明をしたセクションの最後で、絶対値の和や四乗和ではなく二乗和を選ぶのはなぜかという話題に触れました。ここまで読んだ方は答えられるはずです。その理由は、残差が独立な正規分布に従うと仮定したときに尤度が最大となるパラメータが得られるからです。こうして、確率の考え方を使って最小二乗法を正当化することができました。

　では逆に、答えが同じになるなら、なぜわざわざ確率を持ち出して最尤法を定義したのでしょうか？ それは、**モデルの拡張**が可能となるからです。今回は $x$ と $t$ が直線的な関係にあり、残差の分布を正規分布と仮定しました。しかし、例えば $x$ と $t$ が指数関数や対数関数で表される関係だったり、残差の分布に二項分布やポワソン分布を仮定することがあります。そのような場合、もはや最小二乗法では適切にパラメータを推定することができません。一方最尤法を使うと計算は複雑になりますが、確率に基づいた適切な方法でパラメータを推定することができるのです。

　今回の記事はここまでですが、実は最尤法も万能ではありません。次に続く記事では最尤法の限界について触れ、最尤法ではないアプローチとベイズ的なアプローチについて解説したいと思います。

次回URL　→　https://qiita.com/YusukeToda1984/items/39831b75fa86099b80c3
図の描画はRで行いました。スクリプトをGitHubで公開しています　→　https://github.com/YT100100/QiitaArticles/blob/master/20200329_WhatIsBayes_Section1_Fig.R