Examples with SpeciationLikelihood dr.evomodel.speciation.SpeciationLikelihood used on opensource projects

Example 1 with SpeciationLikelihood

use of dr.evomodel.speciation.SpeciationLikelihood in project beast-mcmc by beast-dev.

the class RLYModelTest method randomLocalYuleTester.

private void randomLocalYuleTester(TreeModel treeModel, Parameter I, Parameter b, OperatorSchedule schedule) {
    MCMC mcmc = new MCMC("mcmc1");
    MCMCOptions options = new MCMCOptions(1000000);
    TreeLengthStatistic tls = new TreeLengthStatistic(TL, treeModel);
    TreeHeightStatistic rootHeight = new TreeHeightStatistic(TREE_HEIGHT, treeModel);
    Parameter m = new Parameter.Default("m", 1.0, 0.0, Double.MAX_VALUE);
    SpeciationModel speciationModel = new RandomLocalYuleModel(b, I, m, false, Units.Type.YEARS, 4);
    Likelihood likelihood = new SpeciationLikelihood(treeModel, speciationModel, "randomYule.like");
    ArrayLogFormatter formatter = new ArrayLogFormatter(false);
    MCLogger[] loggers = new MCLogger[2];
    loggers[0] = new MCLogger(formatter, 100, false);
    loggers[0].add(likelihood);
    loggers[0].add(rootHeight);
    loggers[0].add(tls);
    loggers[0].add(I);
    loggers[1] = new MCLogger(new TabDelimitedFormatter(System.out), 100000, false);
    loggers[1].add(likelihood);
    loggers[1].add(rootHeight);
    loggers[1].add(tls);
    loggers[1].add(I);
    mcmc.setShowOperatorAnalysis(true);
    mcmc.init(options, likelihood, schedule, loggers);
    mcmc.run();
    List<Trace> traces = formatter.getTraces();
    ArrayTraceList traceList = new ArrayTraceList("yuleModelTest", traces, 0);
    for (int i = 1; i < traces.size(); i++) {
        traceList.analyseTrace(i);
    }
    TraceCorrelation tlStats = traceList.getCorrelationStatistics(traceList.getTraceIndex("root." + birthRateIndicator));
    System.out.println("mean = " + tlStats.getMean());
    System.out.println("expected mean = 0.5");
    assertExpectation("root." + birthRateIndicator, tlStats, 0.5);
}

Example 2 with SpeciationLikelihood

use of dr.evomodel.speciation.SpeciationLikelihood in project beast-mcmc by beast-dev.

the class YuleModelTest method yuleTester.

//    public void testYuleWithWideExchange() {
//
//        TreeModel treeModel = new TreeModel("treeModel", tree);
// Doesn't compile...
//      yuleTester(treeModel, ExchangeOperatorTest.getWideExchangeSchedule(treeModel));
//    }
private void yuleTester(TreeModel treeModel, OperatorSchedule schedule) {
    MCMC mcmc = new MCMC("mcmc1");
    MCMCOptions options = new MCMCOptions(1000000);
    TreeLengthStatistic tls = new TreeLengthStatistic(TL, treeModel);
    TreeHeightStatistic rootHeight = new TreeHeightStatistic(TREE_HEIGHT, treeModel);
    Parameter b = new Parameter.Default("b", 2.0, 0.0, Double.MAX_VALUE);
    Parameter d = new Parameter.Default("d", 0.0, 0.0, Double.MAX_VALUE);
    SpeciationModel speciationModel = new BirthDeathGernhard08Model(b, d, null, BirthDeathGernhard08Model.TreeType.TIMESONLY, Units.Type.YEARS);
    Likelihood likelihood = new SpeciationLikelihood(treeModel, speciationModel, "yule.like");
    ArrayLogFormatter formatter = new ArrayLogFormatter(false);
    MCLogger[] loggers = new MCLogger[2];
    loggers[0] = new MCLogger(formatter, 100, false);
    loggers[0].add(likelihood);
    loggers[0].add(rootHeight);
    loggers[0].add(tls);
    loggers[1] = new MCLogger(new TabDelimitedFormatter(System.out), 100000, false);
    loggers[1].add(likelihood);
    loggers[1].add(rootHeight);
    loggers[1].add(tls);
    mcmc.setShowOperatorAnalysis(true);
    mcmc.init(options, likelihood, schedule, loggers);
    mcmc.run();
    List<Trace> traces = formatter.getTraces();
    ArrayTraceList traceList = new ArrayTraceList("yuleModelTest", traces, 0);
    for (int i = 1; i < traces.size(); i++) {
        traceList.analyseTrace(i);
    }
    // expectation of root height for 4 tips and lambda = 2
    // rootHeight = 0.541666
    // TL = 1.5
    TraceCorrelation tlStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(TL));
    assertExpectation(TL, tlStats, 1.5);
    TraceCorrelation treeHeightStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(TREE_HEIGHT));
    assertExpectation(TREE_HEIGHT, treeHeightStats, 0.5416666);
}

Example 3 with SpeciationLikelihood

use of dr.evomodel.speciation.SpeciationLikelihood in project beast-mcmc by beast-dev.

the class OperatorAssert method irreducibilityTester.

private void irreducibilityTester(Tree tree, int numLabelledTopologies, int chainLength, int sampleTreeEvery) throws IOException, Importer.ImportException {
    MCMC mcmc = new MCMC("mcmc1");
    MCMCOptions options = new MCMCOptions(chainLength);
    TreeModel treeModel = new TreeModel("treeModel", tree);
    TreeLengthStatistic tls = new TreeLengthStatistic(TL, treeModel);
    TreeHeightStatistic rootHeight = new TreeHeightStatistic(TREE_HEIGHT, treeModel);
    OperatorSchedule schedule = getOperatorSchedule(treeModel);
    Parameter b = new Parameter.Default("b", 2.0, 0.0, Double.MAX_VALUE);
    Parameter d = new Parameter.Default("d", 0.0, 0.0, Double.MAX_VALUE);
    SpeciationModel speciationModel = new BirthDeathGernhard08Model(b, d, null, BirthDeathGernhard08Model.TreeType.UNSCALED, Units.Type.YEARS);
    Likelihood likelihood = new SpeciationLikelihood(treeModel, speciationModel, "yule.like");
    MCLogger[] loggers = new MCLogger[2];
    //        loggers[0] = new MCLogger(new ArrayLogFormatter(false), 100, false);
    //        loggers[0].add(likelihood);
    //        loggers[0].add(rootHeight);
    //        loggers[0].add(tls);
    loggers[0] = new MCLogger(new TabDelimitedFormatter(System.out), 10000, false);
    loggers[0].add(likelihood);
    loggers[0].add(rootHeight);
    loggers[0].add(tls);
    File file = new File("yule.trees");
    file.deleteOnExit();
    FileOutputStream out = new FileOutputStream(file);
    loggers[1] = new TreeLogger(treeModel, new TabDelimitedFormatter(out), sampleTreeEvery, true, true, false);
    mcmc.setShowOperatorAnalysis(true);
    mcmc.init(options, likelihood, schedule, loggers);
    mcmc.run();
    out.flush();
    out.close();
    Set<String> uniqueTrees = new HashSet<String>();
    HashMap<String, Integer> topologies = new HashMap<String, Integer>();
    HashMap<String, HashMap<String, Integer>> treeCounts = new HashMap<String, HashMap<String, Integer>>();
    NexusImporter importer = new NexusImporter(new FileReader(file));
    int sampleSize = 0;
    while (importer.hasTree()) {
        sampleSize++;
        Tree t = importer.importNextTree();
        String uniqueNewick = TreeUtils.uniqueNewick(t, t.getRoot());
        String topology = uniqueNewick.replaceAll("\\w+", "X");
        if (!uniqueTrees.contains(uniqueNewick)) {
            uniqueTrees.add(uniqueNewick);
        }
        HashMap<String, Integer> counts;
        if (topologies.containsKey(topology)) {
            topologies.put(topology, topologies.get(topology) + 1);
            counts = treeCounts.get(topology);
        } else {
            topologies.put(topology, 1);
            counts = new HashMap<String, Integer>();
            treeCounts.put(topology, counts);
        }
        if (counts.containsKey(uniqueNewick)) {
            counts.put(uniqueNewick, counts.get(uniqueNewick) + 1);
        } else {
            counts.put(uniqueNewick, 1);
        }
    }
    TestCase.assertEquals(numLabelledTopologies, uniqueTrees.size());
    TestCase.assertEquals(sampleSize, chainLength / sampleTreeEvery + 1);
    Set<String> keys = topologies.keySet();
    double ep = 1.0 / topologies.size();
    for (String topology : keys) {
        double ap = ((double) topologies.get(topology)) / (sampleSize);
        //          	assertExpectation(ep, ap, sampleSize);
        HashMap<String, Integer> counts = treeCounts.get(topology);
        Set<String> trees = counts.keySet();
        double MSE = 0;
        double ep1 = 1.0 / counts.size();
        for (String t : trees) {
            double ap1 = ((double) counts.get(t)) / (topologies.get(topology));
            //              	assertExpectation(ep1, ap1, topologies.get(topology));
            MSE += (ep1 - ap1) * (ep1 - ap1);
        }
        MSE /= counts.size();
        System.out.println("The Mean Square Error for the topolgy " + topology + " is " + MSE);
    }
}

Example 4 with SpeciationLikelihood

use of dr.evomodel.speciation.SpeciationLikelihood in project beast-mcmc by beast-dev.

the class BirthDeathLikelihoodTest method birthDeathLikelihoodTester.

private void birthDeathLikelihoodTester(Tree tree, double birthRate, double deathRate, double logL) {
    Parameter b = new Parameter.Default("b", birthRate, 0.0, Double.MAX_VALUE);
    Parameter d = new Parameter.Default("d", deathRate, 0.0, Double.MAX_VALUE);
    SpeciationModel speciationModel = new BirthDeathGernhard08Model(b, d, null, BirthDeathGernhard08Model.TreeType.ORIENTED, Units.Type.YEARS);
    Likelihood likelihood = new SpeciationLikelihood(tree, speciationModel, "bd.like");
    assertEquals(logL, likelihood.getLogLikelihood(), 1e-14);
}

Example 5 with SpeciationLikelihood

use of dr.evomodel.speciation.SpeciationLikelihood in project beast-mcmc by beast-dev.

the class BirthDeathSSLikelihoodTest method likelihoodTester.

private void likelihoodTester(Tree tree, double birthRate, double deathRate, Variable<Double> origin, double logL) {
    Variable<Double> b = new Variable.D("b", birthRate);
    Variable<Double> d = new Variable.D("d", deathRate);
    Variable<Double> psi = new Variable.D("psi", this.psi);
    Variable<Double> p = new Variable.D("p", this.p);
    Variable<Double> r = new Variable.D("r", 0.5);
    Variable<Double> fTime = new Variable.D("time", 0.0);
    SpeciationModel speciationModel = new BirthDeathSerialSamplingModel(b, d, psi, p, false, r, true, origin, Units.Type.YEARS);
    Likelihood likelihood = new SpeciationLikelihood(tree, speciationModel, "bdss.like");
    assertEquals(logL, likelihood.getLogLikelihood());
}