自然科学の統計学4章最尤法を読んだその1

はじめに
記事リンク
第4章最尤法
- 4.0 お気持ち
- 4.1 一般線形モデル
続きはまた今度

はじめに

昨日はニコニコ超会議 2018に参加して、バーチャルYouTuberの電脳少女シロを応援してきた。ニコニコ超会議は思ってたより人口密度が低くて快適だったが、どうも過去最高来場者数を達成したらしい。余裕のある会場を取ることは大切だ。

昨日は人生で始めてオタク・イベントに始発参加して、テンションが最高潮に達した結果、キモ・オタク特有の自分語りをドヤ顔で繰り出してしまった。あまりにもキモいので、ぜひみんなにも読んでほしい。そして一緒にシロちゃんを愛でてほしい。

kriver-1.hatenablog.com

さて、最新記事がいつまでも推しへのラブレターなのは流石にちょっと恥ずかしいので、照れ隠しに以前から書かなければならないと思っていた自然科学の統計学の解説記事を書こうと思う。自然科学の統計学はわりと難しい本なので、統計学の基礎知識がないと大変険しい。ただ、その分自分にどこが足りていないかがはっきりするのでとても勉強になる。

僕が今回読んだのは4章と6章で、僕は4章のほうが面白いと思ったので、4章を紹介する。暇だったら6章もそのうち紹介するし、もしかしたら他の章も読んで紹介するかもしれない。そのへんはこの記事の反響次第だ。

記事リンク

次回
kriver-1.hatenablog.com

最終回
kriver-1.hatenablog.com

第4章最尤法

4.0 お気持ち

この章は、最尤法という手法について、最尤法がどんなもので、それがいかに素晴らしいかということを語る章だ。

最尤法は、「たくさんの観測データからデータの生成源を当てるゲーム」の解き方のひとつ。たぶんみんなも無意識のうちに一度はやったことがある手法だと思う。そのくらい内容は単純だ。
ただ、その手法でほんとうに良いのかを考えるのはなかなか難しい。なんたって、そもそも何を以て「良い」とするのかも我々は知らないのだから。

なので、この章ではまず、いくつか「良い」とされる指標を挙げていく。そして最後に、最尤法がそれらの指標をぜんぶ満たしていることを示す。

本の中では語られなかったけども、いやそれみんな知らんやろ、と思った内容については、この 4.0 節と同じように、 .0 を付けた節または項で僕が勝手に説明を加えていく予定だ。
ただし、0以外の数字が付いているからといって、本にも同じ節・項があるとは限らない。より詳細な分別が必要だと判断した場合は、僕の独断と偏見で勝手に節・項を追加したりすることがある。その場合、対応する節や項の番号が本とずれないようには配慮する予定だ。

4.1 一般線形モデル

4.1.0 モデルと最適化

最尤法がどんなものか説明する前に、まず最尤法が使われる場面を説明しなければならない。

最尤法は、「一般線形モデル」を最適化するときによく用いられる道具である。モデルを最適化するというのは少し馴染みにくい表現かもしれないが、次のように考えてほしい。

1. まず、ある観測データの山 ${Y=\{y_1, \cdots, y_N\}}$ に対して、その観測データがどういう分布から来ているのかを推定したいと考える。
これは例えば、次の ${y_{N+1}}$ を推測したいということかもしれないし、 ${\{y_1, \cdots, y_N\}}$ の中で最もよく現れる ${y}$ はどんなものなのか知りたい、ということかもしれないが、何にせよ ${\{y_1, \cdots, y_N\}}$ を解析したいということだ。
問題によっては、観測データが ${\{(x_1, y_1), \cdots, (x_N, y_N)\}}$ と組になっているかもしれない（これは ${x_i}$ が決まれば ${y_i}$ の分布が決まるということだ）が、そのへんはどっちでもかまわない。どっちでもかまわないので、書き方としては後者にそろえておく。後者の ${x}$ を定数にすれば前者と全く同じデータになる。

2. そのために、観測データの山を説明するためのモデルを考える。
モデルというのは、パラメータ ${\theta}$ と説明変数 ${x}$ に対して ${y}$ の分布を与える関数のことだ。我々は ${y}$ がどのような分布から来ているか知りたいので、これは答えの範囲を制限しているということになる。この世に存在するすべての分布を考えることは到底不可能なので、モデルと ${\theta}$ で表せる範囲に答えを絞ってしまおう、ということだ。モデルが複雑すぎるとパラメータを決めるためにデータがたくさん必要になってしまうし、モデルが貧弱すぎるとデータがぜんぜん説明できない。このモデルについては定数だと思って、パラメータを自由に動かすことでデータを表すことを考える。

3. 観測データを用いて、モデルのパラメータを決める。データからできるだけそれっぽいパラメータ ${\theta}$ を決定するということだ。もちろん、サイコロを振ってデタラメに ${\theta}$ を決めたっていいんだけど、それだとせっかくデータがあるのにもったいないので、普通はデータと見比べてできるだけそれっぽい ${\theta}$ を選択する。

この工程3.のことを、「モデルを最適化する」とよく言う。「よく言う」というぼかした表現なのは、僕がそう表現しているだけで、この本には最適化するとは一言も書かれていなかったからである。
たぶん最適化する (optimize) で合ってると思うんだけど、もしかしたら学術用語の使い方を間違えているかもしれない。間違ってる場合は指摘してください。

話を戻して、つまり最尤法は、1. 観測データの山に対して、2. 一般線形モデルというモデルを考えたときに、3. そのモデルのパラメータをいい感じに決める、手法の一つだということだ。

4.1.1 一般線形モデルとは

さて、一般線形モデルというのは、データ ${y}$ の確率分布 ${p=\text{Pr}(y)}$ を何かしら変形した ${f(p)}$ に対して、線形な関数 ${f(p) = \theta X}$ を仮定したモデルである。
すべてのモデルは ${p=g(x; \theta)}$ という形で表すことができるが、その中でも ${p=f^{-1}(\theta X)}$ という形で表せるものを一般線形モデルと呼ぶ、ということである。

変形する関数 $f(z)$ はリンク関数と呼ばれ、式に ${x}$ や ${\theta}$ が含まれていなければたぶんなんだって構わない（全単射とか微分可能とかそういう良さげな性質はたぶん持ってないと困ると思う）。
${\theta X}$ のことを $\theta$ に対して線形、と表現する部分についての説明は省略する。知ってる人は当たり前に知ってる知識だし、知らない人に対してこの部分の説明をするのは結構骨が折れることのように思えるからだ。

気持ちとしては、ただの線形モデルだとつまらないので、変形 ${f}$ を一枚噛ませてもおっけーだよ、としたモデルだと思ってほしい。実際には本質はそこではなく、線形モデルを仮定する時点で値が正規分布に従うことを仮定しているので、それを $y$ について仮定することを避けられるという点が本質なのだが、それは今回の本題ではない。
あくまでも変形はモデルの一部なので、どういう変形にするかは最適化の対象外だ。

この一般線形モデルを最適化する手法の一つが最尤法である。

この教科書では generalized linear model の訳として一般線形モデルを用いているが、general linear model という用語と訳し分けるために、generalized linear model を一般化線形モデルと訳す文化もあるらしい（というかこっちが主流らしい）。
この教科書では、たぶん general linear model の訳としては linear model と同じ線形モデルを使っている。実際、general linear modelと linear model の差はほとんどないし、それに対して general linear model と generalized linear model の差はめちゃくちゃ大きいので、 general linear model に個別の訳語を与えないという姿勢のほうが僕は共感できる。なので、本記事でも generalized linear modelに対して一般線形モデルという訳語を当てた。ググる時は注意してほしい。

僕は詳しくないのでわからないが、たぶん最尤法自体は一般線形モデル以外のモデルにも使える手法だと思う。ただ、今回の内容では、一般線形モデルに対して使うと思って各種の証明をする。実際一般線形モデルはかなり手広い範囲のモデルをカバーしているように見えるし、一般線形モデルより複雑なモデルは数式で扱うのが大変なのでそもそも最尤法なんてやってられない場合が多いんじゃないだろうか（例えばニューラルネットワークとか）。

4.1.2 一般線形モデルの例1

例えば、ラットにある毒物を量 ${x}$ だけ与えたとき、ラットが死んだかどうかを ${y=0,1}$ として（生きているとき ${y=0}$ 、死んでいるとき ${y=1}$ ）、データ点 ${\{(x_1, y_1),\cdots, (x_N, y_N) \}}$ から ${x}$ が与えられたときの ${y}$ の確率分布を考えたいとする。

最初の例にしてはなかなか残酷な話に思えるが、本に載っている例をそのまま載せただけなので勘弁してほしい。

先ほどの手順にしたがって、まずは 2. モデルを考える。ここでは、最初はただの線形モデルを考えてみよう。線形モデルとは一般線形モデルの ${f}$ が恒等写像、つまり何もしない関数になっているモデルのことだ。
すなわち、ラットが死ぬ確率 ${p_i=\text{Pr}(y_i=1|x_i)}$ に対して、 ${p_i = ax_i + b}$
というモデルを考えるということになる。

このモデルのパラメータは ${a}$ と ${b}$ の2つだ。さっきまでパラメータのことを ${\theta}$ と表していたから、パラメータが1つから2つに増えてるじゃないか！と思う人もいるかもしれないが、数学にはベクトルという便利なものがあるので、単に ${\theta=(a, b)}$ だと思ってしまえば何の問題もない。むしろ、 ${\theta}$ と書いていた頃は実際の数字が何個並んでいても構わなかったのだから、パラメータは2個に減ってしまったと言ったっていい。
また、先程 $X$ とぼかして書いていた部分だが、これもいまは $X=(1, x)$ としている（定数項を追加している）。実は線形モデルは $\theta$ に対して線形でさえあればよく、 $x$ についてはどんな関数でも（非線形でも）構わない。というのも、 $x$ をどう変形するかという部分はモデルを構築する時点で決まるものであって、最適化の対象外だからである。

さて、このモデルを3. 最適化する手法はいろいろあるだろうが、ここではその詳細な手順については触れない（本題ではないので）。
ここで言いたいのは、この問題に対して線形モデルを考えるのはあまり適切ではないのではないか、ということだ。

例えば、このモデルを何かしらの方法で最適化したとする。そのときにはきっとすごく最適なパラメータ ${(a^{*}, b^{*})}$ があって、このモデル ${p_i=a^{*}x_i + b^{*}}$ はデータに対してすごく最適な感じになってるんだと思う。
ここで、もし ${x_{N+1}}$ としてすごく大きな値が来たらどうなるだろうか？与える毒物が多いほど死亡率は上がるだろうから、きっと ${a^{*}}$ は正の値だろう。そこにすごく大きな値が掛かるわけだから、 ${p_{N+1}}$ はきっと1より大きな値になってしまう。 ${x}$ が大きくなれば、 ${p}$ はどこまでも大きくなっていく。確率分布を表すモデルとして、これは不合理だ。
${1}$ より大きな確率は ${1}$ に丸める、という自分ルールを定めてもいいが、それはもはや線形モデルとは言えないものになってしまう。

一般線形モデルなら、この問題に対処できる。 ${0}$ から ${1}$ の値を全空間にマップするようなリンク関数を考えてあげればいいからである。

例えば、ロジット関数 ${f(p) = \log {\frac{p}{1-p}} }$ はよさげなリンク関数の一つだ。 ${p}$ が ${0}$ から ${1}$ までの値を取るとき、 ${f(p)}$ は ${-\infty}$ から ${\infty}$ までの値を取る全単射な関数になっている。このようなモデルを仮定することで、「せっかくがんばって最適化したのに、なんだか不合理なモデルになってしまった(;_;)」という事態を回避することができる。

ちなみに、この例は確率分布に対してロジスティックモデルを仮定したのと同値である。ロジスティックモデルは ${p'=\frac{1}{1+\exp(a'x+b')}}$ という式で表されるモデルだが、

$\displaystyle \begin{eqnarray} \log p' &=& -\log \left( 1+\exp(a'x+b') \right) \\ \log (1-p') &=& (a'x+b') - \log \left( 1+\exp(a'x+b') \right) \\ \therefore f(p') &=& \log {\frac{p'}{1-p'}} \\ &=& (-a')x+(-b') \end{eqnarray}$

となるので、ロジット変換の一般線形モデルと $a'=-a, b'=-b$ としたロジスティックモデルは同じものになる。

4.1.3 一般線形モデルの例2

例えば、あるクラスにおいて数学と国語の試験をしたとする。
クラスの人数が $N$ 人なら、観測データは $\{(y^1_1, y^2_1), \cdots, (y^1_N, y^2_N)\}$ のようになるだろう。ここで、 $y^1_i$ が数学の点数を、 $y^2_i$ が国語の点数を表している。
今回は「数学の点数が◯点な人は国語の点数は何点になるか」ということではなく、数学と国語の点数両方に興味があるという設定を考えているので、そういう意味で $x, y$ ではなく両方とも $y$ で表してみた。

わかりやすいのは、数学の点数の確率分布を $p^1$ 、国語の点数の確率分布を $p^2$ として、 $\text{Pr}(y^1, y^2)=p^1(y^1)*p^2(y^2)$ とすることである。これはつまり、数学と国語の点数が独立であると仮定していることになる。
モデルのパラメータは、点数が0点から100点の101通りあるとして、 $p^1(0), \cdots, p^1(100), p^2(0), \cdots, p^2(100)$ の202個になる。

だが、このモデルも実は、見方を変えれば一般線形モデルになっている。
いま、数学の点数が $y^1$ で、国語の点数が $y^2$ である確率 $\text{Pr}(y^1, y^2)$ について、その対数を取ることを考える。すると、
$\displaystyle \begin{eqnarray} \log \text{Pr}(y^1, y^2) &=& \log p^1(y^1) + \log p^2 (y^2) \end{eqnarray}$
となる。ところが、ここで $p^1(0), \cdots, p^1(100), p^2(0), \cdots, p^2(100)$ の代わりに対数を取った $\log p^1(0), \cdots, \log p^1(100), \log p^2(0), \cdots, \log p^2(100)$ をパラメータだと考えてあげれば、右辺はパラメータの線形な演算（ただの足し算）になっているから、これは実は、対数をリンク関数とした一般線形モデルになっていたことがわかる。

この例では説明変数 $x$ が出てこなかったが、前述したように説明変数は別になくても構わない（あるいは、説明変数が0次元ベクトルなんだと考えてもいい）。そうした場合、パラメータのただの足し算になっていればそれは一般線形モデルである。

このように、一般線形モデルの中には様々な有用なモデルが含まれている。したがって、この一般線形モデルを最適化する手法である最尤法もまた、有用な手法である。

次回以降はこの最尤法がどういうもので、どういう意味で良いものなのかをみんなと一緒に確認していきたい。

続きはまた今度

普通に4章を終わらせる予定だったが、まさかの4.1だけでとんでもない長さになってしまった（まだ本のうち4ページ半しか進んでいない）。さすがにこれ以上だらだら続けるとページが重くなってしまうので、一旦このあたりで筆を置きたいと思う。

もし続きが読みたいと思ったら、記事にいいねを付けてくれるとありがたい。
人間は承認欲求を満たしたい生き物なので、承認欲求が得られない記事を書き続けるのはあまり楽しくない。たった1クリックで他人を幸せにできる機会なんてそうそうないので、善行を積む気持ちでその星マークをそっと押してほしい。そう、その星だ。よろしくお願いします。青い鳥の方でもいいよ。

次回
kriver-1.hatenablog.com

最終回
kriver-1.hatenablog.com