richarddmorey / BayesFactor

BayesFactor R package for Bayesian data analysis with common statistical models.
https://richarddmorey.github.io/BayesFactor/
131 stars 48 forks source link

feature request: adding support for one-way table #135

Closed IndrajeetPatil closed 3 years ago

IndrajeetPatil commented 5 years ago

Based on our Twitter discussion: https://twitter.com/patilindrajeets/status/1142763129715208193

IndrajeetPatil commented 3 years ago

@richarddmorey Any possibility that this can be included in the next release, if there is one on the horizon?

richarddmorey commented 3 years ago

I'll take a look at this; wrapping the functions in the BayesFactor objects is a nontrivial job and given everything else going on couldn't be done before Christmas though.

IndrajeetPatil commented 3 years ago

I should also post the code you had sent:

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

# Density function and random generation from the Dirichlet distribution.
rdirichlet<- function (n, alpha) 
{
    l <- length(alpha)
    x <- matrix(stats::rgamma(l * n, alpha), ncol = l, byrow = TRUE)
    sm <- x %*% rep(1, l)
    return(x/as.vector(sm))
}

# Visualize prior density; vis really only works for 3 probabilities
x = 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 = 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)
IndrajeetPatil commented 3 years ago

@richarddmorey Dear Richard, hope you are having good holidays. I just wanted to touch base again about this issue.

Any possibility that the next release will be happening any time soon and if it will contain support for one-way table?

In my package development, I am also relying on the fix implemented in https://github.com/richarddmorey/BayesFactor/issues/149, but currently need to skip all testing on CRAN of this functionality, since the CRAN-version of BayesFactor doesn't contain this fix.

IndrajeetPatil commented 3 years ago

Yikes, I thought I was closing this issue, and not the current one 🙈

Are you still planning to work on this, or do you think this is something better left off for extension packages to implement? If the former, I can re-open it.