[REFACTOR] Parametrization of the exchangeable Markov model

[hsloot]

Summary

The exchangeable Markov model is currently parametrized with the exchangeable intensities $\lambda_i$. However, the relevant quantity for the algorithm are the scaled exchangeable intensities $\binom{d}{I} \lambda_i$. Computationally, it is desirable to calculate the latter directly.

Proposal

Rewrite the algorithm such that the parameter ex_intensity has the interpretation of the scaled exchangeable intensity.

Additional context

Here are three mathematically equivalent ways to calculate the intensity matrix

qmatrix1 <- function(bf, d) {
  ex_intensities <- sapply(1:d, function(i) valueOf(bf, d-i, i))
  out <- matrix(0, nrow = d+1, ncol = d+1)
  for (i in 0:d) {
    if (i < d) {
      for (j in (i+1):d) {
        out[1+i, 1+j] <- rmo:::multiply_binomial_coefficient(
          sum(sapply(
            0:i, 
            function(k) {
              rmo:::multiply_binomial_coefficient(ex_intensities[k + j-i], i, k)
            })), d-i, j-i)
      }
      out[1+i, 1+i] <- -sum(out[1+i, ])
    }
  }

  out
}

qmatrix2 <- function(bf, d) {
  out <- matrix(0, nrow = d+1, ncol = d+1)
  for (i in 0:d) {
    if (i < d) {
      for (j in (i+1):d) {
        out[1+i, 1+j] <- valueOf(bf, d-j, j-i, n=d-i, k=j-i)
      }
      out[1+i, 1+i] <- -valueOf(bf, d-i, n = d-i, k = j-i)
    }
  }

  out
}

qmatrix3 <- function(bf, d) {
  ex_intensities <- sapply(
    1:d, 
    function(i) {
      valueOf(bf, d-i, i, n = d, k = i)
      })
  out <- matrix(0, nrow = d+1, ncol = d+1)
  out[1, -1] <- ex_intensities
  out[1, 1] <- -sum(out[1, -1])
  for (i in 1:d) {
    if (i < d) {
      for (j in (i+1):d) {
        out[1+i, 1+j] <- (d-j+1)/(d-i+1) * out[i, j] + (j+1-i)/(d-i+1) * out[i, 1+j]
      }
      out[1+i, 1+i] <- -sum(out[1+i, ])
    }
  }

  out
}