はしくれエンジニアもどきのメモ

情報系技術・哲学・デザインなどの勉強メモ・備忘録です。

ガンマ分布のパラメータ推定

ガンマ分布のパラメータ推定

ガンマ分布のパラメータ推定メモ. 以下の手法のメモ

  • 標本平均,標本分散からの近似値推定
  • 最尤推定
    • shapeパラメータの近似値推定
    • それを初期値として数値解析的にニュートン・ラフソン法で求める

ガンマ分布

shapeパラメータを$\alpha$, inverse scaleパラメータ$\beta$として, ガンマ分布を$Ga(\alpha, \beta)$で表せるとすると,その密度関数pdfは次のようになる.


\begin{eqnarray}
p(x|\alpha, \beta)
= \mathrm{Gamma}_{x}(\alpha, \beta)
= \frac{\beta^{\alpha}}{\Gamma(\alpha)} x^{\alpha-1}e^{- \beta x }, x > 0
\end{eqnarray}

標本平均,標本分散からの近似値推定

参考: http://web.sfc.keio.ac.jp/~maunz/BS14/BS14-06_GC.pdf

shapeパラメータを$\alpha$, inverse scaleパラメータ$\beta$のガンマ分布を$Ga(\alpha, \beta)$とすると, 平均と分散は,


\begin{eqnarray}
\mu &=& \frac{\alpha}{\beta} \\\\
\sigma^{2} &=& \frac{\alpha}{\beta^{2}} \\\\
\end{eqnarray}

この2式を連立方程式としてみる.第2式を変形すると,


\alpha = \sigma^{2} \beta^{2}

第1式に代入する.


\beta = \frac{\mu}{\sigma^{2}}

これを$\alpha$の式に代入すれば


\alpha = \frac{\mu^{2}}{\sigma^{2}}

まとめると,平均をデータの標本平均,分散をデータの標本分散に置き換えれば 各パラメータの近似が求められる.標本の大きさ$N$のデータを$x_{1}\cdots x_{N}$とすると,


\begin{eqnarray}
\left\{
\begin{array}{l}
\alpha &=& \frac{\mu^{2}}{\sigma^{2}}\\
\beta &=& \frac{\mu}{\sigma^{2}}\\
\mu &\fallingdotseq& \frac{1}{N} \sum_{i=1}^{N} x_{i} \\
\sigma^{2} &\fallingdotseq& \frac{1}{N-1} \sum_{i=1}^{N} (x_{i} - \mu)^{2} \\
\end{array}
\right.
\end{eqnarray}

最尤推定

参考:https://en.wikipedia.org/wiki/Gamma_distribution#Statistical_inference

尤度関数を求める.


\begin{eqnarray}
L(\alpha,\beta )
&=& \prod_{i=1}^{N} f(x_{i}; \alpha,\beta) \\\\
&=& \prod_{i=1}^{N} \frac{\beta^{\alpha}}{\Gamma(\alpha)}
x_{i}^{\alpha-1}e^{- \beta x_{i} } \\\\
\end{eqnarray}

対数尤度関数を求める.


\begin{eqnarray}
\ell(\alpha, \beta)
&=& \ln{(\beta^{N \alpha})} -\ln{(\Gamma(\alpha)^{N}) }
+ \sum{\ln{x_{i}^{\alpha-1}}}
+ \ln{e^{-\beta\sum{x_{i}}}}\\\\
&=& N\alpha \ln{(\beta)} - N\ln{(\Gamma(\alpha)) }
+ (\alpha-1)\sum{\ln{x_{i}}}-\beta\sum{x_{i}}\\\\
\end{eqnarray}

inverse scaleパラメータを推定

inverse scaleパラメータ$\beta$を求める. inverse scaleパラメータ$\beta$を微分すると,


\frac{\partial \ell(\alpha, \beta)}{\partial \beta}
=
N\alpha \frac{1}{\beta} - \sum{x_i}

0とおく.


\begin{eqnarray}
N\alpha \frac{1}{\beta} - \sum{x_i}&=& 0\\\\
\beta \sum{x_i} &=& N\alpha \\\\\
\hat{\beta}
&=& \frac{N\alpha}{\sum{x_i}}
= \frac{1}{\frac{1}{N\alpha}\sum{x_i}}
= \frac{\alpha}{\mu}\\\\
\end{eqnarray}

inverse scaleパラメータ$\beta$はclosed formになる. が,shapeパラメータ$\alpha$は解析解にならないので厳密解を求めるのは難しい.

shapeパラメータ$\alpha$の推定

対数尤度にinverse scale paramterの推定値$\hat{\beta}$を代入する.


\ell (\alpha, \beta )
=N\alpha \ln{\left( \frac{1}{\frac{1}{N\alpha}\sum{x_i}} \right)}
-N\ln{(\Gamma(\alpha)) }
+ (\alpha-1)\sum{\ln{x_{i}}}
- N\alpha

対数尤度をshapeパラメータ$\alpha$で微分する.

ここで,digamma関数が出てくるので次のように定義する.


\psi(\alpha) := \frac{\partial \ln{\Gamma(\alpha)}}{\partial \alpha}

\begin{eqnarray}
\frac{\partial \ell(\alpha, \beta)}{\partial \alpha}
&=& N
\left(
\ln{\left( \frac{1}{\frac{1}{N\alpha}\sum{x_i}} \right)}
+\alpha \frac{\sum{x_i}}{N\alpha}\frac{N}{\sum{x_i}}
\right) - N\psi(\alpha)
+\sum\ln{x_{i}}-N \\\\
&=& N
\left(
\ln{\left( \frac{1}{\frac{1}{N\alpha}\sum{x_i}} \right)}+1 \right) - N\psi(\alpha)
+\sum\ln{x_{i}}-N \\\\
&=& N
\ln{\left( \frac{1}{\frac{1}{N}\sum{x_i}} \right)}
+ N \ln{\alpha} - N\psi(\alpha)
+\sum\ln{x_{i}} +(N-N) \\\\
\end{eqnarray}

\begin{eqnarray}
N
\ln{\left( \frac{1}{\frac{1}{N}\sum{x_i}} \right)}
+ N \ln{\alpha} - N\psi(\alpha)
+\sum\ln{x_{i}} = 0\\\\
\ln{(\alpha)}- \psi(\alpha)
&=& -\ln{\left( \frac{1}{\frac{1}{N}\sum{x_i}} \right)} -\frac{1}{N}\sum\ln{(x_{i})} \\\\
\ln{(\alpha)}- \psi(\alpha)
&=&
\ln{\left( \frac{1}{N}\sum{x_i} \right)} -\frac{1}{N}\sum\ln{(x_{i})} \\\\
\end{eqnarray}

この等式が成立するときの$]alpha$が解になり,closed formがない. なので次の( minkaの近似式 )を利用して解く.


\ln{(\alpha)}- \psi(\alpha)
\approx
\frac{1}{2\alpha}
\left( 1 + \frac{1}{6\alpha+1} \right)

そして観測値$x$から次の定数を定義する.


s:=
\ln \left({\frac {1}{N}}\sum _{i=1}^{N}x_{i}\right)-{\frac {1}{N}}\sum_{i=1}^{N}\ln(x_{i})

minkaの近似式に代入すると次の近似式が成立する.


s \approx
\frac{1}{2\alpha}
\left( 1 + \frac{1}{6\alpha+1} \right)

$\alpha$で解く.2次方程式ができる.


\begin{eqnarray}
2\alpha s &\approx& \left( 1 + \frac{1}{6\alpha+1} \right) \\\\
2\alpha s(6\alpha+1) &\approx& (6\alpha+1) + 1 \\\\
12\alpha^{2}s + (2s-6)\alpha -2 &\approx& 0 \\\\
6s\alpha^{2} + (s-3)\alpha -1 &\approx& 0 \\\\
\end{eqnarray}

解の公式より


\frac{-(s-3)\pm \sqrt{(s-3)^{2} + 24s}}{12s}

2つ解が出てくるが,$s>0$より,$\sqrt{(s-3)^{2} + 24s}>0$そのため,+のみが解になる.


\alpha \approx \frac{-(s-3) + \sqrt{(s-3)^{2} + 24s}}{12s}

これはあくまで近似式なので,精度高くを求める方法として, これを初期値としてニュートン・ラフソン法で解く方法がある.

次の方程式を考える.


f(\alpha)=(\ln{(\alpha)}- \psi(\alpha))-s=0

$(\alpha,0)$を通る接線方程式を考えると更新式は,


\begin{eqnarray}
0-f(\alpha_0)
&=&
f'(\alpha_0)(\alpha-\alpha_0)\\
\alpha
&=&\alpha_0+\frac{-f(\alpha_0)}{f'(\alpha_0)}\\
\alpha&\leftarrow&\alpha_0-\frac{\ln{(\alpha_0)}- \psi(\alpha_0)-s}{\frac{1}{\alpha_0}-\psi'(\alpha)}\\
\end{eqnarray}

$\psi'(\alpha)$はdigamma関数の微分になる.