尺度混合表現

こんにちは，株式会社Nospareの小林です．
前回の記事前回の記事では，正規分布の尺度混合表現について紹介しました．尺度混合表現は，$t$分布などに従う確率変数$Y$を，正の独立な確率変数$V$を条件づけたときに

$$
Y = \mu+\sigma\sqrt{V} Z,\quad Z\sim N(0,1)
$$

となるように表現するものでした．$Y$の周辺確率密度関数$p(y)$は，$V$について積分して

$$
p(y)=\int_0^\infty \phi(y;0,v\sigma^2)p(v)dv,\quad (2)
$$

のように得られます（$\phi(x;\mu,\sigma^2)$は$N(\mu,\sigma^2)$の密度関数，$p(v)$は$V$の密度関数）．

尺度混合表現を利用することによって$Y$の条件付き分布が正規分布になるので，数値計算が効率的・アルゴリズムの一部のステップの適用が簡単になることがあります．一方で尺度混合表現は左右対称な分布しか表現できません．

位置尺度混合表現

前回も少し触れましたように，位置尺度混合表現は例えば

$$
Y = \mu+\beta W + \sigma\sqrt{V} Z,\quad Z\sim N(0,1)
$$

のような形で与えられます．$W$と$V$を所与としたときの$Y$の条件付き分布は$N(\mu+\beta W, V \sigma^2)$です．尺度混合表現と同様に$V$と$W$の分布によって様々な$Y$の周辺分布を表現することができます．$Y$の周辺密度は

$$
p(y)=\int_0^\infty\int_0^\infty\phi(y;\beta w, v\sigma^2)p(v)p(w)dvdw,
$$

となります．$\beta=0$のときは尺度混合表現となり，$Y$の分布は対称になりますが，$\beta\neq0$のときは非対称になり，$\beta$がゼロから離れるにつれて歪度が大きくなります．

便利な位置尺度混合表現

非対称正規分布

正規分布は左右対称な分布ですが，これを非対称化した分布は数多く提案されています．その中でも，Azzalini (1986)の非対称正規分布$SN(\mu,\sigma^2,\alpha)$の密度関数は

$$
f(y)=\frac{2}{\sigma}\phi\left(\frac{y-\mu}{\sigma}\right)\Phi\left(\alpha\left(\frac{y-\mu}{\sigma}\right)\right),
$$

ここで$\Phi(\cdot)$は標準正規分布の分布関数です．$\alpha$は歪度パラメータで$\alpha=0$のときにこの分布は通常の正規分布となります．
下の図は，非対称正規分布の密度関数をプロットしたものです．（期待値や分散は揃えていませんが）$\alpha$の値によって分布の形状が大きく変わることがわかります．

$SN(\mu,\sigma,\alpha)$に対する混合表現は

$$
Y=\mu+\sigma\delta W + \sigma\sqrt{1-\delta^2} Z,\quad Z\sim N(0,1),
$$

ここで$\delta=\alpha/\sqrt{1+\alpha^2}$（$-1\leq \delta\leq 1$）です．混合に使う$W$の分布は切断正規分布$TN_{(0,\infty)}(0,1)$です（尺度部分には混合に使う確率変数は入っていません）．
混合表現が与えられると，その分布からの乱数の生成を簡単に行うことが可能になります．非対称正規分布の場合，

$W$を$TN_{(0,\infty)}(0,1)$から生成する
生成した$W$を所与に$Y$を$N(\mu+\sigma\delta W, \sigma^2(1-\delta^2))$から生成する

$SN(0,1,1)$から$n=1000$個の乱数生成しました．Rのコードは以下のようになります．

n <- 1000
alpha <- 1.0
delta <- alpha/sqrt(1+alpha^2)
w <- abs(rnorm(n))
y <- rnorm(n, mean=delta*w, sd=sqrt(1-delta^2))

下の図はヒストグラムと密度関数を描いています．上記のかんたんな2ステップのアルゴリズムで，ちゃんと乱数が生成できていることがわかります．

非対称ラプラス分布

前回紹介したラプラス分布も非対称のバージョンがあり，位置尺度混合表現を用いて表すことができます．まず，非対称ラプラス分布$AL(\mu,\sigma,\alpha)$の密度関数は

$$
f(y)=\frac{\alpha(1-\alpha)}{\sigma}\exp\left( -\rho_\alpha\left(\frac{y-\mu}{\sigma}\right)\right),
$$

ここで$\rho_\alpha(u)=u(\alpha-I(u<0))$，$0<\alpha<1$は形状パラメータで，$\alpha=0.5$のときにこの分布は対称になります．下の図は異なる$\alpha$の値のもとでの密度関数で，$\alpha$の値が変わることで分布の形状が大きく変わることがわかります．

$AL(0,1,\alpha)$の位置尺度混合表現は次のように与えられます（Kozumi and Kobayashi, 2011）．

$$
Y = \mu+\theta V +\tau\sqrt{\sigma V} Z,\quad Z\sim N(0,1).
$$

ここで，位置と尺度両方に入っている$V$は平均$\sigma$を持つ指数分布$Exp(\sigma)$に従い，

$$
\theta=\frac{1-2\alpha}{\alpha(1-\alpha)},\quad \tau^2=\frac{2}{\alpha(1-\alpha)}
$$

です．

$$
p(y)=\frac{1}{2\sigma} \exp\left(-\frac{|y-\mu|}{\sigma}\right)
$$

で与えられます．ラプラス分布に対する正規尺度混合表現は，$V$に平均が$1$の指数分布$Exp(1)$を仮定することで得られます．

非対称ラプラス分布は分位点回帰（quantile regression）などに用いられる分布で，この混合表現を使うことでギブスサンプラーやEMアルゴリズムの適用が簡単になることが知られています．

$AL(\mu,\sigma,\alpha)$からの乱数の生成方法は先程と同様，次の2ステップからなります．

$V$を$Exp(\sigma)$から生成する
生成した$V$を所与に，$Y$を$N(\mu+\theta V,\tau^2\sigma V)$から生成する

ここでも$AL(0,1,0.75)$から$n=1000$個の乱数を生成しました．Rのコードは以下のとおりです．

p     <- 0.75
theta <- (1-2*p)/(p-p^2)
tau2  <- 2/(p-p^2)
v     <- rexp(n, rate=1)
y <- rnorm(n,theta*v,sd=sqrt(tau2))

非対称t分布

通常の$t$分布が正規分布とガンマ分布の尺度混合で表されるように，上の非対称正規分布もガンマ分布の尺度混合を用いることで非対称$t$分布とすることができます．まず非対称$t$分布$TS(\mu,\sigma,\alpha,\nu)$の密度関数は，

$$
f(y)=\frac{2}{\sigma}t_\nu(y_s)T_{\nu+1}\left(\alpha\sqrt{\frac{\nu+1}{\nu+y_s}}\right)
$$

で，ここで$y_s=(y-\mu)/\sigma$, $t_\nu$と$T_{\nu+1}$はそれぞれ自由度$\nu$の$t$分布の密度関数と分布関数です．$\alpha$は非対称正規分布のときと同じ，非対称性をコントロールするパラメータです（Azzalini and Capitanio, 2003）．

非対称$t$分布の混合表現は，

$$
Y = \mu + \sigma \frac{X}{\sqrt{V}}
$$

で与えられます．ここで，$X\sim SN(0,1,\alpha)$，$V\sim Ga(\nu/2,\nu/2)$です．非対称正規分布の混合表現から，

$$
Y=\mu+\frac{\sigma}{\sqrt{V}}\delta W + \sigma\sqrt{\frac{1-\delta^2}{V}} Z,\quad Z\sim N(0,1),\quad V\sim Ga(\nu/2,\nu/2),\quad W\sim TN_{(0,\infty)}(0,1)
$$

と書けます．

一般化双曲型分布

一般化双曲型線（generalised hyperbolic, GH）分布は左右対称で裾の厚い分布として知られています．GH分布は一般化逆ガウス分布（$GIG(\lambda,\delta,\gamma)$）と呼ばれる分布に従う確率変数との位置尺度混合で表現することができます：

$$
Y = \mu + \beta V + \sqrt{V}Z,\quad V\sim GIG(\lambda,\delta,\gamma),\quad Z\sim N(0,1).
$$

Nakajima and Omori (2012)では，GH分布の特別な場合であるGH非対称$t$分布を用いて，金融時系列データの収益分布のモデリングを行っています．GH非対称$t$分布はGH分布の$\lambda=-\nu/2$，$\delta=\sqrt{\nu}$，$\gamma=0$の場合に対応し，上と同様に

$$
Y = \mu + \beta V + \sqrt{V}Z,\quad V\sim IG(\nu/2,\nu/2),\quad Z\sim N(0,1).
$$

で表されます．

一般化非対称ラプラス分布

混合させる確率変数（上記の$V$や$W$）の入れ方や入れる数によって，様々な柔軟な分布を作り出すことができます．
例えば，Yan and Kottas (2018)は上の非対称ラプラス分布分布の一般化として，次のような混合表現に基づいた確率分布を提案しました．

$$
Y=\mu+\theta V + \alpha S + \tau\sqrt{\sigma V} Z,
$$

ここで，$V\sim Exp(\sigma)$, $S\sim TN_{(0,\infty)}(0,\sigma^2)$, $Z\sim N(0,1)$です．この表現から得られる$Y$の周辺分布の密度関数はかなり複雑な形となるので省きますが，非対称ラプラス分布よりも柔軟な形状を表現できるとされています．

おわりに

統計モデリングにおいて，$t$分布や非対称分布など，文脈に応じて正規分布以外の分布を使用しなければならない場合，使用する分布の混合表現を用いることは有用である場合が多いです．例えば単純に乱数を生成する場合でも，棄却サンプリングなどを使って直接ターゲットの分布からサンプリングするよりも，上で紹介したような混合表現に基づいた生成方法のほうが簡単に実行することができます．

データサイエンス研修，研究に必要なデータ収集やります！
株式会社Nospareではベイズ統計学に限らず統計学の様々な分野を専門とする研究者が所属しており，新たな知見を日々追求しています．統計アドバイザリー，ビジネスデータの分析，データサイエンス研修につきましては弊社までお問い合わせください．また主に社会科学分野における実証研究に必要なデータ収集サービスも行っておりますので，こちらもぜひお問い合わせください．

参考文献

Azzalini, A. (1986). Further results on a class of distributions which includes the normal ones. Statistica, 46, 199-208.
Azzalini, A. and Capitanio, A. (2003). Distributions generated by perturbation of symmetry with emphasis on a multivariate skew t-distribution. Journal of the Royal Statistical Society, Series B, 65, 367–389.
Kozumi, H. and Kobayashi, G. (2011). Gibbs sampling methods for Bayesian quantile regression. Journal of Statistical Computation and Simulation, 81, 1565-1578.
Nakajima, J. and Omori, Y. (2012). Stochastic volatility model with leverage and asymmetrically heavy-tailed error using GH skew Student’s $t$-distribution. Computational Statistics & Data Analysis, 56 3690-3704.
Yan, Y. and Kottas, A. (2017). A New Family of Error Distributions for Bayesian Quantile Regression, arXiv:1701.05666.