Search in sources :

Example 41 with StrictClockBranchRates

use of dr.evomodel.branchratemodel.StrictClockBranchRates in project beast-mcmc by beast-dev.

the class AncestralStateTreeLikelihoodTest method testJointLikelihood.

public void testJointLikelihood() {
    TreeModel treeModel = new DefaultTreeModel("treeModel", tree);
    Sequence[] sequence = new Sequence[3];
    sequence[0] = new Sequence(new Taxon("0"), "A");
    sequence[1] = new Sequence(new Taxon("1"), "C");
    sequence[2] = new Sequence(new Taxon("2"), "C");
    Taxa taxa = new Taxa();
    for (Sequence s : sequence) {
        taxa.addTaxon(s.getTaxon());
    }
    SimpleAlignment alignment = new SimpleAlignment();
    for (Sequence s : sequence) {
        alignment.addSequence(s);
    }
    Parameter mu = new Parameter.Default(1, 1.0);
    Parameter kappa = new Parameter.Default(1, 1.0);
    double[] pi = { 0.25, 0.25, 0.25, 0.25 };
    Parameter freqs = new Parameter.Default(pi);
    FrequencyModel f = new FrequencyModel(Nucleotides.INSTANCE, freqs);
    HKY hky = new HKY(kappa, f);
    AncestralStateTreeLikelihood treeLikelihood = new AncestralStateTreeLikelihood(alignment, treeModel, new GammaSiteModel(hky), new StrictClockBranchRates(mu), false, true, Nucleotides.INSTANCE, "state", false, // useMap = true
    true, false);
    double logLike = treeLikelihood.getLogLikelihood();
    StringBuffer buffer = new StringBuffer();
    TreeUtils.newick(treeModel, treeModel.getRoot(), false, TreeUtils.BranchLengthType.LENGTHS_AS_TIME, null, null, new TreeTraitProvider[] { treeLikelihood }, null, buffer);
    System.out.println(buffer);
    System.out.println("t_CA(2) = " + t(false, 2.0));
    System.out.println("t_CC(1) = " + t(true, 1.0));
    double trueValue = 0.25 * t(false, 2.0) * Math.pow(t(true, 1.0), 3.0);
    assertEquals(logLike, Math.log(trueValue), 1e-6);
}
Also used : FrequencyModel(dr.oldevomodel.substmodel.FrequencyModel) Taxon(dr.evolution.util.Taxon) AncestralStateTreeLikelihood(dr.oldevomodel.treelikelihood.AncestralStateTreeLikelihood) DefaultTreeModel(dr.evomodel.tree.DefaultTreeModel) Sequence(dr.evolution.sequence.Sequence) StrictClockBranchRates(dr.evomodel.branchratemodel.StrictClockBranchRates) TreeModel(dr.evomodel.tree.TreeModel) DefaultTreeModel(dr.evomodel.tree.DefaultTreeModel) Taxa(dr.evolution.util.Taxa) GammaSiteModel(dr.oldevomodel.sitemodel.GammaSiteModel) SimpleAlignment(dr.evolution.alignment.SimpleAlignment) HKY(dr.oldevomodel.substmodel.HKY) Parameter(dr.inference.model.Parameter)

Aggregations

StrictClockBranchRates (dr.evomodel.branchratemodel.StrictClockBranchRates)41 BranchRateModel (dr.evomodel.branchratemodel.BranchRateModel)32 ArrayList (java.util.ArrayList)29 DefaultBranchRateModel (dr.evomodel.branchratemodel.DefaultBranchRateModel)26 Parameter (dr.inference.model.Parameter)26 TreeDataLikelihood (dr.evomodel.treedatalikelihood.TreeDataLikelihood)21 MultivariateElasticModel (dr.evomodel.continuous.MultivariateElasticModel)18 MatrixParameter (dr.inference.model.MatrixParameter)15 ArbitraryBranchRates (dr.evomodel.branchratemodel.ArbitraryBranchRates)10 NewickImporter (dr.evolution.io.NewickImporter)8 DiagonalMatrix (dr.inference.model.DiagonalMatrix)7 CladeNodeModel (dr.evomodel.bigfasttree.constrainedtree.CladeNodeModel)5 ConstrainedTreeBranchLengthProvider (dr.evomodel.bigfasttree.constrainedtree.ConstrainedTreeBranchLengthProvider)5 NodeRef (dr.evolution.tree.NodeRef)4 DefaultTreeModel (dr.evomodel.tree.DefaultTreeModel)4 TreeModel (dr.evomodel.tree.TreeModel)4 Tree (dr.evolution.tree.Tree)3 Taxon (dr.evolution.util.Taxon)3 CladeRef (dr.evomodel.bigfasttree.constrainedtree.CladeRef)3 GammaSiteRateModel (dr.evomodel.siteratemodel.GammaSiteRateModel)3