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

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

標本平均，標本分散からの近似値推定
最尤推定
- 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}$