分布の識別アルゴリズム

本アプリでは、最尤推定法（MLE: Maximum Likelihood Estimation）を用いて各分布のパラメータを推定します。ただし、ガウス混合分布（Gaussian Mixture Distribution）については、EMアルゴリズム（Expectation Maximization）を使用します。

MLEメソッドでは、入力データ $\mathbf{y}$ 、パラメータベクトル $\theta$ について

まず、入力データの確率密度関数の積である分布の尤度関数 $\mathcal{L}_{n}( \theta \,;\mathbf {y} )$ を計算します。
次に、対数尤度関数 $\ell ( \theta \,;\mathbf {y} )$ の偏導関数をゼロにしてパラメータ推定します。
ニュートン・ラフソン法を使用して、ステップ2でこれらのパラメータを求めます。

分布

本アプリで対応している分布は以下の通りです。

分布	パラメータ	確率密度関数（PDF）	累積分布関数（CDF）	逆累積分布関数（Inverse CDF）
正規分布	位置: $\mu$ スケール: $\sigma$	$f(x)=\frac{1}{ \sqrt{2 \pi} \sigma } e^{ -\frac{( x - \mu )^2}{2 \sigma^2} }, \; \sigma > 0$	$P(x) = \Phi \left( \frac{x - \mu}{\sigma} \right)$	$x = \mu + \sigma \Phi^{-1}(p)$
対数正規	位置: $\mu$ スケール: $\sigma$	$f(x)=\frac{1}{ \sqrt{2 \pi} \sigma x } e^{ -\frac{( \text{ln}(x) - \mu )^2}{2 \sigma^2} },$ $x>0 , \; \sigma > 0$	$P(x) = \Phi \left( \frac{\text{ln}(x) - \mu}{\sigma} \right)$	$\text{ln}(x) = \mu + \sigma \Phi^{-1}(p)$
3パラメータ対数正規分布	位置: $\mu$ スケール: $\lambda$ しきい値: $\lambda$	$f(x)=\frac{1}{\sqrt{2 \pi} \sigma (x -\lambda)} e^{ -\frac{( \text{ln}(x -\lambda) - \mu )^2}{2 \sigma^2} },$ $x > \lambda , \; \sigma > 0$	$P(x) = \Phi \left( \frac{\text{ln}(x - \lambda) - \mu}{\sigma} \right)$	$\text{ln}(x - \lambda) = \mu + \sigma \Phi^{-1}(p)$
ガンマ	形状: $\alpha$ スケール: $\beta$	$f(x) = \frac{ x^{\alpha-1} e^{-\frac{x}{\beta}}}{ \beta^{\alpha} \Gamma(\alpha) },$ $x > 0, \; \alpha > 0, \; \beta > 0$	$P(x) = \frac{1}{ \Gamma(\alpha) } \gamma( \alpha, \frac{x}{\beta})$	$x = \beta \gamma^{-1}( \alpha, \Gamma(\alpha) p )$
指数分布	スケール: $\theta$	$f(x) = \frac{ 1 }{ \theta } e^{-\frac{x}{\theta}}, \; x > 0, \; \theta > 0$	$P(x) = 1- e^{ -\frac{x}{\theta} }$	$x = \theta ( - \text{ln}(1-p) )$
2パラメータ指数分布	スケール: $\lambda$ しきい値: $\lambda$	$f(x) = \frac{ 1 }{ \theta } e^{-\frac{x - \lambda}{\theta}}, \; x > \lambda, \; \theta > 0$	$P(x) = 1- e^{ -\frac{x - \lambda}{\theta} }$	$x = \lambda + \theta ( - \text{ln}(1-p) )$
最小極値分布	位置: $\alpha$ スケール: $\beta$	$f(x) = \frac{1}{\beta } e^{ \frac{x- \alpha}{\beta}} e^{ - e^{ \frac{x- \alpha}{\beta} } },$ $\beta > 0$	$P(x) = 1- e^{-e^{\frac{x-\alpha}{\beta}}}$	$x = \alpha + \beta \text{ln}( - \text{ln}(1-p) )$
ワイブル	scale: $\alpha$ shape: $\beta$	$f(x) = \frac{\beta}{\alpha} \left(\frac{ x }{ \alpha }\right)^{\beta - 1} e^{ -(\frac{ x }{ \alpha })^\beta },$ $x > 0, \; \alpha > 0, \; \beta > 0$	$P(x) = 1 - e^{ -(\frac{ x }{ \alpha })^\beta }$	$\text{ln}(x) = \text{ln}(\alpha) + \frac{1}{\beta} \text{ln}( - \text{ln}(1-p) )$
3パラメータワイブル	スケール: $\alpha$ 形状: $\beta$ しきい値: $\lambda$	$f(x) = \frac{\beta}{\alpha} \left(\frac{ x - \lambda }{ \alpha }\right)^{\beta - 1} e^{ -\left(\frac{ x - \lambda }{ \alpha }\right)^\beta },$ $x > \lambda, \; \alpha > 0, \; \beta > 0$	$P(x) = 1 - e^{ -\left(\frac{ x - \lambda }{ \alpha }\right)^\beta }$	$\text{ln}(x - \lambda) = \text{ln}(\alpha) + \frac{1}{\beta} \text{ln}( - \text{ln}(1-p) )$
最大極値分布	位置: $\alpha$ スケール: $\beta$	$f(x) = \frac{1}{\beta } e^{- \frac{x- \alpha}{\beta}} e^{ - e^{ -\frac{x- \alpha}{\beta} } },$ $\beta > 0$	$P(x) = e^{-e^{-\frac{x-\alpha}{\beta}}}$	$x = \alpha - \beta \text{ln}( - \text{ln}(p) )$
ロジスティック	位置: $\mu$ スケール: $\sigma$	$f(x) = \frac{e^{\frac{x- \mu}{\sigma}}}{\sigma(1+e^{\frac{x- \mu}{\sigma}})^2 },\sigma > 0$	$P(x) = \frac{1}{1+e^{-\frac{x- \mu}{\sigma}}}$	$x =\mu+ \sigma\text{ln}(\frac{p}{1-p})$
対数ロジスティック	位置: $\mu$ スケール: $\sigma$	$\begin{array}{l}f(x) = \frac{e^{\frac{\text{ln}(x)- \mu}{\sigma}}}{x\sigma(1+e^{\frac{\text{ln}(x)- \mu}{\sigma}})^2 } , \\ x>0, \sigma > 0\end{array}$	$P(x) = \frac{1}{1+e^{-\frac{\text{ln}(x)- \mu}{\sigma}}}$	$x =e^{\mu+ \sigma\text{ln}(\frac{p}{1-p})}$
折りたたみ正規分布	位置: $\mu$ スケール: $\sigma$	$\begin{array}{l}f(x) = \frac{1}{ \sqrt{ 2 \pi } \sigma} e^{ - \frac{(x - \mu)^2}{ 2 \sigma^2 } } \\ + \frac{1}{ \sqrt{ 2 \pi } \sigma} e^{ - \frac{(x + \mu)^2}{ 2 \sigma^2 } }, x \geqslant 0, \; \sigma > 0\end{array}$	$\begin{matrix}P(x) = \Phi \left( \frac{x + \mu}{\sigma} \right) \\ - \Phi \left( \frac{\mu - x}{\sigma} \right)\end{matrix}$	$x = \text{foldnorminv}(p, \, \mu, \, \sigma)$
レイリー	スケール: $\sigma$	$f(x) = \frac{x}{\sigma^2}e^{-\frac{x^2}{2\sigma^2}},x\geqslant0, \sigma > 0$	$P(x) = 1-e^{-\frac{x^2}{2\sigma^2}}$	$x=\sqrt{2}\sigma \sqrt{-\text{ln}(1-p)}$
混合ガウス	成分数: $g$ 、i番目のコンポーネントに対して重み: $w_i$ 位置: $\mu_i$ スケール: $\sigma_i$	$\begin{array}{l}f(x) = \sum_{i=1}^g \frac{w_i}{ \sqrt{2 \pi} \sigma_i } e^{ -\frac{( x - \mu_i )^2}{2 \sigma_i^2} }, \\ 0 < w_i < 1, \; \sigma > 0\end{array}$	$P(x) = \sum_{i=1}^g w_i \Phi \left( \frac{x - \mu_i}{\sigma_i} \right)$	明示的な解なし。補間によって計算することができます。

ここでΦは標準正規分布のCDF関数、γ は下不完全ガンマ関数です。

Note:

Box-Cox 変換、Johnson変換、Yeo-Johnson変換は、まずデータを変換し、正規分布をフィッティングします。詳しくはデータ変換アルゴリズムを参照して下さい。
正規分布と対数正規分布の尺度パラメータは不偏標準偏差として補正されます。すなわちnをポイント数とするとnの代わりにn-1で割ったものになります。
2パラメータ指数分布におけるしきい値パラメータ λ は以下のように推定されます。

$\lambda = X_{(1)} - \frac{\bar{X} - X_{(1)}}{n-1}$ ,

ここで、X₍₁₎はXの最小値で、 XはXの平均値です。

適合度の評価

Anderson-Darling検定

Anderson-Darling検定の仮説は以下の通りです。

H₀: 入力データは分布に従う。

H₁: 入力データは分布に従わない。

Anderson-Darling統計量

入力データはまずソートされ、各データ点に対する累積分布関数（CDF）値Z_iを計算します。Anderson-Darling統計量 A² を以下のように計算します。

$A^2 = -n -(1/n) \sum_{i=1}^n \left[ (2i-1) \text{ln} (Z_i) + ( 2n+1-i ) \text{ln} (1 - Z_i) \right]$

A²値を使用して適合度を評価できます。値が小さいほどデータは指定された分布に適合しやすくなります。

Anderson-Darling検定のp値

p値は、リファレンス[1]と[2]の表の統計量A²を補間することによって計算されます。p値が棄却値（例えば0.05）より小さい場合、データが指定された分布に適合しないと判断（H₀を棄却）されます。

3パラメータ対数正規分布のようないくつかの分布では、参照がなく、p値が欠落しています。

折りたたみ正規分布、ガウス混合分布及び固定形状パラメータの3パラメータワイブル分布では、モンテカルロ法によりp値を計算します。帰無仮説のもとで、パラメータを推定した分布から 1000個のサンプルを生成します。それぞれのサンプルに対してA²を計算し、値が元データ[4]の Anderson-Darling 統計量以上である回数をカウントします。

$\text{P-value} = \frac{b + 1}{N + 1}$ ,

ここで、bは入力データのA²以上の統計値の数、Nはサンプルの数、ここではN=1000です。

BIC

ガウス混合分布では、BIC値は次のように計算されます。

$\text{BIC} = -2 \, \ell + k \text{ln}(n)$ ,

ここでℓは、推定における尤度の対数、kはパラメータの数、nはポイント数です。

BIC値が小さいほどモデルの適合度が良いと解釈できます。

尤度比検定(LRT)

尤度比検定(LRT)は、2つの統計モデルの適合度を比較するための検定手法です。単純なモデルAと複雑なモデルBがあり、モデルAはモデルBのサブセットです。たとえば、モデルAはワイブル分布で、モデルBは3パラメータワイブル分布です。より複雑なモデル（B）が単純なモデル（A）よりも有意に優れているかどうかを判断します。このアプリでは、尤度比検定は対数正規分布vs 3パラメータ対数正規分布、指数分布 vs 2パラメータ指数分布、ワイブル分布 vs 3パラメータワイブル分布を比較することができます。

尤度比検定の仮説は以下のとおりです。

H₀: モデルBはモデルAと同じ適合度（より複雑なモデルは必要ない）

H₁: モデルBはモデルAよりも有意に優れている

尤度比検定のP値が棄却値（例えば0.05）より小さい場合、H₀を棄却し、複雑なモデルBが単純なモデルAよりも著しく優れていることを示します。

尤度比検定の統計量

$\text{LRT} = -2 \, ( \ell(A) - \ell(B) )$ ,

ここでℓ(A) と ℓ(B) は、それぞれモデルAとモデルBの対数尤度です。

尤度比検定p値

尤度比検定統計量は、自由度dfのカイ二乗分布に従います。dfはモデルBとモデルAのパラメータ数の差です。したがって、尤度比検定p値は次のように計算できます。

$\text{P-val} = \text{chi2cdf}( LRT, 1, 1)$ ,

ここでchi2cdfは上側確率（3番目のパラメータが1の場合上側確率を意味する）であり、𝜒²分布に対するCDF関数です。

パラメータの標準誤差の推定

MLE法でパラメータを推定すると、フィッシャー情報行列(FIM)はヘッセ行列で計算できます。

$[\mathcal{I}(\theta)]_{i,j} = - E\left[ \frac{\partial^2}{\partial \theta_i \, \partial \theta_j} \text{ln}\left( f(Y; \theta) \right) \, \Big| \, \theta \right]$

パラメータの共分散行列は次のように表すことができます。

$C = (n \mathcal{I} )^{-1}$

ここで、nはデータの総数です。

パラメータの標準誤差は共分散行列Cの対角成分の平方根です。

確率プロット

確率プロット内の点は昇順に並べ替えられ、i番目の点については、その確率は中央値ランク(Benardの方法)によって計算されます。

$P = \frac{i - 0.3}{n + 0.4}$

期待パーセンタイルの中央線は、最尤推定（MLE）で求めたパラメータを用いた逆累積分布関数（Inverse CDF）により計算されます。

パーセンタイルの信頼限界の推定では、パーセンタイルの分散は誤差伝播の法則を用いて次のように計算されます。

$\sigma^2_{x_p} = \left( \frac{\partial x_p}{\partial \theta} \right)^T C \left( \frac{\partial x_p}{\partial \theta} \right)$

ここで、x_pは与えられた確率pに対するパーセンタイルであり、(.)^Tはベクトルの転置を表します。

パーセンタイルの信頼限界は以下のように計算できます。

$x_{pL} = x_p - Z_{\alpha} \sigma_{x_p}$

$x_{pU} = x_p + Z_{\alpha} \sigma_{x_p}$

ここで $Z_{\alpha} = \Phi^{-1}(1- \alpha / 2)$

対数正規、ガンマ、指数、ワイブル、対数ロジスティック、折りたたみ正規、レイリー分布など、すべての値が正の分布では、信頼限界は以下のようになります。

$x_{pL} = e^{ \text{ln}(x_p) - Z_{\alpha} \frac{\sigma_{x_p}}{x_p} }$

$x_{pU} = e^{ \text{ln}(x_p) + Z_{\alpha} \frac{\sigma_{x_p}}{x_p} }$

閾値パラメータ λ を持つ確率変数の中には、例えば3パラメータ対数正規分布、 2パラメータ指数分布、 3パラメータワイブル分布などでは、信頼限界は次のようになります。

$\lambda < 0$ の場合