Search in sources :

Example 46 with RealParameter

use of beast.core.parameter.RealParameter in project bacter by tgvaughan.

the class ACGLikelihoodTest method testLikelihoodUsingSimulatedData.

@Test
public void testLikelihoodUsingSimulatedData() throws Exception {
    ConstantPopulation popFunc = new ConstantPopulation();
    popFunc.initByName("popSize", new RealParameter("1.0"));
    Locus locus = new Locus("locus", 10000);
    TaxonSet taxonSet = getTaxonSet(10);
    ConversionGraph acg = new SimulatedACG();
    acg.initByName("rho", 5.0 / locus.getSiteCount(), "delta", 1000.0, "populationModel", popFunc, "locus", locus, "taxonset", taxonSet);
    System.out.println(acg);
    // Site model:
    JukesCantor jc = new JukesCantor();
    jc.initByName();
    SiteModel siteModel = new SiteModel();
    siteModel.initByName("mutationRate", new RealParameter("1"), "substModel", jc);
    // Simulate alignment:
    SimulatedAlignment alignment = new SimulatedAlignment();
    alignment.initByName("acg", acg, "siteModel", siteModel, "outputFileName", "simulated_alignment.nexus", "useNexus", true);
    // Calculate likelihood:
    ACGLikelihood argLikelihood = new ACGLikelihood();
    argLikelihood.initByName("locus", locus, "data", alignment, "tree", acg, "siteModel", siteModel);
    double logP = argLikelihood.calculateLogP();
    // Compare product of likelihoods of "marginal alignments" with
    // likelihood computed using RGL.
    ACGLikelihoodSlow argLikelihoodSlow = new ACGLikelihoodSlow();
    argLikelihoodSlow.initByName("locus", locus, "data", alignment, "tree", acg, "siteModel", siteModel);
    double logPprime = argLikelihoodSlow.calculateLogP();
    double relError = 2.0 * Math.abs(logP - logPprime) / Math.abs(logP + logPprime);
    System.out.format("logP=%g\nlogPprime=%g\nrelError=%g\n", logP, logPprime, relError);
    assertTrue(relError < 1e-13);
}
Also used : ConstantPopulation(beast.evolution.tree.coalescent.ConstantPopulation) RealParameter(beast.core.parameter.RealParameter) SiteModel(beast.evolution.sitemodel.SiteModel) Locus(bacter.Locus) TaxonSet(beast.evolution.alignment.TaxonSet) JukesCantor(beast.evolution.substitutionmodel.JukesCantor) ConversionGraph(bacter.ConversionGraph) Test(org.junit.Test)

Example 47 with RealParameter

use of beast.core.parameter.RealParameter in project bacter by tgvaughan.

the class SkylinePopulationFunctionTest method test5.

@Test
public void test5() throws Exception {
    SkylinePopulationFunction skyline = new SkylinePopulationFunction();
    skyline.initByName("acg", acg, "popSizes", new RealParameter("1.0 1.0 5.0 1.0"), "groupSizes", new IntegerParameter("0"), "nGridPoints", 5, "piecewiseLinear", true);
    for (double t = 0.0; t < 10; t += 0.01) assertTrue(Math.abs(t - skyline.getInverseIntensity(skyline.getIntensity(t))) < 1e-14);
}
Also used : IntegerParameter(beast.core.parameter.IntegerParameter) RealParameter(beast.core.parameter.RealParameter) Test(org.junit.Test)

Example 48 with RealParameter

use of beast.core.parameter.RealParameter in project bacter by tgvaughan.

the class SkylinePopulationFunctionTest method test6.

@Test
public void test6() throws Exception {
    SkylinePopulationFunction skyline = new SkylinePopulationFunction();
    skyline.initByName("acg", acg, "popSizes", new RealParameter("1.0 1.0 5.0 1.0"), "groupSizes", new IntegerParameter("0"), "nGridPoints", 5, "piecewiseLinear", true);
    for (CFEventList.Event cfEvent : acg.getCFEvents()) {
        double t = cfEvent.getHeight();
        assertTrue(Math.abs(t - skyline.getInverseIntensity(skyline.getIntensity(t))) < 1e-14);
    }
}
Also used : IntegerParameter(beast.core.parameter.IntegerParameter) CFEventList(bacter.CFEventList) RealParameter(beast.core.parameter.RealParameter) Test(org.junit.Test)

Example 49 with RealParameter

use of beast.core.parameter.RealParameter in project bacter by tgvaughan.

the class AddRemoveConversionTest method testProbability.

/**
 * Tests whether probability of proposing a conversion lines up with
 * conversion probability found in ACGCoalescent.
 * @throws java.lang.Exception
 */
@Test
public void testProbability() throws Exception {
    ConstantPopulation popFunc = new ConstantPopulation();
    popFunc.initByName("popSize", new RealParameter("1.0"));
    Locus locus = new Locus("locus", 10000);
    TaxonSet taxonSet = getTaxonSet(10);
    SimulatedACG acg = new SimulatedACG();
    acg.initByName("rho", 1.0 / locus.getSiteCount(), "delta", 50.0, "locus", locus, "taxonset", taxonSet, "populationModel", popFunc);
    RealParameter rho = new RealParameter(Double.toString(1.0 / locus.getSiteCount()));
    RealParameter delta = new RealParameter("50.0");
    AddRemoveConversion operator = new AddRemoveConversion();
    operator.initByName("weight", 1.0, "acg", acg, "delta", delta, "populationModel", popFunc);
    ACGCoalescent coal = new ACGCoalescent();
    coal.initByName("tree", acg, "populationModel", popFunc, "rho", rho, "delta", delta);
    double logP1 = 0.0;
    double logP2 = 0.0;
    for (Conversion conv : acg.getConversions(locus)) {
        logP1 += operator.getConversionProb(conv);
        logP2 += coal.calculateConversionLogP(conv);
    }
    System.out.println("logP1 = " + logP1);
    System.out.println("logP2 = " + logP2);
    assertTrue(Math.abs(logP1 - logP2) < 1e-10);
}
Also used : ConstantPopulation(beast.evolution.tree.coalescent.ConstantPopulation) ACGCoalescent(bacter.model.ACGCoalescent) SimulatedACG(bacter.model.SimulatedACG) RealParameter(beast.core.parameter.RealParameter) Locus(bacter.Locus) TaxonSet(beast.evolution.alignment.TaxonSet) Conversion(bacter.Conversion) Test(org.junit.Test)

Example 50 with RealParameter

use of beast.core.parameter.RealParameter in project bacter by tgvaughan.

the class PiecewisePopulationFunction method main.

/**
 * Main method for testing.
 *
 * @param args command line arguments (unused)
 */
public static void main(String[] args) throws Exception {
    String acgString = "[&15,0,1.3905355989030808,31,770,1.597708055397074] " + "[&30,931,2.4351280458424904,36,2486,3.78055549386568] " + "[&15,941,2.0439957300083322,38,2364,6.911056700367016] " + "[&36,1091,4.285505683622974,38,2589,9.867725913197855] " + "((((10:0.5385300170206817,(17:0.116794353049212," + "((3:0.039229346597297564,12:0.039229346597297564)23:0.04582913870888949," + "13:0.08505848530618705)24:0.03173586774302495)26:0.4217356639714697)28:1.8114199763246093," + "((8:0.10883006062265468,2:0.10883006062265468)25:0.556428062025291," + "(6:0.5393311342677402,11:0.5393311342677402)29:0.12592698838020555)31:1.6846918706973453)34:1.4536824928125807," + "(1:0.47184545557390367,14:0.47184545557390367)27:3.331787030583968)37:2.9704369411362554," + "(((15:2.0624287390593707,((16:0.01825347077733299,19:0.01825347077733299)21:0.7668749128372041," + "(7:0.008018731329538273,9:0.008018731329538273)20:0.7771096522849988)32:1.2773003554448337)33:0.7487092404613747," + "4:2.8111379795207454)35:0.1331794525400949,((0:0.0243537216663141," + "5:0.0243537216663141)22:0.5681537100482162,18:0.5925074317145304)30:2.35181000034631)36:3.829751995233287)38:0.0";
    ConversionGraph acg = new ConversionGraph();
    acg.initByName("siteCount", 10000, "fromString", acgString);
    PiecewisePopulationFunction skyline = new PiecewisePopulationFunction();
    skyline.initByName("acg", acg, "popSizes", new RealParameter("1.0 1.0 5.0 1.0 2.0"), "groupSizes", new IntegerParameter("0"), "piecewiseLinear", true);
    try (PrintStream ps = new PrintStream("out.txt")) {
        ps.println("t N intensity intensityInv");
        double t = 0.0;
        while (t < 10) {
            ps.format("%g %g %g %g\n", t, skyline.getPopSize(t), skyline.getIntensity(t), skyline.getInverseIntensity(skyline.getIntensity(t)));
            t += 0.001;
        }
        ps.close();
    }
}
Also used : IntegerParameter(beast.core.parameter.IntegerParameter) PrintStream(java.io.PrintStream) RealParameter(beast.core.parameter.RealParameter) ConversionGraph(bacter.ConversionGraph)

Aggregations

RealParameter (beast.core.parameter.RealParameter)97 Test (org.junit.Test)39 IntegerParameter (beast.core.parameter.IntegerParameter)23 Tree (beast.evolution.tree.Tree)16 State (beast.core.State)14 SCMigrationModel (beast.evolution.tree.SCMigrationModel)13 Alignment (beast.evolution.alignment.Alignment)11 TypeSet (beast.evolution.tree.TypeSet)11 MCMC (beast.core.MCMC)10 SiteModel (beast.evolution.sitemodel.SiteModel)10 Frequencies (beast.evolution.substitutionmodel.Frequencies)10 ConstantPopulation (beast.evolution.tree.coalescent.ConstantPopulation)10 StructuredCoalescentTreeDensity (multitypetree.distributions.StructuredCoalescentTreeDensity)10 TaxonSet (beast.evolution.alignment.TaxonSet)9 MultiTypeTreeStatLogger (multitypetree.util.MultiTypeTreeStatLogger)9 MultiTypeTreeFromNewick (beast.evolution.tree.MultiTypeTreeFromNewick)8 Node (beast.evolution.tree.Node)8 Operator (beast.core.Operator)7 RandomTree (beast.evolution.tree.RandomTree)7 Locus (bacter.Locus)6