get_log_ML_bioassay {perms}

get_log_ML_bioassay

Description

Computes the log marginal likelihood of bioassay data from the log permanents. Given the computed log permanents log_perms, this function computes the log marginal likelihood using the formula (2.3) in [1]. It takes care of repeated trials by adding the appropriate log binomial factor.

Usage

get_log_ML_bioassay(log_perms, successes, trials, n, debug = FALSE)

Arguments

log_perms

A vector length n containing the computed log permanents, where a zero permanent is indicated by a NA value.

successes

A vector of length n contatining the number of successful trials at each level.

trials

A vector of length n containing the number of trials at each level.

n

Sample size.

debug

If TRUE, debug information is printed.

Value

The estimated log marginal likelihood. A NA value is returned if there are no non-zero numbers.

References

[1] Christensen, D (2023). Inference for Bayesian nonparametric models with binary response data via permutation counting. Bayesian Analysis, Advance online publication, DOI: 10.1214/22-BA1353.

Examples

## Dirichlet toy model
library(perms)
set.seed(1996)
n = 500
num_trials = 10
levels = seq(-1, 1, length.out = num_trials)

trials = rep(n %/% num_trials, num_trials)
successes = c(10, 26, 10, 20, 20, 19, 29, 24, 31, 33)

S = 200
alpha = 1.0

get_X = function(S,n,alpha,seed){
  set.seed(seed)
  X = matrix(0, nrow = S, ncol = n)
  for (s in 1:S) {
    X[s,1] = rnorm(1)
    for (i in 2:n) {
      u = runif(1)
      if(u < (alpha/(alpha+i-1))){
        X[s,i] = rnorm(1)
      }else{
        if(i==2){
          X[s,i] = X[s,1]
        }else{
          X[s,i] = sample(X[s, 1:(i-1)],size=1)
        }
      }
      
    }
    
  }
  return(X)
}

seed = 1996
X = get_X(S, n, alpha, seed)
log_perms = get_log_perms_bioassay(X, levels, successes, trials,
           debug=FALSE,parallel = FALSE)
log_ml = get_log_ML_bioassay(log_perms, successes, trials, n)

proportion = sum(!is.na(log_perms)) / S*100

proportion 
log_ml