RigelのR言語メモであーる(R言語だけとは言っていない)

RigelのR言語メモであーる(主にpython)

興味あることや趣味、やったことについて書くよ

最大エントロピー法によるガンマ分布の導出

前回、最大エントロピー法により正規分布を導出しました。
その際の制約条件を少し変えることで、以下のガンマ分布が出てくるそうです。
{ \displaystyle f(x)=\frac{x^{k-1}e^{-x/\theta}}{\Gamma(k)\theta^k} \ \ {\rm for}\ x>0}
今回は、このガンマ分布の導出に挑戦します。

前の記事・twitterのリプライ

前の記事はこちらです。
strawberry-kyon.hateblo.jp
平均と分散を固定してエントロピーを最大化し、正規分布を導出しました。
この記事をアップロードすると、以下の興味深いコメントを頂きました。


ガンマ分布は、平均と対数の平均を固定すると出てくるそうです。
ガンマ分布は、指数分布に従う確率変数の和と覚えていたので衝撃的でした。
ということで、今回は上記の制約条件のもとでエントロピーを最大化してみます。

エントロピーの最大化

変数は前回の記事と同じものを用います。
制約条件は以下の3つとします。
{\begin{eqnarray}
&&\int_{-\infty}^{\infty} p(x)dx=1\\
&&\int_{-\infty}^{\infty} xp(x)dx=\theta_1\\
&&\int_{-\infty}^{\infty} \log x\ p(x)dx=\theta_2
  \end{eqnarray}}
ここで、{ \displaystyle \log x}を扱うので{ \displaystyle p(x)=0 \ \ (x\leq 0)}としておきます(測度とか台とかよく分かっていない)。
このような制約条件のもとでエントロピーを最大化する密度関数{ \displaystyle p(x) }はなんだろうか?という問題です。
さっそく、前回と同様にラグランジュ関数は、
{\begin{eqnarray}
L&=&-\int_{-\infty}^{\infty} p(x)\log p(x)dx+\lambda_1 \left(\int_{-\infty}^{\infty}p(x)dx-1\right)\\
&&\ \ \ +\lambda_2 \left(\int_{-\infty}^{\infty}xp(x)dx-\theta_1\right)+\lambda_3 \left(\int_{-\infty}^{\infty} \log x\ p(x)dx-\theta_2 \right)\\
&=&\int_{-\infty}^{\infty}\left( -p(x)\log p(x) +\lambda_1 p(x) +\lambda_2 xp(x) +\lambda_3 \log x\ p(x) \right)dx-\lambda_1-\lambda_2\theta_1-\lambda_3\theta_2
  \end{eqnarray}}
となります。{ \displaystyle p(x)}偏微分しイコール0と置くことにより、
{\begin{eqnarray}
\frac{\partial L}{\partial {p(x)}}= -\log p(x) -1 +\lambda_1+\lambda_2 x+ \lambda_3 \log x =0
  \end{eqnarray}}
となりまして、{ \displaystyle p(x)}について解くと、
{\begin{eqnarray}
p(x)&=&\exp\left(-1+\lambda_1+\lambda_2 x+\lambda_3 \log x\right)\\
&=&\exp\left(-1+\lambda_1\right)x^{\lambda_3}\exp\left(\lambda_2 x \right)
  \end{eqnarray}}
となります。
係数部分の{ \displaystyle \exp\left(-1+\lambda_1\right)}は正規化定数なのでとりあえず無視しておきます。
ガンマ分布の形に合わせて、定数を
{\begin{eqnarray}
\lambda_2&=&-\frac{1}{\theta}\\
\lambda_3&=&k-1
  \end{eqnarray}}
と置くと、
{\begin{eqnarray}
p(x)&=&\exp\left(-1+\lambda_1\right)x^{\lambda_3}\exp\left(\lambda_2 x \right)\\
&=&\exp\left(-1+\lambda_1\right)x^{k-1}\exp\left(- x/\theta \right)
  \end{eqnarray}}
となります。あとは1つ目の制約条件から正規化定数を求めます。
つまり、
{\begin{eqnarray}
&&\int_{-\infty}^{\infty} p(x)dx\\
=&&\int_{0}^{\infty} \exp\left(-1+\lambda_1\right)x^{k-1}\exp\left(- x/\theta \right) dx\\
=&&\exp\left(-1+\lambda_1\right) \int_{0}^{\infty} x^{k-1}\exp\left(- x/\theta \right) dx\\
  \end{eqnarray}}
ガンマ関数の形に合わせるために、{ \displaystyle x/\theta=y}と変数変換して、
{\begin{eqnarray}
&&\exp\left(-1+\lambda_1\right) \int_{0}^{\infty} x^{k-1}\exp\left(- x/\theta \right) dx\\
=&&\exp\left(-1+\lambda_1\right) \int_{0}^{\infty} (\theta y)^{k-1}\exp\left(-y \right) \theta dy\\
=&&\exp\left(-1+\lambda_1\right) \theta^k \int_{0}^{\infty} y^{k-1}\exp\left(-y \right) dy\\
=&&\exp\left(-1+\lambda_1\right) \theta^k \int_{0}^{\infty} y^{k-1}\exp\left(-y \right) dy\\
=&&\exp\left(-1+\lambda_1\right) \theta^k \Gamma(k)=1
  \end{eqnarray}}
となるので正規化定数は、
{\begin{eqnarray}
\exp\left(-1+\lambda_1\right) =\frac{1}{\Gamma(k) \theta^k}
  \end{eqnarray}}
となります。よって、
{\begin{eqnarray}
p(x)&=&\exp\left(-1+\lambda_1\right)x^{k-1}\exp\left(- x/\theta \right)\\
&=&\frac{x^{k-1}e^{-x/\theta}}{\Gamma(k)\theta^k}\ \ {\rm for}\ x>0
  \end{eqnarray}}
と無事にガンマ分布が出てきました。

まとめ

前回に続き、今回はガンマ分布を導出しました。
パラメータは前回に比べるとズルして決めた感じがありますが、勘弁して下さい。
分布が扱いやすいようにパラメータを決めたと思うと良いですね。
いつか、{ \displaystyle \theta_1,\ \theta_2}{ \displaystyle k,\ \theta}の関係について書きます({ \displaystyle \theta_1}については明らか)。
それから、この辺りをもっと深く理解するためには相対エントロピー、sanovの定理、カノニカル分布、とかを勉強する必要がありそうです。
この記事が誰かの理解の助けになると嬉しいです。
間違いの指摘や質問等ありましたら、お問い合わせフォームやコメントから気軽にお願いします。

(記事のサムネイルはオイラー先生です。ガンマ分布を提案した人が分からなかったため、ガンマ関数を最初に導入したオイラー先生にしました。)