Test for condition on root and on rho

scripts/LikelihoodForConditionOnRootAndRhoTest.R

@@ -0,0 +1,27 @@
+#The formula for calculating tree likelihood is equation (4) in Gavryushkina et al (2014) multiplied by (m+M+k+K)!*q1(x1)/(lambda1*(1-p0hat(x1))^2)
+lambda1<-1.5
+lambda2<-0.5
+mu1<-0.4
+mu2<-0.3
+psi1<-0.3
+psi2<-0.8
+rho<-0.9
+A2<-sqrt((lambda2-mu2-psi2)^2+4*lambda2*psi2)
+B2<-((1-2*(1-rho))*lambda2+mu2+psi2)/A2
+p2from1<-(lambda2+mu2+psi2-A2*(exp(A2)*(1+B2)-(1-B2))/(exp(A2)*(1+B2)+(1-B2)))/(2*lambda2)
+p2from0point8<-(lambda2+mu2+psi2-A2*(exp(0.8*A2)*(1+B2)-(1-B2))/(exp(0.8*A2)*(1+B2)+(1-B2)))/(2*lambda2)
+q2from0point8<-4*exp(0.8*A2)/(exp(0.8*A2)*(1+B2)+(1-B2))^2
+A1<-sqrt((lambda1-mu1-psi1)^2+4*lambda1*psi1)
+B1<-((1-2*p2from1)*lambda1+mu1+psi1)/A1
+q1from2point5<-4*exp(1.5*A1)/(exp(1.5*A1)*(1+B1)+(1-B1))^2
+q2from1<-4*exp(A2)/(exp(A2)*(1+B2)+(1-B2))^2
+A2hat<-sqrt((lambda2-mu2)^2)
+B2hat<-((1-2*(1-rho))*lambda2+mu2)/A2hat
+A1hat<-sqrt((lambda1-mu1)^2)
+p2hatfrom1<-(lambda2+mu2-A2hat*(exp(A2hat)*(1+B2hat)-(1-B2hat))/(exp(A2hat)*(1+B2hat)+(1-B2hat)))/(2*lambda2)
+B1hat<-((1-2*p2hatfrom1)*lambda1+mu1)/A1hat
+p1hatfrom2point5<-(lambda1+mu1-A1hat*(exp(1.5*A1hat)*(1+B1hat)-(1-B1hat))/(exp(1.5*A1hat)*(1+B1hat)+(1-B1hat)))/(2*lambda1)
+logFfromT<-log(2*psi1*psi2*rho*(q1from2point5)^2*p2from0point8) -2*log(1-p1hatfrom2point5)-log(q2from0point8)+2*log(q2from1)
+print(logFfromT, digits=15)
+
+

src/test/beast/evolution/speciation/SABirthDeathSkylineModel.java

@@ -136,6 +136,33 @@ public void testLikelihoodCalculationThreeIntervalsWithRhoConditionOnRhoSampling
         assertEquals(-8.80702177958906, model.calculateTreeLogLikelihood(tree2), 1e-14);
     }
 
+    @Test
+    public void testLikelihoodCalculationTwoIntervalsWithRhoConditionOnRhoSamplingAndRoot() throws Exception {
+
+        ArrayList<String> taxa1 = new ArrayList<String>(Arrays.asList("1", "2", "3"));
+        Tree tree1 = new TreeParser(taxa1, "((1:1.5,3:0.0):1.0,2:1.7):0.0", 1, false);
+
+        BirthDeathSkylineModel model = new BirthDeathSkylineModel();
+        model.setInputValue("tree", tree1);
+        model.setInputValue("birthRate", new RealParameter("1.5 0.5"));
+        model.setInputValue("deathRate", new RealParameter("0.4 0.3"));
+        model.setInputValue("samplingRate", new RealParameter("0.3 0.8"));
+        model.setInputValue("rho", new RealParameter("0.9"));
+        model.setInputValue("removalProbability", new RealParameter("0.0"));
+        model.setInputValue("reverseTimeArrays", "true true true true true");
+        model.setInputValue("birthRateChangeTimes", new RealParameter("1. 0."));
+        model.setInputValue("deathRateChangeTimes", new RealParameter("1. 0."));
+        model.setInputValue("samplingRateChangeTimes", new RealParameter("1. 0."));
+        model.setInputValue("rhoSamplingTimes", new RealParameter("0.0"));
+        model.setInputValue("conditionOnRoot", true);
+        model.setInputValue("conditionOnSurvival", false);
+        model.setInputValue("conditionOnRhoSampling", true);
+        model.initAndValidate();
+
+        // the true value is calculated in R 'scripts/LikelihoodForConditionOnRootAndRhoTest.R'
+        assertEquals(-8.55232499940742, model.calculateTreeLogLikelihood(tree1), 1e-14);
+    }
+
     @Test
     public void testLikelihoodCalculationDiversifiedSampling() throws Exception {