search

BEAST2-Dev / bdsky

Birth Death Serial Skyline Model for BEAST2
GNU General Public License v3.0
2 stars 11 forks source link

logP being positive infinity #11

Closed zhangchicool closed 7 years ago

zhangchicool commented 7 years ago

In rare cases, (say psi is very large, so Ai is very large too, then g() is zero [line 765]), the tree log likelihood becomes positive infinity.
It checks Double.isInfinite(logP)) then returns logP, then logAlpha in MCMC.java is positive infinity too, casing the move being accepted and ending up with prior being Infinity. Should it be if (Double.isInfinite(logP)) return Double.NEGATIVE_INFINITY; so that the move can be rejected properly?

zhangchicool commented 7 years ago

Update: from the conversation at https://github.com/CompEvol/beast2/issues/632, it need to be dealt with in bdsky. One way is to improve numerical instability, another may be if (Double.isInfinite(logP)) return Double.NEGATIVE_INFINITY ?

zhangchicool commented 7 years ago

I found that log( g(int index, double ti, double t) ) is to blame. It first calculates g() then takes the log of g(). If g() is zero due to numerical rounding, then - log( g() ) is infinite. I've changed it to calculate log_g directly, and it solves the positive infinite logP. The unit tests are passed, and I've sent a pull request.