Examples with BooleanLikelihood dr.inference.model.BooleanLikelihood used on opensource projects

Example 1 with BooleanLikelihood

use of dr.inference.model.BooleanLikelihood in project beast-mcmc by beast-dev.

the class TestCalibratedYuleModel method yuleTester.

private void yuleTester(TreeModel treeModel, OperatorSchedule schedule, Parameter brParameter, double S, int chainLength) throws IOException, TreeUtils.MissingTaxonException {
    MCMC mcmc = new MCMC("mcmc1");
    MCMCOptions options = new MCMCOptions(chainLength);
    TreeLengthStatistic tls = new TreeLengthStatistic(TL, treeModel);
    TreeHeightStatistic rootHeight = new TreeHeightStatistic(TREE_HEIGHT, treeModel);
    SpeciationModel speciationModel = new BirthDeathGernhard08Model("yule", brParameter, null, null, BirthDeathGernhard08Model.TreeType.UNSCALED, Units.Type.SUBSTITUTIONS, false);
    Likelihood speciationLikelihood = new SpeciationLikelihood(treeModel, speciationModel, "yule.like");
    Taxa halfTaxa = new Taxa();
    for (int i = 0; i < taxa.getTaxonCount() / 2; i++) {
        halfTaxa.addTaxon(new Taxon("T" + Integer.toString(i)));
    }
    TMRCAStatistic tmrca = new TMRCAStatistic("tmrca(halfTaxa)", treeModel, halfTaxa, false, false);
    DistributionLikelihood logNormalLikelihood = new DistributionLikelihood(new LogNormalDistribution(M, S), // meanInRealSpace="false"
    0);
    logNormalLikelihood.addData(tmrca);
    MonophylyStatistic monophylyStatistic = new MonophylyStatistic("monophyly(halfTaxa)", treeModel, halfTaxa, null);
    BooleanLikelihood booleanLikelihood = new BooleanLikelihood();
    booleanLikelihood.addData(monophylyStatistic);
    //CompoundLikelihood
    List<Likelihood> likelihoods = new ArrayList<Likelihood>();
    likelihoods.add(speciationLikelihood);
    likelihoods.add(logNormalLikelihood);
    likelihoods.add(booleanLikelihood);
    Likelihood prior = new CompoundLikelihood(0, likelihoods);
    prior.setId(CompoundLikelihoodParser.PRIOR);
    ArrayLogFormatter logformatter = new ArrayLogFormatter(false);
    MCLogger[] loggers = new MCLogger[1];
    loggers[0] = new MCLogger(logformatter, (int) (options.getChainLength() / 10000), false);
    loggers[0].add(speciationLikelihood);
    loggers[0].add(rootHeight);
    loggers[0].add(tmrca);
    loggers[0].add(tls);
    loggers[0].add(brParameter);
    mcmc.setShowOperatorAnalysis(false);
    mcmc.init(options, prior, schedule, loggers);
    mcmc.run();
    List<Trace> traces = logformatter.getTraces();
    ArrayTraceList traceList = new ArrayTraceList("yuleModelTest", traces, 1000);
    for (int i = 1; i < traces.size(); i++) {
        traceList.analyseTrace(i);
    }
    NumberFormatter formatter = new NumberFormatter(8);
    TraceCorrelation tlStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(TL));
    TraceCorrelation treeHeightStats = traceList.getCorrelationStatistics(traceList.getTraceIndex("tmrca(halfTaxa)"));
    //        out.write("tmrcaHeight = \t");
    out.write(formatter.format(treeHeightStats.getMean()));
    out.write("\t");
    double expectedNodeHeight = Math.pow(Math.E, (M + (Math.pow(S, 2) / 2)));
    //        out.write("expectation = \t");
    out.write(formatter.format(expectedNodeHeight));
    out.write("\t");
    double error = Math.abs((treeHeightStats.getMean() - expectedNodeHeight) / expectedNodeHeight);
    NumberFormat percentFormatter = NumberFormat.getNumberInstance();
    percentFormatter.setMinimumFractionDigits(5);
    percentFormatter.setMinimumFractionDigits(5);
    //        out.write("error = \t");
    out.write(percentFormatter.format(error));
    out.write("\t");
    //        out.write("tl.ess = \t");
    out.write(Double.toString(tlStats.getESS()));
    System.out.println("tmrcaHeight = " + formatter.format(treeHeightStats.getMean()) + ";  expectation = " + formatter.format(expectedNodeHeight) + ";  error = " + percentFormatter.format(error) + ";  tl.ess = " + tlStats.getESS());
}