diff --git a/src/dr/app/beast/development_parsers.properties b/src/dr/app/beast/development_parsers.properties index 021e964cfa..cd97e0f51d 100644 --- a/src/dr/app/beast/development_parsers.properties +++ b/src/dr/app/beast/development_parsers.properties @@ -1,7 +1,7 @@ # # development_parsers.properties # -# Copyright © 2002-2024 the BEAST Development Team +# Copyright © 2002-2024 the BEAST Development Team # http://beast.community/about # # This file is part of BEAST. diff --git a/src/dr/app/beast/release_parsers.properties b/src/dr/app/beast/release_parsers.properties index 9e1a9ec59f..3aab17d598 100644 --- a/src/dr/app/beast/release_parsers.properties +++ b/src/dr/app/beast/release_parsers.properties @@ -1,7 +1,7 @@ # # release_parsers.properties # -# Copyright © 2002-2024 the BEAST Development Team +# Copyright © 2002-2024 the BEAST Development Team # http://beast.community/about # # This file is part of BEAST. @@ -194,7 +194,6 @@ dr.evomodelxml.branchratemodel.RandomLocalClockModelParser dr.evomodelxml.tree.RLTVLoggerOnTreeParser dr.evomodelxml.branchratemodel.BranchCategoriesParser dr.evomodelxml.branchratemodel.CountableMixtureBranchRatesParser -dr.evomodelxml.branchratemodel.LatentStateBranchRateModelParser dr.evomodelxml.branchratemodel.RelaxedDriftModelParser dr.evomodelxml.branchratemodel.FixedDriftModelParser dr.evomodelxml.branchratemodel.BranchSpecificFixedEffectsParser @@ -204,6 +203,12 @@ dr.evomodelxml.branchratemodel.BranchRateGradientWrtIncrementsParser dr.evomodelxml.branchratemodel.AncestralTraitBranchRatesParser dr.evomodelxml.branchratemodel.RandomEffectsTreeTraitProviderParser +## Latent model +dr.evomodelxml.branchratemodel.latentStateBranchRate.LatentStateBranchRateModelParser +dr.evomodelxml.branchratemodel.latentStateBranchRate.LatentStateBranchRateLengthStatisticParser + + + #COALESCENT dr.evomodelxml.coalescent.demographicmodel.CataclysmicDemographicModelParser dr.evomodelxml.coalescent.demographicmodel.ExpConstExpDemographicModelParser diff --git a/src/dr/app/tools/MarkovRewardDensityCalculator.java b/src/dr/app/tools/MarkovRewardDensityCalculator.java new file mode 100644 index 0000000000..71c657e948 --- /dev/null +++ b/src/dr/app/tools/MarkovRewardDensityCalculator.java @@ -0,0 +1,162 @@ +package dr.app.tools; + +import java.io.IOException; + +import dr.app.util.Arguments; +import dr.inference.markovjumps.SericolaSeriesMarkovReward; +import dr.inference.markovjumps.TwoStateOccupancyMarkovReward; +import dr.inference.markovjumps.TwoStateSericolaSeriesMarkovReward; + +/** + * The goal is to calculate the pdf and cdf of time spent in a given state + * for a continuous-time Markov chain with rewards. + * + * There are several implementations that can achieve this. + * We will compare them here and output the results in a json + * file format. + * @author JT McCrone + */ +public class MarkovRewardDensityCalculator { + private double[] Q ; + private final int dim = 2; + private final double[] r = new double[]{0.0,1.0}; // rewards + + private TwoStateOccupancyMarkovReward twoStateOccupancyMarkovReward; + private TwoStateSericolaSeriesMarkovReward twoStateSericolaSeriesMarkovReward; + private SericolaSeriesMarkovReward sericolaSeriesMarkovReward; + + private MarkovRewardDensityCalculator(double rate,double bias){ + this.Q = new double[]{ + -rate * bias, rate * bias, + rate * (1.0 - bias), -rate * (1.0 - bias) + }; + + this.twoStateOccupancyMarkovReward = new TwoStateOccupancyMarkovReward(Q); + this.twoStateSericolaSeriesMarkovReward = new TwoStateSericolaSeriesMarkovReward(Q, r, dim); + this.sericolaSeriesMarkovReward = new SericolaSeriesMarkovReward(Q, r, dim); + + } + + + void run(double time){ + double s=0.0; + double[] times = new double[100]; + + double pdf; + double cdf; + double conditional; + + System.out.print("{\"data\":["); + + for(int i=0; i getBounds() { + return bounds; + } + + public void addDimension(int index, double value) { + throw new RuntimeException("Not yet implemented."); + } + + public double removeDimension(int index) { + throw new RuntimeException("Not yet implemented."); + } + @Override + public void modelChangedEvent(Model model, Object object, int index) { + fireParameterChangedEvent(); + } + + + @Override + public void modelRestored(Model model) { + // do nothing + } + + private SericolaLatentStateBranchRateModel model = null; + private STATE state; + private Tree tree; + public static String LATENT_STATE_BRANCH_RATE_LENGTH_STATISTIC = "latentStateBranchRateLengthStatistic"; + public enum STATE {REPLICATING, LATENT} + private Bounds bounds = null; + +} diff --git a/src/dr/evomodel/branchratemodel/LatentStateBranchRateModel.java b/src/dr/evomodel/branchratemodel/latentStateBranchRate/LatentStateBranchRateModel.java similarity index 98% rename from src/dr/evomodel/branchratemodel/LatentStateBranchRateModel.java rename to src/dr/evomodel/branchratemodel/latentStateBranchRate/LatentStateBranchRateModel.java index 21a6616baa..31a1cdc5ef 100644 --- a/src/dr/evomodel/branchratemodel/LatentStateBranchRateModel.java +++ b/src/dr/evomodel/branchratemodel/latentStateBranchRate/LatentStateBranchRateModel.java @@ -25,11 +25,13 @@ * */ -package dr.evomodel.branchratemodel; +package dr.evomodel.branchratemodel.latentStateBranchRate; import dr.evolution.tree.NodeRef; import dr.evolution.tree.Tree; import dr.evolution.tree.TreeTrait; +import dr.evomodel.branchratemodel.BranchRateModel; +import dr.evomodel.branchratemodel.CountableBranchCategoryProvider; import dr.evomodel.tree.TreeModel; import dr.evomodel.tree.TreeParameterModel; import dr.inference.markovjumps.TwoStateOccupancyMarkovReward; @@ -50,6 +52,8 @@ * $LastChangedDate$ * $LastChangedRevision$ */ +// Deprecated in favor of SericolaLatentStateBranchRateModel +@Deprecated public class LatentStateBranchRateModel extends AbstractModelLikelihood implements BranchRateModel { public static final String LATENT_STATE_BRANCH_RATE_MODEL = "latentStateBranchRateModel"; diff --git a/src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java b/src/dr/evomodel/branchratemodel/latentStateBranchRate/SericolaLatentStateBranchRateModel.java similarity index 81% rename from src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java rename to src/dr/evomodel/branchratemodel/latentStateBranchRate/SericolaLatentStateBranchRateModel.java index 230e7bee9b..33ef95c21a 100644 --- a/src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java +++ b/src/dr/evomodel/branchratemodel/latentStateBranchRate/SericolaLatentStateBranchRateModel.java @@ -25,15 +25,17 @@ * */ -package dr.evomodel.branchratemodel; +package dr.evomodel.branchratemodel.latentStateBranchRate; import dr.evolution.tree.NodeRef; import dr.evolution.tree.Tree; import dr.evolution.tree.TreeTrait; +import dr.evomodel.branchratemodel.BranchRateModel; +import dr.evomodel.branchratemodel.CountableBranchCategoryProvider; import dr.evomodel.tree.TreeModel; import dr.evomodel.tree.TreeParameterModel; import dr.inference.markovjumps.MarkovReward; -import dr.inference.markovjumps.TwoStateOccupancyMarkovReward; +import dr.inference.markovjumps.TwoStateSericolaSeriesMarkovReward; import dr.inference.model.AbstractModelLikelihood; import dr.inference.model.Model; import dr.inference.model.Parameter; @@ -72,6 +74,7 @@ public class SericolaLatentStateBranchRateModel extends AbstractModelLikelihood private final Parameter latentTransitionRateParameter; private final Parameter latentTransitionFrequencyParameter; private final TreeParameterModel latentStateProportions; + private final TreeParameterModel latentStateIndicators; private final Parameter latentStateProportionParameter; private final CountableBranchCategoryProvider branchCategoryProvider; @@ -79,8 +82,10 @@ public class SericolaLatentStateBranchRateModel extends AbstractModelLikelihood private MarkovReward storedSeries; private boolean likelihoodKnown = false; private boolean storedLikelihoodKnown; + private boolean excludeRoot; private double logLikelihood; private double storedLogLikelihood; + private double epsilon; private double[] branchLikelihoods; private double[] storedbranchLikelihoods; @@ -97,11 +102,16 @@ public SericolaLatentStateBranchRateModel(String name, Parameter latentTransitionRateParameter, Parameter latentTransitionFrequencyParameter, Parameter latentStateProportionParameter, - CountableBranchCategoryProvider branchCategoryProvider) { + Parameter latentStateIndicatorParameter, + CountableBranchCategoryProvider branchCategoryProvider, + double epsilon, + boolean excludeRoot) { super(name); this.tree = treeModel; addModel(tree); + + this.excludeRoot = excludeRoot; this.nonLatentRateModel = nonLatentRateModel; addModel(nonLatentRateModel); @@ -114,12 +124,15 @@ public SericolaLatentStateBranchRateModel(String name, if (branchCategoryProvider == null) { this.latentStateProportions = new TreeParameterModel(tree, latentStateProportionParameter, false, Intent.BRANCH); + this.latentStateIndicators = new TreeParameterModel(tree, latentStateIndicatorParameter, false, Intent.BRANCH); addModel(latentStateProportions); + addModel(latentStateIndicators); this.latentStateProportionParameter = null; this.branchCategoryProvider = null; } else { this.latentStateProportions = null; + this.latentStateIndicators = null; this.branchCategoryProvider = branchCategoryProvider; this.latentStateProportionParameter = latentStateProportionParameter; this.latentStateProportionParameter.setDimension(branchCategoryProvider.getCategoryCount()); @@ -141,8 +154,19 @@ public SericolaLatentStateBranchRateModel(String name, setUpdateAllBranches(); } + this.epsilon=epsilon; } + public SericolaLatentStateBranchRateModel(String name, + TreeModel treeModel, + BranchRateModel nonLatentRateModel, + Parameter latentTransitionRateParameter, + Parameter latentTransitionFrequencyParameter, + Parameter latentStateProportionParameter, + Parameter latentStateIndicatorParameter, + CountableBranchCategoryProvider branchCategoryProvider){ + this(name, treeModel, nonLatentRateModel, latentTransitionRateParameter, latentTransitionFrequencyParameter, latentStateProportionParameter, latentStateIndicatorParameter,branchCategoryProvider, 1E-10,false); + } public SericolaLatentStateBranchRateModel(Parameter rate, Parameter prop) { super(LATENT_STATE_BRANCH_RATE_MODEL); tree = null; @@ -151,6 +175,7 @@ public SericolaLatentStateBranchRateModel(Parameter rate, Parameter prop) { latentTransitionFrequencyParameter = prop; latentStateProportions = null; this.latentStateProportionParameter = null; + this.latentStateIndicators = null; this.branchCategoryProvider = null; } @@ -172,7 +197,8 @@ private static double[] createReward() { private MarkovReward createSeries() { // MarkovReward series = new SericolaSeriesMarkovReward(createLatentInfinitesimalMatrix(), // createReward(), 2); - MarkovReward series = new TwoStateOccupancyMarkovReward(createLatentInfinitesimalMatrix()); + // MarkovReward series = new TwoStateOccupancyMarkovReward(createLatentInfinitesimalMatrix()); + MarkovReward series = new TwoStateSericolaSeriesMarkovReward(createLatentInfinitesimalMatrix(),createReward(),2,epsilon); return series; } @@ -188,7 +214,7 @@ public double getBranchRate(Tree tree, NodeRef node) { public double getLatentProportion(Tree tree, NodeRef node) { if (latentStateProportions != null) { - return latentStateProportions.getNodeValue(tree, node); + return latentStateProportions.getNodeValue(tree, node)*latentStateIndicators.getNodeValue(tree,node); } else { return latentStateProportionParameter.getParameterValue(branchCategoryProvider.getBranchCategory(tree, node)); } @@ -211,7 +237,7 @@ protected void handleModelChangedEvent(Model model, Object object, int index) { } else if (model == nonLatentRateModel) { // rates will change but the latent proportions haven't so the density is unchanged - } else if (model == latentStateProportions) { + } else if (model == latentStateProportions || model ==latentStateIndicators) { likelihoodKnown = false; // argument of density has changed if (index == -1) { @@ -348,24 +374,36 @@ private double calculateLogLikelihood() { double logLike = 0.0; - for (int i = 0; i < tree.getInternalNodeCount(); ++i) { + for (int i = 0; i < tree.getNodeCount(); ++i) { NodeRef node = tree.getNode(i); if (node != tree.getRoot()) { if (updateNeededForNode(tree, node)) { + if(excludeRoot && tree.getParent(node)== tree.getRoot()){ + if(getLatentProportion(tree,node)>0){ + // branchLikelihoods[node.getNumber()]=Double.NEGATIVE_INFINITY; // reject we don't allow latency on these branches + return Double.NEGATIVE_INFINITY; + } + branchLikelihoods[node.getNumber()]=0.0; // if it is zero we don't get credit for it since it couldn't be otherwise + }else{ double branchLength = tree.getBranchLength(node); double latentProportion = getLatentProportion(tree, node); assert(latentProportion < 1.0); + double reward = branchLength * latentProportion; double density = getBranchRewardDensity(reward, branchLength); + if(latentProportion>0){ + density *=branchLength;// jacobian for sampling on 0,1 + } branchLikelihoods[node.getNumber()] = Math.log(density); } + } logLike += branchLikelihoods[node.getNumber()]; // TODO More importantly, MH proposals on [0,1] may be missing a Jacobian for which we should adjust. // TODO This is easy to test and we should do it when sampling appears to work. - } } + } clearUpdateAllBranches(); clearAllCategories(); @@ -390,7 +428,27 @@ public double getBranchRewardDensity(double reward, double branchLength) { // Reward is [0,1], and we want to track time in latent state (= 1). // Therefore all nodes are in state 0 // double joint = series.computePdf(reward, branchLength)[state]; - double joint = series.computePdf(reward, branchLength, 0, 0); + double rate = latentTransitionRateParameter.getParameterValue(0); + double prop = latentTransitionFrequencyParameter.getParameterValue(0); + double noJump = Math.exp(-1 * rate * prop * branchLength); // no jump + + double joint; + + if (reward == 0) { + // joint = series.computeCdf(reward, branchLength, 0, 0); + // no jump. + joint = noJump; + // series.computeCdf(reward, branchLength, 0, 0)- is no jumps using the series + // machinery + assert Math.abs(series.computeCdf(reward, branchLength, 0, 0) - joint) < epsilon + : "0 jumps not calculated correctly: " + + Math.abs(series.computeCdf(reward, branchLength, 0, 0) - joint) + "epsilon: " + epsilon + + " rate: " + rate + " prop: " + prop + " bl: " + branchLength; + // System.out.println(joint + " = " + noJump); + } else { + joint = series.computePdf(reward, branchLength, 0, 0); + } + double marg = series.computeConditionalProbability(branchLength, 0, 0); // TODO Overhead in creating double[] could be saved by changing signature to computePdf @@ -420,7 +478,10 @@ public TreeTrait getTreeTrait(final String key) { return this; } else if (latentStateProportions != null && key.equals(latentStateProportions.getTraitName())) { return latentStateProportions; - } else if (branchCategoryProvider != null && key.equals(branchCategoryProvider.getTraitName())) { + }else if(latentStateIndicators!=null && key.equals(latentStateIndicators.getTraitName())){ + return latentStateIndicators; + } + else if (branchCategoryProvider != null && key.equals(branchCategoryProvider.getTraitName())) { return branchCategoryProvider; } else { throw new IllegalArgumentException("Unrecognised Tree Trait key, " + key); diff --git a/src/dr/evomodelxml/branchratemodel/latentStateBranchRate/LatentStateBranchRateLengthStatisticParser.java b/src/dr/evomodelxml/branchratemodel/latentStateBranchRate/LatentStateBranchRateLengthStatisticParser.java new file mode 100644 index 0000000000..7da535e11f --- /dev/null +++ b/src/dr/evomodelxml/branchratemodel/latentStateBranchRate/LatentStateBranchRateLengthStatisticParser.java @@ -0,0 +1,44 @@ +package dr.evomodelxml.branchratemodel.latentStateBranchRate; + + +import dr.evomodel.branchratemodel.latentStateBranchRate.LatentStateBranchRateLengthStatistic; +import dr.evomodel.branchratemodel.latentStateBranchRate.SericolaLatentStateBranchRateModel; +import dr.evomodel.tree.TreeModel; +import dr.xml.AbstractXMLObjectParser; +import dr.xml.ElementRule; +import dr.xml.XMLObject; +import dr.xml.XMLParseException; +import dr.xml.XMLSyntaxRule; + +public class LatentStateBranchRateLengthStatisticParser extends AbstractXMLObjectParser { + + private static final String STATE="state"; + + public String getParserName() { + return LatentStateBranchRateLengthStatistic.LATENT_STATE_BRANCH_RATE_LENGTH_STATISTIC; + } + + public Object parseXMLObject(XMLObject xo) throws XMLParseException { + SericolaLatentStateBranchRateModel latenBranchRateModel = (SericolaLatentStateBranchRateModel) xo.getChild(SericolaLatentStateBranchRateModel.class); + TreeModel tree = (TreeModel) xo.getChild(TreeModel.class); + LatentStateBranchRateLengthStatistic.STATE state = LatentStateBranchRateLengthStatistic.STATE.valueOf(xo.getAttribute(STATE, LatentStateBranchRateLengthStatistic.STATE.REPLICATING.name()).toUpperCase()); + return new LatentStateBranchRateLengthStatistic(latenBranchRateModel, tree, state); + } + public String getParserDescription() { + return "This element provides a statistic for the length of a tree spent in a replicating or latent state."; + } + + public Class getReturnType() { + return LatentStateBranchRateLengthStatistic.class; + } + + public XMLSyntaxRule[] getSyntaxRules() { + return rules; + } + + private XMLSyntaxRule[] rules = new XMLSyntaxRule[]{ + new ElementRule(SericolaLatentStateBranchRateModel.class, "A branch rate model to provide the rates for the non-latent state"), + new ElementRule(TreeModel.class, "The tree on which this will operate"), + }; + +} diff --git a/src/dr/evomodelxml/branchratemodel/LatentStateBranchRateModelParser.java b/src/dr/evomodelxml/branchratemodel/latentStateBranchRate/LatentStateBranchRateModelParser.java similarity index 71% rename from src/dr/evomodelxml/branchratemodel/LatentStateBranchRateModelParser.java rename to src/dr/evomodelxml/branchratemodel/latentStateBranchRate/LatentStateBranchRateModelParser.java index 95300da423..ee77715145 100644 --- a/src/dr/evomodelxml/branchratemodel/LatentStateBranchRateModelParser.java +++ b/src/dr/evomodelxml/branchratemodel/latentStateBranchRate/LatentStateBranchRateModelParser.java @@ -25,9 +25,10 @@ * */ -package dr.evomodelxml.branchratemodel; +package dr.evomodelxml.branchratemodel.latentStateBranchRate; import dr.evomodel.branchratemodel.*; +import dr.evomodel.branchratemodel.latentStateBranchRate.SericolaLatentStateBranchRateModel; import dr.evomodel.tree.TreeModel; import dr.inference.model.Parameter; import dr.xml.*; @@ -40,6 +41,9 @@ public class LatentStateBranchRateModelParser extends AbstractXMLObjectParser { public static final String LATENT_TRANSITION_RATE = "latentTransitionRate"; public static final String LATENT_TRANSITION_FREQUENCY = "latentTransitionFrequency"; public static final String LATENT_STATE_PROPORTIONS = "latentStateProportions"; + public static final String LATENT_STATE_INDICATORS = "latentStateIndicators"; + public static final String EPSILON = "epsilon"; + public static final String EXCLUDE_ROOT = "excludeRoot"; public String getParserName() { @@ -51,24 +55,30 @@ public Object parseXMLObject(XMLObject xo) throws XMLParseException { TreeModel tree = (TreeModel) xo.getChild(TreeModel.class); Parameter latentTransitionRateParameter = (Parameter) xo.getElementFirstChild(LATENT_TRANSITION_RATE); Parameter latentTransitionFrequencyParameter = (Parameter) xo.getElementFirstChild(LATENT_TRANSITION_FREQUENCY); + boolean excludeRoot = xo.getAttribute("excludeRoot", false); CountableBranchCategoryProvider branchCategoryProvider = (CountableBranchCategoryProvider)xo.getChild(CountableBranchCategoryProvider.class); Parameter latentStateProportionParameter = null; + if (xo.hasChildNamed(LATENT_STATE_PROPORTIONS)) { latentStateProportionParameter = (Parameter) xo.getElementFirstChild(LATENT_STATE_PROPORTIONS); } - + Parameter latentStateIndicatorParameter = null; + if (xo.hasChildNamed(LATENT_STATE_INDICATORS)) { + latentStateIndicatorParameter = (Parameter) xo.getElementFirstChild(LATENT_STATE_INDICATORS); + } + double epsilon = xo.getAttribute(EPSILON, 1E-10); Logger.getLogger("dr.evomodel").info("\nCreating a latent state branch rate model"); - return new LatentStateBranchRateModel(LatentStateBranchRateModel.LATENT_STATE_BRANCH_RATE_MODEL, - tree, nonLatentRateModel, - latentTransitionRateParameter, latentTransitionFrequencyParameter, /* 0/1 CTMC have two parameters */ - latentStateProportionParameter, branchCategoryProvider); -// return new SericolaLatentStateBranchRateModel(SericolaLatentStateBranchRateModel.LATENT_STATE_BRANCH_RATE_MODEL, -// tree, nonLatentRateModel, -// latentTransitionRateParameter, latentTransitionFrequencyParameter, /* 0/1 CTMC have two parameters */ -// latentStateProportionParameter, branchCategoryProvider); + // return new LatentStateBranchRateModel(LatentStateBranchRateModel.LATENT_STATE_BRANCH_RATE_MODEL, + // tree, nonLatentRateModel, + // latentTransitionRateParameter, latentTransitionFrequencyParameter, /* 0/1 CTMC have two parameters */ + // latentStateProportionParameter, branchCategoryProvider); + return new SericolaLatentStateBranchRateModel(SericolaLatentStateBranchRateModel.LATENT_STATE_BRANCH_RATE_MODEL, + tree, nonLatentRateModel, + latentTransitionRateParameter, latentTransitionFrequencyParameter, /* 0/1 CTMC have two parameters */ + latentStateProportionParameter,latentStateIndicatorParameter, branchCategoryProvider,epsilon,excludeRoot); } //************************************************************************ diff --git a/src/dr/inference/markovjumps/TwoStateSericolaSeriesMarkovReward.java b/src/dr/inference/markovjumps/TwoStateSericolaSeriesMarkovReward.java index b8f8a192e9..31830fc93e 100644 --- a/src/dr/inference/markovjumps/TwoStateSericolaSeriesMarkovReward.java +++ b/src/dr/inference/markovjumps/TwoStateSericolaSeriesMarkovReward.java @@ -145,8 +145,18 @@ public double computePdf(double x, double time, int i, int j) { double w = 0.0; final int N = getNfromC() - 1; + + // moving these caluculations out of the loop + final double premultFactor = (Math.log(lambda) + Math.log(time)); + final int h = 1; + final double factor = lambda / (r[h] - r[h - 1]); + double xh = (x - r[h - 1] * time) / ((r[h] - r[h - 1]) * time); + for (int n = 0; n <= N; ++n) { - w += accumulatePdf(x, n, time, uv); + final double premult = Math.exp( + -lambda * time + n * premultFactor - getLnGamma(n + 1.0) + ); + w += accumulatePdf(x, n, time, uv,premult,xh,factor); } // if (DEBUG2) { @@ -257,33 +267,46 @@ private void accumulateCdf(double[][] W, double[] X, int[] H, int n, double time } } - private double accumulatePdf(double x, int n, double time, int uv) { + + private double accumulatePdf(double x, int n, double time, int uv, double premult, double xh,double factor) { double w = 0.0; - final double premult = Math.exp( - -lambda * time + n * (Math.log(lambda) + Math.log(time)) - GammaFunction.lnGamma(n + 1.0) - ); + // final double premult = Math.exp( + // // -lambda * time + n * (Math.log(lambda) + Math.log(time)) - getLnGamma(n + 1.0) + // -lambda * time + n * (Math.log(lambda) + Math.log(time)) - getLnGamma(n + 1.0) + // ); // TODO Make factorial/choose static look-up tables // for (int t = 0; t < X.length; ++t) { // For each time point int h = 1; - final double factor = lambda / (r[h] - r[h - 1]); + // final double factor = lambda / (r[h] - r[h - 1]); - double xh = (x - r[h - 1] * time) / ((r[h] - r[h - 1]) * time); + // double xh = (x - r[h - 1] * time) / ((r[h] - r[h - 1]) * time); // final int dim2 = dim * dim; // double[] inc = new double[dim2]; // W^{\epsilon}(x(i),t,n) double inc = 0.0; + double xhPow = 1.0; + double oneMinusXhPow = Math.pow(1.0 - xh, n); for (int k = 0; k <= n; k++) { - final double binomialCoef = Binomial.choose(n, k) * Math.pow(xh, k) * Math.pow(1.0 - xh, n - k); -// for (int uv = 0; uv < dim2; ++uv) { - inc += binomialCoef * (C(h, n + 1, k + 1)[uv] - C(h, n + 1, k)[uv]); -// } + double binomialCoef = choose(n, k) * xhPow * oneMinusXhPow; + inc += binomialCoef * (C(h, n + 1, k + 1)[uv] - C(h, n + 1, k)[uv]); + // update powers for next iteration + xhPow *= xh; + oneMinusXhPow /= (1.0 - xh); } + +// for (int k = 0; k <= n; k++) { +// final double binomialCoef = choose(n, k) * Math.pow(xh, k) * Math.pow(1.0 - xh, n - k); +// // for (int uv = 0; uv < dim2; ++uv) { +// inc += binomialCoef * (C(h, n + 1, k + 1)[uv] - C(h, n + 1, k)[uv]); +// // } +// } + // for (int uv = 0; uv < dim2; ++uv) { w += factor * premult * inc; // } @@ -471,7 +494,7 @@ private int determineNumberOfSteps(double time, double lambda) { while (Math.abs(sum2 - tolerance2) > epsilon && sum2 < 1.0) { // while (sum2 < tolerance2) { i++; - double logDensity = -lambda * time + i * (Math.log(lambda) + Math.log(time)) - GammaFunction.lnGamma(i + 1); + double logDensity = -lambda * time + i * (Math.log(lambda) + Math.log(time)) -getLnGamma(i + 1); sum2 += Math.exp(logDensity); // sum2 = LogTricks.logSum(sum2, logDensity); if (DEBUG2) { @@ -487,6 +510,25 @@ private int determineNumberOfSteps(double time, double lambda) { return i; } + private double choose(int n, int k){ + if (k < 0 || k > n) return 0; + double lchoose = getLnGamma(n + 1.0) - + getLnGamma(k + 1.0) - getLnGamma(n - k + 1.0); + + return Math.floor(Math.exp(lchoose) + 0.5); + } + // should be the same as calling GammaFunction.lnGamma + private double getLnGamma(double n){ + n = Math.floor(n + 0.5); + if (lnGamma == null || lnGamma.length <= n ) { // fill it up! + lnGamma = new double[(int) (n+50)]; + for (int i = 0; i < (n+50); ++i) { + lnGamma[i] = GammaFunction.lnGamma(i); + } + } + return lnGamma[(int) n]; + } + public String toString() { StringBuilder sb = new StringBuilder(); sb.append("Q: " + new Vector(Q) + "\n"); @@ -528,6 +570,9 @@ public double computeConditionalProbability(double distance, int i, int j) { private final int dim; private final double epsilon; + private double[][] binomialCoefficients; + private double[] lnGamma; + private final EigenSystem eigenSystem; private double maxTime; diff --git a/src/test/dr/evomodel/branchratemodel/latentStateBranchRate/LatentRewardDensityTest.java b/src/test/dr/evomodel/branchratemodel/latentStateBranchRate/LatentRewardDensityTest.java new file mode 100644 index 0000000000..cde5beb195 --- /dev/null +++ b/src/test/dr/evomodel/branchratemodel/latentStateBranchRate/LatentRewardDensityTest.java @@ -0,0 +1,146 @@ +package test.dr.evomodel.branchratemodel.latentStateBranchRate; + +import dr.inference.markovjumps.SericolaSeriesMarkovReward; +import dr.inference.markovjumps.TwoStateOccupancyMarkovReward; +import dr.inference.markovjumps.TwoStateSericolaSeriesMarkovReward; +import junit.framework.TestCase; +/** + * @author JT McCrone + * This class tests the various classes that estimate the reward for a two state + * model. The TwoStateSericolaSeriesMarkovReward is used in the SericolaLatentStateBranchModel. + * TwoStateOccupancyMarkovReward is used in the deprecated LatentStateBranchRateModel + * We compare the cdf, pdf, and conditional probabilities of these classes where applicable. + */ +public class LatentRewardDensityTest extends TestCase{ + + private final int dim = 2; + private final double[] r = new double[]{0.0,1.0}; // rewards + private final double epsilon = 1e-5; + private double[] Q(double rate, double bias){ + return new double[]{ + -rate * bias, rate * bias, + rate * (1.0 - bias), -rate * (1.0 - bias) + }; + } + private final double totalTime = 100.0; + + public void testPdfEqualRates(){ + double rate = 0.05; + double bias = 0.5; + TwoStateOccupancyMarkovReward TSOMR = createTwoStateOccupancyMarkovReward(rate, bias); + TwoStateSericolaSeriesMarkovReward TSSMR = createTwoStateSericolaSeriesMarkovReward(rate, bias); + SericolaSeriesMarkovReward SSMR = createSericolaSeriesMarkovReward(rate, bias); + + // pdf + double[] latentTimes = new double[10]; + + for(int i=0; i