🟠

ベイジアンブートストラップ

2023/10/28に公開

はじめに

今回は大規模計算時代の統計推論10.3のベイズ的ブートストラップについて説明していきます。

$G_1,G_2,...,G_n \text{ i.i.d. } \text{Gam}(1,1)$ のとき

\bm P^*=(G_1,G_2,...,G_n)'\bigg/\sum_{i=1}^nG_i \tag{10.44}

この式の平均と共分散行列を求める。

$\bm P_0=(1/n,...,1/n)'$ とする。

\bm P^*\sim \left[\bm P_0 , \frac{1}{n+1}(\text{diag}(\bm P_0)-\bm P_0 \bm P_0' )\right]\tag{10.45}

これを解く鍵は、ガンマ分布とディリクレ分布の関係を理解すること。

必要な知識を先にまとめておく。

ガンマ分布

X \sim \text{Gam}(a,b) のとき\\ f(x) = \frac{1}{\Gamma(a)b^a}x^{a-1}e^{-x/b} \hspace{5pt},\hspace{5pt} x>0

ディリクレ分布

パラメータを $\bm{\alpha}=(\alpha_1,...,\alpha_n)$ , 確率変数を $\bm x = (x_1,...,x_n)$ とするとき

f(\bm x ; \bm \alpha)=\frac{\Gamma({\sum_{i=1}^n}\alpha_i)}{\prod_{i=1}^n \Gamma(\alpha_i)}\prod_{i=1}^nx_i^{\alpha_i-1} \hspace{5pt},\hspace{5pt} x_i≥0 \hspace{5pt},\hspace{5pt}\sum_{i=1}^nx_i=1

一般的に考えるために次にようなものを考える。記号も変えているので注意。

\begin{align*} X_i &\sim \text{Gam}(\alpha_i,1)\\ Y_i &=\frac{X_i}{\sum_{j=1}^nX_j}\hspace{5pt},\hspace{5pt}\sum_{i=1}^nY_i=1 \end{align*}

こうすると $P^*=(Y_1,...,Y_n)$ である。

データは次のように書ける。

y_i=\frac{x_i}{\sum_{j=1}^nx_j} \hspace{5pt},\hspace{5pt}\sum_{i=1}^ny_i=1

逆変換を考えると、

x_i=y_i\sum_{j=1}^n x_j = y_iz \hspace{5pt},\hspace{5pt}z=\sum_{j=1}^nx_j

ヤコビアンを考える

\begin{align*} J &= \left|\begin{matrix}\partial x_1 / \partial y_1 & \cdots &\partial x_1/ \partial y_{n-1} & \partial x_1 / \partial z\\ & \vdots\\ \partial x_{n-1} / \partial y_1 & \cdots &\partial x_{n-1}/ \partial y_{n-1} & \partial x_{n-1} / \partial z\\ \partial x_n / \partial y_1 & \cdots &\partial x_n/ \partial y_{n-1} & \partial x_n / \partial z \end{matrix} \right| \\ &=\left|\begin{matrix} z & 　 & 　 & & y_1\\ & \ddots & \text{\huge{0}} \\ & \text{\huge{0}} & \ddots & & \vdots \\ & & & z & y_{n-1}\\ -z & \cdots & \cdots & -z & y_n \end{matrix} \right| \\ &= \left|\begin{matrix} z & & & & y_1\\ & \ddots & \text{\huge{0}}& \\ & \text{\huge{0}} & \ddots & & \vdots\\ & & & z & y_{n-1}\\ 0 & &\cdots & 0 & 1 \end{matrix} \right| \\ & = z^{n-1} \end{align*}

$z≥ 0$ より $|J|=z^{n-1}$ である。

次に、 $y_1,...,y_{n-1}$ の同時確率密度関数を計算する。周辺化とガンマ分布を使って次のように計算できる。

\begin{align*} f(y_1,...,y_{n-1}) &= \int_0^{\infty} f(y_1,...,y_{n-1},z) dz\\ &= \int_0^{\infty} \prod_{i=1}^{n}\left\{\frac{1}{\Gamma(\alpha_i)}(y_iz)^{\alpha_i-1} e^{-y_i z} \right\} z^{n-1} dz\\ &= \frac{1}{\Gamma(\alpha_1)\cdots \Gamma(\alpha_n)}\prod_{i=1}^n y_i^{\alpha_i-1} \int_0^{\infty} z^{\sum_{i=1}^n \alpha_i-n }e^{z\sum_{i=1}^ny_i} z^{n-1} dz\\ &= \frac{1}{\Gamma(\alpha_1)\cdots \Gamma(\alpha_n)}\prod_{i=1}^n y_i^{\alpha_i-1} \int_0^{\infty} z^{\sum_{i=1}^n \alpha_i-1}e^{z} dz\\ &= \frac{\Gamma(\sum_{i=1}^n \alpha_i)}{\Gamma(\alpha_1)\cdots \Gamma(\alpha_n)}\prod_{i=1}^n y_i^{\alpha_i-1} \end{align*}

よって $y_1,...,y_{n-1}$ はディリクレ分布に従うことがわかる。

したがって式(10.45)を求めるには、ディリクレ分布の期待値と共分散行列を求めれば良いことがわかった。

期待値は

E[Y_i]=\frac{\alpha_i}{\sum_{i=1}^n\alpha_i} \tag{1}

分散は

Var[Y_i]=\frac{\alpha_i(\sum_{i=1}^n\alpha_i-\alpha_i)}{(\sum_{i=1}^n \alpha_i+1)(\sum_{i=1}^n \alpha_i)^2} \tag{2}

共分散は

Cov[Y_i,Y_j] = -\frac{\alpha_i\alpha_j}{(\sum_{i=1}^n \alpha_i)^2(\sum_{i=1}^n \alpha_i+1)} \tag{3}

$\bm \alpha=(1,...,1)$ のとき

$(1)$ より $\bm P^*=(P_1^*,...,P_n^*)'$ の期待値は

E[P_i^*]= 1/n

よって

E[\bm P^*]=\bm P_0 \tag{4}

$(2),(3)$ より $\bm P^*$ の共分散行列は

\begin{align*} Var[P_i^*]=\frac{n-1}{(n+1)n^2} \hspace{5pt},\hspace{5pt} Cov[P_i^*,P_j^*]= - \frac{1}{(n+1)n^2} \end{align*}

よって

Var[\bm P^*]= \frac{1}{n+1} (\text{diag}(\bm P_0)-\bm P_0 \bm P_0') \tag{5}

よって示すことができた。

ログインするとコメントできます