The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.

Tracing the likelihood calculation of a Gaussian model

Venelin Mitov

2024-03-15

This tutorial shows how the pruning likelihood calculation algorithm works for \(\mathcal{G}_{LInv}\) models. This is an advanced topic, which can be helpful in case you wish to validate another implementation of the algorithm, e.g. in a different programming language, such as Java. For details about the mathematical terms and equations involved in the following calculations, refer to (Mitov et al. 2019).

Example tree and data

library(ape); 
library(PCMBase);

# Non-ultrametric phylogenetic tree of 5 tips in both examples:
treeNewick <- "((5:0.8,4:1.8)7:1.5,(((3:0.8,2:1.6)6:0.7)8:0.6,1:2.6)9:0.9)0;"
tree <- PCMTree(read.tree(text = treeNewick))
# Partitioning the tree in two parts and assign the regimes:
PCMTreeSetPartRegimes(tree, part.regime = c(`6`=2), setPartition = TRUE, inplace = TRUE)

pOrder <- c(PCMTreeGetLabels(tree)[tree$edge[PCMTreePostorder(tree), 2]], "0")

# Trait-data:
X <- cbind(
  c(0.3, NaN, 1.4), 
  c(0.1, NaN, NA), 
  c(0.2, NaN, 1.2), 
  c(NA, 0.2, 0.2), 
  c(NA, 1.2, 0.4))

colnames(X) <- as.character(1:5)

A tree with five tips and two evolutionary regimes

A mixed Gaussian model

model.OU.BM <- MixedGaussian(
  k = nrow(X), 
  modelTypes = c(
    BM = "BM__Omitted_X0__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigma_x__Omitted_Sigmae_x",
    OU = "OU__Omitted_X0__H__Theta__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigma_x__Omitted_Sigmae_x"), 
  mapping = c(2, 1), 
  Sigmae_x = structure(
    0, 
    class = c("MatrixParameter", "_Omitted", 
              description = "upper triangular factor of the non-phylogenetic variance-covariance")))

model.OU.BM <- PCMApplyTransformation(model.OU.BM)
model.OU.BM$X0[] <- c(NA, NA, NA)
model.OU.BM$`1`$H[,,1] <- cbind(
  c(.1, -.7, .6), 
  c(1.3, 2.2, -1.4), 
  c(0.8, 0.2, 0.9))
model.OU.BM$`1`$Theta[] <- c(1.3, -.5, .2)
model.OU.BM$`1`$Sigma_x[,,1] <- cbind(
  c(1, 0, 0), 
  c(1.0, 0.5, 0), 
  c(0.3, -.8, 1))

model.OU.BM$`2`$Sigma_x[,,1] <- cbind(
  c(0.8, 0, 0), 
  c(1, 0.3, 0), 
  c(0.4, 0.5, 0.3))

print(
  PCMTable(model.OU.BM, removeUntransformed = FALSE), 
  xtable = TRUE, type=tableOutputType)


regime type \(X0\) \(H\) \(\Theta\) \(\Sigma_{x}\) \(\Sigma\)
:global: NA \(\begin{bmatrix}{} . \\ . \\ . \\ \end{bmatrix}\)
1 OU \(\begin{bmatrix}{} 0.10 & 1.30 & 0.80 \\ -0.70 & 2.20 & 0.20 \\ 0.60 & -1.40 & 0.90 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.30 \\ -0.50 \\ 0.20 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.00 & 1.00 & 0.30 \\ 0.00 & 0.50 & -0.80 \\ 0.00 & 0.00 & 1.00 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 2.09 & 0.26 & 0.30 \\ 0.26 & 0.89 & -0.80 \\ 0.30 & -0.80 & 1.00 \\ \end{bmatrix}\)
2 BM \(\begin{bmatrix}{} 0.80 & 1.00 & 0.40 \\ 0.00 & 0.30 & 0.50 \\ 0.00 & 0.00 & 0.30 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.80 & 0.50 & 0.12 \\ 0.50 & 0.34 & 0.15 \\ 0.12 & 0.15 & 0.09 \\ \end{bmatrix}\)

Likelihood values for the three example variants

options(digits = 4)
# Variant 1: 
PCMLik(X[, tree$tip.label], tree, model.OU.BM)
## [1] -11.92
## attr(,"X0")
## [1]  9.566 -6.349 15.254
# Variant 2: First we call the function PCMInfo to obtain a meta-information object.
metaI.variant2 <- PCMInfo(X[, tree$tip.label], tree, model.OU.BM)
# Then, we manually change the vector of present coordinates for the root node.
# The pc-matrix is a k x M matrix of logical values, each column corresponding
# to a node. The active coordinates are indicated by the TRUE entries.
# To prevent assigning to the wrong column in the pc-table, we first assign
# the node-labels as column nanmes.
colnames(metaI.variant2$pc) <- PCMTreeGetLabels(tree)
metaI.variant2$pc[, "0"] <- c(TRUE, FALSE, TRUE)
# After the change, the pc-matrix looks like this:
metaI.variant2$pc
##          5     4     3     2     1     0    7     9     8     6
## [1,] FALSE FALSE  TRUE  TRUE  TRUE  TRUE TRUE  TRUE  TRUE  TRUE
## [2,]  TRUE  TRUE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE
## [3,]  TRUE  TRUE  TRUE FALSE  TRUE  TRUE TRUE  TRUE  TRUE  TRUE
# And the log-likelihood value is:
PCMLik(X[, tree$tip.label], tree, model.OU.BM, metaI = metaI.variant2)
## [1] -12.17
## attr(,"X0")
## [1] 10.611    NaN  8.747
# Variant 3: We set all NaN values in X to NA, to indicate that these are
# missing measurements
X3 <- X
X3[is.nan(X3)] <- NA_real_
PCMLik(X3[, tree$tip.label], tree, model.OU.BM)
## [1] -10.71
## attr(,"X0")
## [1]  15.99  18.34 -11.95

Tracing the likelihood calculation using the function PCMLikTrace

Variant 1

traceTable1 <- PCMLikTrace(X[, tree$tip.label], tree, model.OU.BM)
traceTable1[, node:=.I]
setkey(traceTable1, i)

Variant 2

traceTable2 <- PCMLikTrace(
  X[, tree$tip.label], tree, model.OU.BM, metaI = metaI.variant2)
traceTable2[, node:=.I]
setkey(traceTable2, i)

Variant 3

traceTable3 <- PCMLikTrace(X3[, tree$tip.label], tree, model.OU.BM)
traceTable3[, node:=.I]
setkey(traceTable3, i)

A step by step description of the log-likelihood calculation

Step 1: Calculating \(\vec{\omega}\), \(\mathbf{\Phi}\) and \(\mathbf{V}\) for each tip or internal node}

Calculating \(\vec{\omega}\), \(\mathbf{\Phi}\) and \(\mathbf{V}\) for a node in an OU regime

# OU parameters:
H <- model.OU.BM$`1`$H[,,1]
theta <- model.OU.BM$`1`$Theta[,1]
Sigma <- model.OU.BM$`1`$Sigma_x[,,1] %*% t(model.OU.BM$`1`$Sigma_x[,,1])

# Eigenvalues of H: these can be complex numbers
lambda <- eigen(H)$values

# Matrix of eigenvectors of H: again, these can be complex
P <- eigen(H)$vectors
P_1 <- solve(P)

# vectors of active coordinates:
pc <- PCMInfo(X[, tree$tip.label], tree, model.OU.BM)$pc

# length of the branch leading to tip 1 (2.6):
t1 <- PCMTreeDtNodes(tree)[endNodeLab == "1", endTime - startTime]

# active coordinates for tip 1 and its parent:
k1 <- pc[, match("1", PCMTreeGetLabels(tree))]
k9 <- pc[, match("9", PCMTreeGetLabels(tree))]

# k x k matrix formed from the pairs of lambda-values and t1 (see Eq. 19):
LambdaMat <- matrix(0, 3, 3)
for(i in 1:3) 
  for(j in 1:3) 
    LambdaMat[i,j] <- 1/(lambda[i]+lambda[j])*(1-exp(-(lambda[i]+lambda[j])*t1))

# omega, Phi, V for tip 1:
print(omega1 <- (diag(1, 3, 3)[k1, ] - expm::expm(-H*t1)[k1, ]) %*% theta[])
##        [,1]
## [1,] 0.6151
## [2,] 0.3392
print(Phi1 <- expm::expm(-H*t1)[k1, k9])
##          [,1]    [,2]
## [1,]  0.37818 -0.3461
## [2,] -0.06011  0.1261
print(V1 <- (P %*% (LambdaMat * (P_1 %*% Sigma %*% t(P_1))) %*% t(P))[k1, k1])
##             [,1]        [,2]
## [1,]  2.21599+0i -0.07466+0i
## [2,] -0.07466-0i  0.25787+0i

Calculating \(\vec{\omega}\), \(\mathbf{\Phi}\) and \(\mathbf{V}\) for a node in a BM regime

# BM parameter:
Sigma <- model.OU.BM$`2`$Sigma_x[,,1] %*% t(model.OU.BM$`2`$Sigma_x[,,1])

# vectors of active coordinates:
pc <- PCMInfo(X[, tree$tip.label], tree, model.OU.BM)$pc

# length of the branch leading to tip 2 (1.6):
t2 <- PCMTreeDtNodes(tree)[endNodeLab == "2", endTime - startTime]

# active coordinates for tip 1 and its parent:
k2 <- pc[, match("2", PCMTreeGetLabels(tree))]
k6 <- pc[, match("6", PCMTreeGetLabels(tree))]

# omega, Phi, V for tip 1:
print(omega2 <- as.matrix(rep(0, 3)[k2]))
##      [,1]
## [1,]    0
print(Phi2 <- as.matrix(diag(1, 3, 3)[k2, k6]))
##      [,1]
## [1,]    1
## [2,]    0
print(V2 <- as.matrix((t2*Sigma)[k2, k2]))
##      [,1]
## [1,] 2.88

Variant 1

options(digits = 2)
cat(FormatTableAsLatex(
  traceTable1[list(pOrder), list(j, i, t_i, k_i, omega_i, Phi_i, V_i, V_1_i)], 
  type = tableOutputType))









\(j\) \(i\) \(t_i\) \(k_i\) \(\omega_i\) \(\Phi_i\) \(V_i\) \(V^{-1}_i\)
6 2 \(\begin{bmatrix}{} 1.60 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.00 \\ . \\ . \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.00 & . & 0.00 \\ . & . & . \\ . & . & . \\ \end{bmatrix}\) \(\begin{bmatrix}{} 2.88 & . & . \\ . & . & . \\ . & . & . \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.35 & . & . \\ . & . & . \\ . & . & . \\ \end{bmatrix}\)
6 3 \(\begin{bmatrix}{} 0.80 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.00 \\ . \\ 0.00 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.00 & . & 0.00 \\ . & . & . \\ 0.00 & . & 1.00 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.44 & . & 0.10 \\ . & . & . \\ 0.10 & . & 0.07 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.76 & . & -1.02 \\ . & . & . \\ -1.02 & . & 15.24 \\ \end{bmatrix}\)
8 6 \(\begin{bmatrix}{} 0.70 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.00 \\ . \\ 0.00 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.00 & . & 0.00 \\ . & . & . \\ 0.00 & . & 1.00 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.26 & . & 0.08 \\ . & . & . \\ 0.08 & . & 0.06 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.87 & . & -1.16 \\ . & . & . \\ -1.16 & . & 17.42 \\ \end{bmatrix}\)
9 1 \(\begin{bmatrix}{} 2.60 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.62 \\ . \\ 0.34 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.38 & . & -0.35 \\ . & . & . \\ -0.06 & . & 0.13 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 2.22 & . & -0.07 \\ . & . & . \\ -0.07 & . & 0.26 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.46 & . & 0.13 \\ . & . & . \\ 0.13 & . & 3.92 \\ \end{bmatrix}\)
9 8 \(\begin{bmatrix}{} 0.60 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.04 \\ . \\ 0.49 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.89 & . & -0.32 \\ . & . & . \\ -0.17 & . & 0.60 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.00 & . & 0.06 \\ . & . & . \\ 0.06 & . & 0.21 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.01 & . & -0.26 \\ . & . & . \\ -0.26 & . & 4.77 \\ \end{bmatrix}\)
7 4 \(\begin{bmatrix}{} 1.80 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 2 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} . \\ -0.86 \\ 0.44 \\ \end{bmatrix}\) \(\begin{bmatrix}{} . & . & . \\ 0.22 & -0.21 & -0.16 \\ -0.11 & 0.29 & 0.24 \\ \end{bmatrix}\) \(\begin{bmatrix}{} . & . & . \\ . & 0.37 & -0.21 \\ . & -0.21 & 0.25 \\ \end{bmatrix}\) \(\begin{bmatrix}{} . & . & . \\ . & 5.16 & 4.30 \\ . & 4.30 & 7.58 \\ \end{bmatrix}\)
7 5 \(\begin{bmatrix}{} 0.80 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 2 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} . \\ -0.78 \\ 0.53 \\ \end{bmatrix}\) \(\begin{bmatrix}{} . & . & . \\ 0.25 & 0.03 & -0.12 \\ -0.17 & 0.42 & 0.51 \\ \end{bmatrix}\) \(\begin{bmatrix}{} . & . & . \\ . & 0.27 & -0.18 \\ . & -0.18 & 0.22 \\ \end{bmatrix}\) \(\begin{bmatrix}{} . & . & . \\ . & 7.78 & 6.15 \\ . & 6.15 & 9.30 \\ \end{bmatrix}\)
0 9 \(\begin{bmatrix}{} 0.90 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.02 \\ . \\ 0.53 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.81 & -0.61 & -0.39 \\ . & . & . \\ -0.17 & 0.42 & 0.47 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.35 & . & 0.03 \\ . & . & . \\ 0.03 & . & 0.23 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.74 & . & -0.11 \\ . & . & . \\ -0.11 & . & 4.38 \\ \end{bmatrix}\)
0 7 \(\begin{bmatrix}{} 1.50 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ 2 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.22 \\ -0.88 \\ 0.49 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.64 & -0.67 & -0.43 \\ 0.24 & -0.18 & -0.16 \\ -0.13 & 0.34 & 0.30 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.82 & 0.45 & -0.02 \\ 0.45 & 0.34 & -0.20 \\ -0.02 & -0.20 & 0.25 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.29 & -3.11 & -2.44 \\ -3.11 & 13.06 & 10.44 \\ -2.44 & 10.44 & 12.43 \\ \end{bmatrix}\)
ND 0 ND \(\begin{bmatrix}{} 1 \\ 2 \\ 3 \\ \end{bmatrix}\) ND ND ND ND

Variant 2

cat(FormatTableAsLatex(
  traceTable2[list(pOrder), list(j, i, t_i, k_i, omega_i, Phi_i, V_i, V_1_i)], 
  type = tableOutputType))









\(j\) \(i\) \(t_i\) \(k_i\) \(\omega_i\) \(\Phi_i\) \(V_i\) \(V^{-1}_i\)
6 2 \(\begin{bmatrix}{} 1.60 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.00 \\ . \\ . \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.00 & . & 0.00 \\ . & . & . \\ . & . & . \\ \end{bmatrix}\) \(\begin{bmatrix}{} 2.88 & . & . \\ . & . & . \\ . & . & . \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.35 & . & . \\ . & . & . \\ . & . & . \\ \end{bmatrix}\)
6 3 \(\begin{bmatrix}{} 0.80 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.00 \\ . \\ 0.00 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.00 & . & 0.00 \\ . & . & . \\ 0.00 & . & 1.00 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.44 & . & 0.10 \\ . & . & . \\ 0.10 & . & 0.07 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.76 & . & -1.02 \\ . & . & . \\ -1.02 & . & 15.24 \\ \end{bmatrix}\)
8 6 \(\begin{bmatrix}{} 0.70 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.00 \\ . \\ 0.00 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.00 & . & 0.00 \\ . & . & . \\ 0.00 & . & 1.00 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.26 & . & 0.08 \\ . & . & . \\ 0.08 & . & 0.06 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.87 & . & -1.16 \\ . & . & . \\ -1.16 & . & 17.42 \\ \end{bmatrix}\)
9 1 \(\begin{bmatrix}{} 2.60 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.62 \\ . \\ 0.34 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.38 & . & -0.35 \\ . & . & . \\ -0.06 & . & 0.13 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 2.22 & . & -0.07 \\ . & . & . \\ -0.07 & . & 0.26 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.46 & . & 0.13 \\ . & . & . \\ 0.13 & . & 3.92 \\ \end{bmatrix}\)
9 8 \(\begin{bmatrix}{} 0.60 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.04 \\ . \\ 0.49 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.89 & . & -0.32 \\ . & . & . \\ -0.17 & . & 0.60 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.00 & . & 0.06 \\ . & . & . \\ 0.06 & . & 0.21 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.01 & . & -0.26 \\ . & . & . \\ -0.26 & . & 4.77 \\ \end{bmatrix}\)
7 4 \(\begin{bmatrix}{} 1.80 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 2 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} . \\ -0.86 \\ 0.44 \\ \end{bmatrix}\) \(\begin{bmatrix}{} . & . & . \\ 0.22 & -0.21 & -0.16 \\ -0.11 & 0.29 & 0.24 \\ \end{bmatrix}\) \(\begin{bmatrix}{} . & . & . \\ . & 0.37 & -0.21 \\ . & -0.21 & 0.25 \\ \end{bmatrix}\) \(\begin{bmatrix}{} . & . & . \\ . & 5.16 & 4.30 \\ . & 4.30 & 7.58 \\ \end{bmatrix}\)
7 5 \(\begin{bmatrix}{} 0.80 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 2 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} . \\ -0.78 \\ 0.53 \\ \end{bmatrix}\) \(\begin{bmatrix}{} . & . & . \\ 0.25 & 0.03 & -0.12 \\ -0.17 & 0.42 & 0.51 \\ \end{bmatrix}\) \(\begin{bmatrix}{} . & . & . \\ . & 0.27 & -0.18 \\ . & -0.18 & 0.22 \\ \end{bmatrix}\) \(\begin{bmatrix}{} . & . & . \\ . & 7.78 & 6.15 \\ . & 6.15 & 9.30 \\ \end{bmatrix}\)
0 9 \(\begin{bmatrix}{} 0.90 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.02 \\ . \\ 0.53 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.81 & . & -0.39 \\ . & . & . \\ -0.17 & . & 0.47 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.35 & . & 0.03 \\ . & . & . \\ 0.03 & . & 0.23 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.74 & . & -0.11 \\ . & . & . \\ -0.11 & . & 4.38 \\ \end{bmatrix}\)
0 7 \(\begin{bmatrix}{} 1.50 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ 2 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.22 \\ -0.88 \\ 0.49 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.64 & . & -0.43 \\ 0.24 & . & -0.16 \\ -0.13 & . & 0.30 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.82 & 0.45 & -0.02 \\ 0.45 & 0.34 & -0.20 \\ -0.02 & -0.20 & 0.25 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.29 & -3.11 & -2.44 \\ -3.11 & 13.06 & 10.44 \\ -2.44 & 10.44 & 12.43 \\ \end{bmatrix}\)
ND 0 ND \(\begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix}\) ND ND ND ND

Variant 3

options(digits = 2)
cat(FormatTableAsLatex(
  traceTable3[list(pOrder), list(j, i, t_i, k_i, omega_i, Phi_i, V_i, V_1_i)], 
  type = tableOutputType))









\(j\) \(i\) \(t_i\) \(k_i\) \(\omega_i\) \(\Phi_i\) \(V_i\) \(V^{-1}_i\)
6 2 \(\begin{bmatrix}{} 1.60 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.00 \\ . \\ . \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.00 & 0.00 & 0.00 \\ . & . & . \\ . & . & . \\ \end{bmatrix}\) \(\begin{bmatrix}{} 2.88 & . & . \\ . & . & . \\ . & . & . \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.35 & . & . \\ . & . & . \\ . & . & . \\ \end{bmatrix}\)
6 3 \(\begin{bmatrix}{} 0.80 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.00 \\ . \\ 0.00 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.00 & 0.00 & 0.00 \\ . & . & . \\ 0.00 & 0.00 & 1.00 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.44 & . & 0.10 \\ . & . & . \\ 0.10 & . & 0.07 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.76 & . & -1.02 \\ . & . & . \\ -1.02 & . & 15.24 \\ \end{bmatrix}\)
8 6 \(\begin{bmatrix}{} 0.70 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ 2 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.00 \\ 0.00 \\ 0.00 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.00 & 0.00 & 0.00 \\ 0.00 & 1.00 & 0.00 \\ 0.00 & 0.00 & 1.00 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.26 & 0.35 & 0.08 \\ 0.35 & 0.24 & 0.10 \\ 0.08 & 0.10 & 0.06 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 2.23 & -7.44 & 9.42 \\ -7.44 & 40.67 & -57.87 \\ 9.42 & -57.87 & 99.76 \\ \end{bmatrix}\)
9 1 \(\begin{bmatrix}{} 2.60 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.62 \\ . \\ 0.34 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.38 & -0.52 & -0.35 \\ . & . & . \\ -0.06 & 0.17 & 0.13 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 2.22 & . & -0.07 \\ . & . & . \\ -0.07 & . & 0.26 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.46 & . & 0.13 \\ . & . & . \\ 0.13 & . & 3.92 \\ \end{bmatrix}\)
9 8 \(\begin{bmatrix}{} 0.60 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ 2 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.04 \\ -0.69 \\ 0.49 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.89 & -0.51 & -0.32 \\ 0.23 & 0.17 & -0.10 \\ -0.17 & 0.40 & 0.60 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.00 & 0.19 & 0.06 \\ 0.19 & 0.24 & -0.17 \\ 0.06 & -0.17 & 0.21 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 2.09 & -4.60 & -4.16 \\ -4.60 & 19.50 & 16.55 \\ -4.16 & 16.55 & 18.82 \\ \end{bmatrix}\)
7 4 \(\begin{bmatrix}{} 1.80 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 2 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} . \\ -0.86 \\ 0.44 \\ \end{bmatrix}\) \(\begin{bmatrix}{} . & . & . \\ 0.22 & -0.21 & -0.16 \\ -0.11 & 0.29 & 0.24 \\ \end{bmatrix}\) \(\begin{bmatrix}{} . & . & . \\ . & 0.37 & -0.21 \\ . & -0.21 & 0.25 \\ \end{bmatrix}\) \(\begin{bmatrix}{} . & . & . \\ . & 5.16 & 4.30 \\ . & 4.30 & 7.58 \\ \end{bmatrix}\)
7 5 \(\begin{bmatrix}{} 0.80 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 2 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} . \\ -0.78 \\ 0.53 \\ \end{bmatrix}\) \(\begin{bmatrix}{} . & . & . \\ 0.25 & 0.03 & -0.12 \\ -0.17 & 0.42 & 0.51 \\ \end{bmatrix}\) \(\begin{bmatrix}{} . & . & . \\ . & 0.27 & -0.18 \\ . & -0.18 & 0.22 \\ \end{bmatrix}\) \(\begin{bmatrix}{} . & . & . \\ . & 7.78 & 6.15 \\ . & 6.15 & 9.30 \\ \end{bmatrix}\)
0 9 \(\begin{bmatrix}{} 0.90 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ 2 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.02 \\ -0.81 \\ 0.53 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.81 & -0.61 & -0.39 \\ 0.25 & -0.02 & -0.13 \\ -0.17 & 0.42 & 0.47 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.35 & 0.29 & 0.03 \\ 0.29 & 0.28 & -0.18 \\ 0.03 & -0.18 & 0.23 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.60 & -3.66 & -3.14 \\ -3.66 & 15.62 & 12.92 \\ -3.14 & 12.92 & 15.08 \\ \end{bmatrix}\)
0 7 \(\begin{bmatrix}{} 1.50 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ 2 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.22 \\ -0.88 \\ 0.49 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.64 & -0.67 & -0.43 \\ 0.24 & -0.18 & -0.16 \\ -0.13 & 0.34 & 0.30 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.82 & 0.45 & -0.02 \\ 0.45 & 0.34 & -0.20 \\ -0.02 & -0.20 & 0.25 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.29 & -3.11 & -2.44 \\ -3.11 & 13.06 & 10.44 \\ -2.44 & 10.44 & 12.43 \\ \end{bmatrix}\)
ND 0 ND \(\begin{bmatrix}{} 1 \\ 2 \\ 3 \\ \end{bmatrix}\) ND ND ND ND

Step 2: Calculating \(\mathbf{A}\), \(\vec{b}\), \(\mathbf{C}\), \(\vec{d}\), \(\mathbf{E}\) and \(f\) for tips and internal nodes

# For tip 1. We directly apply Eq. 2, Thm 1:
# We can safely use the real part of V1 (all imaginary parts are 0):
print(V1)
##           [,1]      [,2]
## [1,]  2.216+0i -0.075+0i
## [2,] -0.075-0i  0.258+0i
V1 <- Re(V1)
V1_1 <- solve(V1)

print(A1 <- -0.5*V1_1)
##        [,1]   [,2]
## [1,] -0.228 -0.066
## [2,] -0.066 -1.958
print(E1 <- t(Phi1) %*% V1_1)
##       [,1]  [,2]
## [1,]  0.16 -0.19
## [2,] -0.14  0.45
print(b1 <- V1_1 %*% omega1)
##      [,1]
## [1,] 0.33
## [2,] 1.41
print(C1 <- -0.5 * E1 %*% Phi1)
##        [,1]   [,2]
## [1,] -0.037  0.040
## [2,]  0.040 -0.053
print(d1 <- -E1 %*% omega1)
##        [,1]
## [1,] -0.038
## [2,] -0.065
print(f1 <- -0.5 * (t(omega1) %*% V1_1 %*% omega1 + sum(k1)*log(2*pi) + log(det(V1))))
##      [,1]
## [1,] -1.9

Variant 1

cat(FormatTableAsLatex(
  traceTable1[list(pOrder), list(j, i, k_i, A_i, b_i, C_i, d_i, E_i, f_i)], 
  type = tableOutputType))









\(j\) \(i\) \(k_i\) \(A_i\) \(b_i\) \(C_i\) \(d_i\) \(E_i\) \(f_i\)
6 2 \(\begin{bmatrix}{} 1 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.17 & . & . \\ . & . & . \\ . & . & . \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.00 \\ . \\ . \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.17 & . & -0.00 \\ . & . & . \\ -0.00 & . & -0.00 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.00 \\ . \\ -0.00 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.35 & . & . \\ . & . & . \\ 0.00 & . & . \\ \end{bmatrix}\) \(\begin{bmatrix}{} -1.45 \\ \end{bmatrix}\)
6 3 \(\begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.38 & . & 0.51 \\ . & . & . \\ 0.51 & . & -7.62 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.00 \\ . \\ 0.00 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.38 & . & 0.51 \\ . & . & . \\ 0.51 & . & -7.62 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.00 \\ . \\ 0.00 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.76 & . & -1.02 \\ . & . & . \\ -1.02 & . & 15.24 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.66 \\ \end{bmatrix}\)
8 6 \(\begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.44 & . & 0.58 \\ . & . & . \\ 0.58 & . & -8.71 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.00 \\ . \\ 0.00 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.44 & . & 0.58 \\ . & . & . \\ 0.58 & . & -8.71 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.00 \\ . \\ 0.00 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.87 & . & -1.16 \\ . & . & . \\ -1.16 & . & 17.42 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.52 \\ \end{bmatrix}\)
9 1 \(\begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.23 & . & -0.07 \\ . & . & . \\ -0.07 & . & -1.96 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.33 \\ . \\ 1.41 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.04 & . & 0.04 \\ . & . & . \\ 0.04 & . & -0.05 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.04 \\ . \\ -0.07 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.16 & . & -0.19 \\ . & . & . \\ -0.14 & . & 0.45 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -1.89 \\ \end{bmatrix}\)
9 8 \(\begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.51 & . & 0.13 \\ . & . & . \\ 0.13 & . & -2.39 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.17 \\ . \\ 2.37 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.50 & . & 0.46 \\ . & . & . \\ 0.46 & . & -0.96 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.54 \\ . \\ -1.48 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.94 & . & -1.02 \\ . & . & . \\ -0.48 & . & 2.95 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -1.65 \\ \end{bmatrix}\)
7 4 \(\begin{bmatrix}{} 2 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} . & . & . \\ . & -2.58 & -2.15 \\ . & -2.15 & -3.79 \\ \end{bmatrix}\) \(\begin{bmatrix}{} . \\ -2.54 \\ -0.35 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.07 & 0.05 & 0.04 \\ 0.05 & -0.17 & -0.14 \\ 0.04 & -0.14 & -0.11 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.53 \\ -0.43 \\ -0.32 \\ \end{bmatrix}\) \(\begin{bmatrix}{} . & 0.67 & 0.11 \\ . & 0.17 & 1.31 \\ . & 0.19 & 1.10 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -1.34 \\ \end{bmatrix}\)
7 5 \(\begin{bmatrix}{} 2 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} . & . & . \\ . & -3.89 & -3.07 \\ . & -3.07 & -4.65 \\ \end{bmatrix}\) \(\begin{bmatrix}{} . \\ -2.85 \\ 0.10 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.12 & -0.00 & 0.07 \\ -0.00 & -0.89 & -0.87 \\ 0.07 & -0.87 & -0.89 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.73 \\ 0.05 \\ -0.40 \\ \end{bmatrix}\) \(\begin{bmatrix}{} . & 0.88 & -0.06 \\ . & 2.81 & 4.07 \\ . & 2.19 & 4.00 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -1.21 \\ \end{bmatrix}\)
0 9 \(\begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.37 & . & 0.06 \\ . & . & . \\ 0.06 & . & -2.19 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.04 \\ . \\ 2.34 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.32 & 0.36 & 0.32 \\ 0.36 & -0.55 & -0.55 \\ 0.32 & -0.55 & -0.57 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.43 \\ -1.00 \\ -1.12 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.62 & . & -0.83 \\ -0.50 & . & 1.89 \\ -0.34 & . & 2.11 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -1.87 \\ \end{bmatrix}\)
0 7 \(\begin{bmatrix}{} 1 \\ 2 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.64 & 1.55 & 1.22 \\ 1.55 & -6.53 & -5.22 \\ 1.22 & -5.22 & -6.22 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.82 \\ -7.04 \\ -3.63 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.15 & 0.23 & 0.17 \\ 0.23 & -0.76 & -0.58 \\ 0.17 & -0.58 & -0.44 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.06 \\ 1.16 \\ 0.75 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.40 & -0.22 & -0.70 \\ -1.13 & 3.25 & 3.98 \\ -0.79 & 2.38 & 3.09 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -3.47 \\ \end{bmatrix}\)
ND 0 \(\begin{bmatrix}{} 1 \\ 2 \\ 3 \\ \end{bmatrix}\) ND ND ND ND ND ND

Variant 2

cat(FormatTableAsLatex(
  traceTable2[list(pOrder), list(j, i, k_i, A_i, b_i, C_i, d_i, E_i, f_i)], 
  type = tableOutputType))









\(j\) \(i\) \(k_i\) \(A_i\) \(b_i\) \(C_i\) \(d_i\) \(E_i\) \(f_i\)
6 2 \(\begin{bmatrix}{} 1 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.17 & . & . \\ . & . & . \\ . & . & . \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.00 \\ . \\ . \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.17 & . & -0.00 \\ . & . & . \\ -0.00 & . & -0.00 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.00 \\ . \\ -0.00 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.35 & . & . \\ . & . & . \\ 0.00 & . & . \\ \end{bmatrix}\) \(\begin{bmatrix}{} -1.45 \\ \end{bmatrix}\)
6 3 \(\begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.38 & . & 0.51 \\ . & . & . \\ 0.51 & . & -7.62 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.00 \\ . \\ 0.00 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.38 & . & 0.51 \\ . & . & . \\ 0.51 & . & -7.62 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.00 \\ . \\ 0.00 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.76 & . & -1.02 \\ . & . & . \\ -1.02 & . & 15.24 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.66 \\ \end{bmatrix}\)
8 6 \(\begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.44 & . & 0.58 \\ . & . & . \\ 0.58 & . & -8.71 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.00 \\ . \\ 0.00 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.44 & . & 0.58 \\ . & . & . \\ 0.58 & . & -8.71 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.00 \\ . \\ 0.00 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.87 & . & -1.16 \\ . & . & . \\ -1.16 & . & 17.42 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.52 \\ \end{bmatrix}\)
9 1 \(\begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.23 & . & -0.07 \\ . & . & . \\ -0.07 & . & -1.96 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.33 \\ . \\ 1.41 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.04 & . & 0.04 \\ . & . & . \\ 0.04 & . & -0.05 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.04 \\ . \\ -0.07 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.16 & . & -0.19 \\ . & . & . \\ -0.14 & . & 0.45 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -1.89 \\ \end{bmatrix}\)
9 8 \(\begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.51 & . & 0.13 \\ . & . & . \\ 0.13 & . & -2.39 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.17 \\ . \\ 2.37 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.50 & . & 0.46 \\ . & . & . \\ 0.46 & . & -0.96 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.54 \\ . \\ -1.48 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.94 & . & -1.02 \\ . & . & . \\ -0.48 & . & 2.95 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -1.65 \\ \end{bmatrix}\)
7 4 \(\begin{bmatrix}{} 2 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} . & . & . \\ . & -2.58 & -2.15 \\ . & -2.15 & -3.79 \\ \end{bmatrix}\) \(\begin{bmatrix}{} . \\ -2.54 \\ -0.35 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.07 & 0.05 & 0.04 \\ 0.05 & -0.17 & -0.14 \\ 0.04 & -0.14 & -0.11 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.53 \\ -0.43 \\ -0.32 \\ \end{bmatrix}\) \(\begin{bmatrix}{} . & 0.67 & 0.11 \\ . & 0.17 & 1.31 \\ . & 0.19 & 1.10 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -1.34 \\ \end{bmatrix}\)
7 5 \(\begin{bmatrix}{} 2 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} . & . & . \\ . & -3.89 & -3.07 \\ . & -3.07 & -4.65 \\ \end{bmatrix}\) \(\begin{bmatrix}{} . \\ -2.85 \\ 0.10 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.12 & -0.00 & 0.07 \\ -0.00 & -0.89 & -0.87 \\ 0.07 & -0.87 & -0.89 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.73 \\ 0.05 \\ -0.40 \\ \end{bmatrix}\) \(\begin{bmatrix}{} . & 0.88 & -0.06 \\ . & 2.81 & 4.07 \\ . & 2.19 & 4.00 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -1.21 \\ \end{bmatrix}\)
0 9 \(\begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.37 & . & 0.06 \\ . & . & . \\ 0.06 & . & -2.19 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.04 \\ . \\ 2.34 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.32 & . & 0.32 \\ . & . & . \\ 0.32 & . & -0.57 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.43 \\ . \\ -1.12 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.62 & . & -0.83 \\ . & . & . \\ -0.34 & . & 2.11 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -1.87 \\ \end{bmatrix}\)
0 7 \(\begin{bmatrix}{} 1 \\ 2 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.64 & 1.55 & 1.22 \\ 1.55 & -6.53 & -5.22 \\ 1.22 & -5.22 & -6.22 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.82 \\ -7.04 \\ -3.63 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.15 & . & 0.17 \\ . & . & . \\ 0.17 & . & -0.44 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.06 \\ . \\ 0.75 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.40 & -0.22 & -0.70 \\ . & . & . \\ -0.79 & 2.38 & 3.09 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -3.47 \\ \end{bmatrix}\)
ND 0 \(\begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix}\) ND ND ND ND ND ND

Variant 3

cat(FormatTableAsLatex(
  traceTable3[list(pOrder), list(j, i, k_i, A_i, b_i, C_i, d_i, E_i, f_i)], 
  type = tableOutputType))









\(j\) \(i\) \(k_i\) \(A_i\) \(b_i\) \(C_i\) \(d_i\) \(E_i\) \(f_i\)
6 2 \(\begin{bmatrix}{} 1 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.17 & . & . \\ . & . & . \\ . & . & . \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.00 \\ . \\ . \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.17 & -0.00 & -0.00 \\ -0.00 & -0.00 & -0.00 \\ -0.00 & -0.00 & -0.00 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.00 \\ -0.00 \\ -0.00 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.35 & . & . \\ 0.00 & . & . \\ 0.00 & . & . \\ \end{bmatrix}\) \(\begin{bmatrix}{} -1.45 \\ \end{bmatrix}\)
6 3 \(\begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.38 & . & 0.51 \\ . & . & . \\ 0.51 & . & -7.62 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.00 \\ . \\ 0.00 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.38 & -0.00 & 0.51 \\ -0.00 & -0.00 & -0.00 \\ 0.51 & -0.00 & -7.62 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.00 \\ 0.00 \\ 0.00 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.76 & . & -1.02 \\ 0.00 & . & 0.00 \\ -1.02 & . & 15.24 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.66 \\ \end{bmatrix}\)
8 6 \(\begin{bmatrix}{} 1 \\ 2 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -1.12 & 3.72 & -4.71 \\ 3.72 & -20.34 & 28.94 \\ -4.71 & 28.94 & -49.88 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.00 \\ 0.00 \\ 0.00 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -1.12 & 3.72 & -4.71 \\ 3.72 & -20.34 & 28.94 \\ -4.71 & 28.94 & -49.88 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.00 \\ 0.00 \\ 0.00 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 2.23 & -7.44 & 9.42 \\ -7.44 & 40.67 & -57.87 \\ 9.42 & -57.87 & 99.76 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.41 \\ \end{bmatrix}\)
9 1 \(\begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.23 & . & -0.07 \\ . & . & . \\ -0.07 & . & -1.96 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.33 \\ . \\ 1.41 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.04 & 0.06 & 0.04 \\ 0.06 & -0.11 & -0.08 \\ 0.04 & -0.08 & -0.05 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.04 \\ -0.07 \\ -0.07 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.16 & . & -0.19 \\ -0.22 & . & 0.61 \\ -0.14 & . & 0.45 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -1.89 \\ \end{bmatrix}\)
9 8 \(\begin{bmatrix}{} 1 \\ 2 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -1.05 & 2.30 & 2.08 \\ 2.30 & -9.75 & -8.28 \\ 2.08 & -8.28 & -9.41 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.04 \\ -5.12 \\ -1.98 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.65 & 1.19 & 1.04 \\ 1.19 & -4.36 & -3.66 \\ 1.04 & -3.66 & -3.26 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.08 \\ 2.17 \\ 1.01 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.50 & -2.38 & -3.04 \\ -3.49 & 12.17 & 12.36 \\ -2.71 & 9.46 & 10.98 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -1.75 \\ \end{bmatrix}\)
7 4 \(\begin{bmatrix}{} 2 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} . & . & . \\ . & -2.58 & -2.15 \\ . & -2.15 & -3.79 \\ \end{bmatrix}\) \(\begin{bmatrix}{} . \\ -2.54 \\ -0.35 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.07 & 0.05 & 0.04 \\ 0.05 & -0.17 & -0.14 \\ 0.04 & -0.14 & -0.11 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.53 \\ -0.43 \\ -0.32 \\ \end{bmatrix}\) \(\begin{bmatrix}{} . & 0.67 & 0.11 \\ . & 0.17 & 1.31 \\ . & 0.19 & 1.10 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -1.34 \\ \end{bmatrix}\)
7 5 \(\begin{bmatrix}{} 2 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} . & . & . \\ . & -3.89 & -3.07 \\ . & -3.07 & -4.65 \\ \end{bmatrix}\) \(\begin{bmatrix}{} . \\ -2.85 \\ 0.10 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.12 & -0.00 & 0.07 \\ -0.00 & -0.89 & -0.87 \\ 0.07 & -0.87 & -0.89 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.73 \\ 0.05 \\ -0.40 \\ \end{bmatrix}\) \(\begin{bmatrix}{} . & 0.88 & -0.06 \\ . & 2.81 & 4.07 \\ . & 2.19 & 4.00 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -1.21 \\ \end{bmatrix}\)
0 9 \(\begin{bmatrix}{} 1 \\ 2 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.80 & 1.83 & 1.57 \\ 1.83 & -7.81 & -6.46 \\ 1.57 & -6.46 & -7.54 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.34 \\ -5.90 \\ -2.54 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.36 & 0.64 & 0.52 \\ 0.64 & -2.26 & -1.83 \\ 0.52 & -1.83 & -1.53 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.01 \\ 1.76 \\ 0.96 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.89 & -1.17 & -1.80 \\ -2.22 & 7.31 & 7.94 \\ -1.63 & 5.50 & 6.66 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -2.53 \\ \end{bmatrix}\)
0 7 \(\begin{bmatrix}{} 1 \\ 2 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.64 & 1.55 & 1.22 \\ 1.55 & -6.53 & -5.22 \\ 1.22 & -5.22 & -6.22 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.82 \\ -7.04 \\ -3.63 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.15 & 0.23 & 0.17 \\ 0.23 & -0.76 & -0.58 \\ 0.17 & -0.58 & -0.44 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.06 \\ 1.16 \\ 0.75 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.40 & -0.22 & -0.70 \\ -1.13 & 3.25 & 3.98 \\ -0.79 & 2.38 & 3.09 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -3.47 \\ \end{bmatrix}\)
ND 0 \(\begin{bmatrix}{} 1 \\ 2 \\ 3 \\ \end{bmatrix}\) ND ND ND ND ND ND

Step 3: Calculating \(\mathbf{L}\), \(\vec{m}\), \(r\) for the internal nodes and the root

# For tip 2 with parent node 6, we use the following terms stored in Table S5:
A2 <- matrix(-0.17)
b2 <- 0.0
C2 <- rbind(c(-0.17, 0), 
            c(0, 0))
d2 <- c(0.0, 0.0)
E2 <- matrix(c(0.35, 0), nrow = 2, ncol = 1)
f2 <- -1.45
k2 <- 1

# Now we apply Eq. S3:
print(L62 <- C2)
##       [,1] [,2]
## [1,] -0.17    0
## [2,]  0.00    0
print(m62 <- d2 + E2 %*% X[k2, "2", drop = FALSE])
##          2
## [1,] 0.035
## [2,] 0.000
print(r62 <- t(X[k2, "2", drop = FALSE]) %*% A2 %*% X[k2, "2", drop = FALSE] + 
        t(X[k2, "2", drop = FALSE]) %*% b2 + f2)
##      2
## 2 -1.5
# For tip 3 with parent node 6, applying Eq. S3, we obtain (see Table S8):
L63 <- rbind(c(-0.38, 0.51),
             c(0.51, -7.62))
m63 <- c(-1.07, 18.09)
r63 <- -11.41

# Now, we sum the terms L6i, m6i and r6i over all daughters of 6 (i) to obtain:
print(L6 <- L62 + L63)
##       [,1]  [,2]
## [1,] -0.55  0.51
## [2,]  0.51 -7.62
print(m6 <- m62 + m63)
##       2
## [1,] -1
## [2,] 18
print(r6 <- r62 + r63)
##     2
## 2 -13
# Using Eq. S3, we obtain L86, m86, r86, using the values for A,b,C,d,E,f in Table S5:
A6 <- rbind(c(-0.44, 0.58),
            c(0.58, -8.71))
b6 <- c(0.0, 0.0)
C6 <- rbind(c(-0.44, 0.58),
            c(0.58, -8.71))
d6 <- c(0.0, 0.0)
E6 <- rbind(c(0.87, -1.16),
            c(-1.16, 17.42))
f6 <- -0.52
k6 <- c(1, 3)

print(L86 <- C6 - (1/4)*E6 %*% solve(A6 + L6) %*% t(E6))
##       [,1]  [,2]
## [1,] -0.25  0.27
## [2,]  0.27 -4.06
print(m86 <- d6 - (1/2)*E6 %*% solve(A6 + L6) %*% (b6+m6))
##          2
## [1,] -0.57
## [2,]  9.65
print(r86 <- f6+r6+(length(k6)/2)*log(2*pi) - 
        (1/2)*log(det(-2*(A6+L6))) -
        (1/4)*t(b6+m6) %*% solve(A6+L6) %*% (b6+m6))
##      2
## 2 -8.6
# Because 8 is a singleton node, we immediately obtain L8, m8, r8:
L8 <- L86; m8 <- m86; r8 <- r86;

Variant 1

options(digits = 3)
cat(FormatTableAsLatex(
  traceTable1[
    list(pOrder), 
    list(
      j, i, X_i, k_i, 
      L_i, m_i, r_i,
      `L_{ji}`, `m_{ji}`, `r_{ji}`)], 
  
  type = tableOutputType))









\(j\) \(i\) \(X_i\) \(k_i\) \(L_i\) \(m_i\) \(r_i\) \(L_{ji}\) \(m_{ji}\) \(r_{ji}\)
6 2 \(\begin{bmatrix}{} 0.1 \\ NaN \\ NA \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ \end{bmatrix}\) ND ND ND \(\begin{bmatrix}{} -0.174 & . & -0.000 \\ . & . & . \\ -0.000 & . & -0.000 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.035 \\ . \\ 0.000 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -1.450 \\ \end{bmatrix}\)
6 3 \(\begin{bmatrix}{} 0.2 \\ NaN \\ 1.2 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix}\) ND ND ND \(\begin{bmatrix}{} -0.381 & . & 0.508 \\ . & . & . \\ 0.508 & . & -7.622 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -1.067 \\ . \\ 18.089 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -11.405 \\ \end{bmatrix}\)
8 6 \(\begin{bmatrix}{} NA \\ NaN \\ NA \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.555 & . & 0.508 \\ . & . & . \\ 0.508 & . & -7.622 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -1.032 \\ . \\ 18.089 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -12.855 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.243 & . & 0.271 \\ . & . & . \\ 0.271 & . & -4.065 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.568 \\ . \\ 9.648 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -8.571 \\ \end{bmatrix}\)
9 1 \(\begin{bmatrix}{} 0.3 \\ NaN \\ 1.4 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix}\) ND ND ND \(\begin{bmatrix}{} -0.037 & . & 0.040 \\ . & . & . \\ 0.040 & . & -0.053 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.249 \\ . \\ 0.520 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -3.735 \\ \end{bmatrix}\)
9 8 \(\begin{bmatrix}{} NA \\ NaN \\ NA \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.243 & . & 0.271 \\ . & . & . \\ 0.271 & . & -4.065 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.568 \\ . \\ 9.648 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -8.571 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.195 & . & 0.250 \\ . & . & . \\ 0.250 & . & -0.594 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.402 \\ . \\ 1.267 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -4.248 \\ \end{bmatrix}\)
7 4 \(\begin{bmatrix}{} NA \\ 0.2 \\ 0.2 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 2 \\ 3 \\ \end{bmatrix}\) ND ND ND \(\begin{bmatrix}{} -0.068 & 0.054 & 0.040 \\ 0.054 & -0.174 & -0.141 \\ 0.040 & -0.141 & -0.114 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.683 \\ -0.138 \\ -0.066 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -2.349 \\ \end{bmatrix}\)
7 5 \(\begin{bmatrix}{} NA \\ 1.2 \\ 0.4 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 2 \\ 3 \\ \end{bmatrix}\) ND ND ND \(\begin{bmatrix}{} -0.115 & -0.001 & 0.070 \\ -0.001 & -0.893 & -0.869 \\ 0.070 & -0.869 & -0.890 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.762 \\ 5.043 \\ 3.834 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -13.880 \\ \end{bmatrix}\)
0 9 \(\begin{bmatrix}{} NA \\ NaN \\ NA \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.232 & . & 0.291 \\ . & . & . \\ 0.291 & . & -0.647 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.651 \\ . \\ 1.786 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -7.983 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.140 & 0.162 & 0.143 \\ 0.162 & -0.200 & -0.183 \\ 0.143 & -0.183 & -0.170 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.262 \\ 0.422 \\ 0.429 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -7.430 \\ \end{bmatrix}\)
0 7 \(\begin{bmatrix}{} NA \\ NA \\ NA \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ 2 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.183 & 0.053 & 0.110 \\ 0.053 & -1.066 & -1.010 \\ 0.110 & -1.010 & -1.004 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 2.444 \\ 4.906 \\ 3.769 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -16.230 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.052 & 0.052 & 0.035 \\ 0.052 & -0.113 & -0.082 \\ 0.035 & -0.082 & -0.060 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.222 \\ -0.396 \\ -0.174 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -10.947 \\ \end{bmatrix}\)
ND 0 \(\begin{bmatrix}{} NA \\ NA \\ NA \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ 2 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.192 & 0.214 & 0.178 \\ 0.214 & -0.313 & -0.265 \\ 0.178 & -0.265 & -0.230 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.960 \\ 0.026 \\ 0.255 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -18.377 \\ \end{bmatrix}\) ND ND ND

Variant 2

cat(FormatTableAsLatex(
  traceTable2[
    list(pOrder), 
    list(
      j, i, X_i, k_i, 
      L_i, m_i, r_i,
      `L_{ji}`, `m_{ji}`, `r_{ji}`)], 
  
  type = tableOutputType))









\(j\) \(i\) \(X_i\) \(k_i\) \(L_i\) \(m_i\) \(r_i\) \(L_{ji}\) \(m_{ji}\) \(r_{ji}\)
6 2 \(\begin{bmatrix}{} 0.1 \\ NaN \\ NA \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ \end{bmatrix}\) ND ND ND \(\begin{bmatrix}{} -0.174 & . & -0.000 \\ . & . & . \\ -0.000 & . & -0.000 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.035 \\ . \\ 0.000 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -1.450 \\ \end{bmatrix}\)
6 3 \(\begin{bmatrix}{} 0.2 \\ NaN \\ 1.2 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix}\) ND ND ND \(\begin{bmatrix}{} -0.381 & . & 0.508 \\ . & . & . \\ 0.508 & . & -7.622 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -1.067 \\ . \\ 18.089 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -11.405 \\ \end{bmatrix}\)
8 6 \(\begin{bmatrix}{} NA \\ NaN \\ NA \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.555 & . & 0.508 \\ . & . & . \\ 0.508 & . & -7.622 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -1.032 \\ . \\ 18.089 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -12.855 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.243 & . & 0.271 \\ . & . & . \\ 0.271 & . & -4.065 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.568 \\ . \\ 9.648 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -8.571 \\ \end{bmatrix}\)
9 1 \(\begin{bmatrix}{} 0.3 \\ NaN \\ 1.4 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix}\) ND ND ND \(\begin{bmatrix}{} -0.037 & . & 0.040 \\ . & . & . \\ 0.040 & . & -0.053 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.249 \\ . \\ 0.520 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -3.735 \\ \end{bmatrix}\)
9 8 \(\begin{bmatrix}{} NA \\ NaN \\ NA \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.243 & . & 0.271 \\ . & . & . \\ 0.271 & . & -4.065 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.568 \\ . \\ 9.648 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -8.571 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.195 & . & 0.250 \\ . & . & . \\ 0.250 & . & -0.594 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.402 \\ . \\ 1.267 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -4.248 \\ \end{bmatrix}\)
7 4 \(\begin{bmatrix}{} NA \\ 0.2 \\ 0.2 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 2 \\ 3 \\ \end{bmatrix}\) ND ND ND \(\begin{bmatrix}{} -0.068 & 0.054 & 0.040 \\ 0.054 & -0.174 & -0.141 \\ 0.040 & -0.141 & -0.114 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.683 \\ -0.138 \\ -0.066 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -2.349 \\ \end{bmatrix}\)
7 5 \(\begin{bmatrix}{} NA \\ 1.2 \\ 0.4 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 2 \\ 3 \\ \end{bmatrix}\) ND ND ND \(\begin{bmatrix}{} -0.115 & -0.001 & 0.070 \\ -0.001 & -0.893 & -0.869 \\ 0.070 & -0.869 & -0.890 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.762 \\ 5.043 \\ 3.834 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -13.880 \\ \end{bmatrix}\)
0 9 \(\begin{bmatrix}{} NA \\ NaN \\ NA \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.232 & . & 0.291 \\ . & . & . \\ 0.291 & . & -0.647 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.651 \\ . \\ 1.786 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -7.983 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.140 & . & 0.143 \\ . & . & . \\ 0.143 & . & -0.170 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.262 \\ . \\ 0.429 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -7.430 \\ \end{bmatrix}\)
0 7 \(\begin{bmatrix}{} NA \\ NA \\ NA \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ 2 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.183 & 0.053 & 0.110 \\ 0.053 & -1.066 & -1.010 \\ 0.110 & -1.010 & -1.004 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 2.444 \\ 4.906 \\ 3.769 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -16.230 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.052 & 0.053 & 0.035 \\ 0.053 & -1.066 & -1.010 \\ 0.035 & -1.010 & -0.060 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.222 \\ 4.906 \\ -0.174 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -10.947 \\ \end{bmatrix}\)
ND 0 \(\begin{bmatrix}{} NA \\ NaN \\ NA \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.192 & . & 0.178 \\ . & . & . \\ 0.178 & . & -0.230 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.960 \\ . \\ 0.255 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -18.377 \\ \end{bmatrix}\) ND ND ND

Variant 3

cat(FormatTableAsLatex(
  traceTable3[
    list(pOrder), 
    list(
      j, i, X_i, k_i, 
      L_i, m_i, r_i, 
      `L_{ji}`, `m_{ji}`, `r_{ji}`)], 
  
  type = tableOutputType))









\(j\) \(i\) \(X_i\) \(k_i\) \(L_i\) \(m_i\) \(r_i\) \(L_{ji}\) \(m_{ji}\) \(r_{ji}\)
6 2 \(\begin{bmatrix}{} 0.1 \\ NA \\ NA \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ \end{bmatrix}\) ND ND ND \(\begin{bmatrix}{} -0.174 & -0.000 & -0.000 \\ -0.000 & -0.000 & -0.000 \\ -0.000 & -0.000 & -0.000 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.035 \\ 0.000 \\ 0.000 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -1.450 \\ \end{bmatrix}\)
6 3 \(\begin{bmatrix}{} 0.2 \\ NA \\ 1.2 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix}\) ND ND ND \(\begin{bmatrix}{} -0.381 & -0.000 & 0.508 \\ -0.000 & -0.000 & -0.000 \\ 0.508 & -0.000 & -7.622 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -1.067 \\ 0.000 \\ 18.089 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -11.405 \\ \end{bmatrix}\)
8 6 \(\begin{bmatrix}{} NA \\ NA \\ NA \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ 2 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.555 & 0.000 & 0.508 \\ 0.000 & 0.000 & 0.000 \\ 0.508 & 0.000 & -7.622 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -1.032 \\ 0.000 \\ 18.089 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -12.855 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.243 & -0.000 & 0.271 \\ -0.000 & 0.000 & -0.000 \\ 0.271 & -0.000 & -4.065 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.568 \\ -0.000 \\ 9.648 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -8.571 \\ \end{bmatrix}\)
9 1 \(\begin{bmatrix}{} 0.3 \\ NA \\ 1.4 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix}\) ND ND ND \(\begin{bmatrix}{} -0.037 & 0.059 & 0.040 \\ 0.059 & -0.109 & -0.076 \\ 0.040 & -0.076 & -0.053 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.249 \\ 0.711 \\ 0.520 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -3.735 \\ \end{bmatrix}\)
9 8 \(\begin{bmatrix}{} NA \\ NA \\ NA \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ 2 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.243 & -0.000 & 0.271 \\ -0.000 & 0.000 & -0.000 \\ 0.271 & -0.000 & -4.065 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.568 \\ -0.000 \\ 9.648 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -8.571 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.195 & 0.213 & 0.250 \\ 0.213 & -0.318 & -0.426 \\ 0.250 & -0.426 & -0.594 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.402 \\ 0.860 \\ 1.267 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -4.248 \\ \end{bmatrix}\)
7 4 \(\begin{bmatrix}{} NA \\ 0.2 \\ 0.2 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 2 \\ 3 \\ \end{bmatrix}\) ND ND ND \(\begin{bmatrix}{} -0.068 & 0.054 & 0.040 \\ 0.054 & -0.174 & -0.141 \\ 0.040 & -0.141 & -0.114 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.683 \\ -0.138 \\ -0.066 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -2.349 \\ \end{bmatrix}\)
7 5 \(\begin{bmatrix}{} NA \\ 1.2 \\ 0.4 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 2 \\ 3 \\ \end{bmatrix}\) ND ND ND \(\begin{bmatrix}{} -0.115 & -0.001 & 0.070 \\ -0.001 & -0.893 & -0.869 \\ 0.070 & -0.869 & -0.890 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.762 \\ 5.043 \\ 3.834 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -13.880 \\ \end{bmatrix}\)
0 9 \(\begin{bmatrix}{} NA \\ NA \\ NA \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ 2 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.232 & 0.272 & 0.291 \\ 0.272 & -0.427 & -0.502 \\ 0.291 & -0.502 & -0.647 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.651 \\ 1.571 \\ 1.786 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -7.983 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.089 & 0.140 & 0.106 \\ 0.140 & -0.258 & -0.204 \\ 0.106 & -0.204 & -0.163 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.394 \\ 1.044 \\ 0.833 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -8.380 \\ \end{bmatrix}\)
0 7 \(\begin{bmatrix}{} NA \\ NA \\ NA \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ 2 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.183 & 0.053 & 0.110 \\ 0.053 & -1.066 & -1.010 \\ 0.110 & -1.010 & -1.004 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 2.444 \\ 4.906 \\ 3.769 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -16.230 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.052 & 0.052 & 0.035 \\ 0.052 & -0.113 & -0.082 \\ 0.035 & -0.082 & -0.060 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1.222 \\ -0.396 \\ -0.174 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -10.947 \\ \end{bmatrix}\)
ND 0 \(\begin{bmatrix}{} NA \\ NA \\ NA \\ \end{bmatrix}\) \(\begin{bmatrix}{} 1 \\ 2 \\ 3 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -0.141 & 0.192 & 0.141 \\ 0.192 & -0.372 & -0.286 \\ 0.141 & -0.286 & -0.223 \\ \end{bmatrix}\) \(\begin{bmatrix}{} 0.827 \\ 0.648 \\ 0.659 \\ \end{bmatrix}\) \(\begin{bmatrix}{} -19.328 \\ \end{bmatrix}\) ND ND ND

Step 4: Calculating the log-likelihood value using \(\mathbf{L}_{0}\), \(\vec{m}_{0}\), and \(r_{0}\).

# Variant 1.
# Copy the values of L0, m0 and r0 from Table S8:
L0 <- rbind(c(-0.192, 0.214, 0.178),
            c(0.214, -0.313, -0.265),
            c(0.178, -0.265, -0.230))
m0 <- c(0.96, 0.026, 0.255)
r0 <- -18.377

# Use Eq. S2 to estimate the optimal X0:
print(t(x0Hat <- -0.5*solve(L0) %*% m0))
##      [,1]  [,2] [,3]
## [1,] 9.64 -6.25 15.2
# Use Eq. S1 to calculate the log-likelihood:
print(ll0 <- t(x0Hat) %*% L0 %*% x0Hat + t(x0Hat) %*% m0 + r0)
##       [,1]
## [1,] -11.9
# Variant 2.
# Copy the values of L0, m0 and r0 from Table S9:
L0 <- rbind(c(-0.192, 0.178),
            c(0.178, -0.230))
m0 <- c(0.96, 0.255)
r0 <- -18.377

# Use Eq. S2 to estimate the optimal X0:
print(t(x0Hat <- -0.5*solve(L0) %*% m0))
##      [,1] [,2]
## [1,] 10.7 8.81
# Use Eq. S1 to calculate the log-likelihood:
print(ll0 <- t(x0Hat) %*% L0 %*% x0Hat + t(x0Hat) %*% m0 + r0)
##       [,1]
## [1,] -12.1
# Variant 3.
# The function PCMLikTrace generates a data.table with the values of 
# omega, Phi, V, A, b, C, d, E, f, L, m, r. 
traceTable3 <- PCMLikTrace(X3[, tree$tip.label], tree, model.OU.BM)
# The column i corresponds to the node label in the tree as depicted on Fig. 1:
setkey(traceTable3, i)

options(digits = 4)
# Variant 3.
# Copy the values of L0, m0 and r0 from the traceTable object (these values have
# the maximal double floating point precision):
print(L0 <- traceTable3[list("0")][["L_i"]][[1]])
##         [,1]    [,2]    [,3]
## [1,] -0.1408  0.1919  0.1407
## [2,]  0.1919 -0.3715 -0.2862
## [3,]  0.1407 -0.2862 -0.2233
print(m0 <- traceTable3[list("0")][["m_i"]][[1]])
## [1] 0.8273 0.6484 0.6590
print(r0 <- traceTable3[list("0")][["r_i"]][[1]])
## [1] -19.33
# Notice the exact match with the values for variant 3 reported in Fig. S5:
print(t(x0Hat <- -0.5*solve(L0) %*% m0))
##       [,1]  [,2]   [,3]
## [1,] 15.99 18.34 -11.95
print(ll0 <- t(x0Hat) %*% L0 %*% x0Hat + t(x0Hat) %*% m0 + r0)
##        [,1]
## [1,] -10.71

References

Mitov, Venelin, Krzysztof Bartoszek, Georgios Asimomitis, and Tanja Stadler. 2019. “Fast likelihood calculation for multivariate Gaussian phylogenetic models with shifts.” Theor. Popul. Biol., December. https://doi.org/10.1016/j.tpb.2019.11.005.

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.