Hi all, Just a note on this topic, which Rich may be able to help clarify if he sees this. It can be shown that Rich's equation which uses log(abs(...)) is exactly equivalent to the regular version of the Nee equation. However, for an ultrametric tree, the ML of lambda will always be greater than mu. As I understand it, the advantage of the approach Diversitree uses is primarily in the context of using a mcmc chain to estimate the parameters
myLikelihoodFunction <- make.bd(myTree, sampling.f = 1.0) p <- starting.point.bd(myTree) pars_init <- c(p[1], p[2]) myprior <- make.prior.exponential(1 / (2 * (p[1] - p[2]))) runmcmc <- mcmc(myLikelihoodFunction, x.init=pars_init, w=rep(1,2), prior=myprior) Allowing the chain to sample from parameter space where lambda < mu allows for better mixing of the mcmc. cheers, matt On Thu, Apr 19, 2012 at 9:29 PM, Emmanuel Paradis <emmanuel.para...@ird.fr>wrote: > Hi Simon, > > This is fixed and you can find the new version on ape's SVN. Another thing > I fixed: it could happen that the calculation of the Hessian-based SEs > fails: the error is now caught and the MLEs are returned with SEs set to NA > (the profile-likelihood CIs are still computed of course). > > Matt is right that the distribution of branching times is "weird" with a > coalescent tree, but an answer is still expected from birthdeath(). > > Regarding the transformation on mu and lambda, the present version imposes > r (= lambda - mu) < 0 and/or a (= mu/lambda) > 1 because the likelihood > computation uses log(r) and log(1 - a). This also avoids most warnings now. > diversitree uses log(abs(... which, at first sight, is strange to me. But > there may be some justification for this. Rich may comment on it? > > Best, > > Emmanuel > > Simon Greenhill wrote on 20/04/2012 06:52: > > Morning all, >> >> Whenever I try to use ape's (v 3.0-2) birthdeath function e.g. >> >> birthdeath(rcoal(sample(10:**100, 1))) >>> >> >> I get the following error: >> >> Error in while (foo(up)< 0) up<- up + inc : >> missing value where TRUE/FALSE needed >> >> The warnings all appear to be this sort of thing: >> >> 1: In log(1 - a) : NaNs produced >> 2: In log(exp(r * x[2:N]) - a) : NaNs produced >> 3: In nlm(function(p) dev(p[1], p[2]), c(0.1, 0.2), hessian = TRUE) : >> NA/Inf replaced by maximum positive value >> >> and seem to be generated from this line: >> out<- nlm(function(p) dev(p[1], p[2]), c(0.1, 0.2), hessian = TRUE) >> >> >> Is there a fix/workaround for this? >> >> Kind regards, >> >> Simon >> >> ______________________________**_________________ >> R-sig-phylo mailing list >> R-sig-phylo@r-project.org >> https://stat.ethz.ch/mailman/**listinfo/r-sig-phylo<https://stat.ethz.ch/mailman/listinfo/r-sig-phylo> >> >> > -- > Emmanuel Paradis > IRD, Jakarta, Indonesia > http://ape.mpl.ird.fr/ > > > ______________________________**_________________ > R-sig-phylo mailing list > R-sig-phylo@r-project.org > https://stat.ethz.ch/mailman/**listinfo/r-sig-phylo<https://stat.ethz.ch/mailman/listinfo/r-sig-phylo> > [[alternative HTML version deleted]] _______________________________________________ R-sig-phylo mailing list R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
