Bayes factor for one-way table

bf_one_way

# 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)