diff --git a/README.md b/README.md
index be19bca..0b14931 100644
--- a/README.md
+++ b/README.md
@@ -4,6 +4,12 @@
 
 ## Citation
 
-- Birth-death serial skyline: [stadler2013](http://www.pnas.org/content/110/1/228.full)	
-- Birth-death serial sampling: [stadler2013](http://www.pnas.org/content/110/1/228.full)
-- Birth-death incomplete sampling (no ψ): [stadler2009](http://www.ncbi.nlm.nih.gov/pubmed/19631666)
+Please cite [Stadler2013](http://www.pnas.org/content/110/1/228.full) for the following models: 
+
+- Birth-death skyline serial; 	
+- Birth-death skyline contemporary;
+- Birth-death skyline contemporary BDS (birth, death, sampling rate);
+- Birth-death skyline multi-rho.
+
+If sampled ancestors are used, please also cite: 
+[Gavryushkina2014](https://doi.org/10.1371/journal.pcbi.1003919) 
diff --git a/scripts/LikelihoodForConditionOnRootAndRhoTest.R b/scripts/LikelihoodForConditionOnRootAndRhoTest.R
new file mode 100644
index 0000000..05978b8
--- /dev/null
+++ b/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)
+
+              
diff --git a/src/beast/evolution/speciation/BirthDeathSkylineModel.java b/src/beast/evolution/speciation/BirthDeathSkylineModel.java
index ccea9d2..b4e49aa 100644
--- a/src/beast/evolution/speciation/BirthDeathSkylineModel.java
+++ b/src/beast/evolution/speciation/BirthDeathSkylineModel.java
@@ -1017,7 +1017,7 @@ public double calculateTreeLogLikelihood(TreeInterface tree) {
                 if (temp == 1)
                     return Double.NEGATIVE_INFINITY;
                 if (conditionOnRootInput.get()) {
-                    temp = log_q(index, times[index], x0) - 2 * Math.log(1 - temp);
+                    temp = log_q(index, times[index], x0) - 2 * Math.log(1 - temp)- Math.log(birth[index]);
                 } else {
                     temp = log_q(index, times[index], x0) - Math.log(1 - temp);
                 }
diff --git a/src/test/beast/evolution/speciation/SABirthDeathSkylineModel.java b/src/test/beast/evolution/speciation/SABirthDeathSkylineModel.java
index 4ba8f6c..0c5117e 100644
--- a/src/test/beast/evolution/speciation/SABirthDeathSkylineModel.java
+++ b/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 {
 
diff --git a/version.xml b/version.xml
index e969a73..bacf970 100644
--- a/version.xml
+++ b/version.xml
@@ -1,3 +1,3 @@
-<addon name='BDSKY' version='1.4.2'>
+<addon name='BDSKY' version='1.4.3'>
 	<depends on='beast2' atleast='2.5.1'/>
 </addon>