Try simulating the errors just once,

`ERR ~ MVN(nrep, 0, Sigma)`

, resulting in an`nrep*nsnp`

matrix. Let`mexp.zm`

be a matrix with`exp.zm`

replicated in each row (`nrep`

times). Then`zstar = exp.zm + ERR`

and we do not have to call`rmvnorm::mvtnorm`

for every SNP considered causal.The corresponding code from the

`corrcov`

function is,

```
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)
}
```

This method possibly makes the results dependent between different SNPs (only one ERR matrix is used - they have the same errors), but this should be ok for large numbers of simulations.

Check that the original method,

`zstar = rmvnorm(nrep, exp.zm, Sigma)`

, is equivalent to this new method,`zstar = exp.zm + ERR`

, where`ERR ~ N(0, Sigma)`

.