giappo / mbd Goto Github PK
View Code? Open in Web Editor NEWMacro Evolutionary Likelihood models: Multiple Birth Death model.
License: GNU General Public License v3.0
Macro Evolutionary Likelihood models: Multiple Birth Death model.
License: GNU General Public License v3.0
The conditional probability should be something like:
1 - exp(-nu * deltaT)
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
A declarative, efficient, and flexible JavaScript library for building user interfaces.
๐ Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. ๐๐๐
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google โค๏ธ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.