CRABS (version 1.2.0)

sample.congruence.class: 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
)

Value

A named list with congruent rates.

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

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