sample.congruence.class {CRABS}

Stochastic exploration of congruent models.

Description

Stochastic exploration of congruent models.

Usage

sample.congruence.class(
  model,
  num.samples,
  rate.type = "both",
  sample.speciation.rates = NULL,
  sample.extinction.rates = NULL,
  sample.joint.rates = NULL
)

Arguments

model

the reference model, an object of class "CRABS"

num.samples

The number of samples to be drawn

rate.type

either "extinction", "speciation", "both" or "joint"

sample.speciation.rates

a function that when called returns a speciation rate function

sample.extinction.rates

a function that when called returns a extinction rate function

sample.joint.rates

a function that when called returns a list with a speciation rate function and an extinction rate function

Value

A named list with congruent rates.

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, primates_ebd[["time"]])

# Sampling extinction rates

extinction_rate_samples <- function(){
   res <- sample.basic.models(times = times, 
                              rate0 = 0.05, 
                              model = "MRF", 
                              MRF.type = "HSMRF", 
                              fc.mean = 2.0, 
                              min.rate = 0.0, 
                              max.rate = 1.0)
   return(res)
} 

samples <- sample.congruence.class(model, 
                                   num.samples = 8,
                                   rate.type = "extinction",
                                   sample.extinction.rates = extinction_rate_samples)

samples

# Jointly sampling speciation and extinction rates

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