Sampling Multivariate Gaussian distribution with a Gibbs Sampler

Computing conditional distributions

Let’s consider the following matrix \(\Sigma\)

\[\Sigma=\left(\begin{array}{ccc} 2 & 1 & 0 \\ 1 & 2 & 1 \\ 0 & 1 & 1 \end{array}\right)\]

Sigma=matrix(c(2,1,0,1,2,1,0,1,2),3,3)
Sigma
##      [,1] [,2] [,3]
## [1,]    2    1    0
## [2,]    1    2    1
## [3,]    0    1    2
eigen(Sigma)$values
## [1] 3.4142136 2.0000000 0.5857864

\(\Sigma\) is then definite positive definite matrix and our aim is to simulate the mulitivariate Gaussian distribution with zero mean and covariance matrix \(\Sigma\)

\[(\theta_1,\theta_2,\theta_3)\sim\mathcal{N}_3(0,\Sigma)\]

We will use here the Gibbs sampler and we need then to compute the conditional distributions : \(\theta_1\mid\theta_2,\theta_3\), \(\theta_2\mid\theta_1,\theta_3\) and \(\theta_3\mid\theta_1,\theta_2\).

Let us recall the famous Theorem on the conditional distribution of Multivariate Gaussian Distributions:

Theorem Let \(X\sim\mathcal{N}(\mu,\Sigma)\) be a Gaussian random vector in \(\mathbb{R}^p\). \(A\subset \{1,\ldots,p\}\), \(X_A\) be subvector of \(X\) indexed by \(A\), and \(X_{\overline{A}}\) be subvector of \(X\) indexed by \(\overline{A}=\{1,\ldots,p\}\setminus A\). Then, if we decompose \(\Sigma\) and \(\mu\) as follows \[\mu=\left(\begin{array}{c}\mu_A\\ \mu_{\overline{A}}\end{array}\right)\] and \[ \Sigma=\left(\begin{array}{cc} \Sigma_A & \Sigma_{A \overline{A}} \\ \Sigma_{\overline{A}A} & \Sigma_{\overline{A}} \\ \end{array}\right)\] the random vector \(X_A\mid X_{\overline{A}}=x_{\overline{A}}\) has a Gaussian distribution such that \[\left(X_A\mid X_{\overline{A}}=x_{\overline{A}}\right)\sim \mathcal{N}(\mu_{A\mid \overline{A}},\Sigma_{A\mid \overline{A}})\] where

\[\mu_{A\mid \overline{A}}=\mu_A+ \Sigma_{A \overline{A}}\Sigma_{\overline{A}}^{-1}(x_{\overline{A}}-\mu_{\overline{A}})\]

and \[ \Sigma_{A\mid \overline{A}}=\Sigma_A- \Sigma_{A \overline{A}} \Sigma_{\overline{A}}^{-1} \Sigma_{\overline{A}A}.\]

Let us then applied the previous theorem to compute the conditional distribution of \(\theta_1\mid \theta_2,\theta_3\). Let us consider \(A=\{1\}\) and \(\overline{A}=\{2,3\}\). Then

> sigmaAAb=2-t(c(1,0))%*%solve(matrix(c(2,1,1,2),2,2))%*%c(1,0)
> sigmaAAb
##          [,1]
## [1,] 1.333333
> t(c(1,0))%*%solve(matrix(c(2,1,1,2),2,2))
##           [,1]       [,2]
## [1,] 0.6666667 -0.3333333

Then \[\mu_{1\mid 2,3}=.6666\theta_2-.3333 \theta_3\] Finaly we conclude that \[\theta_1\mid \theta_2,\theta_3\sim\mathcal{N}(.6666\,\theta_2-.3333\, \theta_3,\sqrt{1.333333})\]

Let’s do the same for \(\theta_2\mid \theta_1,\theta_3\) and \(\theta_3\mid \theta_1,\theta_2\)

  • for \(\theta_2\mid \theta_1,\theta_3\), \(A=\{2\}\) and \(\overline{A}=\{1,3\}\).
> sigmaAAb=2-t(c(1,1))%*%solve(matrix(c(2,0,0,2),2,2))%*%c(1,1)
> sigmaAAb
##      [,1]
## [1,]    1
> t(c(1,1))%*%solve(matrix(c(2,0,0,2),2,2))
##      [,1] [,2]
## [1,]  0.5  0.5

Then \[\mu_{2\mid 1,3}=.5\theta_1+.5 \theta_3\] Finaly we conclude that \[\theta_2\mid \theta_1,\theta_3\sim\mathcal{N}(.5\,\theta_1+.5\, \theta_3,1)\]

  • for \(\theta_3\mid \theta_1,\theta_2\), \(A=\{3\}\) and \(\overline{A}=\{1,2\}\).
> sigmaAAb=2-t(c(0,1))%*%solve(matrix(c(2,1,1,2),2,2))%*%c(0,1)
> sigmaAAb
##          [,1]
## [1,] 1.333333
> t(c(0,1))%*%solve(matrix(c(2,1,1,2),2,2))
##            [,1]      [,2]
## [1,] -0.3333333 0.6666667

Then \[\mu_{3\mid 1,2}=-.3333 \,\theta_1+ .6666\, \theta_2\] Finaly we conclude that \[\theta_3\mid \theta_1,\theta_2\sim\mathcal{N}(-.3333 \,\theta_1+ .6666\, \theta_2,\sqrt{1.333333})\]

R implementation

Here it’s the R code that will allow use to obtain a sample of \((\theta_1,\theta_2,\theta_3)\) using the Gibbs sampler:

> set.seed(55679)
> theta0=rnorm(3)
> N=1000
> theta=matrix(NA,ncol = 3,nrow = N)
> theta=rbind(theta0,theta)
> for(i in 1:N){
+   theta1<-rnorm(1,mean = .6666667*theta[i,2]-.33333*theta[i,3],sd = sqrt(1.3333))
+   theta2<-rnorm(1,mean = .5*theta1+.5*theta[i,3],sd = 1)
+   theta3<-rnorm(1,mean = -.33333*theta1+.66666*theta2,sd = sqrt(1.3333))
+   theta[i+1,]=c(theta1,theta2,theta3)
+ }

Let’s check some characteristics of the distribution of the sample

> colnames(theta)=c("theta1","theta2","theta3")
> theta=as.data.frame(theta[501:N,])
> m1=lm(theta1~theta2+theta3,data=theta)
> summary(m1)
## 
## Call:
## lm(formula = theta1 ~ theta2 + theta3, data = theta)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.3642 -0.7880  0.0400  0.6766  3.8589 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.02611    0.05199   0.502    0.616    
## theta2       0.62595    0.04470  14.002   <2e-16 ***
## theta3      -0.37615    0.04263  -8.823   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.16 on 497 degrees of freedom
## Multiple R-squared:  0.287,  Adjusted R-squared:  0.2841 
## F-statistic:   100 on 2 and 497 DF,  p-value: < 2.2e-16
> m1$coef
## (Intercept)      theta2      theta3 
##  0.02611218  0.62594589 -0.37614945
> m2=lm(theta2~theta1+theta3,data=theta)
> summary(m2)
## 
## Call:
## lm(formula = theta2 ~ theta1 + theta3, data = theta)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.5278 -0.6562  0.0188  0.6822  2.3690 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.002279   0.044186  -0.052    0.959    
## theta1       0.451920   0.032276  14.002   <2e-16 ***
## theta3       0.516927   0.031307  16.511   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9856 on 497 degrees of freedom
## Multiple R-squared:  0.4674, Adjusted R-squared:  0.4653 
## F-statistic: 218.1 on 2 and 497 DF,  p-value: < 2.2e-16
> m2$coef
##  (Intercept)       theta1       theta3 
## -0.002279313  0.451920201  0.516927367
> m3=lm(theta3~theta1+theta2,data=theta)
> summary(m3)
## 
## Call:
## lm(formula = theta3 ~ theta1 + theta2, data = theta)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.5586 -0.7140  0.0092  0.7020  3.6074 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.05641    0.05081  -1.110    0.267    
## theta1      -0.36001    0.04080  -8.823   <2e-16 ***
## theta2       0.68526    0.04150  16.511   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.135 on 497 degrees of freedom
## Multiple R-squared:  0.3579, Adjusted R-squared:  0.3553 
## F-statistic: 138.5 on 2 and 497 DF,  p-value: < 2.2e-16
> m3$coef
## (Intercept)      theta1      theta2 
## -0.05640736 -0.36000971  0.68526468