giappo / mbd Goto Github PK

The conditional probability should be something like:
1 - exp(-nu * deltaT)

Hi @Giappo,

I wanted to use a reference to a manuscript that uses MDB (?razzo), so I looked for a CITATION file in the inst folder. There is none yet.

Therefore, I'd like to request a CITATION file (tip: peek at babette)

Would be great!

Sure I'll help 👍

Once we have decided which conditioning probability should be used, we should FORTRANize that computation.

(note to self: PMB = Pure Multiple Birth)

Add this as a test...

For some parameter choices (e.g. very high numbers, see test-mbd_ml:R147 (pars = c(60, 50, 100, 0.5)). The integration yields negative probabilities. This is best exposed considering the result of the conditioning probability, calculated with "calculate_conditional_prob".

For example:
function_names <- sls_get_function_names(
models = models
)
model_names <- sls_get_model_names(
function_names = function_names,
verbose = verbose
)

Sometimes when running mbd_loglik (and, therefore, mbd_ml) some values of C and D are less than zero. This is exposed by flags put in the code. It probably comes from negative results in the integration process.
Depends on #69

The function prob_cond is still to slow and not accurate. It must be improved!

There are several ways to do that:
1: Start multiple chains from different points and choose the ending results with the greatest likelihood;
2: Use a kind of simulated annealing to thermally move when a local maximum is reached;
3: ...

This code

mbd::mbd_loglik(
  pars = c(0.2, 0.15, 2.0, 0.1),
  brts = c(1, 2, 2, 2, 3),
  soc = 2, # Crown age
  cond = 1  # Non-extinction
)

produces the warning

number of items to replace is not a multiple of replacement length

This should not be.

I put a sketch of a bug-report here, will improve later:

richel@sonic:~/data$ ./do_all 

razzo_project_20190815/
./scripts/10_create_n_mb_species_file.sh: line 20: module: command not found
Error: sum(time_intervals) == max(abs(brts)) is not TRUE
Execution halted

...

razzo_project_20190911_unfinished/
./scripts/10_create_n_mb_species_file.sh: line 20: module: command not found
Error: sum(time_intervals) == max(abs(brts)) is not TRUE
Execution halted

Verify the assert is in mbd:

richel@sonic:~/GitHubs$ egrep -R "sum.time_intervals. == max.abs.brts.." --include=*.R
mbd/R/mbd_loglik_utils.R:  testit::assert(sum(time_intervals) == max(abs(brts)))

I predict the assert is too strict and adding a little tolerance would fix this.

But again, this is just a (temporary) sketch

The maximum likelihood estimation of the MBD package fails on (Single)BD trees.

# Simulate a BD tree
set.seed(10)
lambda <- 0.3
mu <- 0.1
phylogeny <- mbd_sim_checked(
  mbd_params = create_mbd_params(lambda = lambda, mu = mu, nu = 0.0, q = 0.0),
  crown_age = 2,
  conditioned_on = "non_extinction"
)$tes

# Maximum likelihood of BD tree
ml_est <- mbd_calc_max_lik(
  branching_times = ape::branching.times(phylogeny),
  init_param_values = create_mbd_params(
    lambda = lambda,
    mu = mu,
    nu = 0.0,
    q = 0.0
  ),
  estimated_params = create_mbd_params_selector(lambda = TRUE, mu = TRUE),
  fixed_params = create_mbd_params_selector(nu = TRUE, q = TRUE),
  init_n_species = 2,
  n_missing_species = 0,
  conditioned_on = "non_extinction"
)
expect_true(ml_est$lambda >= 0)

View this YouTube video or else you can send me the code, and I will put it here. Awesome plan?

@Giappo : mbd uses the packages 'hyperA' and 'rsexpm'. These are not CRAN packages. Where did you get these from?

I suggest to rename the package to 'mbd' (from 'MBD') for multiple reasons: (1) it follows the Hadley Wickham standard (e.g. ggplot2, tidyr, etc), (2) @lintr-bot will love us for that,

@Giappo: agree?

If yes, I will do a replace all.

You can do it changing the method deSolve::ode uses to do that.
Methods that are worth a try are "ode45", "lsodes", ...

calculate_conditional_prob yields prob > 1 for this setting:
prob <- calculate_conditional_prob(
brts = c(1, 0.5),
pars = c(0.359969011, 0.001163674, 2.051024209, 0.101152488) ,
cond = 1,
lx = 300,
n_0 = 2,
tips_interval = c(n_0 * (cond > 0), Inf),
methode = "lsodes",
abstol = 1e-16,
reltol = 1e-10
)
prob is 1.0022