data("primates_ebd")
set.seed(123)
l <- approxfun(primates_ebd[["time"]], primates_ebd[["lambda"]])
mu <- approxfun(primates_ebd[["time"]], primates_ebd[["mu"]])
times <- primates_ebd[["time"]]
model <- create.model(l, mu, times)
sample.joint.rates <- function(n) {
sample.basic.models.joint(times = times,
p.delta = model$p.delta,
beta.param = c(0.5,0.3),
lambda0 = l(0.0),
mu0.median = mu(0.0))
}
joint.samples <- sample.congruence.class(model = model,
num.samples = 100,
rate.type = "joint",
sample.joint.rates = sample.joint.rates)
full.plot.regularity.thresholds(joint.samples)
Run the code above in your browser using DataLab