一般化ガンマ分布の数式の学習記録

【学習動機】

生存時間解析のパラメトリックモデルを見ていると、一般化ガンマ分布というものが現れるのだが、生存時間解析の和書や日本語資料を読んでもよく分からない。式の形も分かりにくいし、ソフトウェアを見てると2種類の確率密度関数があって、何が何か分からない。今回は、数式関連で英語で調べて分かったことをまとめておく。次回はシミュレーションして推定の様子を見る。

【学習内容】

対数正規分布とワイブル分布を包含する一般化として、Stacy(1962)でGeneralized Gamma Distribution(GGD)が出てきた。しかし、最尤法が上手く収束しない問題があった。とくに、対数正規分布に近いパラメータのところ(\(k \rightarrow \infty\)あたり)で上手くいかなかった。

GGDに従う確率変数の対数の分布を調べたのがPrentice(1974)で、こちら側からパラメータ化すると、推定も安定したようだ。こういう経緯があるため、こちらのパラメータ化したものGGDはGeneralized log-gamma distributionとも呼ばれる。

Rの{flexsurv}パッケージで言うと、Stacy(1962)のパラメータ化したGGDがgengamma.origに該当して、Prentice(1974)のものがgengammaに該当する。後者の方が、より一般化されていて、広い分布がとれる({flexsurv}のgengammaの\(Q<0\)の部分。後で説明する式で言う\(\lambda\))。

Stacy(1962)の定式化では、確率変数\(T\)が一般化ガンマ分布に従う(\(T\)～GG( \(\alpha\) , τ , \(k\) ))とき、確率密度関数は、

f:id:moratoriamuo271:20210405224218p:plain

と表される。一般化ガンマ分布は特殊形として、指数分布、ワイブル分布、ガンマ分布、対数正規分布を含む。

\(k= τ =1のときは指数分布\)

\[ f = \frac{1}{\alpha}(\frac{t}{\alpha})^0 exp[-(\frac{t}{\alpha})] =\frac{1}{\alpha} exp[-(\frac{t}{\alpha})] \]

\(k=1のときはワイブル分布\)

\[ f = \frac{τ}{\alpha}(\frac{t}{\alpha})^{τ -1} exp[-(\frac{t}{\alpha})^τ] \]

τ =1のときはガンマ分布

f:id:moratoriamuo271:20210405224327p:plain

\(k \rightarrow \infty\)のときは対数正規分布となる

この式は未確認。

パラメータを\((\alpha,τ,k)\)を以下のように、\((μ,σ,\lambda)\)と書き換える。

\[ μ = log(\alpha) + \frac{1}{τ} log(k) \] \[ σ = \frac{1}{τ \sqrt{k}} \] \[ \lambda = \frac{1}{\sqrt{k}} \]

として、\(T\)の対数をとって変数変換やらなんやらで操作して、元に戻すとPrentice(1974)の定式化した確率密度関数が得られる。この\(\lambda\)はRの{flexsurv}パッケージの\(Q\)に該当する。

\[ f(t;\lambda,σ,μ)= \left\{ \begin{array}{ll} \frac{|\lambda|}{σ t \Gamma(\lambda^{-2})} (\lambda^{-2})^{\lambda^{-2}} exp[\lambda^{-2}[(\frac{log(t)-μ}{σ})\lambda - exp[(\frac{log(t)-μ}{σ})\lambda]]] & (\lambda \neq 0) \\ \frac{1}{\sqrt{2\pi}σ t}exp[-\frac{(log(t)-μ)^2}{2σ^2}] & (\lambda=0) \end{array} \right. \]

指数分布となるのは\(\lambda=σ=1\)のとき、ワイブル分布となるのは\(\lambda=1\)のとき、ガンマ分布となるのは\(\lambda=σ\)のとき、対数正規分布となるのは\(\lambda=0\)のとき。