Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Bayes factor for one-way table
# Prior setup
prior_concentration = 10
# null
p0 = c(0.3, 0.3, 0.4)
# fake data, null
N = 100
y = rmultinom(n = 1,prob = p0, size = N )
# Visualize prior density; vis really only works for 3 probabilities
x = MCMCpack::rdirichlet(10000, prior_concentration * p0)
plot(x[,1], x[,2], ylim = c(0,1), xlim = c(0,1), pch = 19, col = rgb(0,0,1,.1))
abline(h = p0[1], v = p0[2])
# (log) prob of data under null
pr_y_h0 = dmultinom(x = y, prob = p0, log = TRUE)
# Estimate log prob of data under null with Monte Carlo
M = 100000
p1s = MCMCpack::rdirichlet(M, prior_concentration * p0)
tmp_pr_h1 = sapply(1:M, function(i) dmultinom(x = y, prob = p1s[i,], log = TRUE))
pr_y_h1 = BayesFactor::logMeanExpLogs(tmp_pr_h1)
bf_10 = exp(pr_y_h1 - pr_y_h0)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.