グンベル極値分布

グンベル極値分布 (Gumbel extreme value distribution) は連続確率分布のひとつである．分布名のグンベルとはドイツの数学者 Emil Julius Gumbel に由来する．また，日本語ではガンベル極値分布とも発音されることがある．グンベル極値分布は，タイプ1，タイプ2 (フレシェ型) およびタイプ3 (ワイブル型) からなる一般極値分布 (Generalized extreme value distribution) のタイプ1に相当するものである．単にグンベル分布と呼ばれることもある．グンベル極値分布は，互いに独立で同一の何らかの分布に従う標本における最大値を集めてそれらを新たな確率変数Xとみなしたとき，その確率変数Xの分布が漸近的に従う確率分布である．ある期間においてある大きな値が出現する確率を予測するのに利用される．例えば，ある期間における大きなマグニチュードを持つ地震が起こる確率を予測することや大雨による洪水が発生する確率の予測をすること等がそれにあたる．また，本分布はバイオインフォマティクス研究にとっても重要な確率分布のひとつである．BLASTをはじめとする配列検索法におけるアラインメントスコアの分布はグンベル極値分布にてモデル化される．パラメーターは位置パラメーター μ (>0) および尺度パラメーター θ (>0) であり，グンベル極値分布は Gb(μ, θ) にて略記される．確率密度関数は以下で与えられる．

\begin{eqnarray*}f(x)=\frac{1}{\theta}e^{-\frac{x-\mu}{\theta}}e^{-e^{-\frac{x-\mu}{\theta}}}\tag{1}\end{eqnarray*}

累積分布関数は以下で与えられる．位置パラメータの μ は分布の最頻値 (モード) を指す．

\begin{eqnarray*}F(x)=e^{-e^{-\frac{x-\mu}{\theta}}}\tag{2}\end{eqnarray*}

確率変数Xの範囲は以下である．

\begin{eqnarray*}-\infty < x < \infty\tag{3}\end{eqnarray*}

モーメント母関数は $t < \frac{1}{\theta}$ の範囲で定義され以下の式で与えられる．式中の $\Gamma$ はガンマ関数を指す．

\begin{eqnarray*}M_X(t)=e^{\mu t}\Gamma(1-\theta t)\ \ \ \ \left(t < \frac{1}{\theta}\right)\tag{4}\end{eqnarray*}

期待値は以下で与えられる値である．ここで，γ はオイラーの定数を指す．

\begin{eqnarray*}E(X)=\mu+\gamma\theta\tag{5}\end{eqnarray*}

分散は以下の式で与えられる．

\begin{eqnarray*}V(X)=\frac{\theta^2\pi^2}{6}\tag{6}\end{eqnarray*}

期待値の導出

期待値は定義式に従って以下のように求める．

\begin{eqnarray*}E(X)&=&\int_{-\infty}^{\infty}xf(x)dx\\&=&\int_{-\infty}^{\infty}x\frac{1}{\theta}e^{-\frac{x-\mu}{\theta}}e^{-e^{-\frac{x-\mu}{\theta}}}dx\\\tag{7}\end{eqnarray*}

ここで以下の変換を考える．

\begin{eqnarray*}z=e^{-\frac{x-\mu}{\theta}}\tag{8}\end{eqnarray*}

これは同時に以下の関係も表す．

\begin{eqnarray*}x=\mu-\theta\ln z\tag{9}\end{eqnarray*}

\begin{eqnarray*}dx=-\frac{\theta}{z}dz\tag{10}\end{eqnarray*}

以上の関係を用いて上式を変換すると以下のようになる．また，積分範囲も変わる．

\begin{eqnarray*}(上式)&=&-\int_{0}^{\infty}(\mu-\theta\ln z)\frac{1}{\theta}ze^{-z}\left(-\frac{\theta}{z}\right)dz\\&=&\int_{0}^{\infty}(\mu-\theta\ln z)e^{-z}dz\\&=&\mu\int_{0}^{\infty}e^{-z}dz-\theta\int_{0}^{\infty}e^{-z}\ln zdz\\&=&\mu+\theta\left(-\int_{0}^{\infty}e^{-z}\ln zdz\right)\tag{11}\end{eqnarray*}

ここで，上式の括弧内はオイラーの定数を表すものである．

\begin{eqnarray*}\gamma=-\int_{0}^{\infty}e^{-z}\ln zdz\tag{12}\end{eqnarray*}

よって，以上の関係を用いることで期待値は以下のように表される．

\begin{eqnarray*}E(X)=\mu+\gamma\theta\tag{13}\end{eqnarray*}

一般極値分布との関係

以下の確率密度関数を持つ確率分布を一般極値分布という．パラメーターは位置パラメーター μ，尺度パラメーター θ および形状パラメーター γ の3つであり，GEV(μ, σ, γ) によって略記される．

\begin{eqnarray*}f(x)=\frac{1}{\theta}\left[1+\gamma\left(\frac{x-\mu}{\theta}\right)\right]^{-\left(\frac{1}{\gamma}+1\right)}e^{-\displaystyle \left[1+\gamma\left(\frac{x-\mu}{\theta}\right)\right]^{-\frac{1}{\gamma}}}\tag{14}\end{eqnarray*}

グンベル極値分布は一般極値分布において，γ → 0 であるときに相当する．すなわち，確率変数Xが以下のようにグンベル極値分布に従うとき，

\begin{eqnarray*}X\sim Gb(\mu,\theta)\tag{15}\end{eqnarray*}

この確率変数Xは同時に以下の一般極値分布に従う．

\begin{eqnarray*}X\sim GEV(\mu,\theta,0)\tag{16}\end{eqnarray*}

もし，γ > 0 の場合はタイプ2のフレシェ型，γ < 0 の場合はタイプ3のワイブル型となる．

ワイブル分布との関係

ワイブル分布 W(a, b) の確率密度関数 g(y) は以下で与えられる．

\begin{eqnarray*}g(y)=\frac{by^{b-1}}{a^b}e^{-\left(\displaystyle \frac{y}{a}\right)^b}\tag{17}\end{eqnarray*}

ここで以下のような変数変換を考える．

\begin{eqnarray*}x=-b\ln\frac{y}{a}\tag{18}\end{eqnarray*}

これは同時に以下の式も意味する．

\begin{eqnarray*}y=ae^{-\displaystyle \left(\frac{x}{b}\right)}\tag{19}\end{eqnarray*}

以上の関係を用いて確率変数の変換を行うと変換後の確率密度関数 f(x) は以下のようになる．

\begin{eqnarray*}f(x)=g(y)\frac{dy}{dx}=e^{-x}e^{-e^{-x}}\tag{20}\end{eqnarray*}

このようにして得られた確率密度関数 f(x) はまさに標準グンベル極値分布 Gb(0, 1) に他ならない．このような背景から，グンベル極値分布はしばしば対数ワイブル分布と呼ばれることがある．

配列類似性検索法における p-value の導出

バイオインフォマティクス分野におけるBLASTをはじめとする配列類似性検索法において，得られたアラインメントスコアの統計的有意性はグンベル極値分布を利用して求められる．すなわち，BLAST (等の配列類似性検索法) では検索で得られたアラインメントスコアをグンベル極値分布に従う確率変数とみなしている．そもそも，グンベル極値分布は何らかの分布に従う標本の最大値を集めてそれらを確率変数としたときに，その確率変数が従う確率分布であった．BLASTにおいても同様に考えることができる．BLASTではクエリ配列に対してデータベース中のサブジェクト配列と可能な限りのアラインメントを試行し，その結果得られた最大スコアに対してその統計的有意性を議論するものであるので，それらのスコアを集めて確率変数としたときに，その確率変数が従う確率分布はグンベル極値分布となるのは当然である．得られたアラインメントスコアの統計的有意性の論理の流れは以下のようになる．

検索によってアラインメントスコアSが得られた．このSはグンベル極値分布に従う．
グンベル極値分布の累積分布関数を利用することで，このS以上のスコアが得られる確率が計算できる．
S以上のスコアが得られる確率が極めて低いなら，Sが得られた事実は統計的に有意であると判断する．

以上のS以上のスコアが (偶然に) 得られる確率のことを p-value (p値) といい P(S) で表す．P(S) は以下のように求める．ここで，Altschul の論文や NCBI のウェブサイトにおいてはグンベル極値分布を表す場合，位置パラメーター μ を u に，尺度パラメーター θ を 1/λ にて置換して表現する．これに従うと確率密度関数 f(x) および累積分布関数 F(x) は以下のように表される．

\begin{eqnarray*}f(x)=\lambda e^{-\lambda(x-u)}e^{-e^{-\lambda(x-u)}}\tag{21}\end{eqnarray*}

\begin{eqnarray*}F(x)=e^{-e^{-\lambda(x-u)}}\tag{22}\end{eqnarray*}

以上の表記を用いて，まず，得られたアラインメントスコアが S であるときその S より低いスコアが得られる確率 P₀(S) は以下のように求められる．これは単に，負の無限大から S の範囲における確率密度関数 f(x) の積分であり，それは累積分布関数 F(x) と定義されるものである．

\begin{eqnarray*}P_0(S)&=&\int_{-\infty}^{S}f(t)dt\\&=&F(S)\\&=&e^{-e^{-\lambda(S-u)}}\tag{23}\end{eqnarray*}

一方で，S 以上のスコアが偶然に得られる確率 P(S) は上の関係を用いて以下のように計算される．

\begin{eqnarray*}P(S)&=&1-P_0(S)\\&=&1-e^{-e^{-\lambda(S-u)}}\tag{24}\end{eqnarray*}

ここで，λ はアミノ酸スコアテーブル (アミノ酸とアミノ酸を並べるときに用いる20行*20列の400要素からなる得点表) から求められる値である．これは，s_ij をそのスコア，p_i および p_j を20種類のアミノ酸のそれぞれの存在比 (データベース中の存在比でも良いし，アミノ酸スコアテーブル自体から逆算することもできる) としたとき，以下の式を満たす値となる．

\begin{eqnarray*}\sum_{i=1}^{20}\sum_{j=i}^{20}p_ip_je^{\lambda s_{ij}}=1\tag{25}\end{eqnarray*}

また，u は以下の式で与えられる値である．ここで，K はアミノ酸の頻度 p_i およびアミノ酸スコアテーブル s_ij から計算される定数であり，m はクエリの残基長，n はデータベースの残基長 (データベース中の全配列をひとつに繋げたときの長さ) である．

\begin{eqnarray*}u=\frac{1}{\lambda}\ln Kmn\tag{26}\end{eqnarray*}

以上の関係を用いることで上の P(S) は以下のように変形される．

\begin{eqnarray*}P(S)&=&1-e^{-e^{-\lambda(S-u)}}\\&=&1-e^{-e^{-\lambda(S-\frac{1}{\lambda}\ln Kmn)}}\\&=&1-e^{-e^{-(\lambda S+\ln Kmn)}}\\&=&1-e^{-Kmne^{-\lambda S}} \tag{27}\end{eqnarray*}

以上が求める p-value である．実際には，この p値が大体 10^-4 より小さい場合，検出されたサブジェクト配列はクエリ配列と相同であるとみなす場合が多い．医学・薬学・生物学分野で一般に用いられる p-value よりも随分厳しい．

[式の簡略化1] ここまでで，p値の導出は終了となるが，計算の節約ために，上の p値は近似を用いることでさらに簡略化することができる．グンベル極値分布の確率密度関数は，X軸の値が大きいとき，指数関数的に減少する．よって，S>2 の範囲では P(S) は以下のように近似される．

\begin{eqnarray*}P(S)\approx Kmne^{-\lambda S}\tag{28}\end{eqnarray*}

上の値は p-value の近似式であるが，これは E-value (クエリ配列をデータベースに対して検索したときに，得られたスコアが S 以上の値になるアラインメントの本数の期待値) の値と等しい．得られたスコア S が大きい範囲では p-value と E-value は一致する．

\begin{eqnarray*}P(S)\approx E(S)\ \ \ \ (S>2)\tag{29}\end{eqnarray*}

ただし，本来の p-value と E-value の関係はそれぞれを P(S) および E(S) で表したとき，以下のように与えられる．

\begin{eqnarray*}P(S)=1-e^{-E(S)}\tag{30}\end{eqnarray*}

[式の簡略化2] この他にも，スコア S を正規化することで式の簡略化が可能である．正規化後のスコアを S_N とすると，変換式は以下のように表される．λ および K は事前に求まる値なのでこのような変換が可能となる．ウェブ上では自然対数の真数を Kmn としている場合も多く存在するが，NCBI の公式サイトには以下の変換式が記されている．

\begin{eqnarray*}S_N=\lambda S-\ln K\tag{31}\end{eqnarray*}

この関係式と，式27の3行目の関係を用いることで P(S_N) は以下のように表される．

\begin{eqnarray*}P(S_N)=1-e^{-mne^{-S_N}}\tag{32}\end{eqnarray*}

このときの E-value は以下のようになる．

\begin{eqnarray*}E(S_N)=mne^{-S_N}\tag{33}\end{eqnarray*}

[式の簡略化3] また，以上とほぼ同じ変換ではあるが，以下のような変換式を用いてスコア S を正規化することで式の簡略化が可能となる．この正規化後のスコア S' は bit score とよばれる値である (一方でスコア S は raw score とよばれる)．

\begin{eqnarray*}S'=\frac{\lambda S-\ln K}{\ln 2}\tag{34}\end{eqnarray*}

この関係式と上の式27の3行目の関係を用いることで P(S') は以下のように表される．

\begin{eqnarray*}P(S')=1-e^{-mn2^{-S'}}\tag{35}\end{eqnarray*}

このときの E-value は以下のようになる．

\begin{eqnarray*}E(S')=mn2^{-S'}\tag{36}\end{eqnarray*}

正規分布

グンベル極値分布