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

n<-100
b<-seq(from=0.05, to=1, by=0.05)
pr.true<-c(0.1, 0.2, 0.5)

update<-matrix(data=NA, ncol=length(pr.true), nrow=length(b))
for(j in 1:length(b)){
  
  sig.true<-c()
  sig.false<-c()
  for(i in 1:1000){
    
    x<-runif(n)
    y<-b[j]*x+rnorm(n, mean=0, sd=1)
    sig.true[i]<-ifelse(summary(lm(y~x))$coefficients[2,4]<0.05, 1, 0)
    
    x<-runif(n)
    y<-0*x+rnorm(n, mean=0, sd=1)
    lm(y~x)
    sig.false[i]<-ifelse(summary(lm(y~x))$coefficients[2,4]<0.05, 1, 0)
    
  }
  
  pr.sig.true<-sum(sig.true)/1000
  pr.sig.false<-sum(sig.false)/1000
  
  for(k in 1:length(pr.true)){
    update[j,k]<-(1-((pr.sig.true*pr.true[k])/(pr.sig.true*pr.true[k]+pr.sig.false*(1-pr.true[k]))))
  }

}
  
plot(update[,1]~b, ylim=c(0,1), xlim=c(0,1), xlab=expression(paste("true ", beta, " coefficient")), ylab=expression(paste("Pr(", beta," = 0 | ", hat(beta), " is statistically significant)")), main=expression(paste("Updated Probability that ", beta," = 0")))
lines(predict(loess(update[,1]~b))~b)
points(update[,2]~b, col="red")
lines(predict(loess(update[,2]~b))~b, col="red")
points(update[,3]~b, col="blue")
lines(predict(loess(update[,3]~b))~b, col="blue")
legend("topright", y.intersp=1.4, col=c("black", "red", "blue"), lty=c(1,1,1), legend=c(expression(paste("prior Pr(", beta, " = 0) = 0.9   ")), expression(paste("prior Pr(", beta, " = 0) = 0.8   ")), expression(paste("prior Pr(", beta, " = 0) = 0.5   "))))