Search in sources :

Example 1 with AncestralStateTreeLikelihood

use of dr.oldevomodel.treelikelihood.AncestralStateTreeLikelihood in project beast-mcmc by beast-dev.

the class AncestralStateTreeLikelihoodTest method testJointLikelihood.

public void testJointLikelihood() {
    TreeModel treeModel = new TreeModel("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) Sequence(dr.evolution.sequence.Sequence) StrictClockBranchRates(dr.evomodel.branchratemodel.StrictClockBranchRates) TreeModel(dr.evomodel.tree.TreeModel) 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

SimpleAlignment (dr.evolution.alignment.SimpleAlignment)1 Sequence (dr.evolution.sequence.Sequence)1 Taxa (dr.evolution.util.Taxa)1 Taxon (dr.evolution.util.Taxon)1 StrictClockBranchRates (dr.evomodel.branchratemodel.StrictClockBranchRates)1 TreeModel (dr.evomodel.tree.TreeModel)1 Parameter (dr.inference.model.Parameter)1 GammaSiteModel (dr.oldevomodel.sitemodel.GammaSiteModel)1 FrequencyModel (dr.oldevomodel.substmodel.FrequencyModel)1 HKY (dr.oldevomodel.substmodel.HKY)1 AncestralStateTreeLikelihood (dr.oldevomodel.treelikelihood.AncestralStateTreeLikelihood)1