bnormal=function(N,rho){
rho2=sqrt(1-rho*rho)
k=1	
theta1=0
theta2=0
while(k<=N){
theta1=rnorm(1,rho*theta2,rho2)
theta2=rnorm(1,rho*theta1,rho2)
write.table(matrix(c(k,theta1,theta2),ncol=3),"bivariatenormal.txt",sep="\t",append=T,row.names=F,col.names=F)
k=k+1
}
}

bnormal(1000,0.5)  

gen.theta1=read.table("bivariatenormal.txt")
dim(gen.theta1)
head(gen.theta1)

par(mfrow=c(2,1))
plot(gen.theta1[,1],gen.theta1[,2],xlab="Iteration",ylab="Theta1")
plot(gen.theta1[,1],gen.theta1[,3],xlab="Iteration",ylab="Theta2")

par(mfrow=c(2,1))
plot(gen.theta1[,1],gen.theta1[,2],xlab="Iteration",ylab="Theta1",type='n')
lines(gen.theta1[,1],gen.theta1[,2])
plot(gen.theta1[,1],gen.theta1[,3],xlab="Iteration",ylab="Theta2",type='n')
lines(gen.theta1[,1],gen.theta1[,3])

var(gen.theta1[,2:3])   # expect variances to be close to 1, covariances close to rho

colMeans(gen.theta1[,2:3])   # expect means to be close to 0


par(mfrow=c(2,1))
hist(gen.theta1[,2])
hist(gen.theta1[,3])

bnormal(10000,0.5)  
bnormal(20000,0.5)  




bnormal2=function(N,rho){
rho2=sqrt(1-rho*rho)
k=1	
theta1=100
theta2=100
while(k<=N){
theta1=rnorm(1,rho*theta2,rho2)
theta2=rnorm(1,rho*theta1,rho2)
write.table(matrix(c(k,theta1,theta2),ncol=3),"bnrho05.txt",sep="\t",append=T,row.names=F,col.names=F)
k=k+1
}
}

bnormal2(10000,0.5)


gen2.theta=read.table("bnrho05.txt")
dim(gen2.theta)

par(mfrow=c(2,1))
plot(gen2.theta[,1],gen2.theta[,2],xlab="Iteration",ylab="Theta1",type='n')
lines(gen2.theta[,1],gen2.theta[,2])
plot(gen2.theta[,1],gen2.theta[,3],xlab="Iteration",ylab="Theta2",type='n')
lines(gen2.theta[,1],gen2.theta[,3])

par(mfrow=c(2,1))
plot(gen2.theta[1:20,1],gen2.theta[1:20,2],xlab="Iteration",ylab="Theta1",type='n')
lines(gen2.theta[1:20,1],gen2.theta[1:20,2])
plot(gen2.theta[1:20,1],gen2.theta[1:20,3],xlab="Iteration",ylab="Theta2",type='n')
lines(gen2.theta[1:20,1],gen2.theta[1:20,3])


par(mfrow=c(2,1))
plot(gen2.theta[-c(1:20),1],gen2.theta[-c(1:20),2],xlab="Iteration",ylab="Theta1",type='n')
lines(gen2.theta[-c(1:20),1],gen2.theta[-c(1:20),2])
plot(gen2.theta[-c(1:20),1],gen2.theta[-c(1:20),3],xlab="Iteration",ylab="Theta2",type='n')
lines(gen2.theta[-c(1:20),1],gen2.theta[-c(1:20),3])


# BOA

install.packages('boa')
library(boa) 
boa.menu()

### don't run ####
boa.menu(recover=T) # if you have troubleshooting
#### end of don't run ####


Options
Import Files 
Working Directory

File
Import Data
Flat ASCII file
bnrho05

Plot
Descriptive
Trace
Autocorrelations

Analysis
Descriptive Statistics
Summary Statistics
Autocorrelations
Correlation Matrix
Convergence Diagnostics
Geweke  # there is evidence against convergence when p-value is less than 0.05 or greater than 0.95.
# Heidelberger and welch if the stationary test fails, chain needs to be run longer for     
#  convergence purposes. If the halfwidth test fails, chain might    
# be run longer to increase the accuracy in estimating posterior 
# estimate 
Raftery and lewis	# dependence factor greater than 5 implies convergence problem

Data Management
Chains
Subset
#Press enter until you see "Specify iterations[all]", then
100:10000

File
Exit BOA