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

はじめに

　この記事は、統計学を勉強し始めた人に向けて書いています。最尤法やベイズという用語は聞いたけれども、具体的にイメージがしにくいという人に向けて、それらの手法の一面を解説する内容です。

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

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

統計モデルの基本

　統計モデルについて考えるために、まずは簡単な例を通じて必要最低限の道具を揃えていきます。今、あなたは植物の大きさと水やりの量の関係を調べる実験をしているとします。鉢植えに種を蒔いて、毎日決まった量の水をやり、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のようになりました（このデータはサンプルです）。

　Fig 1 を見ると、$x$ の増加に伴って $t$ も増加する直線的な関係があるように見えます。そこで次にやることは、$x$ と $t$ の関係性についての仮説を数式化すること、すなわちモデルを考えることです。今回あなたは、水の量と植物体の重さの間に以下のような関係があるという仮説を立てました。

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$ はそれぞれ直線の切片と傾きに対応します。

　先ほども触れた通り、どんなにうまいモデルを作っても現象を完璧に表現することはできません。水やりの量と植物体の関係性は厳密には直線的ではないでしょうし、計測の誤差なども影響するでしょう。そこで、モデルによる推定値 $y$ と計測された値 $t_i$ の差を $e_i$ で表します。

t_i = y_i + e_i \tag{2}

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

ここまでのまとめ

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

• $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$ は残差で、推定値とデータの差を表します。

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

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

最小二乗法

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

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

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

\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$ は計算で導き出すことができ、その結果は黒の点線で示してあります。

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

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

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

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

• $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には正規分布、ラプラス分布、一様分布の３つの分布を描きました。

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

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

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

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$ は以下のように表されます。

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$ がどれほどの大きさになるのか、事前には分からないので、データから推定することになります。

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

最尤法

　さっそく、導入された確率を使ってパラメータを推定してみましょう。先の議論によって、$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$ の確率分布を図にしてみます。データが生成される確率の高い部分は青色、低い部分は白色で表しています。

　この場合、$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点のデータ全体が生成される確率密度を計算するには、各データ点の確率密度を全てかけ合わせれば良いと積の法則から分かります。数式で表すとこうなります。

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)$ です。

　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)$ はどちらを採用するのがよさそうでしょうか。ここまでの議論してきた通り、確率の考え方を導入したことで、データ点が生成される確率密度を計算することができました。今回は確率密度の高い後者の方が、データ点を表す上でもっともらしい値と言えます。

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

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

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}$ でした。

最尤法と最小二乗法

　ここまでで、パラメータを推定する手法として最小二乗法と最尤法を紹介してきました。サンプルデータに対して最小二乗法を用いてモデルパラメータを推定した結果が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