Analysis of the Gamma prior parameter in the Hierarchical Dawid-Skene model

hierarchical-ds-model-prior-analysis.R

# From rater
softmax <- function(x) {
  exp(x - logsumexp(x))
}

logsumexp <- function(x) {
  y <- max(x)
  y + log(sum(exp(x - y)))
}

sample_theta <- function(K, on_diag) {
  mu <- matrix(0, ncol = K, nrow = K)
  diag(mu) <- on_diag

  gamma <- matrix(ncol = K, nrow = K)
  for (i in seq_len(K)) {
    for (j in seq_len(K)) {
      gamma[i, j] <- rnorm(1, mean = mu[i, j], 1)
    }
  }

  theta <- matrix(nrow = K, ncol = K)
  for (i in seq_len(K)) {
    theta[i, ] <- softmax(gamma[i, ])
  }
  theta
}

reps <- replicate(10000, sample_theta(K = 2, on_diag = 0.5), simplify = FALSE)

theta_1_1 <- sapply(reps, function(x) x[1, 1])
theta_1_2 <- sapply(reps, function(x) x[1, 2])
theta_2_1 <- sapply(reps, function(x) x[2, 1])
theta_2_2 <- sapply(reps, function(x) x[2, 2])

hist(theta_1_1)
hist(theta_1_2)
hist(theta_2_1)
hist(theta_2_2)