スキップしてメイン コンテンツに移動

最尤推定

Introduction


今日は最尤推定について加工と思います。これは統計的推定でよく使われる手法です。最尤推定の例も書こうと思います。初めに尤度の説明をし、そのあとで最尤推定の説明をします。

概要


  • 尤度
  • 最尤推定
  • 最尤推定の問題点

尤度

前提条件から得られる観察データを考えます。この時、えられた観測データに対して前提条件が尤もらしい条件であるかの値を尤度といいます。

なにをゆっているのかわからない人がほとんどだと思います。。。
尤度の例を扱っていきます。

コインを投げることを考えます。このコインは確率Pで表、確率1-Pで裏を出すコインだとします。
例えば、100回コインを投げたとき、全て表だったとします。この時このコインが表が出る確率はかなり1に近いことが予想されます。

ではもし、表が出る確率PがP=0.5だとします。この時、表が100回連続で出る確率は$0.5^{100} = 7.88860e-31$になります。あり得ない確率ですね。これがP=0.5としたときのもっともらしさです。つまり、あまり現実的ではないということです。

もしP=0.99とするとき、100回とも表が出る確率は$0.99^{100} = 0.3666....$となります。つまり、P=0.99としたときの尤度は0.36くらいということです。よって、P=0.5よりかは現実見があることになります。まだまだ低い数字ではありますが。


観測データである、100回表が出るという事象を固定したとき、尤度はPを変数としたP(100回表|P)を尤度関数と呼びます。この関数の値を尤度と呼びます。


尤度が高いほうが尤もらしい値、つまり理にかなっているなと感じることができる値ということになります。

例えば、先ほどの例でゆうと、

P=0.5としたときの尤度は7.88860e-31でした。P=0.99としたときの尤度は0.3666でした。よってP=0.5より、P=0.99のほうが尤もらしい自然な値ということになります。


最尤推定

最尤推定とは得られた観測データからデータが依存している分布のパラメーターを推測するための手法です。
最尤推定では尤度を最大化して、最も尤もらしいパラメーターを求めます。

確率密度関数$f$と$X_1,X_2,...,X_n$が$f$に従う確率変数とします。 $$X_1,X_2, ..., X_n \sim f$$

この時、$X_1,X_2,..,X_n$が同時におこる確率は
$$\Pi_{i=1}^{N} P(X_i)$$
いわゆる同時確率です。

よって、尤度関数を次のように定義します。
$$L(\theta) = f(x_1,x_2,...,x_n|\theta)$$

この時、
$$\theta^{\star} \in \arg_{\theta} \max L(\theta)$$
$\theta$最尤推定量といいます。

そして、
$$\frac{\partial}{\partial \theta} \log L(\theta)$$
これを尤度方程式といいます。

なぜ、$\log$がいきなり登場しているのかは後の最尤推定の例で説明します。


最尤推定の例

$x_1,x_2,...,x_n \in {0,1}$について考えます。 $\forall i \in {1,2,..,n}$について、$x_i = 1$とするとき、i回目のコイン投げは表とします。$x_i$とするとき、コインは裏になったとします。

この時、尤度関数は
$$L(\theta) = P(x_1,x_2,...,x_n|\theta) = \Pi_{i=1}^{n} \theta^{x_i} (1-\theta)^{1-x_i}$$

$\forall i \in {1,2,..,n},  \sim p(k;\theta) = \theta^k (1-\theta)^{1-k} ~~~~\textrm{for} k \in {0,1}$

コインはベルヌーイ分布に従うので、このような形になります。
ここで、 $\theta$表が出る確率とします。

$L(\theta)$を$\theta$について最大化したいのですが、微分がかなり難しい形になっています。なぜなら$L(\theta)$について線形でないからです。つまり、$\theta$について掛け算の形になっていることが微分を難しくしています。

この問題を解決するために、$\log$を使います。
$\log$は単調増加関数なので$L(\theta)$と$\log L(\theta)$の局所解は変わりません。

よって、$\log L(\theta)$を最大化します。

\begin{eqnarray*}
\log L(\theta) &=& \log \Pi_{i=1}^{n} \theta^{x_i} (1-\theta)^{1-x_i} \\
&=& \sum_{i=1}^N \log \theta^{x_i} + \log (1-\theta)^{1-x_i} \\
&=& \sum_{i=1}^N  x_i \log \theta + (1-x_i)\log(1-\theta)
\end{eqnarray*}

これの微分は

\begin{eqnarray*}
\frac{\partial}{\partial \theta} \log L(\theta) &=& 0 \\
\frac{\partial}{\partial \theta} \sum_{i=1}^N x_i \log \theta + (1-x_i) \log (1-\theta) &=& 0 \\
\sum_{i=1}^N \frac{x_i}{\theta} - \frac{1-x_i}{1-\theta} &=& 0 \\
\frac{1}{\theta} \sum_{i=1}^N x_i - \frac{1}{1-\theta_i} \sum_{i=1}^N (1-x_i) &=& 0 \\
(1-\theta) \sum_{i=1}^N x_i - \theta \sum_{i=1}^N 1-x_i &=& 0 \\
\sum_{i=1}^N x_i - \theta \sum_{i=1}^N x_i - \theta \sum_{i=1}^N 1 + \theta \sum_{i=1}^N x_i &=& 0 \\
\sum_{i=1}^N x_i - n \theta &=& 0 \\
\theta &=& \frac{\sum_{i=1}^N }{n} \\
\end{eqnarray*}

この最適解は$x_1,x_2,..,x_n$の平均を表していることがわかります。
もし、コインが100回表、裏が0回だった時、最尤推定により、得た$\theta$の値は$\theta=1$
また、コインが50回表、裏が50回であれば、最尤推定により、得られた$\theta$の値は$\theta = 0.5$

最尤推定の問題点

最尤推定には問題点もあります。例えば、先ほどみたとおり、100回表がでて、裏が0解の時、$\theta=1$とするのが尤もらしいという結果が最尤推定から得られました。

ただ、もし3回表がでて、0回裏が出たとき、この時の最尤推定量も$\theta=1$となってします。しかし、3回表が出たからと言って、コインの表が出る確率が1と考えるのはあまりに危険すぎます。


つまり、観測データが少ない場合は最尤推定の結果はあまり信用できません。


Reference

https://ja.wikipedia.org/wiki/%E5%B0%A4%E5%BA%A6%E9%96%A2%E6%95%B0

コメント

  1. 対数尤度をさらに変形していくと確率モデルの分布をデータの出現頻度の分布(経験分布)に近づけるって話にももっていけますよねー。

    返信削除
    返信
    1. コメントありがとうございます!
      すいません。勉強不足で、よろしければご教授していただけないでしょうか?

      削除
  2. このコメントは投稿者によって削除されました。

    返信削除

コメントを投稿

このブログの人気の投稿

Implementation of Robbins monro

Robbins monro の実装 sorry, this page is Japanese only.   今回はRobbins monro の実装をしてみました。 Robbins monroは確率勾配降下法の学習率を入りテーション回数の逆数で割っていくものです。 使っているprogram言語はpython 3です。osはwindowsです。(macほしい...) アルゴリズム 確率勾配降下方とは目的関数の最適解を求めるアルゴリズムです。目的関数をf(X)とすると、手順は以下のようになっています。 初期学習率$n_0$を決めます。訓練データDを用意します。この訓練データは複数の初期値の集まりです。 訓練データから一つ初期値をランダムに取り出し、これを$x_0$とし、最初の予測値とします。 次の式に現在の予測値$x_0$を代入し、新たな予測値$x_{n+1}$を得ます。$$x_{n+1} = x_{n} - \frac{n_0}{n} grad f(X_n)$$ 収束して入れば4へ、収束していなければ2で得られた値$x{n+1}$を新たに$x_n$としてもう一度2を行う。 訓練データを一周していなければ2へ、一周していれば各初期値から得られた解の中から目的関数を最も小さくするものを選ぶ。   実装例 以下の目的関数を最小化させてみましょう。 $$f(x,y) = (x-2)^2 + (y-3)^2 $$ コマンドラインでpythonを実行していきます。 予想通り、(2,3)という解を導き出してくれました。目的関数が簡単だったので、初期値をどの値でとってもばっちり正解にたどり着いてくれました。 CODE 以下にRobbins monroの関数だけ置いておきます。 こちら にすべてのコードを載せています。 def Robbins_monro(function,grad,number_variable_gradient): init_learning_rate = 1.5 stepsize = 1000 init_value = np.array([range(-1000,1020,20) for i in range(number_v...

カーネルSVMの実装

Introduction English ver 今日はカーネルSVMを実装しました。 最適化には内点法を実装しました。 この記事では実装のみ書いていますが、理論編もそのうちに書きたいと思います。書き次第、リンクを貼っておきます。 # 第一弾書きました。 SVMの理論 part 1 僕のcomputerはwindowsでOSもwindowsです。 実装はPython3で行います。 概要 カーネルとは データセットについて 実装の結果 カーネルとは カーネルとは非線形問題を解くための手法です。 直線で分類できないようなデータをある写像をもとで、直線で分類できるような分布に変形することです。次の画像を見たほうがわかりやすいかもしれません。 ⇓↓⇓↓⇓↓⇓↓⇓↓⇓↓⇓↓⇓↓⇓↓⇓↓⇓↓⇓↓⇓↓⇓↓⇓↓⇓↓⇓↓⇓↓⇓↓⇓↓⇓↓ このように変形することをカーネルトリックといいます。 変形されたデータは写像\(\phi\)によって\(\phi(x)\)と表現されます。 \[\phi:x -> \phi(x)\] しかし、SVMではカーネル関数は次のように表されることが多いです。 \[K(x,y)=\phi(x)^T \phi(y)\] なぜなら、SVMでは \(\phi(x)^T \phi\)のみを計算するからです。 最もよく使われているカーネルはRBFとpolynomialです。 RBF \[K(x,y) = \exp(-\gamma ||x-y||^2)\] polynomial \[K(x,y) = (x^Ty+c)^p\] \(K(x,y) = x^Ty\)とすればこれは線形SVMになります。 データセット 今日使うデータセットはnp.random.randnによって生成された正規分布に従う乱数です。 と このデータを作ったコードは ここ に置いてあります。 Implementation 今回はRBFカーネルを使いました。 \(\gamma = 0.5\) 正規化係数50とすると次のようになりました。 次に二つ目のデータセットに対して、RBF kernelを使いました。 \(\gamma = 0.5\) 正規化係数を100とすると、次のよ...

Graph theory in Python

Introduction 日本語 ver Today, I write the networkx module for graph theory in Python. I write about basis of Gragh theory. Please look at this page . This post is not written Graph theory, but I want to share module in Python for Graph theory. Networkx Python3 have Networkx module. If you have python2, you should install by pip pip install networkx Make Graph by Networkx first, import networkx module and make instance. import networkx as nx import matplotlib.pyplot as plt G = nx.Graph() next, add node and edge G.add_node(1) # add Multiple nodes G.add_nodes_from([2,3,4]) G.add_edge(1,2) # add Multiple edges G.add_edges_from([(3,4),(1,2),(4,6)]) next, plot G nx.draw(G) plt.show() networkx have many function. I will write new learning about networkx in this post. Reference https://qiita.com/kzm4269/items/081ff2fdb8a6b0a6112f http://akiniwa.hatenablog.jp/entry/2013/05/12/012459