add test for ML estimate of bf

KlausVigo · Dec 3, 2024 · 4aef7fc · 4aef7fc
1 parent 1704f7a
commit 4aef7fc
Showing 1 changed file with 13 additions and 0 deletions.
diff --git a/inst/tinytest/test_pml.R b/inst/tinytest/test_pml.R
@@ -67,6 +67,7 @@ expect_error(pml_bb(dat, model="GTR", method="tipdated"))
     weights <- as.vector(1000 * exp(fit_T$siteLik))
     dat_tmp <- dat
     attr(dat_tmp, "weight") <- weights
+    # empirical
     fit0 <- pml(treeU1, dat_tmp)
     fit.bf <- optim.pml(fit0, optEdge=FALSE, optBf = TRUE,
                         control = pml.control(epsilon=1e-10, trace=0))
@@ -75,6 +76,18 @@ expect_error(pml_bb(dat, model="GTR", method="tipdated"))
     fit.F81 <- pml_bb(dat_tmp, model="F81", control = pml.control(trace=0))
     expect_equal(logLik(fit.F81), logLik(pml(treeU1, dat_tmp, bf=bf)))
     expect_equal(bf, fit.F81$bf, tolerance=5e-4)
+    #estimate
+    fit0 <- pml(treeU1, dat_tmp)
+    fit.bf <- optim.pml(fit0, optEdge=FALSE, optBf = TRUE,
+                        control = pml.control(epsilon=1e-10, trace=0,
+                                              statefreq="estimate"))
+    expect_equal(logLik(fit.bf), logLik(pml(treeU1, dat_tmp, bf=bf)))
+    expect_equal(bf, fit.bf$bf, tolerance=5e-4)
+    fit.F81 <- pml_bb(dat_tmp, model="F81",
+                      control = pml.control(trace=0, statefreq="estimate"))
+    expect_equal(logLik(fit.F81), logLik(pml(treeU1, dat_tmp, bf=bf)))
+    expect_equal(bf, fit.F81$bf, tolerance=5e-4)
+
 # test Q optimisation
     Q <- c(6:1)
     fit_T <- pml(treeU1, dat, Q=Q)