## Bibby(Mi) Zhou, STAT203 Presentation code#####
install.packages("BayesSummaryStatLM")
install.packages("mvnfast")
library("coda")
library(BayesSummaryStatLM)
library(mvnfast)
###### simulate data#####
# SIGMA MATRIX
S <- diag(10)
mu<-rep(0,10)
S[upper.tri(S, diag=TRUE)] 
for(i in 1:10){
  for(j in 1:10){
    if(S[i,j]==0){
      S[i,j]=0.2
    }
  }
}

set.seed(123)
N=10000
a<-rep(1,N)
X<-mvrnorm(n=N,mu,S)
e<-rnorm(n=N,0,1)
b0<-0.4623
b1<--0.8638
b2<-1.479
b3<--0.5139
b4<-0.4335
b5<-1.7971
b6<--0.0874
b7<--0.3138
b8<--0.8459
b9<-1.4213
b10<-1.6152
B<-matrix(c(b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10),11,1)
X<-cbind(a,X)
Y<-X%*%B+e
mydata<-cbind(Y,X[,-1])
write.csv(x=mydata,file = "mydata2.csv",row.names = F)


####### input data########
sim.regress.data <- read.regress.data.ff(filename ='mydata2.csv', predictor.cols = c(2:11), 
                                  response.col =1, first.rows = 10, next.rows = 10)
xtx<-sim.regress.data$xtx
yty<-sim.regress.data$yty
xty<-sim.regress.data$xty

##### prior #############

beta.prior.2 <- list(type="mvnorm.known", mean.mu =
                       rep(0.0, dim(xtx)[1]), cov.C = diag(1.0, dim(xtx)[1]))
sigmasq.prior.1 <- list(type = "inverse.gamma",
                        inverse.gamma.a = 1, inverse.gamma.b = 1, sigmasq.init = 1)

#### regression ########
sim.beta.sigmasq.out <- bayes.regress(data.values =
                                        sim.regress.data, beta.prior = beta.prior.2, sigmasq.prior
                                      = sigmasq.prior.1, Tsamp.out = 10000, zero.intercept =
                                        FALSE)

#### furtehr check on b1 and sigma #########
plot(mcmc(sim.beta.sigmasq.out$beta[,2])) #no burn-in
plot(mcmc(sim.beta.sigmasq.out$beta[500:10000,2])) # burnin first 500
summary(mcmc(sim.beta.sigmasq.out$beta[,2]))
summary(mcmc(sim.beta.sigmasq.out$beta[500:10000,2]))

plot(mcmc(sim.beta.sigmasq.out$sigmasq))
summary(mcmc(sim.beta.sigmasq.out$sigmasq))