use of dr.evolution.tree.MultivariateTraitTree in project beast-mcmc by beast-dev.
the class MultivariateTraitUtils method computeTreeTraitMean.
public static double[] computeTreeTraitMean(FullyConjugateMultivariateTraitLikelihood trait, double[] rootValue, boolean conditionOnRoot) {
double[] root = trait.getPriorMean();
if (conditionOnRoot) {
System.err.println("WARNING: Not yet fully implemented (conditioning on root in simulator)");
//root = new double[root.length];
root = rootValue;
}
final int nTaxa = trait.getTreeModel().getExternalNodeCount();
double[] mean = new double[root.length * nTaxa];
for (int i = 0; i < nTaxa; ++i) {
System.arraycopy(root, 0, mean, i * root.length, root.length);
}
if (trait.driftModels != null) {
MultivariateTraitTree myTreeModel = trait.getTreeModel();
for (int i = 0; i < nTaxa; ++i) {
double[] shiftContribution = getShiftContributionToMean(myTreeModel.getExternalNode(i), trait);
for (int j = 0; j < trait.dimTrait; ++j) {
mean[i * trait.dimTrait + j] = mean[i * trait.dimTrait + j] + shiftContribution[j];
}
}
}
return mean;
}
use of dr.evolution.tree.MultivariateTraitTree in project beast-mcmc by beast-dev.
the class MultivariateTraitUtils method computeTreeVariance.
public static double[][] computeTreeVariance(FullyConjugateMultivariateTraitLikelihood trait, boolean conditionOnRoot) {
MultivariateTraitTree treeModel = trait.getTreeModel();
final int tipCount = treeModel.getExternalNodeCount();
double[][] variance = new double[tipCount][tipCount];
for (int i = 0; i < tipCount; i++) {
// Fill in diagonal
double marginalTime = trait.getRescaledLengthToRoot(treeModel.getExternalNode(i));
variance[i][i] = marginalTime;
for (int j = i + 1; j < tipCount; j++) {
NodeRef mrca = findMRCA(trait, i, j);
if (DEBUG) {
System.err.println(trait.getTreeModel().getRoot().getNumber());
System.err.print("Taxa pair: " + i + " : " + j + " (" + mrca.getNumber() + ") = ");
}
double length = trait.getRescaledLengthToRoot(mrca);
if (DEBUG) {
System.err.println(length);
}
variance[i][j] = length;
}
}
// Make symmetric
for (int i = 0; i < tipCount; i++) {
for (int j = i + 1; j < tipCount; j++) {
variance[j][i] = variance[i][j];
}
}
if (DEBUG) {
System.err.println("");
System.err.println("New tree conditional variance:\n" + new Matrix(variance));
}
if (!conditionOnRoot) {
double priorSampleSize = trait.getPriorSampleSize();
for (int i = 0; i < tipCount; ++i) {
for (int j = 0; j < tipCount; ++j) {
variance[i][j] += 1.0 / priorSampleSize;
}
}
if (DEBUG) {
System.err.println("");
System.err.println("New tree unconditional variance:\n" + new Matrix(variance));
}
}
return variance;
}
Aggregations