Nothing
## sim.coalescent.R (2015-06-29)
## Coalescent Simulation and Visualisation
## Copyright 2014-2015 Emmanuel Paradis
## This file is part of the R-package `coalescentMCMC'.
## See the file ../COPYING for licensing issues.
sim.coalescent <-
function(n = 5, TIME = 50, growth.rate = NULL, N.0 = 50, N.final = 20,
col.lin = "grey", col.coal = "blue", pch = NULL, ...)
{
if (n > N.final) stop("sample size smaller than final population size")
if (is.null(TIME)) TIME <- log(N.final/N.0) / growth.rate
if (is.null(growth.rate)) growth.rate <- log(N.final/N.0) / TIME
if (is.null(N.0)) N.0 <- N.final/exp(growth.rate * TIME)
if (is.null(N.final)) N.final <- N.0 * exp(growth.rate * TIME)
N <- round(N.0 * exp(growth.rate * 1:TIME))
x <- numeric(0)
for (j in 1:TIME) x <- c(x, (-N[j]/2 + 0.5):(N[j]/2 - 0.5))
y <- rep(1:TIME, N)
x1 <- x2 <- y1 <- y2 <- numeric(0)
for (j in TIME:2) {
anc <- sort(sample(1:N[j - 1], N[j], replace = TRUE))
x1 <- c(x1, x[which(y == j)])
y1 <- c(y1, rep(j, N[j]))
x2 <- c(x2, x[which(y == j - 1)][anc])
y2 <- c(y2, rep(j - 1, N[j]))
}
plot(x, y, type = "n", xlab = "", ylab = "", frame = FALSE,
xaxt = "n", yaxt = "n")
segments(x1, y1, x2, y2, col = col.lin)
dat <- sample(x[which(y == TIME)], n)
for (j in TIME:2) {
for (i in seq_len(n)) {
index <- which(x1 == dat[i] & y1 == j)
segments(x1[index], j, x2[index], j - 1, lwd = 2, col = col.coal)
dat[i] <- x2[index]
}
}
if (!is.null(pch))
points(x, y, pch = pch, ...)
}
proba.coalescent <- function(t, N = 1e4, n = 2, exact = TRUE)
{
n <- round(n)
if (exact) {
q <- prod(1 - 1:(n - 1)/N)
p <- 1 - q
} else {
p <- (n * (n - 1)/2)/N
q <- 1 - p
}
t <- round(t)
z <- t <= 0
res <- cbind(t, ifelse(z, 0, p * q^(t - 1)))
colnames(res) <- c("t (generations)", "Pr(coal)")
res
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.