対数ロジスティック分布
はじめに
統計検定一級（理工学）の為、「一般化線形モデル入門 原著第２版（Annette J. Dobson著)」を使って生存時間解析の勉強をしていたら、「対数ロジスティック分布」という聞いた事があるようでよく知らない確率分布が出現。
ネットで調べてみるも日本語のサイトが全くと言っていいほど見つからず、結局英語版WikipediaやらMATLAB公式サイトの記事等を使って手計算するしかなかった。
以下、備忘録として途中式を中心に記録を残しておきます。
対数ロジスティック分布（log-logistic distribution）とは
対数を取ると「ロジスティック分布」に従うような分布を「対数ロジスティック分布」という。確率変数が非負の値を取り、かつ密度関数の裾が対数正規分布よりも厚いことから、生存時間解析の分野で使われることが多い（らしい）
ロジスティック分布（logistic distribution）
正規分布と同じく釣鐘型の密度関数を持つ連続型確率分布。正規分布よりも裾が厚い（平均から離れても下がりにくい）。
パラメータ(μ,s)はそれぞれ位置 (location) パラメータと尺度(scale)パラメータと呼ばれ、期待値、分散はそれぞれ
\begin{align}
E[X] &= \mu \\
V[Y] &= \frac{\pi ^2 s^2}{3}
\end{align}
である。また、累積分布関数は
\begin{equation}
F_X(x) = \frac{1}{1 + e ^ {-\frac{(x-\mu)}{s}}}
\end{equation}
である。
以下のグラフは標準正規分布と標準ロジスティック分布（(μ,s) = (0,1)となるようなロジスティック分布）の確率密度関数を比較したものである。
累積分布関数(cumulative distribution function)
上で挙げたロジスティック分布の累積分布関数から計算することが出来る。
Xがロジスティック分布に従い、Yが対数ロジスティック分布に従うとすると
$$
log(Y) \overset{\mathrm{dist}}{=} X \Leftrightarrow Y \overset{\mathrm{dist}}{=} exp(X)
$$
であることから、Yの累積分布関数は
\begin{align}
F_Y(y) &= P(Y \leq y) \\
&= P(e^X \leq y) \\
&= P(X \leq log(y)) \\
&= F_X(log y) \\
&= \frac{1}{1 + exp\{ -\frac{(logy - \mu)}{s}\}} \\
&= \frac{1}{1 + exp\{(logy^{-\frac{1}{s}} + \frac{\mu}{s}\}} \\
&= \frac{1}{1 + y^{-\frac{1}{s}} e^{\frac{\mu}{s}}} \\
&= \frac{1}{1 + \left( \frac{e^{\mu}}{y} \right)^{\frac{1}{s}}} \\
&= \frac{1}{1 + \left( \frac{y}{e^{\mu}} \right)^{-\frac{1}{s}}} \tag{1}
\end{align}
となる。このままでも良いが、WikipediaやTadikamalla, Pandu R. (1980)では
$$\mu = log(\alpha), s = \frac{1}{\beta}$$
としてパラメータを変換した形で書いている。つまり累積分布関数は
\begin{align}
F_Y(y) &= \frac{1}{1 + \left( \frac{y}{\alpha} \right)^{-\beta} } \tag{2} \\
&= \frac{y^{\beta}}{y^{\beta} + \alpha^{\beta}} \tag{3}
\end{align}
と表すことが出来る。以下、この形式を定式として進める。
確率密度関数(probability density function)
上で得られた累積分布関数(2)式をyで微分すればよい。つまり、
\begin{align}
f_Y(y) &= - \frac{ \left( 1 + \left( \frac{y}{\alpha} \right)^{-\beta} \right)^{'} } {\left( 1 + \left( \frac{y}{\alpha} \right)^{-\beta} \right)^2 } \\
&= \frac{ \left( \frac{\beta}{\alpha} \right) \left( \frac{y}{\alpha} \right)^{-\beta - 1}} {\left( 1 + \left( \frac{y}{\alpha} \right)^{-\beta} \right)^2 } \\
&= \frac{ \left( \frac{\beta}{\alpha} \right) \left( \frac{y}{\alpha} \right)^{-\beta - 1} \left( \frac{y}{\alpha} \right) ^ {2 \beta}} {\left( 1 + \left( \frac{y}{\alpha} \right)^{-\beta} \right)^2 \left( \frac{y}{\alpha} \right) ^ {2 \beta}} \\
&= \frac{\left( \frac{\beta}{\alpha} \right) \left( \frac{y}{\alpha} \right)^{\beta - 1}} {\left( 1 + \left( \frac{y}{\alpha} \right)^{\beta} \right)^2 } \tag{4}
\end{align}
生存関数（survivor function）
生存時間解析における「生存関数（survivor function）」とは、Yを生存時間とした時、Yがある時間yを超える確率を意味している。つまり、
\begin{align}
S_Y(y) &= P(Y \geq y)\\
&= 1 - F_Y(y)
\end{align}
と定義される。
従って、対数ロジスティック分布の生存関数は(3)式より
\begin{align}
S_Y(y) &= 1 - \frac{y^{\beta}}{y^{\beta} + \alpha^{\beta}} \\
&= \frac{\alpha^{\beta}}{y^{\beta} + \alpha^{\beta}} \\
&= \frac{1}{1 + \left( \frac{y}{\alpha} \right) ^{\beta}} \tag{5}
\end{align}
となる。
ハザード関数（hazard function）
生存時間解析における「ハザード関数（hazard function）」とは、Yを生存時間とした時、時間yまでは生存しているという条件のもとでyとy + δyの間に死亡する確率の極限値（δy→０）を意味している。
すなわち、
\begin{align}
h(y) &= \frac{f(y)}{S(y)} \tag{6} \\
\end{align}
と定義される。また、$$f(y) = -\frac{dS(y)}{dy}$$であることから、
\begin{align}
h(y) &= -\frac{S(y)^{'}}{S(y)} \\
&= -\frac{d}{dy} (log[S(y)]) \tag{7}
\end{align}
という表し方も出来る。
対数ロジスティック分布のハザード関数は、(6)(5)(4)式より
\begin{align}
h_Y(y) &= \frac{f_Y(y)}{S_Y(y)} \\
&= \dfrac{\frac{\left( \frac{\beta}{\alpha} \right) \left( \frac{y}{\alpha} \right)^{\beta - 1}} {\left( 1 + \left( \frac{y}{\alpha} \right)^{\beta} \right)^2 }}{\frac{1}{1 + \left( \frac{y}{\alpha} \right) ^{\beta}}} \\
&= \frac{\left( \frac{\beta}{\alpha} \right) \left( \frac{y}{\alpha} \right)^{\beta - 1}} {1 + \left( \frac{y}{\alpha} \right)^{\beta}}
\end{align}
となる。
k次モーメント
k次のモーメントは、上で求めた確率密度関数を用いて以下のように計算できる。
\begin{align}
E[Y^k] &= \int^\infty_0 y^k f(y) dy \\
&= \int^\infty_0 y^k \frac{\left( \frac{\beta}{\alpha} \right) \left( \frac{y}{\alpha} \right)^{\beta - 1}} {\left( 1 + \left( \frac{y}{\alpha} \right)^{\beta} \right)^2 } dy \\
\end{align}
$$ここで、t = \left( \frac{y}{\alpha} \right)^{\beta} とすると、$$
$$dt = \left( \frac{\beta}{\alpha} \right) \left( \frac{y}{\alpha} \right)^{\beta -1} dy$$
また自然数kを用いて
\begin{align}
t^{\frac{1}{\beta}} &= \frac{y}{\alpha}より、\\
y^k &= (\alpha t^{\frac{1}{\beta}})^k
\end{align}
であることから、
\begin{align}
E[Y^k] &= \int^\infty_0 (\alpha t^{\frac{1}{\beta}})^k \frac{1}{(t + 1)^2} dt \\
&= \alpha^k \int^\infty_0 \frac{t^{\frac{k}{\beta}}}{(t + 1)^2} dt \\
&= \alpha^k \int^\infty_0 \frac{t^{(1 + \frac{k}{\beta}) - 1}}{(t + 1)^{(1 + \frac{k}{\beta})+(1 - \frac{k}{\beta})}} dt \\
&= \alpha^k B(1 - \frac{k}{\beta}, 1 + \frac{k}{\beta}) \tag{8}
\end{align}
が得られる。ただし、$$B(・,・)$$はベータ関数。
（Wikipediaだとここから更に
$$ E[Y^k] = \alpha^k \frac{\frac{k \pi}{\beta}}{sin(\frac{k \pi}{\beta})}$$
という風に変形していますが仕組みがよく分かっていません... 分かる方は是非コメントで教えて下さい）
これにより、平均、分散はそれぞれ
\begin{align}
E[Y] &= \alpha B(1 - \frac{1}{\beta}, 1 + \frac{1}{\beta})\\
V(Y) &= E[Y^2] - E[Y] \\
&= \alpha^2 B(1 - \frac{2}{\beta}, 1 + \frac{2}{\beta}) - \alpha B(1 - \frac{1}{\beta}, 1 + \frac{1}{\beta}) \\
\end{align}
で求めることが出来る。
参考文献
・Annette J. Dobson 著、田中豊、森川敏彦、山中竹春、冨田誠 訳「一般化線形モデル入門 原著第２版」
https://www.amazon.co.jp/%E4%B8%80%E8%88%AC%E5%8C%96%E7%B7%9A%E5%BD%A2%E3%83%A2%E3%83%87%E3%83%AB%E5%85%A5%E9%96%80-%E5%8E%9F%E8%91%97%E7%AC%AC2%E7%89%88-Annette-J-Dobson/dp/4320018672
・Tadikamalla, Pandu R. (1980), A Look at the Burr and Related Distributions, International Statistical Review, 48 (3)
・Log-logistic distribution - Wikipedia
・対数ロジスティック分布 - MATLAB & Simulink - MathWorks 日本
