文章标题

The Seventh Homework

标签（空格分隔）： statistical_computation

---

1.problem

Complete the example “Bivariate Normal Covariance Matrix” on the page 93 of the textbook .

2.solution

step a:
a1. Draw Σ from the current estimate of the observed posterior.
a2. Generate the latent data.

step b:
Set the posterior density of Σ equal to the mixture of inverted Wishart distribution.

###the original plotf=function(x)(1-x^2)^4.5/(1.25-x^2)^8s=integrate(f,-1,1)$valuex=seq(-1,1,length=100000)plot(x,p=f(x)/s,type="l")##using the data in the bookx=matrix(c(1,1,-1,-1,2,2,-2,-2,0,0,0,0,1,-1,1,-1,0,0,0,0,2,2,-2,-2),2,12,byrow=T)n=12m=6400##data augmentationlibrary(MCMCpack)A=x%*%t(x)v=12MM=matrix(0,2*m,12)MMn=matrix(0,2*m,12)for(l in 1:m){MM[c(2*l-1,2*l),]=x}for(h in 1:15){for(j in 1:m){k=sample.int(m,1,replace=T)A=MM[c(2*k-1,2*k),]%*%t(MM[c(2*k-1,2*k),])sigma=riwish(v,A)for(i in 5:8)x[2,i]=rnorm(1,sigma[1,2]/sigma[1,1]*x[1,i],sqrt(det(sigma)/sigma[1,1]))for(i in 9:12)x[1,i]=rnorm(1,sigma[1,2]/sigma[2,2]*x[2,i],sqrt(det(sigma)/sigma[2,2]))MMn[c(2*j-1,2*j),]=x}MM=MMn}##plot the densityp=c()for(i in 1:10000){k=sample.int(m,1,replace=T)A=MM[c(2*k-1,2*k),]%*%t(MM[c(2*k-1,2*k),])sigma=riwish(v,A)p[i]=sigma[1,2]/sqrt(sigma[1,1]*sigma[2,2])}points(density(p),type="l",col="red")

Then we can obtain the following gragh:

The red line is the result of the augmentation algorithm, the black line is the real posterior of ρ.