多変量t分布の導出とpythonによる乱数生成の実装

はじめに

多変量t分布について調べていたので，そのまとめです．

（この記事を書き上げた瞬間，データを吹っ飛ばして書き直しました...死にたい...）

多変量t分布の導出

定義

ググって一番最初にみつかるという理由から，wikipediaのMultivariate t-distributionのノーテーションにに従います．
多変量t分布の確率密度関数は，以下で定義されます．

$$
p(\textbf{x}) = \frac{\Gamma\Bigl[{\frac{\nu+p}{2}}\Bigr]}{\Gamma\Bigl[{\frac{\nu}{2}}\Bigr] \nu^{\frac{p}{2}}\pi^{\frac{p}{2}} \det(\Sigma)^{\frac{1}{2}}} \Bigl[ 1 + \frac{1}{\nu} (\textbf{x} - {\bf \mu})^{T} {\bf \Sigma}^{-1}(\textbf{x} - {\bf \mu})\bigr]^{\frac{-(\nu+p)}{2}}. \tag{1}
$$

理論

ポイント

そもそもt分布とはなんぞやという話になりますが，t分布とは，正規分布の分散を増減させるパラメーターを考え，それにガンマ分布を仮定し，結局は積分消去したもの，と言えます．
ここで，ガンマ分布ではなく逆ガンマ分布，カイ2乗分布を仮定することもありますが，少しの変更で全て同値になりそうです．（試してない）

これによって，分散の値がよくわからないときに，その平均的な振る舞いを見ることができます．

導出

二つの確率変数 $V \sim Gam(\alpha, \beta)$，$X \sim N(\textbf{0}, {\bf \Sigma})$ を考えます．$X$は$d$次元正規分布とします．
この二つの確率変数を，関数 $f,g$ によって，$Y,U$に変換します．
$$
Y = f(V,X) = \mu + \frac{1}{\sqrt{V}}X, \tag{2}
$$
$$
U = g(V,X) = g(V)= V. \tag{3}
$$

ここで，$g$は$V$の恒等関数であり，$\mu \in \mathbb{ R }^{d}$は平均を表すd次元ベクトル，$V$は$X$と独立と仮定します．

この$f(\cdot)$は，$V$を固定すると直線（平面）の関数であり，1対1対応の関数です．従って，逆関数$f^{-1}(\cdot)$が存在し，
$$
f^{-1}(Y) = X = \sqrt{U}(Y - \mu), \tag{4}
$$
となります．

すでに分布の分かっている$V,X$から，$Y$の分布を求めます．確率変数は$U$と$Y$の2つあるので，$U,Y$の同時分布を考えて，$U$について積分消去することで$Y$の分布を求めます．

Pr(u,y) = Pr(U < u, Y < y) = Pr(g(V) < u,f(X) < y) \\
= Pr(V < g^{-1}(u), X < f^{-1}(y)) \\
= Pr(V < u, X<\sqrt{u}(y - \mu)) = Pr(V < u) Pr(X < \sqrt{u}(y - \mu)). \tag{5}

ここで，$V$と$X$は独立の仮定，
$$
Pr(v,x) = Pr(v)Pr(x), \tag{6}
$$
を使いました．

$p_{V},p_{X}$をそれぞれ$V$と$X$の確率密度関数とすると，
$$
Pr(u,y) = Pr(V < u) Pr(X < \sqrt{u}(y - \mu)) = \int_{-\infty}^{u} \int_{-\infty}^{y} p_{V}(u)p_{X} (\sqrt{u}(y - \mu)) |\det(J)| dudy. \tag{7}
$$

ここで，$|\det(J)|$は変数変換によって生じたヤコビアンであり，ヤコビ行列$J$ （$(d+1)*(d+1)$行列）は，

J = \begin{pmatrix} \frac{\partial x}{\partial y^{T}} & \frac{\partial x}{\partial u} \\ \frac{\partial v}{\partial y^{T}} & \frac{\partial v}{\partial u} \end{pmatrix} =
 \begin{pmatrix} \sqrt{u}I_{d} & \frac{1}{2\sqrt{u}(y-\mu)} \\ \textbf{0} & 1 \end{pmatrix}. \tag{8}

ヤコビ行列$J$が上三角行列なので，その行列式は対角要素の積となり，
$$
\det(J) = u^{\frac{d}{2}}. \tag{9}
$$
　　
$(7)$の積分の中身は，$Pr(u,y)$の確率密度関数$p_{U,Y}(u,y)$であるため，これを求めます．
$V \sim Gam(\alpha, \beta)$，$X \sim N(\textbf{0}, {\bf \Sigma})$より，
$$
p_{V}(v) = \frac{\beta^{\alpha}}{\Gamma(\alpha)} v^{\alpha-1} \exp{( - \beta v)}, \tag{10}
$$
$$
p_{X}(\textbf{x}) = \frac{1}{(2\pi)^{\frac{d}{2}}} \frac{1}{\sqrt{\det(\Sigma)}} \exp{\Bigl( - \frac{1}{2} \textbf{x}^{T} \Sigma^{-1} \textbf{x} \Bigr)}, \tag{11}
$$
従って，
$$
p_{V}(u) = \frac{\beta^{\alpha}}{\Gamma(\alpha)} u^{\alpha-1} \exp{( - \beta u)} = Gam(\alpha, \beta), \tag{12}
$$

p_{X}(\sqrt{u}(\textbf{y} - \mu)) = \frac{1}{(2\pi)^{\frac{d}{2}}} \frac{1}{\sqrt{\det(\Sigma)}} \exp{\Bigl( - \frac{1}{2} \sqrt{u}(\textbf{y} - \mu)^{T} \Sigma^{-1} \sqrt{u}(\textbf{y} - \mu) \Bigr)} \\
= \frac{1}{(2\pi)^{\frac{d}{2}}} u^{-\frac{d}{2}} u^{\frac{d}{2}} \frac{1}{\sqrt{\det(\Sigma)}} \exp{\Bigl( - \frac{1}{2} (\textbf{y} - \mu)^{T} (u \Sigma^{-1}) (\textbf{y} - \mu) \Bigr)} \\
=  u^{-\frac{d}{2}} \frac{1}{(2\pi)^{\frac{d}{2}}}  \frac{1}{\sqrt{\det(\frac{\Sigma}{u})}} \exp{\Bigl( - \frac{1}{2} (\textbf{y} - \mu)^{T} (u \Sigma^{-1}) (\textbf{y} - \mu) \Bigr)} = u^{-\frac{d}{2}} N(\mu, \frac{\Sigma}{u}). \tag{13}

$Pr(u,y)$の確率密度関数$p_{U,Y}(u,y)$は，
$$
p_{U,Y}(u,y) = Gam(\alpha, \beta) u^{-\frac{d}{2}} N(\mu, \frac{\Sigma}{u}) |\det(J)| = Gam(\alpha, \beta) N(\mu, \frac{\Sigma}{u}) u^{-\frac{d}{2}}u^{\frac{d}{2}} = Gam(\alpha, \beta) N(\mu, \frac{\Sigma}{u}). \tag{14}
$$
（ここで，ガンマ分布の台は$u \in [0,\infty), d>0$より，絶対値は外してよい．）

$p_{U,Y}(u,y)$を$u$について積分することによって，$Y$の確率密度関数$p_{Y}$を得ます．
$$
p_{Y} = \int_{0}^{\infty} p_{U,Y}(u,y)du = \int_{0}^{\infty} Gam(\alpha, \beta) N(\mu, \frac{\Sigma}{u})du \tag{15}
$$

ここで，記法の簡略化のために，正規分布のマハラノビスの距離を，$\Delta^{2} = (\textbf{y} - \mu)^{T}\Sigma^{-1}(\textbf{y} - \mu)$と表記します．

$$
p_{Y}(\textbf{y}) = \int_{0}^{\infty} \frac{\beta^{\alpha}}{\Gamma(\alpha)} u^{\alpha-1} \exp{( - \beta u)} u^{-\frac{d}{2}} \frac{1}{(2\pi)^{\frac{d}{2}}} \frac{1}{\sqrt{\det(\frac{\Sigma}{u})}} \exp{\Bigl( - \frac{u}{2} \Delta^{2} \Bigr)} du. \tag{16}
$$
$u$に依存した部分のみを抜き出し，
$$
p_{Y}(\textbf{y}) \propto \int_{0}^{\infty} u^{\alpha-1} \exp{( - \beta u)} u^{\frac{d}{2}} \exp{\Bigl( - \frac{u}{2} \Delta^{2} \Bigr)} du \
= \int_{0}^{\infty} u^{\frac{d}{2}+\alpha-1} \exp{\Bigl( - u( \beta + \frac{\Delta^{2}}{2} ) \Bigr)} du. \tag{17}
$$

ここで，$z = u(\beta + \frac{\Delta^{2}}{2})$とおくと，$u = z(\beta + \frac{\Delta^{2}}{2})^{-1}$，また，$du = (\beta + \frac{\Delta^{2}}{2})^{-1} dz$．$u=0$のとき，$z=0$，$u \to \infty$で$z \to \infty$．

$$
\int_{0}^{\infty} z^{\frac{d}{2} + \alpha -1} ( \beta + \frac{\Delta^{2}}{2} )^{-\frac{d}{2} -\alpha +1} \exp(-z) (\beta + \frac{\Delta^{2}}{2})^{-1} dz \
= ( \beta + \frac{\Delta^{2}}{2} )^{-\frac{d}{2} -\alpha} \int_{0}^{\infty} z^{(\frac{d}{2} + \alpha) -1} \exp(-z) dz.
$$
ガンマ関数の定義より，
$$
p_{Y}(\textbf{y}) \propto \Bigl( \beta + \frac{\Delta^{2}}{2} \Bigr)^{-\frac{d}{2} -\alpha} \Gamma \Bigl( \frac{d}{2} + \alpha \Bigr) \tag{18}
$$

$(17)$の$u$に依存しない項も考慮して，

p_{Y}(\textbf{y}) = \frac{\beta^{\alpha}}{\Gamma(\alpha)} \frac{1}{(2\pi)^{\frac{d}{2}}}  \frac{1}{\sqrt{\det(\Sigma)}} \Bigl( \beta + \frac{\Delta^{2}}{2} \Bigr)^{-\frac{d}{2} -\alpha}  \\
= \frac{\Gamma \Bigl( \frac{d}{2} + \alpha \Bigr)}{\Gamma(\alpha)}\frac{\beta^{\alpha}}{(2\pi)^{\frac{d}{2}}} \frac{1}{\sqrt{\det(\Sigma)}} \Bigl( \beta + \frac{\Delta^{2}}{2} \Bigr)^{-\frac{d}{2} -\alpha} \Bigl( \frac{1}{\beta} \Bigr)^{\frac{d}{2} + \alpha} \beta^{-\frac{d}{2}-\alpha} \\
= \frac{\Gamma \Bigl( \frac{d}{2} + \alpha \Bigr)}{\Gamma(\alpha)}\frac{1}{(2 \beta \pi)^{\frac{d}{2}}} \frac{1}{\sqrt{\det(\Sigma)}} \Bigl( 1 + \frac{\Delta^{2}}{2 \beta} \Bigr)^{-\frac{d}{2} -\alpha}. \tag{19}

ここで，$\alpha = \beta = \frac{\nu}{2}$，$p = d$とおくと，
$$
p_{Y}(\textbf{y}) = \frac{\Gamma \Bigl( \frac{\nu+p}{2} \Bigr)}{\Gamma \Bigl( \frac{\nu}{2} \Bigr)}\frac{1}{(\nu \pi)^{\frac{p}{2}}} \frac{1}{\sqrt{\det(\Sigma)}} \Bigl( 1 + \frac{\Delta^{2}}{\nu} \Bigr)^{-\frac{(\nu+p)}{2}}, \tag{20}
$$
となり，$(1)$のt分布の定義と一致します．$\Box$

実践

Pythonによる乱数発生の実装

$(2)$式がそのまま乱数の生成方法なので，これをpythonで実装します．

# -*- coding: utf-8 -*-
#!/usr/bin/python

import numpy as np
import matplotlib.pyplot as plt

#パラメーターをセットする．
N = 5000
d = 10
df = 4

mean = np.zeros(d)
cov = np.eye(d)

#ガンマ分布からの乱数生成
V = np.random.gamma(shape=df/2., scale=df/2.,size=N)

#多変量正規分布からの乱数生成
X = np.random.multivariate_normal(mean=mean, cov=cov, size=N)

#Xを変換する．
V = V.reshape((N,1))
denom = np.sqrt(np.tile(V,d))
Y = mean + X / denom


#比較用の多変量正規分布からの乱数を発生させる．
#（疑問:比較用はこんな感じでいいのだろうか？）
cov_ = ((df - 2) / float(df)) * np.corrcoef(Y.T)
X_ = np.random.multivariate_normal(mean=mean, cov=cov_, size=N)

#プロットする．
fig = plt.figure()
ax = fig.add_subplot(111)
ax.hist(X_[:,0], color="blue", normed=True, label="Normal dist", bins=30)
ax.hist(Y[:,0], color="red", normed=True, label="t dist", bins=30, alpha=0.5)
ax.set_title("Comparison of normal and t distribution")
ax.legend()
fig.show()

fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(X_[:,0], X_[:,1],color="blue",label="Normal dist")
ax.scatter(Y[:,0], Y[:,1],color="red",label="t dist",alpha=0.5)
ax.set_title("scatter plot of t and normal dist")
ax.legend()
fig.show()

分散も大きいしt分布っぽいかな？
（一次元ずつ以外の比較もしたかったけど，書き直して力尽きたのでこのくらいで...）

最後に

雰囲気数学マンなので間違っていたら教えてください．
時間があればいろいろな分布のつながりを追ってみたいなあ．