Java tutorial
/* * BirthDeathGernhard08Model.java * * Copyright (c) 2002-2015 Alexei Drummond, Andrew Rambaut and Marc Suchard * * This file is part of BEAST. * See the NOTICE file distributed with this work for additional * information regarding copyright ownership and licensing. * * BEAST is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as * published by the Free Software Foundation; either version 2 * of the License, or (at your option) any later version. * * BEAST is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU Lesser General Public License for more details. * * You should have received a copy of the GNU Lesser General Public * License along with BEAST; if not, write to the * Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, * Boston, MA 02110-1301 USA */ package dr.evomodel.speciation; import dr.evolution.tree.NodeRef; import dr.evolution.tree.Tree; import dr.evomodelxml.speciation.BirthDeathModelParser; import dr.inference.model.Parameter; import dr.util.Author; import dr.util.Citable; import dr.util.Citation; import dr.util.CommonCitations; import java.util.ArrayList; import java.util.Collections; import java.util.List; import static org.apache.commons.math.special.Gamma.logGamma; /** * Birth Death model based on Gernhard 2008 "The conditioned reconstructed process" * Journal of Theoretical Biology Volume 253, Issue 4, 21 August 2008, Pages 769-778 * doi:10.1016/j.jtbi.2008.04.005 (http://dx.doi.org/10.1016/j.jtbi.2008.04.005) * <p/> * This derivation conditions directly on fixed N taxa. * <p/> * The inference is directly on b-d (strictly positive) and d/b (constrained in [0,1)) * <p/> * Vefified using simulated trees generated by Klass tree sample. (http://www.klaashartmann.com/treesample/) * <p/> * Sampling proportion not verified via simulation. Proportion set by default to 1, an assignment which makes the expressions * identical to the expressions before the change. * * @author Joseph Heled * Date: 24/02/2008 */ public class BirthDeathGernhard08Model extends UltrametricSpeciationModel implements Citable { public enum TreeType { UNSCALED, // no coefficient TIMESONLY, // n! ORIENTED, // n LABELED, // 2^(n-1)/(n-1)! (conditional on root: 2^(n-1)/n!(n-1) ) } public static final String BIRTH_DEATH_MODEL = BirthDeathModelParser.BIRTH_DEATH_MODEL; /* * mu/lambda * * null means default (0), or pure birth (Yule) */ private Parameter relativeDeathRateParameter; /** * lambda - mu */ private Parameter birthDiffRateParameter; private Parameter sampleProbability; private TreeType type; private boolean conditionalOnRoot; /** * rho * */ public BirthDeathGernhard08Model(Parameter birthDiffRateParameter, Parameter relativeDeathRateParameter, Parameter sampleProbability, TreeType type, Type units) { this(BIRTH_DEATH_MODEL, birthDiffRateParameter, relativeDeathRateParameter, sampleProbability, type, units, false); } public BirthDeathGernhard08Model(String modelName, Parameter birthDiffRateParameter, Parameter relativeDeathRateParameter, Parameter sampleProbability, TreeType type, Type units, boolean conditionalOnRoot) { super(modelName, units); this.birthDiffRateParameter = birthDiffRateParameter; addVariable(birthDiffRateParameter); birthDiffRateParameter.addBounds(new Parameter.DefaultBounds(Double.POSITIVE_INFINITY, 0.0, 1)); this.relativeDeathRateParameter = relativeDeathRateParameter; if (relativeDeathRateParameter != null) { addVariable(relativeDeathRateParameter); relativeDeathRateParameter.addBounds(new Parameter.DefaultBounds(1.0, 0.0, 1)); } this.sampleProbability = sampleProbability; if (sampleProbability != null) { addVariable(sampleProbability); sampleProbability.addBounds(new Parameter.DefaultBounds(1.0, 0.0, 1)); } this.conditionalOnRoot = conditionalOnRoot; if (conditionalOnRoot && sampleProbability != null) { throw new IllegalArgumentException( "Not supported: birth death prior conditional on root with sampling probability."); } this.type = type; } @Override public boolean isYule() { // Yule only return (relativeDeathRateParameter == null && sampleProbability == null && !conditionalOnRoot); } @Override public double getMarginal(Tree tree, CalibrationPoints calibration) { // Yule only return calibration.getCorrection(tree, getR()); } public double getR() { return birthDiffRateParameter.getParameterValue(0); } public double getA() { return relativeDeathRateParameter != null ? relativeDeathRateParameter.getParameterValue(0) : 0; } public double getRho() { return sampleProbability != null ? sampleProbability.getParameterValue(0) : 1.0; } private double logCoeff(int taxonCount) { switch (type) { case UNSCALED: break; case TIMESONLY: return logGamma(taxonCount + 1); case ORIENTED: return Math.log(taxonCount); case LABELED: { final double two2nm1 = (taxonCount - 1) * Math.log(2.0); if (!conditionalOnRoot) { return two2nm1 - logGamma(taxonCount); } else { return two2nm1 - Math.log(taxonCount - 1) - logGamma(taxonCount + 1); } } } return 0.0; } public double logTreeProbability(int taxonCount) { double c1 = logCoeff(taxonCount); if (!conditionalOnRoot) { c1 += (taxonCount - 1) * Math.log(getR() * getRho()) + taxonCount * Math.log(1 - getA()); } return c1; } public double logNodeProbability(Tree tree, NodeRef node) { final double height = tree.getNodeHeight(node); final double r = getR(); final double mrh = -r * height; final double a = getA(); if (!conditionalOnRoot) { final double rho = getRho(); final double z = Math.log(rho + ((1 - rho) - a) * Math.exp(mrh)); double l = -2 * z + mrh; if (tree.getRoot() == node) { l += mrh - z; } return l; } else { double l; if (tree.getRoot() != node) { final double z = Math.log(1 - a * Math.exp(mrh)); l = -2 * z + mrh; } else { // Root dependent coefficient from each internal node final double ca = 1 - a; final double emrh = Math.exp(-mrh); if (emrh != 1.0) { l = (tree.getTaxonCount() - 2) * Math.log(r * ca * (1 + ca / (emrh - 1))); } else { // use exp(x)-1 = x for x near 0 l = (tree.getTaxonCount() - 2) * Math.log(ca * (r + ca / height)); } } return l; } } public boolean includeExternalNodesInLikelihoodCalculation() { return false; } @Override public Citation.Category getCategory() { return Citation.Category.TREE_PRIORS; } @Override public String getDescription() { return "Gernhard 2008 Birth Death Tree Model"; } @Override public List<Citation> getCitations() { return Collections.singletonList( new Citation(new Author[] { new Author("T", "Gernhard"), }, "The conditioned reconstructed process", 2008, "Journal of Theoretical Biology", 253, 769, 778, "10.1016/j.jtbi.2008.04.005")); } }