rm(list=ls())
set.seed(123812)

library(snow)
library(MCMCpack)
cl<-makeCluster(10, type="SOCK")
clusterEvalQ(cl, library(MCMCpack))

test.sim<-function(b, n, reps, B0, c0=0.001, d0=0.001){
  
  sig.true<-c()
  sig.false<-c()
  
  for(i in 1:reps){
    x<-runif(n)
    y<-b*x+rnorm(n, mean=0, sd=1)
    bayes.reg<-summary(MCMCregress(y~x, b0=c(0,0), B0=B0, c0=c0, d0=d0))
    sig.true[i]<-ifelse(sign(bayes.reg$quantiles[2,1])==sign(bayes.reg$quantiles[2,5]), 1, 0)
    
    x<-runif(n)
    y<-0*x+rnorm(n, mean=0, sd=1)
    bayes.reg.2<-summary(MCMCregress(y~x, b0=c(0,0), B0=B0, c0=c0, d0=d0))
    sig.false[i]<-ifelse(sign(bayes.reg.2$quantiles[2,1])==sign(bayes.reg.2$quantiles[2,5]), 1, 0)
  }
  
  return(cbind(sig.true, sig.false))
}

clusterExport(cl, "test.sim")
clusterSetupRNG(cl, seed=328979)

x<-runif(100)
y<-1*x+rnorm(100, mean=0, sd=1)

b.est<-summary(lm(y~x))$coefficients[,1]

result<-clusterCall(cl, test.sim, b=b.est, n=100, reps=200, B0=0)
result.mat<-do.call("rbind", result)

pr.sig.true<-sum(result.mat[,1])/2000
pr.sig.false<-sum(result.mat[,2])/2000

pr.false<-seq(from=0.01, to=0.99, by=0.01)
update.line<-c()
for(k in 1:length(pr.false)){
  update.line[k]<-((pr.sig.false*pr.false[k])/(pr.sig.true*(1-pr.false[k])+pr.sig.false*pr.false[k]))
}

par(mar=c(5,5,4,2))
plot(update.line~pr.false, type="l", ylim=c(0,1), xlim=c(0,1), lwd=2, main=c("Discrimination power of 95% CIs that exclude zero"), xlab=expression(paste("Overall Pr(", beta, " = 0)")), ylab=expression(paste("Pr(", beta, "= 0 | CI excludes zero)")))
abline(0,1,lty=2)
dev.copy2eps(file="information_b1_100.eps")



clusterSetupRNG(cl, seed=328979)
sigma.test<-c(seq(0.1, 1, 0.1), seq(1.5, 20, 0.5))
update.line<-c()

pb<-txtProgressBar(min=0, max=length(sigma.test), style=3)
for(k in 1:length(sigma.test)){
  
  setTxtProgressBar(pb, value=k)
  result<-clusterCall(cl, test.sim, b=b.est, n=100, reps=200, B0=sigma.test[k])
  result.mat<-do.call("rbind", result)
  
  pr.sig.true<-sum(result.mat[,1])/2000
  pr.sig.false<-sum(result.mat[,2])/2000
  
  update.line[k]<-(pr.sig.false)/(pr.sig.true+pr.sig.false)
}
close(pb)

par(mar=c(5,5,4,2))
plot(update.line~sigma.test, lwd=2, xlim=c(0, 20), main=c("Discrimination power of 95% CIs that exclude zero"), xlab=expression(paste("Std. deviation of prior (", sigma, ")")), ylab=expression(paste("Pr( ", beta, " = 0 | 95% CI excludes 0)")))

lines(predict(loess(update.line~sigma.test), deg=1)~sigma.test)

dev.copy2eps(file="information_prior_b1_100.eps")



stopCluster(cl)
cl<-makeCluster(2, type="SOCK")
clusterEvalQ(cl, library(MCMCpack))

set.seed(324985)

x<-runif(1000)
y<-1*x+rnorm(1000, mean=0, sd=1)

b.est<-summary(lm(y~x))$coefficients[,1]

clusterSetupRNG(cl, seed=213981)
result<-clusterCall(cl, test.sim, b=b.est, n=1000, reps=1000, B0=0)
result.mat<-do.call("rbind", result)

pr.sig.true<-sum(result.mat[,1])/2000
pr.sig.false<-sum(result.mat[,2])/2000

pr.false<-seq(from=0.01, to=0.99, by=0.01)
update.line<-c()
for(k in 1:length(pr.false)){
  update.line[k]<-((pr.sig.false*pr.false[k])/(pr.sig.true*(1-pr.false[k])+pr.sig.false*pr.false[k]))
}

par(mar=c(5,5,4,2))
plot(update.line~pr.false, type="l", ylim=c(0,1), xlim=c(0,1), lwd=2, main=c("Posterior Beliefs", "based on 95% CIs that exclude zero"), xlab=expression(paste("Prior Pr(", beta, " = 0)")), ylab=expression(paste("Updated Pr(", beta, "= 0 | results)")))
abline(0,1,lty=2)
dev.copy2eps(file="information_b1_1000.eps")






set.seed(2342)

x<-runif(100)
y<-0*x+rnorm(100, mean=0, sd=1)

b.est<-summary(lm(y~x))$coefficients[,1]

clusterSetupRNG(cl, seed=432830)
result<-clusterCall(cl, test.sim, b=b.est, n=100, reps=1000, B0=0)
result.mat<-do.call("rbind", result)

pr.sig.true<-sum(result.mat[,1])/2000
pr.sig.false<-sum(result.mat[,2])/2000

pr.false<-seq(from=0.01, to=0.99, by=0.01)
update.line<-c()
for(k in 1:length(pr.false)){
  update.line[k]<-((pr.sig.false*pr.false[k])/(pr.sig.true*(1-pr.false[k])+pr.sig.false*pr.false[k]))
}

par(mar=c(5,5,4,2))
plot(update.line~pr.false, type="l", ylim=c(0,1), xlim=c(0,1), lwd=2, main=c("Posterior Beliefs", "based on 95% CIs that exclude zero"), xlab=expression(paste("Prior Pr(", beta, " = 0)")), ylab=expression(paste("Updated Pr(", beta, "= 0 | results)")))
abline(0,1,lty=2)
dev.copy2eps(file="information_b0_100.eps")