- tmp <- eigen(Q * phy$edge.length[i], symmetric = FALSE)
- P1 <- tmp$vectors %*% diag(exp(tmp$values)) %*% solve(tmp$vectors)
- tmp <- eigen(Q * phy$edge.length[j], symmetric = FALSE)
- P2 <- tmp$vectors %*% diag(exp(tmp$values)) %*% solve(tmp$vectors)
- liks[anc, ] <- P1 %*% liks[des1, ] * P2 %*% liks[des2, ]
+ v.l <- matexpo(Q * phy$edge.length[i]) %*% liks[des1, ]
+ v.r <- matexpo(Q * phy$edge.length[j]) %*% liks[des2, ]
+ v <- v.l * v.r
+ comp[anc] <- sum(v)
+ liks[anc, ] <- v/comp[anc]