Friday, November 25, 2011

Improvement to anc.Bayes()

I have already identified a major area in which I can easily improve the computation time of anc.Bayes(), my new function for Bayesian MCMC ancestral state estimation. Specifically, after having thought myself quite clever for using dmnorm() in the mnormt library (see post here), I realized that this was actually extremely computationally inefficient - especially for large matrices. This is because by using dmnorm() I'm essentially forced R to invert our VCV matrix (which is of order number of tips + number of nodes) every generation of the MCMC. This is totally unnecessary and I can easily create a custom function to compute the MVN probability density in which I avoid inverting the phylogenetic VCV matrix on each call. This function is below (it still requires that we invert C once), and I will post a new version of anc.Bayes() tomorrow.

likelihood<-function(C,invC,x,sig2,a,y){
  logLik<--(c(x,y)-a)%*%invC%*%(c(x,y)-a)/(2*sig2) -nrow(invC)*log(2*pi)/2 -determinant(sig2*C,log=T)$modulus[1]/2
  return(logLik)
}

No comments:

Post a Comment