X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=R%2FrTrait.R;h=b87dd6f8a60ff023271acfb1339c36fdc53d190e;hb=bfaeca35ec130110327a3bb7a1f0fe3b66076a95;hp=8cc48e17c9e45296c082b33ca33f29e01c3f3296;hpb=06c3113db74a7cfa54c15a6f18163cd9b2c1f6db;p=ape.git diff --git a/R/rTrait.R b/R/rTrait.R index 8cc48e1..b87dd6f 100644 --- a/R/rTrait.R +++ b/R/rTrait.R @@ -1,4 +1,4 @@ -## rTrait.R (2010-05-06) +## rTrait.R (2010-12-24) ## Trait Evolution @@ -10,7 +10,7 @@ rTraitDisc <- function(phy, model = "ER", k = if (is.matrix(model)) ncol(model) else 2, rate = 0.1, states = LETTERS[1:k], freq = rep(1/k, k), - ancestor = FALSE, root.value = 1) + ancestor = FALSE, root.value = 1, ...) { if (is.null(phy$edge.length)) stop("tree has no branch length") @@ -18,7 +18,7 @@ rTraitDisc <- stop("at least one branch length negative") if (is.character(model)) { - switch(model, "ER" = { + switch(toupper(model), "ER" = { if (length(rate) != 1) stop("`rate' must have one element") Q <- matrix(rate, k, k) @@ -56,7 +56,7 @@ rTraitDisc <- if (is.function(model)) { environment(model) <- environment() # to find 'k' - for (i in N:1) x[des[i]] <- model(x[anc[i]], el[i]) + for (i in N:1) x[des[i]] <- model(x[anc[i]], el[i], ...) } else { freq <- rep(freq, each = k) Q <- Q * freq @@ -81,8 +81,8 @@ rTraitDisc <- } rTraitCont <- - function(phy, model = "BM", sigma = 0.1, alpha = 1, - theta = 0, ancestor = FALSE, root.value = 0) + function(phy, model = "BM", sigma = 0.1, alpha = 1, theta = 0, + ancestor = FALSE, root.value = 0, linear = TRUE, ...) { if (is.null(phy$edge.length)) stop("tree has no branch length") @@ -102,9 +102,9 @@ rTraitCont <- if (is.function(model)) { environment(model) <- environment() - for (i in N:1) x[des[i]] <- model(x[anc[i]], el[i]) + for (i in N:1) x[des[i]] <- model(x[anc[i]], el[i], ...) } else { - model <- pmatch(model, c("BM", "OU")) + model <- pmatch(toupper(model), c("BM", "OU")) if (length(sigma) == 1) sigma <- rep(sigma, N) else if (length(sigma) != N) stop("'sigma' must have one or Nedge(phy) elements") @@ -115,6 +115,7 @@ rTraitCont <- if (length(theta) == 1) theta <- rep(theta, N) else if (length(theta) != N) stop("'theta' must have one or Nedge(phy) elements") + if (!linear) model <- model + 1L } .C("rTraitCont", as.integer(model), as.integer(N), as.integer(anc - 1L), as.integer(des - 1L), el, sigma, alpha, theta, x, DUP = FALSE, NAOK = TRUE, PACKAGE = "ape") }