good.multinomial.test.R

Jim — Nov 9, 2013, 3:45 PM

# Good testing

# test of equiprobability of a multinomial vector
# reference, Good (1967), JRSS

# observe vector of frequencies
# y = (#1, #2, ..., #K)
# compute the log bayes factor in support of equiprobability

# inputs to function are vector of frequencies y
# and Dirichlet smoothing parameter a

log.bf <- function(a, y){
  ldirich <- function(a) {
    val <- sum(lgamma(a)) - lgamma(sum(a))
    return(val)
  }
  K <- length(y)
  N <- sum(y)
  ldirich(y + a) - ldirich(rep(a, K)) + N * log(K)
}

# roll "fair?" dice 200 times

# rolls <- sample(6, size=200, replace=TRUE)
# y <- tabulate(rolls, 6)

y <- c(43, 24, 37, 24, 31, 41)

# compare with chi-squared test of equiprobability

chisq.test(y)

    Chi-squared test for given probabilities

data:  y
X-squared = 10.36, df = 5, p-value = 0.06565

# the log Bayes factor against equiprobability with a = 10

log.bf(10, y)

[1] 0.3794

# show the sensitivity of the log Bayes factor with respect
# to the smoothing parameter a

vector.a <- seq(5, 200, length=100)
plot(vector.a,
     sapply(vector.a, log.bf, y), 
     type="l",
     xlab="Smoothing Parameter a",
     ylab="Log Bayes Factor",
     main="Test of Equiprobability")