diff --git a/src/dr/app/beast/development_parsers.properties b/src/dr/app/beast/development_parsers.properties index 9be4f1febc..bd9d7f2d5a 100644 --- a/src/dr/app/beast/development_parsers.properties +++ b/src/dr/app/beast/development_parsers.properties @@ -125,6 +125,8 @@ dr.evomodelxml.coalescent.PowerLawGrowthModelParser dr.evomodelxml.coalescent.PeakAndDeclineModelParser dr.evomodelxml.coalescent.AsymptoticGrowthModelParser +# INTEGRATED COALESCENT +dr.evomodelxml.coalescent.IGCoalescentLikelihoodParser # UNIFORM INTERNAL NODE HEIGHT PRIOR dr.evomodelxml.operators.FunkyPriorMixerOperatorParser diff --git a/src/dr/evomodel/coalescent/IGCoalescentLikelihood.java b/src/dr/evomodel/coalescent/IGCoalescentLikelihood.java new file mode 100644 index 0000000000..24b7e86325 --- /dev/null +++ b/src/dr/evomodel/coalescent/IGCoalescentLikelihood.java @@ -0,0 +1,137 @@ +/* + * CoalescentLikelihood.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.coalescent; + +import dr.evolution.coalescent.DemographicFunction; +import dr.evolution.coalescent.IntervalList; +import dr.evolution.coalescent.IntervalType; +import dr.evolution.tree.Tree; +import dr.evolution.tree.TreeUtils; +import dr.evolution.util.TaxonList; +import dr.evolution.util.Units; +import dr.evomodelxml.coalescent.CoalescentLikelihoodParser; +import dr.math.Binomial; +import dr.math.LogTricks; + +import java.util.List; + + +/** + * A likelihood function for the coalescent. Takes a tree and a demographic model. + * + * Parts of this class were derived from C++ code provided by Oliver Pybus. + * + * @version $Id: NewCoalescentLikelihood.java,v 1.6 2005/05/24 20:25:57 rambaut Exp $ + * + * @author Andrew Rambaut + * @author Alexei Drummond + * @author Luiz Max Carvalho + */ +public final class IGCoalescentLikelihood extends AbstractCoalescentLikelihood implements Units { + + + + // PUBLIC STUFF + public IGCoalescentLikelihood(Tree tree, + TaxonList includeSubtree, + List excludeSubtrees, + double alpha, double beta) throws TreeUtils.MissingTaxonException { + + super(CoalescentLikelihoodParser.COALESCENT_LIKELIHOOD, tree, includeSubtree, excludeSubtrees); + + this.alpha = alpha; + this.beta = beta; + } + + + /** + * Calculates the log likelihood of this set of coalescent intervals, + * given doubles alpha and beta + */ + public double calculateLogLikelihood() { + double lnL = calculateLogLikelihood(getIntervals(), alpha, beta); + return lnL; + } + + public static double calculateLogLikelihood(IntervalList intervals, double alpha, double beta) { + + double logL = 0.0; + double ldens = 0.0; + final int n = intervals.getIntervalCount(); + for (int i = 0; i < n; i++) { + + final double duration = intervals.getInterval(i); + final int lineageCount = intervals.getLineageCount(i); + final double kChoose2 = lineageCount*(lineageCount-1.0)/2.0; + + if (intervals.getIntervalType(i) == IntervalType.COALESCENT) { + ldens = marginalLogDensity(duration, kChoose2, alpha, beta); + logL += ldens; + }else { + // sampling interval + if(lineageCount == 1.0 || duration == 0.0) { + ldens = 0.0; + logL += ldens; + }else { + ldens = invGammaNoCoalescentlogProb(duration, kChoose2, alpha, beta); + logL += ldens; + } + } + // System.err.println("i=" + i + " k=" + lineageCount + " delta= " + duration + " isCoal=" + intervals.getIntervalType(i) + " ldens= " + ldens + "\n"); + } +// System.err.println("logLik=" + logL + "\n"); + return logL; + } + + public static double marginalLogDensity(double t, double coeff, double alpha, double beta) { + double lnum = Math.log(coeff) + Math.log(alpha) + alpha * Math.log(beta); + double ldenom = (alpha + 1) * Math.log(coeff * t + beta); + return(lnum - ldenom); + } + + public static double invGammaNoCoalescentlogProb(double x, double coeff, double alpha, double beta) { + //Gives (log) Pr(no coalescence | time = x, coeff, alpha, beta) + // Pr(no coal | time =x, coeff, Ne) = exp(-x * coeff/Ne). Then integrate against prior on Ne. + double logL = alpha * ( Math.log(beta)-Math.log(coeff*x + beta) ); + return logL; + } + + @Override + public Type getUnits() { + // TODO Auto-generated method stub + return null; + } + + + @Override + public void setUnits(Type units) { + // TODO Auto-generated method stub + } + + // PRIVATE STUFF + private double alpha; + private double beta; +} \ No newline at end of file diff --git a/src/dr/evomodelxml/coalescent/IGCoalescentLikelihoodParser.java b/src/dr/evomodelxml/coalescent/IGCoalescentLikelihoodParser.java new file mode 100644 index 0000000000..c695cb4044 --- /dev/null +++ b/src/dr/evomodelxml/coalescent/IGCoalescentLikelihoodParser.java @@ -0,0 +1,165 @@ +/* + * IGCoalescentLikelihoodParser.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.evomodelxml.coalescent; + +import dr.evolution.tree.TreeUtils; +import dr.evolution.util.Taxa; +import dr.evolution.util.TaxonList; +import dr.evomodel.coalescent.IGCoalescentLikelihood; +import dr.evomodel.coalescent.DemographicModel; +import dr.evomodel.coalescent.MultiLociTreeSet; +import dr.evomodel.coalescent.OldAbstractCoalescentLikelihood; +import dr.evomodel.tree.TreeModel; +import dr.xml.*; + +import java.util.ArrayList; +import java.util.List; + +/** + */ +public class IGCoalescentLikelihoodParser extends AbstractXMLObjectParser { + + public static final String IGCOALESCENT_LIKELIHOOD = "coalescentLikelihoodIG"; + public static final String MODEL = "model"; + public static final String POPULATION_TREE = "populationTree"; + public static final String POPULATION_FACTOR = "factor"; + + public static final String ALPHA = "alpha"; + public static final String BETA = "beta"; + + public static final String INCLUDE = "include"; + public static final String EXCLUDE = "exclude"; + + public String getParserName() { + return IGCOALESCENT_LIKELIHOOD; + } + + public Object parseXMLObject(XMLObject xo) throws XMLParseException { + + XMLObject cxo = xo.getChild(MODEL); + DemographicModel demoModel = null; + + double alpha = xo.getDoubleAttribute(ALPHA); + double beta = xo.getDoubleAttribute(BETA); + + List trees = new ArrayList(); + List popFactors = new ArrayList(); + MultiLociTreeSet treesSet = demoModel instanceof MultiLociTreeSet ? (MultiLociTreeSet) demoModel : null; + + for (int k = 0; k < xo.getChildCount(); ++k) { + final Object child = xo.getChild(k); + if (child instanceof XMLObject) { + cxo = (XMLObject) child; + if (cxo.getName().equals(POPULATION_TREE)) { + final TreeModel t = (TreeModel) cxo.getChild(TreeModel.class); + assert t != null; + trees.add(t); + + popFactors.add(cxo.getAttribute(POPULATION_FACTOR, 1.0)); + } + } +// in the future we may have arbitrary multi-loci element +// else if( child instanceof MultiLociTreeSet ) { +// treesSet = (MultiLociTreeSet)child; +// } + } + + TreeModel treeModel = null; + if (trees.size() == 1 && popFactors.get(0) == 1.0) { + treeModel = trees.get(0); + } else if (trees.size() > 1) { + treesSet = new MultiLociTreeSet.Default(trees, popFactors); + } else if (!(trees.size() == 0 && treesSet != null)) { + throw new XMLParseException("Incorrectly constructed likelihood element"); + } + + TaxonList includeSubtree = null; + + if (xo.hasChildNamed(INCLUDE)) { + includeSubtree = (TaxonList) xo.getElementFirstChild(INCLUDE); + } + + List excludeSubtrees = new ArrayList(); + + if (xo.hasChildNamed(EXCLUDE)) { + cxo = xo.getChild(EXCLUDE); + for (int i = 0; i < cxo.getChildCount(); i++) { + excludeSubtrees.add((TaxonList) cxo.getChild(i)); + } + } + + if (treeModel != null) { + try { + return new IGCoalescentLikelihood(treeModel, includeSubtree, excludeSubtrees, alpha, beta); + } catch (TreeUtils.MissingTaxonException mte) { + throw new XMLParseException("treeModel missing a taxon from taxon list in " + getParserName() + " element"); + } + } else { + if (includeSubtree != null || excludeSubtrees.size() > 0) { + throw new XMLParseException("Include/Exclude taxa not supported for multi locus sets"); + } + // Use old code for multi locus sets. + // This is a little unfortunate but the current code is using AbstractCoalescentLikelihood as + // a base - and modifing it will probsbly result in a bigger mess. + return new OldAbstractCoalescentLikelihood(treesSet, demoModel); + } + } + + //************************************************************************ + // AbstractXMLObjectParser implementation + //************************************************************************ + + public String getParserDescription() { + return "This element represents the likelihood of the tree given the demographic function."; + } + + public Class getReturnType() { + return IGCoalescentLikelihood.class; + } + + public XMLSyntaxRule[] getSyntaxRules() { + return rules; + } + + private final XMLSyntaxRule[] rules = { + new ElementRule(MODEL, new XMLSyntaxRule[]{ + new ElementRule(DemographicModel.class) + }, "The demographic model which describes the effective population size over time"), + + new ElementRule(POPULATION_TREE, new XMLSyntaxRule[]{ + AttributeRule.newDoubleRule(POPULATION_FACTOR, true), + new ElementRule(TreeModel.class) + }, "Tree(s) to compute likelihood for", 0, Integer.MAX_VALUE), + + new ElementRule(INCLUDE, new XMLSyntaxRule[]{ + new ElementRule(Taxa.class) + }, "An optional subset of taxa on which to calculate the likelihood (should be monophyletic)", true), + + new ElementRule(EXCLUDE, new XMLSyntaxRule[]{ + new ElementRule(Taxa.class, 1, Integer.MAX_VALUE) + }, "One or more subsets of taxa which should be excluded from calculate the likelihood (should be monophyletic)", true) + }; +} diff --git a/src/dr/math/distributions/InverseGammaDistribution.java b/src/dr/math/distributions/InverseGammaDistribution.java index 6ad4bc9737..3cc145051e 100644 --- a/src/dr/math/distributions/InverseGammaDistribution.java +++ b/src/dr/math/distributions/InverseGammaDistribution.java @@ -41,13 +41,13 @@ public class InverseGammaDistribution implements Distribution { private double shape, scale; private final double factor; - private final double logFacor; + private final double logFactor; public InverseGammaDistribution(double shape, double scale) { this.shape = shape; this.scale = scale; this.factor = Math.pow(scale, shape) / Math.exp(GammaFunction.lnGamma(shape)); - this.logFacor = shape * Math.log(scale) - GammaFunction.lnGamma(shape); + this.logFactor = shape * Math.log(scale) - GammaFunction.lnGamma(shape); } public double getShape() { @@ -71,7 +71,7 @@ public double pdf(double x) { } public double logPdf(double x) { - return logPdf(x, shape, scale, logFacor); + return logPdf(x, shape, scale, logFactor); } public double cdf(double x) { @@ -138,11 +138,11 @@ public static double pdf(double x, double shape, double scale, double factor) { * @param scale scale parameter * @return log pdf value */ - public static double logPdf(double x, double shape, double scale, double factor) { + public static double logPdf(double x, double shape, double scale, double lfactor) { if (x <= 0) return Double.NEGATIVE_INFINITY; - return factor + shape*Math.log(scale) - (shape + 1)*Math.log(x) - (scale/x) - GammaFunction.lnGamma(shape); + return lfactor + shape*Math.log(scale) - (shape + 1)*Math.log(x) - (scale/x) - GammaFunction.lnGamma(shape); } /**