use of dr.oldevomodel.sitemodel.GammaSiteModel 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);
}
use of dr.oldevomodel.sitemodel.GammaSiteModel in project beast-mcmc by beast-dev.
the class LikelihoodTest method testLikelihoodGTRI.
public void testLikelihoodGTRI() {
System.out.println("\nTest Likelihood using GTRI:");
// Sub model
Parameter freqs = new Parameter.Default(alignment.getStateFrequencies());
FrequencyModel f = new FrequencyModel(Nucleotides.INSTANCE, freqs);
Variable<Double> rateACValue = new Parameter.Default(GTRParser.A_TO_C, 1.0, 1.0E-8, Double.POSITIVE_INFINITY);
Variable<Double> rateAGValue = new Parameter.Default(GTRParser.A_TO_G, 1.0, 1.0E-8, Double.POSITIVE_INFINITY);
Variable<Double> rateATValue = new Parameter.Default(GTRParser.A_TO_T, 1.0, 1.0E-8, Double.POSITIVE_INFINITY);
Variable<Double> rateCGValue = new Parameter.Default(GTRParser.C_TO_G, 1.0, 1.0E-8, Double.POSITIVE_INFINITY);
Variable<Double> rateCTValue = new Parameter.Default(GTRParser.C_TO_T, 1.0, 1.0E-8, Double.POSITIVE_INFINITY);
Variable<Double> rateGTValue = new Parameter.Default(GTRParser.G_TO_T, 1.0, 1.0E-8, Double.POSITIVE_INFINITY);
GTR gtr = new GTR(rateACValue, rateAGValue, rateATValue, rateCGValue, rateCTValue, rateGTValue, f);
//siteModel
Parameter mu = new Parameter.Default(GammaSiteModelParser.MUTATION_RATE, 1.0, 0, Double.POSITIVE_INFINITY);
Parameter invar = new Parameter.Default(GammaSiteModelParser.PROPORTION_INVARIANT, 0.5, 0, 1.0);
GammaSiteModel siteModel = new GammaSiteModel(gtr, mu, null, 4, invar);
//treeLikelihood
SitePatterns patterns = new SitePatterns(alignment, null, 0, -1, 1, true);
TreeLikelihood treeLikelihood = new TreeLikelihood(patterns, treeModel, siteModel, null, null, false, false, true, false, false);
assertEquals("treeLikelihoodGTRI", format.format(-1948.84175), format.format(treeLikelihood.getLogLikelihood()));
}
use of dr.oldevomodel.sitemodel.GammaSiteModel in project beast-mcmc by beast-dev.
the class LikelihoodTest method testLikelihoodGTRG.
public void testLikelihoodGTRG() {
System.out.println("\nTest Likelihood using GTRG:");
// Sub model
Parameter freqs = new Parameter.Default(alignment.getStateFrequencies());
FrequencyModel f = new FrequencyModel(Nucleotides.INSTANCE, freqs);
Variable<Double> rateACValue = new Parameter.Default(GTRParser.A_TO_C, 1.0, 1.0E-8, Double.POSITIVE_INFINITY);
Variable<Double> rateAGValue = new Parameter.Default(GTRParser.A_TO_G, 1.0, 1.0E-8, Double.POSITIVE_INFINITY);
Variable<Double> rateATValue = new Parameter.Default(GTRParser.A_TO_T, 1.0, 1.0E-8, Double.POSITIVE_INFINITY);
Variable<Double> rateCGValue = new Parameter.Default(GTRParser.C_TO_G, 1.0, 1.0E-8, Double.POSITIVE_INFINITY);
Variable<Double> rateCTValue = new Parameter.Default(GTRParser.C_TO_T, 1.0, 1.0E-8, Double.POSITIVE_INFINITY);
Variable<Double> rateGTValue = new Parameter.Default(GTRParser.G_TO_T, 1.0, 1.0E-8, Double.POSITIVE_INFINITY);
GTR gtr = new GTR(rateACValue, rateAGValue, rateATValue, rateCGValue, rateCTValue, rateGTValue, f);
//siteModel
Parameter mu = new Parameter.Default(GammaSiteModelParser.MUTATION_RATE, 1.0, 0, Double.POSITIVE_INFINITY);
Parameter shape = new Parameter.Default(GammaSiteModelParser.GAMMA_SHAPE, 0.5, 0, 100.0);
GammaSiteModel siteModel = new GammaSiteModel(gtr, mu, shape, 4, null);
//treeLikelihood
SitePatterns patterns = new SitePatterns(alignment, null, 0, -1, 1, true);
TreeLikelihood treeLikelihood = new TreeLikelihood(patterns, treeModel, siteModel, null, null, false, false, true, false, false);
assertEquals("treeLikelihoodGTRG", format.format(-1949.03601), format.format(treeLikelihood.getLogLikelihood()));
}
use of dr.oldevomodel.sitemodel.GammaSiteModel in project beast-mcmc by beast-dev.
the class LikelihoodTest method testLikelihoodHKY85GI.
public void testLikelihoodHKY85GI() {
System.out.println("\nTest Likelihood using HKY85GI:");
// Sub model
Parameter freqs = new Parameter.Default(alignment.getStateFrequencies());
Parameter kappa = new Parameter.Default(HKYParser.KAPPA, 39.464538, 0, 100);
FrequencyModel f = new FrequencyModel(Nucleotides.INSTANCE, freqs);
HKY hky = new HKY(kappa, f);
//siteModel
Parameter mu = new Parameter.Default(GammaSiteModelParser.MUTATION_RATE, 1.0, 0, Double.POSITIVE_INFINITY);
Parameter shape = new Parameter.Default(GammaSiteModelParser.GAMMA_SHAPE, 0.587649, 0, 1000.0);
Parameter invar = new Parameter.Default(GammaSiteModelParser.PROPORTION_INVARIANT, 0.486548, 0, 1.0);
GammaSiteModel siteModel = new GammaSiteModel(hky, mu, shape, 4, invar);
//treeLikelihood
SitePatterns patterns = new SitePatterns(alignment, null, 0, -1, 1, true);
TreeLikelihood treeLikelihood = new TreeLikelihood(patterns, treeModel, siteModel, null, null, false, false, true, false, false);
assertEquals("treeLikelihoodHKY85GI", format.format(-1789.63923), format.format(treeLikelihood.getLogLikelihood()));
}
use of dr.oldevomodel.sitemodel.GammaSiteModel in project beast-mcmc by beast-dev.
the class LikelihoodTest method testLikelihoodGTRGI.
public void testLikelihoodGTRGI() {
System.out.println("\nTest Likelihood using GTRGI:");
// Sub model
Parameter freqs = new Parameter.Default(alignment.getStateFrequencies());
FrequencyModel f = new FrequencyModel(Nucleotides.INSTANCE, freqs);
Variable<Double> rateACValue = new Parameter.Default(GTRParser.A_TO_C, 1.0, 1.0E-8, Double.POSITIVE_INFINITY);
Variable<Double> rateAGValue = new Parameter.Default(GTRParser.A_TO_G, 1.0, 1.0E-8, Double.POSITIVE_INFINITY);
Variable<Double> rateATValue = new Parameter.Default(GTRParser.A_TO_T, 1.0, 1.0E-8, Double.POSITIVE_INFINITY);
Variable<Double> rateCGValue = new Parameter.Default(GTRParser.C_TO_G, 1.0, 1.0E-8, Double.POSITIVE_INFINITY);
Variable<Double> rateCTValue = new Parameter.Default(GTRParser.C_TO_T, 1.0, 1.0E-8, Double.POSITIVE_INFINITY);
Variable<Double> rateGTValue = new Parameter.Default(GTRParser.G_TO_T, 1.0, 1.0E-8, Double.POSITIVE_INFINITY);
GTR gtr = new GTR(rateACValue, rateAGValue, rateATValue, rateCGValue, rateCTValue, rateGTValue, f);
//siteModel
Parameter mu = new Parameter.Default(GammaSiteModelParser.MUTATION_RATE, 1.0, 0, Double.POSITIVE_INFINITY);
Parameter shape = new Parameter.Default(GammaSiteModelParser.GAMMA_SHAPE, 0.5, 0, 100.0);
Parameter invar = new Parameter.Default(GammaSiteModelParser.PROPORTION_INVARIANT, 0.5, 0, 1.0);
GammaSiteModel siteModel = new GammaSiteModel(gtr, mu, shape, 4, invar);
//treeLikelihood
SitePatterns patterns = new SitePatterns(alignment, null, 0, -1, 1, true);
TreeLikelihood treeLikelihood = new TreeLikelihood(patterns, treeModel, siteModel, null, null, false, false, true, false, false);
assertEquals("treeLikelihoodGTRGI", format.format(-1947.58294), format.format(treeLikelihood.getLogLikelihood()));
}
Aggregations