# Simulates BM evolution more quickly. # A trend can be simulated by mu!=0. # mu=0 is standard BM; mu<0 downward trend; mu>0 upward trend. # Written by Liam J. Revell 2011 fastBM<-function(tree,a=0,mu=0,sig2=1,internal=FALSE){ n<-length(tree$tip) # first simulate changes along each branch x<-rnorm(n=length(tree$edge.length),mean=mu*tree$edge.length,sd=sqrt(sig2*tree$edge.length)) # now add them up y<-matrix(0,nrow(tree$edge),ncol(tree$edge)) for(i in 1:length(x)){ if(tree$edge[i,1]==(n+1)) y[i,1]<-a else y[i,1]<-y[match(tree$edge[i,1],tree$edge[,2]),2] y[i,2]<-y[i,1]+x[i] } rm(x); x<-c(y[1,1],y[,2]) names(x)<-c(n+1,tree$edge[,2]) x<-x[as.character(1:(n+tree$Nnode))] names(x)[1:n]<-tree$tip.label # return simulated data if(internal==TRUE) return(x) # include internal nodes else return(x[1:length(tree$tip.label)]) # tip nodes only }