任意の確率分布に従う標本生成：Markov連鎖 Gibbs Sampling

　Gibbs Samplingは、Markov連鎖のMonteCarlo法を用いて、多次元標本の次元毎に標本作成を行う方法で、複雑な確率分布に従う標本作成において、強力なツールである。
互いに依存し合う確率変数$${{\bf x}=(x_1,\cdots,x_d)}$$が、それぞれの平均値$${{\bf \mu}=(\mu_1,\cdots,\mu_d)}$$を持ち、分散共分散行列が
$${\begin{pmatrix}\sigma_{1} &\sigma_{12} & \cdots & \sigma_{1d} \\ \sigma_{21} &\sigma_{2} & \cdots & \sigma_{2d} \\ \vdots & \vdots & \ddots & \vdots \\ \sigma_{1d} &\sigma_{2d} & \cdots & \sigma_{d} \\ \end{pmatrix} }$$
で与えられる$${P({\bf x})}$$の確率分布に従うとする。
この$${d}$$次元の確率変数$${{\bf x}}$$の生成に、Gibbs Samplingを適用すると、
ある$${t}$$における$${x_j}$$の値を$${x_{t,j}}$$と表せば、Markov連鎖を使用することから、これは、$${x_{t-1,j}}$$と、$${j}$$以外の$${x_i, i\neq j}$$に依存することになる。$${x_j}$$前にすでに$${t}$$での更新を終えている変数をそのまま使い、$${x_j}$$は、
$${Pr(x_{t,j}|x_{t,1},\cdots,x_{t,j-1},x_{t-1,j+1},\cdots,x_{t-1,d})}$$
に従いサンプリングされる。
これを繰り返せば、最終的な$${{\bf x}}$$の分布は、$${Pr({\bf x})}$$に定常となる。
　これを簡単化のため、$${x,y}$$の2変数でpythonで実装してみる。
　$${x,y}$$はそれぞれ、平均値$${\mu_x}$$、$${\mu_y}$$を持ち、分散巨分散行列は、
$${\begin{pmatrix}\sigma_{x} &\sigma_{xy} \\ \sigma_{xy} &\sigma_{y} \\ \end{pmatrix} }$$
から、これを用いて、確率分布は
$${P(x,y)=\displaystyle{\frac{1}{2\sigma_x\sigma_y\sqrt{1-\rho_{xy}^2}} }}$$
$${\displaystyle{ \cdot \exp\Big[ - \frac{1}{2(1-\rho_{xy}^2)} \Big(\frac{(x-\mu_x)^2}{\sigma^2}-2\rho_{xy}\frac{(x-\mu_x)(y-\mu_y)}{\sigma_x\sigma_y} + \frac{(y-\mu_y)^2}{\sigma_y^2} \Big)\Big]}}$$
と、与えられているとする。ここで、$${\rho_{xy}}$$は$${x}$$と$${y}$$の相関係数であり、$${\rho_{xy}=\displaystyle{\frac{\sigma_{xy}}{\sigma_x\sigma_y}}}$$である。
互いの条件付き確率は、
$${P(x|y)=\mathcal{N}\Big(\mu_x+\rho_{xy}\frac{\sigma_x}{\sigma_y}(y-\mu_y),(1-\rho_{xy}^2)\sigma_x^2\Big)}$$
で与えられるとし、$${P(x|y)=P(y|x)}$$である。
これは他変量正規分布に他ならない。
この条件付き確率は、２次元と仮定して、以下のように実装される。

def Pxy(x,mean,cov,j):
    #P(x|y)
    #x:return
    #y:given condition
    #Note now j = 0 or 1
    mu = mean[j] + cov[0, j] / cov[j, j] * (x[1-j] - mean[1-j])
    sigma = cov[j, j] - cov[0, j] / cov[1-j, 1-j] * cov[0, j]
    return np.random.normal(mu, sigma)

これを用いたGibbsSampling の関数の実装は以下となる。

def GibbsSampler( mean, cov,itrnum):
    x=rand.multivariate_normal(mean,np.diag(np.diag(cov)))
    
    gs= np.zeros([itrnum,len(mean)])

    for i in range(itrnum):
        for j in range(len(mean)):
            xj = Pxy(x,mean,cov,j)
            x[j]=xj
            gs[i,j]= xj

    return gs

初期値の$${{\bf x}}$$を平均値と分散からの正規分布で与えている。
これを、numpy.randomのmultivariate_normalと比較したのが以下の結果である。

mean = (1, 2)
cov = np.array([[1, 0.7], [0.7, 1]])
gs=GibbsSampler(mean,cov,randnum)
print(gs[:,0].mean(),gs[:,1].mean(),np.cov(gs.T))
plt.plot(gs[:, 0], gs[:, 1], '.', alpha=1.0,label='GS')
x=rand.multivariate_normal(mean, cov, 800_000)
print(x[:,0].mean(),x[:,1].mean(),np.cov(x.T))
plt.plot(x[:, 0], x[:, 1], '.', alpha=0.08,label='numpy.random')
plt.grid()
plt.legend()
plt.show();