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\)
> 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)\]
> 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
implementationHere 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