#################################
#
##################################

# cancel all variables from the memory
rm(list=ls())

# set working directory (change this)
setwd("C:/Users/amele1/Dropbox/BayesSocNet/programs/june2012/new_routines_out2012/new_routines_out2012/recip_example_n200/")


# import data (need to change the file name)
x<-read.table("parbayes_mh_recip.txt")
y<-read.table("parbayes_exch_recip_R1000.txt")
z<-read.table("parbayes_exch_recip_R5000.txt")
w<-read.table("parbayes_exch_recip_R10000.txt")
r<-read.table("parbayes_exch_recip_R50000.txt")
s<-read.table("parbayes_exch_recip_R100000.txt")
t<-read.table("parbayes_exch_recip_R1mil.txt")
q<-read.table("parbayes_exch_recip_R10mil.txt")

# plot time series of mh vs exchange
par(mfrow=c(3,2))
plot(x[,1],col="blue", type="l", ylim=c(-1.5,-2.5))
plot(y[,1],col="red",type="l", ylim=c(-1.5,-2.5))
plot(z[,1],col="green",type="l", ylim=c(-1.5,-2.5))
plot(w[,1],col="purple",type="l", ylim=c(-1.5,-2.5))
plot(r[,1],col="black",type="l", ylim=c(-1.5,-2.5))
plot(s[,1],col="cyan",type="l", ylim=c(-1.5,-2.5))

par(mfrow=c(3,2))
plot(x[,2],col="blue", type="l", ylim=c(-1,1))
plot(y[,2],col="red",type="l", ylim=c(-1,1))
plot(z[,2],col="green",type="l", ylim=c(-1,1))
plot(w[,2],col="purple",type="l", ylim=c(-1,1))
plot(r[,2],col="black",type="l", ylim=c(-1,1))
plot(s[,2],col="cyan",type="l", ylim=c(-1,1))






# plot contour maps
library(KernSmooth)
xkern<-bkde2D(x,c(.1,.01))
ykern<-bkde2D(y,c(.1,.01))
zkern<-bkde2D(z,c(.1,.01))
wkern<-bkde2D(w,c(.1,.01))
rkern<-bkde2D(r,c(.1,.01))
skern<-bkde2D(s,c(.1,.01))
tkern<-bkde2D(t,c(.1,.01))
qkern<-bkde2D(q,c(.1,.01))



contour(xkern$x1,xkern$x2,xkern$fhat, lty=1, lwd=2, col="black", xlim=c(-2.2,-1.6),ylim=c(0.0,0.6))
contour(zkern$x1,zkern$x2,zkern$fhat, add=TRUE, lty=2, lwd=2, col="red")
contour(skern$x1,skern$x2,skern$fhat, add=TRUE, lty=3, lwd=2, col="blue")


contour(xkern$x1,xkern$x2,xkern$fhat, lty=1, lwd=2, col="black", xlim=c(-2.2,-1.6),ylim=c(0.0,0.6), levels=seq(0,14,2))
contour(zkern$x1,zkern$x2,zkern$fhat, add=TRUE, lty=2, lwd=2, col="red", levels=seq(0,14,2))
contour(skern$x1,skern$x2,skern$fhat, add=TRUE, lty=3, lwd=2, col="blue", levels=seq(0,14,2))



contour(xkern$x1,xkern$x2,xkern$fhat)
contour(zkern$x1,zkern$x2,zkern$fhat, add=TRUE, col="red")
contour(skern$x1,skern$x2,skern$fhat, add=TRUE, lty=2, col="green")

points(mean(z[,1]),mean(z[,2]), cex=1, pch=17, col="red")
points(mean(x[,1]),mean(x[,2]), cex=1, pch=19, col ="black")
points(mean(s[,1]),mean(s[,2]), cex=1, pch=15, col ="purple")




contour(xkern$x1,xkern$x2,xkern$fhat, lty=1, lwd=1, col="black", xlim=c(-2.05,-1.75),ylim=c(0.1,0.45), levels=seq(0,14,1))
contour(zkern$x1,zkern$x2,zkern$fhat, add=TRUE, lty=2, lwd=1, col="red", levels=seq(0,14,1))
contour(skern$x1,skern$x2,skern$fhat, add=TRUE, lty=3, lwd=1, col="blue", levels=seq(0,14,1))
points(mean(z[,1]),mean(z[,2]), cex=1, pch=17, col="red")
points(mean(x[,1]),mean(x[,2]), cex=1, pch=19, col ="black")
points(mean(s[,1]),mean(s[,2]), cex=1, pch=15, col ="purple")




jpeg("prova.jpeg")
contour(xkern$x1,xkern$x2,xkern$fhat, lty=1, lwd=2, col="black", xlim=c(-2.05,-1.75),ylim=c(0.1,0.45), levels=seq(0,14,2))
contour(zkern$x1,zkern$x2,zkern$fhat, add=TRUE, lty=2, lwd=2, col="red", levels=seq(0,14,2))
points(mean(z[,1]),mean(z[,2]), cex=1, pch=17, col="red")
points(mean(x[,1]),mean(x[,2]), cex=1, pch=19, col ="black")
dev.off()

pdf("prova.pdf")
contour(xkern$x1,xkern$x2,xkern$fhat, lty=1, lwd=2, col="black", xlim=c(-2.05,-1.75),ylim=c(0.1,0.45), levels=seq(0,14,2))
contour(zkern$x1,zkern$x2,zkern$fhat, add=TRUE, lty=2, lwd=2, col="red", levels=seq(0,14,2))
points(mean(z[,1]),mean(z[,2]), cex=1, pch=17, col="red")
points(mean(x[,1]),mean(x[,2]), cex=1, pch=19, col ="black")
dev.off()

pdf("prova.pdf")
contour(xkern$x1,xkern$x2,xkern$fhat, lty=1, lwd=2, col="black", xlim=c(-2.05,-1.75),ylim=c(0.1,0.45), levels=seq(0,14,2))
contour(zkern$x1,zkern$x2,zkern$fhat, add=TRUE, lty=2, lwd=2, col="red", levels=seq(0,14,2))
points(mean(z[,1]),mean(z[,2]), cex=1, pch=24, lwd=2, col="red")
points(mean(x[,1]),mean(x[,2]), cex=1, pch=25, lwd=2, col ="black")
dev.off()




pdf("prova.pdf")
contour(xkern$x1,xkern$x2,xkern$fhat, lty=1, lwd=2, col="black", xlim=c(-2.05,-1.75),ylim=c(0.1,0.45), levels=seq(0,14,3))
contour(skern$x1,skern$x2,skern$fhat, add=TRUE, lty=2, lwd=2, col="red", levels=seq(0,14,2))
points(mean(s[,1]),mean(s[,2]), cex=.7, pch=1, lwd=2, col="red")
points(mean(x[,1]),mean(x[,2]), cex=.7, pch=2, lwd=2, col ="black")
legend("topright", c("R = 100000", "exact"), text.col=c("red","black"), lty=c(2,1), lwd=1, col=c("red","black") , cex=0.75)

dev.off()



contour(xkern$x1,xkern$x2,xkern$fhat, lwd=2)
contour(skern$x1,skern$x2,skern$fhat, add=TRUE, lty=2, lwd=2, col="red")



# plot autocorrelations of mh vs exchange
par(mfrow=c(4,2))
acf(x[,1])
acf(y[,1])
acf(z[,1])
acf(w[,1])
acf(r[,1])
acf(s[,1])
acf(t[,1])
acf(q[,1])



par(mfrow=c(3,2))
acf(x[,2])
acf(y[,2])
acf(z[,2])
acf(w[,2])
acf(r[,2])
acf(s[,2])



# plot histograms of mh vs exchange
par(mfrow=c(3,2))
brk<-seq(-1.5,-2.5,-0.005)
hist(x[,1], breaks=brk, col="blue", xlim=c(-1.5,-2.5), ylim=c(0,800))
hist(y[,1], breaks=brk, col="red", xlim=c(-1.5,-2.5), ylim=c(0,800))
hist(z[,1], breaks=brk, col="green", xlim=c(-1.5,-2.5), ylim=c(0,800))
hist(w[,1], breaks=brk, col="purple", xlim=c(-1.5,-2.5), ylim=c(0,800))
hist(r[,1], breaks=brk, col="black", xlim=c(-1.5,-2.5), ylim=c(0,800))
hist(s[,1], breaks=brk, col="cyan", xlim=c(-1.5,-2.5), ylim=c(0,800))

par(mfrow=c(3,2))
brk<-seq(-1,1,0.005)
hist(x[,2], breaks=brk, col="blue", xlim=c(-1,1), ylim=c(0,800))
hist(y[,2], breaks=brk, col="red", xlim=c(-1,1), ylim=c(0,800))
hist(z[,2], breaks=brk, col="green", xlim=c(-1,1), ylim=c(0,800))
hist(w[,2], breaks=brk, col="purple", xlim=c(-1,1), ylim=c(0,800))
hist(r[,2], breaks=brk, col="black", xlim=c(-1,1), ylim=c(0,800))
hist(s[,2], breaks=brk, col="cyan", xlim=c(-1,1), ylim=c(0,800))



# kernel density
par(mfrow=c(1,1))
plot(density(x[,1],kernel="biweight"),col="blue", xlim=c(-1.5,-2.5))
lines(density(y[,1],kernel="biweight"),col="red")
lines(density(z[,1],kernel="biweight"),col="green")
lines(density(w[,1],kernel="biweight"),col="purple")
lines(density(r[,1],kernel="biweight"),col="black")
lines(density(s[,1],kernel="biweight"),col="cyan")
lines(density(t[,1],kernel="biweight"),col="yellow")
lines(density(q[,1],kernel="biweight"),col="darkgreen")


# kernel density
par(mfrow=c(1,1))
plot(density(x[,2],kernel="biweight"),col="blue", xlim=c(-1,1))
lines(density(y[,2],kernel="biweight"),col="red")
lines(density(z[,2],kernel="biweight"),col="green")
lines(density(w[,2],kernel="biweight"),col="purple")
lines(density(r[,2],kernel="biweight"),col="black")
lines(density(s[,2],kernel="biweight"),col="cyan")
lines(density(t[,2],kernel="biweight"),col="yellow")
lines(density(q[,2],kernel="biweight"),col="darkgreen")




# kolmogorov smirnov
ks.test(x[,1],y[,1])
ks.test(x[,1],z[,1])
ks.test(x[,1],w[,1])
ks.test(x[,1],r[,1])
ks.test(x[,1],s[,1])
ks.test(x[,1],t[,1])
ks.test(x[,1],q[,1])




# kolmogorov smirnov
ks.test(x[,2],y[,2])
ks.test(x[,2],z[,2])
ks.test(x[,2],w[,2])
ks.test(x[,2],r[,2])
ks.test(x[,2],s[,2])
ks.test(x[,2],t[,2])
ks.test(x[,2],q[,2])




# mean
mean(x[,1])
mean(y[,1])
mean(z[,1])
mean(w[,1])
mean(r[,1])
mean(s[,1])
mean(t[,1])
mean(q[,1])



# sd
sd(x[,1])
sd(y[,1])
sd(z[,1])
sd(w[,1])
sd(r[,1])
sd(s[,1])
sd(t[,1])
sd(q[,1])





# mean
mean(x[,2])
mean(y[,2])
mean(z[,2])
mean(w[,2])
mean(r[,2])
mean(s[,2])

# sd
sd(x[,2])
sd(y[,2])
sd(z[,2])
sd(w[,2])
sd(r[,2])
sd(s[,2])

