sample.basic.models.joint {CRABS}

Jointly samples speciation and extinction trajectories through time, with noise.

Description

Jointly samples speciation and extinction trajectories through time, with noise.

Usage

sample.basic.models.joint(
  times,
  p.delta,
  lambda0,
  mu0 = NULL,
  MRF.type = "HSMRF",
  beta.param = c(0.3, 0.3),
  mu0.median = 0.1,
  mu0.logsd = 1.17481,
  mrf.sd.scale = 1,
  min.lambda = 0,
  min.mu = 0,
  max.lambda = 10,
  max.mu = 10,
  min.p = -0.05,
  max.p = 1.05
)

Arguments

times

the time knots

p.delta

The The pulled diversification rate function (measured in time before present).

lambda0

The speciation rate at present.

mu0

The extinction rate at present, otherwise drawn randomly.

MRF.type

"HSMRF" or "GMRF", type for stochastic noise.

beta.param

Parameters of the Beta distribution used for

mu0.median

When not specified, extinction rate at present is drawn from a lognormal distribution with this median.

mu0.logsd

When not specified, extinction rate at present is drawn from a lognormal distribution with this sd

mrf.sd.scale

scale the sd of the mrf process up or down. defaults to 1.0

min.lambda

The minimum speciation rate (rescaling done after after drawing rates).

min.mu

The minimum extinction rate (rescaling done after after drawing rates).

max.lambda

The maximum speciation rate (rescaling done after after drawing rates).

max.mu

The maximum extinction rate (rescaling done after after drawing rates).

min.p

The lower bound of parameter p's trajectory.

max.p

The upper bound of parameter p's trajectory.

Value

Speciation or extinction rate at a number of timepoints.

Examples

data("primates_ebd")

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 = 40, 
                                         rate.type = "joint", 
                                         sample.joint.rates = sample.joint.rates)

joint.samples