#Computing the Bayes factor as on pages 152-153 of the slides
n=100
Bnumerator=((.088-.1078)^2+.00875)^(-n+.5)

Bh=function(xi){((2*xi-.088-.1078)^2+.00875)^(-n+.5)}
Brnorm=rnorm(10^5)
Bdenominator=mean(Bh(Brnorm))
plot(1:10^5,cumsum(Bh(Brnorm))/(1:10^5),type="l")
Brnorm=rnorm(10^5)
# impact of monte carlo variation
lines(1:10^5,cumsum(Bh(Brnorm))/(1:10^5),col="tomato")
Brnorm=rnorm(10^5)
lines(1:10^5,cumsum(Bh(Brnorm))/(1:10^5),col="green")
Bnumerator/Bdenominator

# checking for proper simulation
hist(Brnorm,nclass=500,prob=T)
curve(dnorm,add=T,col="sienna",lwd=5)

# estimate of the variance of the MC estimator:
var(Bh(Brnorm))/10^5
# we get infinity - so correct because numbers are too big
# thus change the scales
exp(log(var(Bh(Brnorm)/10^(202))/10^5)+404*log(10)/2)
# result 8.050638e+199

mean(Bh(Brnorm))
#[1] 1.937489e+202

#Given these figures, I can trust the estimate up to the 3rd digit (202-199)

#Now redo by importance sampling using as importance density a T with 3 df:
Brt=rt(10^4,df=3)

Bnewh=function(xi){Bh(xi)*dnorm(xi)/dt(xi,df=3)}

Bdenominator=mean(Bnewh(Brt))
plot(1:10^4,cumsum(Bnewh(Brt))/(1:10^4),type="l")

Brt=rt(10^4,df=3)
lines(1:10^4,cumsum(Bnewh(Brt))/(1:10^4),type="l",add=T,col="tomato")

Brt=rt(10^4,df=3)
lines(1:10^4,cumsum(Bnewh(Brt))/(1:10^4),type="l",add=T,col="green")