# BF for comparing H0 to H1 variance known (see page 128 Lee)
BF.z = function(n, z, lambda=1) {
  BF = sqrt((n + lambda)/lambda)*exp(-.5*(z^2*n/(n+lambda)))
  prob.0= 1/(1 + 1/BF)
  return(list(BF=BF, prob.0=prob.0))
  
}

# BF for comparing H0 to H1 variance unknown using the reference prior on phi
BF.t = function(n, t, lambda=1) {
  d = n-1
  BF = sqrt((n + lambda)/lambda)*((t^2*lambda/(n+lambda) + d)/(t^2 + d))^((d+1)/2)
  prob.0 = 1/(1/BF + 1)
return(list(BF=BF, prob.0=prob.0))
}

BF.c = function(n, t, var=F) {
  if (!var) {
  BF = 1/integrate(f=function(x)
            { dgamma(x, 1/2, rate=1/2)/BF.t(n=n, t=t, lambda=x)[[1]]},
            0, Inf)$value}
  else { BF = 1/integrate(f=function(x)
            { dgamma(x, 1/2, rate=1/2)/BF.z(n=n, z=t, lambda=x)[[1]]},
            0, Inf)$value}
prob.0=1/(1/BF + 1)
return(list(BF=BF, prob.0=prob.0))
}

BF.lb= function(n, t) {
  p = 2*pt(-abs(t), n-1)
  BF = -exp(1)*p*log(p)
  prob.0 = 1/(1 + 1/BF )
  return(list(BF=BF, prob.0=prob.0, pvalue=p))
}


# data from Ch 5 p 140
diff = c(1.2, 2.4, 1.3, 1.3, 0., 1,1.9, .9, 4.6, 1.4)
t.test(diff)
#Lindley: reject H0 if mu0 is outside HPD region

BF.t(10, 4.128)
BF.c(10, 4.128)
BF.lb(10, 4.128)