#############################################################################################################
#
#  Exercise 2.26 (page 79) – Bayes factor sensitivity to vague/improper priors.
#
##############################################################################################################
#
# Author : Hedibert Freitas Lopes
#          Graduate School of Business
#          University of Chicago
#          5807 South Woodlawn Avenue
#          Chicago, Illinois, 60637
#          Email : hlopes@ChicagoGSB.edu
#
###############################################################################################################
set.seed(129096)

n     = 10
mu    = 0
x     = rnorm(n,mu,0.1)
mu0   = 0
th0   = 0
tau2s = seq(0.01,5,by=0.01)
sig2  = 1
bf = NULL
for (tau2 in tau2s)
   bf = c(bf,sum(dnorm(x,mu0,sqrt(1+tau2),log=T))-sum(dnorm(x,th0,sqrt(1+sig2),log=T)))
plot(tau2s,exp(bf),type="l",ylab="",xlab="tau2/sig2",main="Bayes Factor \n sig2=1")
segments(1,0,1,1,col=2)
segments(0,1,1,1,col=2)

# Sample of size n=10 from the standard normal
set.seed(129096)
n     = 10
x     = c(-0.411,-0.778,0.437,0.723,1.670,0.404,0.444,0.280,-0.345,1.384)
tau2  = 1.0
alpha = 5
sig20 = 0.8
sig20 = 0.5
stdev = sqrt(sig20*alpha/(alpha-2))
tau2s = seq(0.01,5,by=0.01)
bf = NULL
for (tau2 in tau2s)
  bf = c(bf,sum(dnorm(x,0,sqrt(1+tau2),log=T))-sum(dt(x/stdev,alpha,log=T))-n*log(stdev))
plot(tau2s/sig2,exp(bf),type="l",ylab="",xlab="tau2/sig2",main="Bayes Factor \n sig2=0.8")
exp(bf)

tau2s = c(0.0216,1.208,2.806)
bf = NULL
for (tau2 in tau2s)
  bf = c(bf,sum(dnorm(x,0,sqrt(1+tau2),log=T))-sum(dt(x/stdev,alpha,log=T))-n*log(stdev))
cbind(tau2s,exp(bf))

# Sample of size n=10 from the standard normal
n     = 10
x     = c(-0.411,-0.778,0.437,0.723,1.670,0.404,0.444,0.280,-0.345,1.384)
tau2  = 1.0
alpha = 5
bf = NULL
for (i in 1:100){
  tau2  = i/100*0.2
  sig20 = i/100*0.2
  stdev = sqrt(sig20*alpha/(alpha-2))
  bf = c(bf,sum(dnorm(x,0,sqrt(1+tau2),log=T))-sum(dt(x/stdev,alpha,log=T))-n*log(stdev))
}
bf