Skip to content

Instantly share code, notes, and snippets.

@EtzAlex
Last active June 19, 2016 04:15
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save EtzAlex/5f28d0e0917dab8b00304a2c3fe02b15 to your computer and use it in GitHub Desktop.
Save EtzAlex/5f28d0e0917dab8b00304a2c3fe02b15 to your computer and use it in GitHub Desktop.
This is the code for Alexander Etz's blog at http://wp.me/p4sgtg-SQ
#This is the code for Alexander Etz's blog at http://wp.me/p4sgtg-SQ
#Code to make the figure and find max oracle BF for normal data
maxL <- function(mean=1.96,se=1,h0=0){
L1 <-dnorm(mean,mean,se)
L2 <-dnorm(mean,h0,se)
Ratio <- L1/L2
curve(dnorm(x,mean,se), xlim = c(-2*mean,2.5*mean), ylab = "Likelihood", xlab = "Population mean", las=1,
main = "Likelihood function for the mean", lwd = 3)
points(mean, L1, cex = 2, pch = 21, bg = "cyan")
points(h0, L2, cex = 2, pch = 21, bg = "cyan")
lines(c(mean, h0), c(L1, L1), lwd = 3, lty = 2, col = "cyan")
lines(c(h0, h0), c(L1, L2), lwd = 3, lty = 2, col = "cyan")
return(Ratio) ## Returns the Bayes factor for oracle hypothesis vs null
}
#Code to find the max subtle oracle prior BF for normal data
t = 1.96 #Critical value corresponding to the p value
maxBF = 1/(sqrt(exp(1))*t*exp(-t^2/2)) #Formula from Berger and Sellke 1987
#multiply by 2*(1-p/2) if using a one-sided prior
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment