new image.DNAbin()

[ape.git] / R / cherry.R

diff --git a/R/cherry.R b/R/cherry.R
index 7b952b87e74cbaa06b55ebe97ad6224c8f4719c0..347465cfa094a085bdbbbce7a536e9135ea957e7 100644 (file)
--- a/R/cherry.R
+++ b/R/cherry.R
@@ -1,18 +1,18 @@
-## cherry.R (2006-10-03)
+## cherry.R (2009-05-10)
 
 ##     Number of Cherries and Null Models of Trees
 
-## Copyright 2002-2006 Emmanuel Paradis
+## Copyright 2002-2009 Emmanuel Paradis
 
 ## This file is part of the R-package `ape'.
 ## See the file ../COPYING for licensing issues.
 
 cherry <- function(phy)
 {
-    if (class(phy) != "phylo") stop("object \"phy\" is not of class \"phylo\"")
+    if (!inherits(phy, "phylo")) stop('object "phy" is not of class "phylo"')
     n <- length(phy$tip.label)
     nb.node <- phy$Nnode
-    if (nb.node != n - 1) stop("\"phy\" is not fully dichotomous")
+    if (nb.node != n - 1) stop('"phy" is not fully dichotomous')
     if (n < 4) stop("not enough tips in your phylogeny for this analysis")
     cherry <- sum(tabulate(phy$edge[, 1][phy$edge[, 2] <= n]) == 2)
     small.n <- n < 20
@@ -39,16 +39,14 @@ cherry <- function(phy)
 
 f.cherry.yule <- function(n, k)
 {
-    P <- if (k == 0 || k > floor(n/2)) 0 else if (n == 4) if (k == 1) 2/3 else if (k == 2) 1/3 else 0
-             else (1 - 2*(k - 1)/(n - 1))*f.cherry.yule(n - 1, k - 1) +
-              2*k/(n - 1)*f.cherry.yule(n - 1, k)
-    P
+    if (k == 0 || k > floor(n/2)) 0 else if (n == 4) if (k == 1) 2/3 else if (k == 2) 1/3 else 0
+    else (1 - 2*(k - 1)/(n - 1))*f.cherry.yule(n - 1, k - 1) +
+        2*k/(n - 1)*f.cherry.yule(n - 1, k)
 }
 
 f.cherry.uniform <- function(n, k)
 {
-    P <- if (k == 0 || k > floor(n/2)) 0 else if (n == 4) if (k == 1) 4/5 else if (k == 2) 1/5 else 0
-        else if (k == 1) 0 else (gamma(n + 1)*gamma(n - 2 + 1)*gamma(n - 4 + 1) * 2^(n-2*k)) /
-            (gamma(n - 2*k + 1)*gamma(2*n - 4 + 1)*gamma(k + 1)*gamma(k - 2 + 1))
-    P
+    if (k == 0 || k > floor(n/2)) 0 else if (n == 4) if (k == 1) 4/5 else if (k == 2) 1/5 else 0
+    else if (k == 1) 0 else (gamma(n + 1)*gamma(n - 2 + 1)*gamma(n - 4 + 1) * 2^(n-2*k)) /
+        (gamma(n - 2*k + 1)*gamma(2*n - 4 + 1)*gamma(k + 1)*gamma(k - 2 + 1))
 }