連結ガウス分布とは、 正規分布 (ガウス分布) を区間で切断し、隙間を持たせてX軸方向に並べ、累積分布関数のグラフを徐々に変化させて、連続になるように繋げた確率分布です。 あるいは、それを拡張した、複数の正規分布のグラフが横方向や縦方向にグラデーション的に徐々に移り変わっていくような確率分布です。
連結ガウス分布を使えば、どんなクォンタイルを指定しても、それを極めて小さい誤差で再現する分布を構成できます。
中央値が等しく、分散が異なる正規分布を連結させることによって、 どんな経路 (*) を指定しても、 累積分布関数がその経路のすべての点を極めて小さい誤差で通過する確率分布を構成できるからです。
(*) 経路 : ここでは、任意のクォンタイル (分位点) の列の意。
……と書くと、なんとなく凄そうな気がするかも知れませんが、 連結ガウス分布の確率密度関数は、一般に不連続で、いびつで、自然界には絶対に存在しなさそうな分布になります。 連結する隙間の部分に、凄いひずみを作るからです。
たとえば、こんな形の不連続な分布です。
ただし、クォンタイルが8点以下の場合は (たいていの場合、クォンタイルは3点か5点、あるいは7点でしょう)、 確率密度関数が連続な、自然界に存在しうるような分布モデルを構成することもできます。 設定 (type1.type, continuous / symmetric / v.grad) によって、
- 2つの正規分布の平均 (type1.type = 1, continuous)、
- 横方向 (X軸方向) にグラデーション的に正規分布の混合比率が変化する分布 (type1.type = 2, continuous)、
- 確率密度関数の中央が鋭利に尖ったり、逆に凹んだりしている、左右対称な分布 (type1.type = 2, symmetric)、
- 縦方向 (Y軸方向) にグラデーション的に正規分布の混合比率が変化する分布 (type1.type = 3)
- 縦横両方向にグラデーション的に正規分布の混合比率が変化する分布 (type1.type = 4)
のような連続分布の分布モデルが作れます。たとえば、こんな形の分布が作れます。
これらの分布モデルは、正規分布に従わない、歪んだ分布や、裾野の広い分布のモデルなどに使えると思われます。
ともかく、このパッケージはそういった確率分布を構成するためのリファレンスクラスのパッケージです。
連結ガウス分布は混合ガウス分布とは少し違います。 混合ガウス分布は、複雑な分布のデータを複数の正規分布の混合と考え、1つのデータを複数の正規分布に従うデータに分離することを目的とした、クラスタリングの手法と考えられます。 それに対して、連結ガウス分布は、複雑な分布のデータをあるがままの一団のデータとして捉えて、解析したり、他のデータと比較したりすることを目的と考えています。
連結ガウス分布は以下のような目的に使うために、このパッケージの作者が考案しました。
-
不連続な連結ガウス分布
正規分布に従わない、分布モデルが不明なデータのクォンタイルが与えられているときに、リサンプルしたランダムデータが欲しい場合
-
連続な連結ガウス分布
無限個の正規分布の混合モデルを考える場合
このパッケージの作者は、とある医療・人体関係のデータを扱っていて、上記のようなデータに遭遇したので、このようなモデルを考案しました。
# Install devtools from CRAN install.packages( "devtools" ) # Then use devtools::install_github( "user/repository" ) to install cgd package from GitHub devtools::install_github( "Kimitsuna-Goblin/cgd" )
> library( cgd ) # ライブラリを読み込みます
> a <- CGD$new() # 連結ガウス分布クラスのオブジェクトを生成します
>
> ########################################################################
> # 経路 (クォンタイルの列) をトレースする方法
> ########################################################################
> # CGD$set.waypoints() : 経路 (クォンタイルの列) を指定し、
> # 指定されたすべての点を通過する累積分布関数を持つ連結ガウス分布を構成します。
> # p は確率、q はクォンタイルのX座標値です。
> a$set.waypoints(
+ data.frame(
+ p = c( 0.1, 0.3, 0.5, 0.6, 0.7 ),
+ q = c( qnorm( c( 0.1, 0.3, 0.5, 0.6 ), 0, 1 ), 0.5 ) ) )
> NULL
>
> # trace.q() : クォンタイルの列をトレースする連結ガウス分布クラスオブジェクトを生成します
> # CGD$set.waypoints() と同じです (new しなくて良い分、こちらの方が簡単です)
>
> a <- trace.q(
+ data.frame(
+ p = c( 0.1, 0.3, 0.5, 0.6, 0.7 ),
+ q = c( qnorm( c( 0.1, 0.3, 0.5, 0.6 ), 0, 1 ), 0.5 ) ) )
>
> # CGD$d() : X座標を指定して、確率密度を返します
> dev.new(); plot.new()
> plot( seq( -3, 3, 0.01 ), a$d( seq( -3, 3, 0.01 ) ), type = "l" )
>
> # CGD$p() : X座標を指定して、確率を返します
> dev.new(); plot.new()
> plot( seq( -3, 3, 0.01 ), a$p( seq( -3, 3, 0.01 ) ), type = "l" )
>
> # CGD$q() : 確率を指定して、クォンタイル (X座標) を返します
> # 確率が同一となるX座標が、ある区間内に無限に存在し、
> # 一意に定まらない場合は、該当区間の中点の座標を返します
> dev.new(); plot.new()
> plot( seq( 0, 1, 0.01 ), a$q( seq( 0, 1, 0.01 ) ), type = "l" )
>
> # CGD$r() : ランダムサンプルを生成します (高速化は全然やってません)
> dev.new(); plot.new()
> sample <- a$r( 1000 )
> hist( sample )
>
> a # オブジェクトの内部構造を表示します
Reference class object of class "CGD"
Field "kind.index":
[1] 16
Field "kind":
[1] "Discontinuous Connected Gaussian Distribution"
Field "median":
[1] 0
Field "intervals":
[[1]]
Reference class object of class "CGDInterval"
Field "mean":
[1] 0
Field "sd":
[1] 1
Field "q.ind":
[1] -Inf 0.2533471
Field "q.conn.prev":
[1] -Inf -Inf
Field "q.conn.next":
[1] 0.2533471 0.5244005
Field "p.ind":
[1] 0.0 0.6
Field "p.conn.prev":
[1] 0 0
Field "p.conn.next":
[1] 0.6 0.7
[[2]]
Reference class object of class "CGDInterval"
Field "mean":
[1] 0
Field "sd":
[1] 0.9534697
Field "q.ind":
[1] 0.5 Inf
Field "q.conn.prev":
[1] 0.2415588 0.5000000
Field "q.conn.next":
[1] Inf Inf
Field "p.ind":
[1] 0.7 1.0
Field "p.conn.prev":
[1] 0.6 0.7
Field "p.conn.next":
[1] 1 1
Field "type1.type":
[1] 1
Field "mean":
[1] -0.01719301
Field "sd":
[1] 0.9772327
Field "lsd":
[1] 0.9863357
Field "usd":
[1] 0.968044
Field "lsd.abs.error":
[1] 0
Field "usd.abs.error":
[1] 0
> # 各要素の意味はマニュアルまたはソースファイルのコメントを参照してください
>
> ########################################################################
> # 度数分布を近似する方法
> ########################################################################
> # 度数分布のデータを用意します
> x <- seq( -2, 2, 0.2 )
> freq <- c( 15164, 22923, 25134, 27631, 37239, 40464,
+ 47126, 79469, 109966, 118241, 111333, 78674,
+ 46921, 41026, 36975, 27403, 25493, 22838,
+ 14992, 11468, 9174 )
> total <- sum( freq )
>
> # CGD$nls.freq() : 与えられた度数分布に近くなるように連続分布を構成します
> # 引数に、構成したい分布の種類を指定してください
> a$nls.freq( x, freq, total, kind = "Median-Equaled Sigma-Differed Vertical Gradational Distribution" )
Nonlinear regression model
model: d ~ dp.t3(x, c(mean, mean, mean), c(sqrt.sd.1^2, sqrt.sd.2^2, sqrt.sd.1^2), f.t3.d)
data: list(d = get.d(x, freq, total), x = x)
mean sqrt.sd.1 sqrt.sd.2
-0.1992 1.0198 0.8050
residual sum-of-squares: 0.02285
Number of iterations to convergence: 30
Achieved convergence tolerance: 9.922e-06
>
> # a <- nls.freq() のような書き方もできます (new しなくて良い分、こちらの方が簡単です)
> a <- nls.freq( x, freq, total, kind = "Median-Equaled Sigma-Differed Vertical Gradational Distribution" )
Nonlinear regression model
model: d ~ dp.t3(x, c(mean, mean, mean), c(sqrt.sd.1^2, sqrt.sd.2^2, sqrt.sd.1^2), f.t3.d)
data: list(d = get.d(x, freq, total), x = x)
mean sqrt.sd.1 sqrt.sd.2
-0.1992 1.0198 0.8050
residual sum-of-squares: 0.02285
Number of iterations to convergence: 30
Achieved convergence tolerance: 9.922e-06
>
> # CGD$d() : X座標を指定して、確率密度を返します
> dev.new(); plot.new()
> plot( seq( -3, 3, 0.01 ), a$d( seq( -3, 3, 0.01 ) ), type = "l" )
>
> # CGD$tex() : 得られた分布の累積分布関数と確率密度関数の表現式を TeX 形式で表示します
> a$tex()
\begin{align}
\Psi(x) &= \Phi_1(x) - \dfrac{1}{\sqrt{2}} \Phi^\ast_1(x) + \dfrac{1}{\sqrt{2}} \Phi^\ast_2(x),\\
g(x) &= \left( 1 - \dfrac{f_1(x)}{f_1(\mu_1)} \right) f_1(x) + \dfrac{f_2(x)}{f_2(\mu_2)} f_2(x),\\
\\
\Phi_i(x) &= \dfrac{1}{\sqrt{2 \pi \sigma_i^2}} \int_{-\infty}^{x}
\exp \left( -\dfrac{(t - \mu_i)^2}{2 \sigma_i^2} \right) dt,\\
\Phi^\ast_i(x) &= \dfrac{1}{\sqrt{2 \pi \left( \begin{array}{c}
\dfrac{\sigma_i}{\sqrt{2}} \end{array} \right)^2}}
\int_{-\infty}^{x} \exp \left( \begin{array}{c} -\dfrac{(t - \mu_i)^2}{2
\left( \begin{array}{c} \dfrac{\sigma_i}{\sqrt{2}} \end{array} \right)^2} \end{array} \right) dt,\\
f_i(x) &= \dfrac{1}{\sqrt{2 \pi \sigma_i^2}} \exp \left( -\dfrac{(x - \mu_i)^2}{2 \sigma_i^2} \right),\\
\\
& \begin{array}{l}
\mu_1 = -0.199217, & \sigma_1 = 1.039974,\\
\mu_2 = \mu_1, & \sigma_2 = 0.648027.
\end{array}
\end{align}
>
> # nls.freq.all() : パッケージでサポートされているすべての連続分布の中から、
> # 与えられた度数分布に最も近い分布を探します (少し時間がかかります)
> result <- nls.freq.all( x, freq, total )
>
> # cgd:::kinds : サポートされている分布の種類です ( 1 ～ 15 は連続分布、 16 は不連続分布です)
> cgd:::kinds
[1] "Normal Distribution"
[2] "Mean of Median-Equaled Sigma-Differed 2 Normal Distributions"
[3] "Symmetric Horizontal Gradational Distribution"
[4] "Median-Differed Sigma-Equaled Horizontal Gradational Distribution"
[5] "Median-Equaled Sigma-Differed Horizontal Gradational Distribution"
[6] "Median-Differed Sigma-Differed Horizontal Gradational Distribution"
[7] "Median-Differed Sigma-Equaled Vertical Gradational Distribution"
[8] "Median-Equaled Sigma-Differed Vertical Gradational Distribution"
[9] "Median-Differed Sigma-Differed Vertical Gradational Distribution"
[10] "3-Median-Differed Sigma-Equaled Vertical Gradational Distribution"
[11] "Median-Equaled 3-Sigma-Differed Vertical Gradational Distribution"
[12] "3-Median-Differed 3-Sigma-Differed Vertical Gradational Distribution"
[13] "Median-Differed Sigma-Equaled Vertical-Horizontal Gradational Distribution"
[14] "Median-Equaled Sigma-Differed Vertical-Horizontal Gradational Distribution"
[15] "Median-Differed Sigma-Differed Vertical-Horizontal Gradational Distribution"
[16] "Discontinuous Connected Gaussian Distribution"
>
> # nls.freq.all()$cor : 各モデルと度数分布との相対係数です
> result$cor
[1] 0.9286582 0.9878450 0.9277467 0.8812164 0.9551502 0.8499192 0.9551271
[8] 0.9867571 0.9867073 0.9945393 0.9867147 0.9538990 0.8879331 0.9732565
[15] 0.9964684
>
> # nls.freq.all()$best.cor : 各モデルと度数分布との相対係数のうち、最大の値です
> result$best.cor
[1] 0.9964684
>
> # nls.freq.all()$best : 最も度数分布に近い結果が得られたモデルです
> # CGD$kind : その分布モデルの種類を表す文字列です
> result$best$kind
[1] "Median-Differed Sigma-Differed Vertical-Horizontal Gradational Distribution"
>
> # CGD$d() : X座標を指定して、確率密度を返します
> dev.new(); plot.new()
> plot( seq( -3, 3, 0.01 ), result$best$d( seq( -3, 3, 0.01 ) ), type = "l" )
一般に、 (不連続な) 連結ガウス分布は次のようにして構成することができる。
- まず、確率が取りうる値の全体である
$[0, 1]$ 区間の部分閉区間を有限個 ($n$ 個) 持つ集合$P$ を、閉区間同士が互いに重なり合わないようにして用意する。 ただし、閉区間のどれか1つに必ず 0 が含まれ、どれか1つに必ず 1 が含まれるようにする。
つまり、
を満たすような、閉区間の集合
各要素
それぞれの閉区間
ここで、予め、有限の
この際、もし、
- 閉区間の集合
$P$ の要素数と等しい、$n$ 個の正規分布からなる集合$N = \left\lbrace N_1, N_2, \cdots, N_n \right\rbrace$ を用意して、 その要素の正規分布$N_i$ $( 1 \leq i \leq n )$ に、 閉区間$P_i$ をそれぞれ1つずつ割り当てる。
このとき、集合
- 任意の2つの
$N_i, N_j ( i \neq j )$ について、$N_i, N_j$ に閉区間$P_i, P_j$ が割り当てられているとき、$N_i, N_j$ の累積分布関数$\Phi_i( x ), \Phi_j( x )$ について、 それぞれの区間に対する定義域$\lbrace x \mid \Phi_i( x ) \in P_i \rbrace$ と$\lbrace x \mid \Phi_j( x ) \in P_j \rbrace$ が、互いに重ならないこと。 逆関数を使って言い換えれば、$\left\lbrace \Phi_i^{-1}( p ) \mid p \in P_i \right\rbrace ∩ \left\lbrace \Phi_j^{-1}( p ) \mid p \in P_j \right\rbrace = \varnothing$ であること。 - 経路
$W = \left\lbrace ( x_j, p_j ) \right\rbrace$ が与えられており、 経路の点の確率$p_j$ が閉区間$P_i$ に含まれているときは、 その区間が割り当てられている確率分布$N_i$ の累積分布関数$\Phi_i$ が、必ずその経路の点$( x_j, p_j )$ を通過すること。 すなわち、$p_j \in P_i \Rightarrow \Phi_i( x_j ) = p_j$ が成り立つこと。
なお、集合
ある1つの閉区間
さて、これまでの工程によって、 与えられた経路の点をすべて通過する、累積分布関数の不連続な断片が構成できたので、 あとは、それらの断片を連続になるように接続すればよい。 そこで、最終的な工程として、次のように断片の接続を行う。
以下、記述を簡単にするため、
閉区間
- 開区間
$Q_i = (b_i, a_{i+1})$ $( 1 \leq i \leq n - 1 )$ に対して、$\Psi_i( \beta_i ) = b_i, \Psi_i( \alpha_{i+1} ) = a_{i+1}$ を満たすような、 区間$[\beta_i, \alpha_{i+1}]$ において単調増加する連続関数$\Psi_i(x)$ を適当に定める。 注意点として、$x$ が$\beta_i < x < \alpha_{i+1}$ の範囲を動くとき、$\Psi_i(x)$ の値は必ず$Q_i$ の範囲の中に収まらなければならず、決して$Q_i$ の範囲外に出てはならない。
以上のように関数
を累積分布関数とする確率分布が定義できる。 以上のように構成した累積分布関数を持つ確率分布を 「連結ガウス分布」 と呼ぶ。
これまで述べてきた閉区間
以上、任意の経路の点を通過する、不連続な連結ガウス分布の構成方法を述べたが、
独立区間を
経路をトレースするような (つまり、累積分布関数が経路の点をすべて通過するような) 不連続な連結ガウス分布を構成する場合、 本パッケージでは、以下のように独立区間を構成する。
-
原則として、経路の点
$(x_i, a_i)$ の確率$a_i$ に対し、1点のみからなる区間$[a_i, a_i]$ を独立区間の確率とする。 ただし、以下の2つの場合を例外とする。-
経路の点のうち、確率が 0 または 1 でない最初の点
$( x_1, b_1 )$ および 最後の点$( x_n, a_n )$ $( 0 < b_1, a_n < 1 )$ の確率$b_1, a_n$ に対しては、 確率 0 および 1 を含む区間$[0, b_1], [a_n, 1]$ を独立区間の確率とする。 -
2つ以上の連続する経路の点が、 1つの同じ正規分布 (ただし、中央値は経路の点
$( x, 0.5 )$ の$x$ の値とする) の累積分布関数上の点となる場合は、 それらの点の確率をまとめて1つの独立区間の確率に含める。 ただし、独立区間の確率の上限と下限は、必ず経路に含まれる点の確率か、 または 0 または 1 のいずれかとする。
-
独立区間の定義域を負担する正規分布は、以下のように構成する。
-
正規分布の中央値は、経路で与えられた点
$( x, 0.5 )$ の$x$ の値とする。 -
最初の独立区間の確率
$[0, b_1]$ の定義域$[-\infty, \beta_1]$ に対しては、 累積分布関数が経路の点$( \beta_1, b_1 )$ を通るような正規分布を採用し、 それ以外の独立区間の確率$[a_i, b_i]$ の定義域$[\alpha_i, \beta_i]$ に対しては、 累積分布関数が経路の点$( \alpha_i, a_i )$ を通るような正規分布を採用する (その累積分布関数は必ず点$( \beta_i, b_i )$ も通る)。
ここで、ある正規分布
点
で表される。この式の右辺を
となる。ただし、
という式が得られる。ここに、標準偏差
この式から、正規分布の平均値 (すなわち中央値) と、 それ以外の任意の1点のX座標とその点における確率を与えれば、 それらの条件を満たす正規分布の標準偏差が得られる。 したがって、上に述べたような正規分布は容易に見つけることができる。
本パッケージでは、この式を使って正規分布を構成し、
連結ガウス分布の累積分布関数
不連続な連結ガウス分布の接続関数は、 以下のように、type 1, type 2, type 3a, type 3b の4つの場合に分けて構成する (オプション type1.type = 1 または 2 の場合。オプションの詳細については後述)。
以下、独立区間の定義域
また、接続区間の確率
- 接続区間の定義域の上限
$\alpha_{i+1}$ が$m$ より小さく、標準偏差が$\sigma_i < \sigma_{i + 1}$ の場合 - 接続区間の定義域の下限
$\beta_i$ が$m$ より大きく、標準偏差が$\sigma_i > \sigma_{i + 1}$ の場合
分布の山側の標準偏差が、裾側よりも大きい場合である。
この場合は、すべての
本パッケージでは、 type1.type というオプションによって、以下の 表1 のように接続関数を構成する。 このオプションの名前は「 type 1 の接続区間における計算方法が異なる」という意味で付けられた。
表中の式で、
| オプション | 接続関数 |
通過できる経路の点の個数 | 確率密度関数の連続性 |
|---|---|---|---|
| type1.type = 0 | 接続関数なし (独立区間の両端の2点を通る正規分布を横に並べただけの分布) |
任意 (2個以上、上限なし) |
不連続 |
| type1.type = 1 |
|
任意 (1個以上、上限なし) |
不連続 |
| type1.type = 2 |
|
任意 (1個以上、上限なし) |
不連続 |
- 接続区間の定義域の上限
$\alpha_{i+1}$ が$m$ より小さく、標準偏差が$\sigma_i \geq \sigma_{i + 1}$ の場合 - 接続区間の定義域の下限
$\beta_i$ が$m$ より大きく、標準偏差が$\sigma_i \leq \sigma_{i + 1}$ の場合
分布の山側の標準偏差が、裾側よりも小さい場合である。
この場合は、一部の
本パッケージでは、以下のように
このとき、確率密度関数
となる。このとき、
-
$\beta_i \leq \mu \leq \alpha_{i+1}$ であり、標準偏差が$\sigma_i < \sigma_{i + 1}$ の場合
この場合は
このとき、確率密度関数
となる。
type1.type = 1, continuous = TRUE のときは、この type 3 の計算を拡張して、2つの正規分布の平均とする。
-
$\beta_i \leq \mu \leq \alpha_{i+1}$ であり、標準偏差が$\sigma_i > \sigma_{i + 1}$ の場合
この場合は type 3a と対称的である。
すなわち、
このとき、確率密度関数
となる。
連結ガウス分布の独立区間の確率を
したがって、
不連続分布の type 1 の接続関数
以下、連続な連結ガウス分布の累積分布関数を添字のない
連続な連結ガウス分布の場合、与えられた経路の点を通過させるには、
経路の点
を満たせばよい。 この連立方程式は nleqslv パッケージを使って解くことができる。
ただし、連立方程式の式の数は、変数の個数を超えられないので、通過できる経路の点の個数には限界がある。 ここでの変数は、連結ガウス分布を構成する正規分布の平均値と標準偏差である。 つまり、正規分布の平均値と標準偏差の個数を超える点を通過できる分布は、一般に構成できない。 また、経路の点が多くなくても、経路の点の並び方によっては、連立方程式の実解が得られず、分布が構成できないことがある。
本パッケージでは、 type1.type オプションによって、以下の 表2 のような分布を構成できる。 オプションによって、通過できる経路の点の条件が異なるので、条件に応じたオプションを選ぶとよい。
表中の式で、
| オプション | 累積分布関数 |
通過できる経路の点の個数 | 独立区間の確率 | 確率密度関数の連続性 |
|---|---|---|---|---|
| type1.type = 1, continuous = TRUE (or symmetric = TRUE) (2つの正規分布の平均) |
where |
3点 |
|
連続 ( |
| type1.type = 2, continuous = TRUE (横方向グラデーション) |
|
3～4点 |
|
連続 ( |
| type1.type = 2, symmetric = TRUE (左右対称) |
|
|
|
連続 |
| type1.type = 3 (下側の分布が左右で異なる、 歪んだ縦方向グラデーション) |
|
3～6点 |
あるいは (ただし、条件によっては、 |
連続 ( |
| type1.type = 3, v.grad = TRUE (2つの正規分布の縦方向グラデーション) |
|
3～4点 |
あるいは |
連続 ( |
| type1.type = 4 (縦横グラデーション) |
|
5～8点 |
|
連続 ( |
なお、 type1.type = 1, continuous = TRUE (or symmetric = TRUE) (2つの正規分布の平均) は、正確には不連続分布の type 1 の接続関数の拡張ではなく、 type 2/3a/3b の接続関数の拡張であるが、便宜上、このようなオプションにて実装した。
v1.2.0 にあった type1.type = 3, symmetric = TRUE (縦方向グラデーションの旧バージョン) のオプションは廃止され、 v1.3.8 以降、 type1.type = 3, v.grad = TRUE に上位互換的に変更された。
type1.type = 4 (縦横グラデーション) は正規分布の連結ではなく、2つの連結ガウス分布 (2つの正規分布の縦方向グラデーション) を横方向グラデーションで連結した分布である。
著作者が思いつきで考えて作った確率分布なので、特に外部資料はありません。
ある程度統計をやってたら誰でも思いつきそうな確率分布だし、どっかにあるんじゃね？という気がするので、見つけたら教えてください。
Kimitsuna-Goblin (浦 公統; Ura Kimitsuna)