よく「ほげ分布を仮定してパラメータを最尤推定しました」というのを聞くと思うんですが、それってどういうことなんだろうと戸惑ったことがあったので、どういうことなのか考えたことを書きました。
※ 全体を通して統計的推定の点推定の話です。ベイズ推定には触れません。
参考文献
以下の書籍の3章を参考にしていますが解釈の誤りは私に帰属します。お気付きの点があればご指摘いただけますと幸いです。
日本統計学会公式認定 統計検定1級対応 統計学 | 二宮嘉行, 大西俊郎, 小林 景, 椎名 洋, 笛田 薫, 田中研太郎, 岡田謙介, 大屋幸輔, 廣瀬英雄, 折笠秀樹, 日本統計学会, 竹村彰通, 岩崎学 |本 | 通販 | Amazon
「分布のパラメータを推定する」ってどういうことなの
データを分析するときに、データが何らかの確率密度分布から生成されたものと考えて、どんな確率密度分布から生成されたんだろうと推測するということがあると思います。そのやり方の一つに、パラメトリックな確率密度分布（＝いくつかのパラメータで分布の形を調整できる仕様になっている確率密度分布）を何か決め打って、そのパラメータの値を推定するというものがあると思います。例えば、「ある学校の生徒たちの身長の測定結果が正規分布にしたがうと仮定して、そのパラメータである $\mu$ と $\sigma$ の値を推定する」という具合にです。
ただ、別に生徒たちの身長の測定結果は確率分布から生成されているわけではないと思います。もしかしたらこの世界はメタ世界の人にプログラミングされていて、生徒たちの身長の測定結果が確率分布から生成されているということもあるかもしれませんが、正規分布かはわからないし、正規分布だったとしてパラメータの真の値はメタ世界の人にコンタクトしない限りわからないです。というのも、手元にある有限個のサイズのデータから真の分布を知ることはできないからです。なので、分布のパラメータの真の値を得る方法は正直ないです。
でも、世の中にはデータから分布のパラメータを得ようとしている人たちがたくさんいると思います。そういう人たちは、「自分はこういうやり方でパラメータの推定値を得る」と勝手に決めていると思います。「尤度を最大にする（最尤推定）」とか「誤差の2乗和を最小にする（最小2乗推定）」とか決めていることが多いと思います。「推定する」というといかにもパラメータの真の値に近い値を得ようと努力しているようですが（実際みんなそうしようと努力しているとは思いますが）、どちらかというとここでの「推定する」とは「値を求めるやり方を勝手に決めて、そのやり方で値を求めること」といった方が感覚に合っている気がします。
「それってそんなに勝手なことなの」と思うかもしれません。確かに、パラメータの真の値に近い値を得ようとするのに「尤度を大きくする」「誤差を小さくする」などというのはそんなに自分勝手なやり方とは思われません。ただそう感じられるのは、尤度がより大きくなるような／誤差がより小さくなるようなパラメータにするほど分布の形とデータの広がり方が合致するように感じられる例を頭の中に思い描いているからだと思います。例えば次のような場合を考えてみてください。ある学校では相撲がとても盛んで、全校生徒の実に約3割ほどの生徒たちが相撲部に所属しています。相撲部の生徒たちは強靭な肉体をしているので、全校生徒の体重測定結果のヒストグラムを描くと、相撲部ではないと思われる生徒たちがつくる山と、相撲部と思われる生徒たちがつくる山と、2つの山があるヒストグラムになります。そこであなたはこの学校の生徒たちの体重測定結果を2つの正規分布の混合分布と決め打つことにし、それぞれの $\mu$ と $\sigma$ を推定することにしました。推定のやり方は、最尤推定という方法をよく聞く気がするので、最尤推定をすることにしました。しかしあなたはすぐ気付くはずです。2つの正規分布のうちどちらか一方の $\mu$ を、誰か一人の生徒（相撲部でも相撲部でなくても構いません）の体重にぴったり合わせて、その正規分布をどんどん尖らせれば（ $\sigma$ を小さくすれば）、それだけで全データに対する尤度がどこまでも大きくなることに。でもそれは、なんだかあなたの推定したかった分布の姿ではありません。しかし、「尤度を大きくする」というルールだとこのような尖った分布の方が「よりよいパラメータの推定値」ということになってしまいます。実はこれは尤度関数が最大値をもたない例（つまり、そもそも最尤推定ができない例）ですが、尤度関数が最大値をもつ場合においても、最尤推定は推定のやり方として望ましい性質を満たさないことがあります。では「誤差を小さくする」という推定のやり方はどうでしょうか。これこそ問題がある方針とは感じられません。しかしこの素朴な方針においても、「誤差の分布」がある条件を満たしていなければ、推定のやり方が望ましい性質を満たさなかったり、見せかけの分布の姿が見えてしまったりすることがあります。
なので、あるパラメータの「最尤推定値」などといったとき、それは「データに対する尤度が最大になるような値」という意味でしかなく、最尤推定値であれば絶対的によい推定値であるとかは全くないのです。真の分布を知り得ない以上、自分ルールで決めたよい推定値を求めるしかなく、その意味では推定というのは勝手なものにしかなりません。
「でも上に挙げたのは特殊な例でしょ」と思うかもしれません。どんなデータをどんなやり方で推定することが多いかは人それぞれなので何が特殊かは全然わからないですが、特殊な例なのかもしれません。しかし、真の分布に近づくのに絶対よいやり方というものがなく、常に何か勝手なやり方をしているのだというのは念頭に置いておいた方がいいと思います。
※ この記事では扱いませんが、統計的推定には点推定や区間推定のほかに仮説検定というものがあって、その枠組みでは「データがある確率分布から生成されていることを否定できるか」を判定しようとします。が、 やはり真の分布というものは知り得ないので、「データが正規分布にしたがうか否か」という検定ひとつをとっても、数多くのやり方があり、それぞれに得手不得手があり、真の分布を判定することの難しさをみることができると思います。
※ 上に挙げた「混合正規分布が最尤推定できない」という問題を回避する方法に「EMアルゴリズム」というものがあります。興味がある方はEMアルゴリズムが一体どのように考えて上の問題を回避したのか調べてみてください。
※ 実用的には、データの分布推定は中間目標にすぎず、推定した分布を利用して予測をしたり、異常検知をしたり、その他データに対して何らかの知見を得ることなどを目指すことも多いと思います。その場合はそちらの最終目標があるので、自分勝手なルールではない推定のよさの基準が存在するとは思います。
といっても何か「こういう推定の方がよい」とかないの
「分布のパラメータを推定する」とは勝手なことだと書いてしまいましたが、じゃあどんなやり方を選べばいいのか困ってしまうと思います。ある学校の生徒たちの身長の測定結果にあるパラメトリックな分布を仮定して、AさんとBさんとCさんがそれぞれ異なる推定のやり方で異なるパラメータ推定値を出したとします。みんな口々に「自分のやり方では自分の値が一番いい推定値だ」と主張します。誰の推定値が一番よいものなのか全然わかりません。
一度冷静になって、どうしてこんなことになってしまったのか考えてみます。誰の推定値が一番よいものなのかわからないのは、パラメータの真の値を知り得ないからです。ではなぜパラメータの真の値を知り得ないのかというと、手元に有限個のサイズのデータしかないからです。…なら、データのサイズを無限にとばしたら、真のパラメータが得られるはずです。そこで、次のような検証方法を思いつきます。既にパラメータがわかっている分布からデータを生成して、AさんとBさんとCさんの推定のやり方によってパラメータの推定値を出してもらいます。次に、データをもっと追加して、また推定値を出してもらいます。もっともっとデータを追加して、繰り返します。推定値が真のパラメータに収束していけば、想定通りです。もし推定値が真のパラメータとは違う値に収束したり、いつまでも収束しなかったりしたら、ちょっと推定のやり方として心配だと思います。「データサイズを大きくするほど真のパラメータに近づいてほしい」という当然の期待を満たしてくれないからです。この期待を満たしてくれるとき推定量（※）が一致性をもつといいます（厳密には、データサイズ $n$ を $n \to \infty$ にしたときに推定量 $\hat{\theta}(X)$ が真のパラメータ $\theta$ に確率収束するならば、推定量 $\hat{\theta}(X)$ が一致性をもつといいます）。
※ 「推定量」といったとき、「データ $X$ を代入すると推定値を返す関数」といった意味があります。「最尤推定量」というと「データに対してそのデータに対する尤度関数を最大にするパラメータ値を対応させる関数」で、「最尤推定値」は実際に手元のデータを代入して得た値というイメージです。ただ「推定量」と「推定値」の用語の使い分けはそれほど厳密ではないと参考文献55ページにあります。
ひとつ「こういう性質を満たす推定のやり方がよいのでは」という基準が得られました。さらに考えてみます。一致性は「データのサイズを無限にとばしたら」という状況での推定値のふるまいへの期待でした。でも、実際データのサイズは有限個しかないので、パラメータの推定値は真の値からずれると思います。真の値より大きい方にずれるかもしれないし、小さい方にずれるかもしれません。…こんなことを考えてみます。既にパラメータがわかっている分布からデータを例えば100個だけ生成して、AさんとBさんとCさんの推定のやり方によってパラメータの推定値を出してもらいます。次にその100個は捨ててしまって、改めて100個のデータを生成して、また推定値を出してもらいます。またその100個も捨てて、何度も繰り返します。このとき、推定値が真のパラメータから大きい方にずれたり小さい方にずれたりするのは仕方ない気がしますが、やたら大きい方にばかりずれるとか、やたら小さい方にばかりずれるとかだと、「なんだか偏りのある推定のやり方だなあ」ということになってしまうと思います。データを再生成する度に推定値がずれるのは仕方ないにしても、それをたくさん繰り返した推定値たちの平均値を取ったらパラメータの真の値になっていてほしいような気がします。これが満たされるのはパラメータの推定量の期待値が真の値に一致するときで（※）、このとき推定量が不偏性をもつといいます。
※ 何度も繰り返した推定値たちの平均値がパラメータの真の値になるかを確かめようとして、実際に何度も推定を繰り返させたらAさんもBさんもCさんも怒ってしまうと思いますが、その必要はありません。推定量 $\hat{\theta}(X)$ をデータ $X=(X_1, X_2, \cdots, X_n)$ の式で表して、ありうるデータ全体に対して期待値 $E[\hat{\theta}(X)]$ を取ればすみます。この期待値は当然、真の分布上での期待値です。
さっきの段落には少しまやかしがありました。それに気付いたあなたは「推定値が偏らないことってそんなに大事？ 大事なのは推定値が真のパラメータに近いかどうかでしょ？」と思うかもしれません。実際、ある分布のパラメータ推定においては、不偏推定量の中で真のパラメータとの平均2乗誤差が最も小さい推定量よりも、不偏性をもたないある推定量の方が真のパラメータとの平均2乗誤差が小さくなります（参考文献70ページコラム）。不偏性にこだわってかえって誤差が大きくなるかもしれないのでは目的を見失っているような気がします。それでもしばしば不偏性が重視されるのには理由があります。あるときあなたは仕事である分布のパラメータの推定量を考えることになりました。あなたの上司は「真のパラメータとの平均2乗誤差が最小になるような推定量を考えるように」と勝手なことをいいました。あなたはとりあえず真のパラメータ $\theta$ と推定量 $\hat{\theta}$ との平均2乗誤差を書き出してみました。次のようになりました。
E_{\theta}[(\hat{\theta} - \theta)^2 ] = E_{\theta}[\hat{\theta}^2] - 2\theta E_{\theta}[\hat{\theta}] + \theta^2 = V[\hat{\theta}] + E_{\theta}[\hat{\theta}]^2 - 2\theta E_{\theta}[\hat{\theta}] + \theta^2 = V[\hat{\theta}] + \left( E_{\theta}[\hat{\theta}] - \theta \right)^2
つまり、平均2乗誤差は推定量の分散 $V[\hat{\theta}]$ と「パラメータの推定量の期待値と真の値の差」の2乗 $\left( E_{\theta}[\hat{\theta}] - \theta \right)^2$ の和であることがわかりました。このうち、$V[\hat{\theta}]$ の方は分布のパラメータを推定しようという状況では絶対にゼロになりません（もし $\hat{\theta}$ の分散がゼロだったら、もう生成されたデータがどんなであろうと、すでに $\hat{\theta}$ は決まっていることになってしまいます）。他方、$E_{\theta}[\hat{\theta}] - \theta$ は上手く推定量を設計すればゼロにすることができるかもしれません。もちろん $V[\hat{\theta}] + \left( E_{\theta}[\hat{\theta}] - \theta \right)^2$ 全体を最小化するような推定量 $\hat{\theta}$ を考えるのが理想的ですが、それはとても難しいので、$E_{\theta}[\hat{\theta}] - \theta = 0$ であるような推定量 $\hat{\theta}$ の中で（つまり、不偏推定量の中で）分散 $V[\hat{\theta}]$ が最小になるものを探すことにしました。しかも、よく調べると、$\hat{\theta}$ が不偏推定量のときは分散 $V[\hat{\theta}]$ の理論的な下限値がわかっていることがわかりました（つまりここでは、そのまま平均2乗誤差の理論的下限値になります）。これは朗報で、なぜならあなたの上司は「この推定量で本当に平均2乗誤差が最小なの？ もっと小さくできるんじゃないの？」とかしつこくいいそうですが、理論的な下限値を達成していさえすれば「いえ、不偏推定量の誤差はこれが理論下限値です」などと返せば上司を退けることができそうだからです。不偏推定量であってこの理論的な下限値を達成しているとき、この推定量は有効性をもつといいます。
「不偏性と有効性については大人の事情じゃないか」と思うかもしれません。大人は大変なのです。限られた時間でまあまあな答えをみつけて、それをきちんと説明しなければなりません。そのとき、対象を不偏推定量にしぼって有効推定量を目指すことは役に立ちます。逆に、他の方法で誤差を議論できるならば不偏性と有効性は必ずしも要りません。でもそれは一般にとても難しいです。真のパラメータはやはり知り得ないからです。
クラメール・ラオの下限
先ほど、$\hat{\theta}(X)$ が不偏推定量の場合はその分散（不偏推定量なのでイコール真のパラメータとの平均2乗誤差）の理論下限値（クラメール・ラオの下限といいます）がわかっているという話がありました。その不等式の証明自体は参考文献の76～77ページにありますが、「分散はこれより小さくならない」という下限値がわかっているなんてちょっと不思議な気がします。また、証明にスコア関数（対数尤度のパラメータに関する偏微分）を用いるのだというのもなんだか天から降ってきたような気がします。なので、ちょっと一から考えてみたいと思います。この節の話はこの記事の他の節に輪をかけて厳密さを欠きますし、全然すっきりもしません。蛇足です。
といっても、分散の下限値を一から考えるのはやっぱり大変です。ただ、分散の不等式というと、2つの確率変数の共分散の不等式 ${\rm Cov(A, B)} \leqq \sqrt{V(A)} \sqrt{V(B)}$ を思い出します（これ自体はシュワルツの不等式の積分系から導かれます）。これの $B$ の方に $\hat{\theta}(X)$ を当てはめて式変形すれば、
$$ \frac{{\rm Cov}(A, \hat{\theta}(X))^2}{V(A)} \leqq V(\hat{\theta}(X))$$ こうなって「下限値」の形はできるので、あとは左辺が教科書と同じになってほしいということになります。では確率変数 $A$ をどのように選べばいいのでしょうか。明らかに $A$ は、定数とか $\hat{\theta}(X)$ と独立な変数とかであってはならないです。もしそうなら $\hat{\theta}(X)$ との共分散はゼロとなり、上の不等式は「分散はゼロ以上である」という、全く何もうれしくない不等式にしかなりません。さらに、クラメール・ラオの下限は「不偏推定量をどのようにとっても、分散はこれより小さくならない」という下限なので、上式の左辺は推定量 $\hat{\theta}(X)$ を含んでいてはならないし、であれば $A$ も ${\rm Cov}(A, \hat{\theta}(X))$ も推定量を含んでいてほしくないです（というのは厳密ではないですが、少なくともそうであればうれしいので、そういうことにします）。となると $A$ が依存できるのが $X$ と $\theta$ くらいしかないので、$A(X, \theta)$ とかきます。これより、上式の左辺の分子の共分散は以下のようになります。ここで、簡単のため $A(X, \theta)$ は平均がゼロになるようにシフトさせておくものとします。
$${\rm Cov}(A(X, \theta), \hat{\theta}(X)) = \int A(x, \theta)\left( \hat{\theta}(x) - \theta \right) f(x; \theta) dx = \int A(x, \theta) \hat{\theta}(x) f(x; \theta) dx - \theta \int A(x, \theta) f(x; \theta) dx = \int A(x, \theta) \hat{\theta}(x) f(x; \theta) dx$$ これをみると、$\hat{\theta}(X)$ の積分が含まれています。しかし、共分散に $\hat{\theta}(X)$ は残ってはなりません。いま、$\hat{\theta}(X)$ は不偏推定量なので、期待値さえとってしまえば真のパラメータ $\theta$ になります。上式の最右辺の積分から $A(x,\theta)$ が外に出て行ってほしいです。しかし、外に出て行ける、つまり $A(x,\theta)=A(\theta)$ ならば、結局それは定数 $0$ ということになってしまいます。
であれば、$A(x,\theta)$ は、$f(x;\theta)$ を割り消して、かつ、もう一度 $f(x;\theta)$ を出してくれるものであってほしいです。といっても $\displaystyle A(x,\theta) = \frac{f(x;\theta)}{f(x;\theta)}$ であってはなりません。それでは定数になってしまいす。 $\displaystyle A(x,\theta) = \frac{D \left(f(x;\theta) \right)}{f(x;\theta)}$ であって、$D$ が積分 $\displaystyle \int dx$ と交換する演算子であってほしいです。それなら共分散は $D (\theta)$ になります。結局 $D$ は $\theta$ に関する偏微分くらいしかなく、 $ \displaystyle D=\frac{\partial}{\partial \theta}$ となり、$\displaystyle A(x,\theta) = \frac{\frac{\partial }{\partial \theta}f(x;\theta)}{f(x;\theta)} = \frac{\partial }{\partial \theta} \log f(x;\theta)$ であり（この対数尤度のパラメータに関する偏微分をスコア関数とよびます）、上の共分散は1であり、$\hat{\theta}(X)$ の分散の下限値は以下であったことがわかります（とてもふわふわした話で申し訳ありませんでした）。
$$ \left[ V\left( \frac{\partial }{\partial \theta} \log f(x;\theta) \right) \right]^{-1} \leqq V(\hat{\theta}(X))$$
推定のやり方にはどんなものがあるの
参考文献にあるものについてダイジェストだけかきます。勝手にかいている内容もあります。「適当な条件」が肝心ですが、参考文献や、その他の文献を参照してください。
- 最尤推定量： 「尤度関数を最大にすればいいのでは」という考えのもと、尤度関数を最大にするパラメータを推定値とする。尤度関数が最大値をもたない場合は推定値を得ることができない。適当な条件下では一致性をもつ。一般に不偏性をもたない。
- モーメント推定量： 「モーメントを一致させればいいのでは」という考えのもと、パラメータをモーメントの式で表し標本モーメントを代入したものを推定値とする（標本平均を平均パラメータの推定値とするのは、モーメント推定量の最も簡単な例であるといえる）。モーメントが定義されない分布については推定値を得ることができない。一致性をもつ。一般に不偏性をもたない（標本分散が不偏性をもたないことはあまりに有名である）。
- 最小2乗推定量： 「誤差の2乗和を最小化すればいいのでは」という考えのもと、回帰モデルの誤差の2乗和を最小にするパラメータを推定値とする。真の値が発散する点がなければ推定値は得られるが、誤差系列自体が何らかの構造をもっている場合にはデータ間に依存関係がないのにゼロでない回帰係数が推定されてしまう場合がある（Ex. 見せかけの回帰）。適当な条件下で一致性と不偏性をもち、かつ、データの被説明変数の線形和で表せる推定量のうちでは最小の真のパラメータとの平均2乗誤差を達成する（最良線形不偏推定量；BLUE）。
おまけ
ガンマ分布という連続確率分布があると思います。その分布関数は $0$ 以上の実数の値を取る $x$ に対して以下のように定義されていると思います（確率分布の台が $x \in [0, \infty)$ ）。
f(x ; \alpha, \beta) = \frac{\beta^\alpha}{\Gamma(\alpha)} x^{\alpha - 1} \exp(-\beta x)
$\alpha$ と $\beta$ はガンマ分布の形を決めるパラメータで $\alpha > 0, \; \beta > 0$ であり、$\alpha$ の方は形状パラメータなどとよばれていると思います（もう一方のパラメータは $\beta$ ではなくて $1/\beta$ や $\alpha/\beta$ とする流儀もあると思います）。ガンマ分布についてちょっと心配になるのが、$x=0$ で確率分布がぶつりと終わっているところです。分布関数の $x \to 0$ の極限を取ると $\alpha$ の値に応じて以下のようになります。形状パラメータ $\alpha$ が $0 < \alpha < 1$ のとき、$x \to 0$ でこの分布関数は発散することがわかります（式の中に $x^{\alpha - 1}$ というのがあるのでそうなるだろうなとは思います）。
f(x ; \alpha, \beta) \underset{x \to +0}{\longrightarrow} \left\{ \begin{array}{ll} 0 & (1 < \alpha) \\ \beta & (\alpha = 1) \\ +\infty & (0 < \alpha < 1) \end{array} \right.
ここで、手元のデータにガンマ分布を仮定して $\alpha$ と $\beta$ を推定することを考えます。$\alpha$ の値によっては分布関数が $x=0$ で発散するので、最尤推定をするのは心配な気がします。でも、発散するのは $x=0$ においてのみなので、手元のデータに 0 が混ざっていなければ少なくとも計算上は大丈夫かもしれません。ただあいにく、手元のデータには 0 が多く混ざっています。もし 0 が混ざっていないとしても、データの中に小さな値があるとなんだか心配です。尤度がいきなり振り切れて、終了すると思います。
どうしてこんなことになってしまったのだろうと困惑します。手元のデータの広がりに近いガンマ分布というものは何かしらあるような気がします。でも、尤度はいきなり振り切れてしまいます。手元のデータに混ざっている 0 は何か悪いのでしょうか。…と考えると、どちらかというと尤度が悪い気がしてきました。尤度というのは各データに対応する確率密度の積の形をしていますが、確率密度をそのまま掛け合わせているのがよくないと思います。確率密度は積分しなければ確率になりません。$x=0$ での密度をそのまま取ろうとするから inf とかいうことになるのであって、微小区間でも積分すればそうならなかったはずです。なのに尤度は「データがそのパラメータから生成されている尤もらしさ」とかよくいわれている気がします。尤もらしさとは何なのでしょうか。なんだか何もわからなくなってきました。
ちなみに参考文献の64ページに、ガンマ分布にはモーメント法が有効であるとかいてありました。