1. Speeding up functions

ERR = mvtnorm:::rmvnorm(nrep, rep(0, ncol(Sigma)), Sigma)
pp_ERR = function(Zj, nrep, Sigma) {
    exp.zm = Zj %*% Sigma
    mexp.zm = matrix(exp.zm, nrep, length(Zj), byrow = TRUE)  # matrix of Zj replicated in each row
    zstar = mexp.zm + ERR
    bf = 0.5 * (log(1 - r) + (r * zstar^2))
    denom = coloc:::logsum(bf)  
    pp.tmp = exp(bf - denom)  
    pp.tmp/rowSums(pp.tmp)
}