Search in sources :

Example 1 with MCMCOptions

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

the class StrictClockTest method testStrictClock.

public void testStrictClock() throws Exception {
    Parameter popSize = new Parameter.Default(ConstantPopulationModelParser.POPULATION_SIZE, 380.0, 0, 38000.0);
    ConstantPopulationModel constantModel = createRandomInitialTree(popSize);
    CoalescentLikelihood coalescent = new CoalescentLikelihood(treeModel, null, new ArrayList<TaxonList>(), constantModel);
    coalescent.setId("coalescent");
    // clock model
    Parameter rateParameter = new Parameter.Default(StrictClockBranchRates.RATE, 2.3E-5, 0, 100.0);
    StrictClockBranchRates branchRateModel = new StrictClockBranchRates(rateParameter);
    // Sub model
    Parameter freqs = new Parameter.Default(alignment.getStateFrequencies());
    Parameter kappa = new Parameter.Default(HKYParser.KAPPA, 1.0, 0, 100.0);
    FrequencyModel f = new FrequencyModel(Nucleotides.INSTANCE, freqs);
    HKY hky = new HKY(kappa, f);
    //siteModel
    GammaSiteModel siteModel = new GammaSiteModel(hky);
    Parameter mu = new Parameter.Default(GammaSiteModelParser.MUTATION_RATE, 1.0, 0, Double.POSITIVE_INFINITY);
    siteModel.setMutationRateParameter(mu);
    //treeLikelihood
    SitePatterns patterns = new SitePatterns(alignment, null, 0, -1, 1, true);
    TreeLikelihood treeLikelihood = new TreeLikelihood(patterns, treeModel, siteModel, branchRateModel, null, false, false, true, false, false);
    treeLikelihood.setId(TreeLikelihoodParser.TREE_LIKELIHOOD);
    // Operators
    OperatorSchedule schedule = new SimpleOperatorSchedule();
    MCMCOperator operator = new ScaleOperator(kappa, 0.75);
    operator.setWeight(1.0);
    schedule.addOperator(operator);
    operator = new ScaleOperator(rateParameter, 0.75);
    operator.setWeight(3.0);
    schedule.addOperator(operator);
    Parameter allInternalHeights = treeModel.createNodeHeightsParameter(true, true, false);
    operator = new UpDownOperator(new Scalable[] { new Scalable.Default(rateParameter) }, new Scalable[] { new Scalable.Default(allInternalHeights) }, 0.75, 3.0, CoercionMode.COERCION_ON);
    schedule.addOperator(operator);
    operator = new ScaleOperator(popSize, 0.75);
    operator.setWeight(3.0);
    schedule.addOperator(operator);
    Parameter rootHeight = treeModel.getRootHeightParameter();
    rootHeight.setId(TREE_HEIGHT);
    operator = new ScaleOperator(rootHeight, 0.75);
    operator.setWeight(3.0);
    schedule.addOperator(operator);
    Parameter internalHeights = treeModel.createNodeHeightsParameter(false, true, false);
    operator = new UniformOperator(internalHeights, 30.0);
    schedule.addOperator(operator);
    operator = new SubtreeSlideOperator(treeModel, 15.0, 1.0, true, false, false, false, CoercionMode.COERCION_ON);
    schedule.addOperator(operator);
    operator = new ExchangeOperator(ExchangeOperator.NARROW, treeModel, 15.0);
    //        operator.doOperation();
    schedule.addOperator(operator);
    operator = new ExchangeOperator(ExchangeOperator.WIDE, treeModel, 3.0);
    //        operator.doOperation();
    schedule.addOperator(operator);
    operator = new WilsonBalding(treeModel, 3.0);
    //        operator.doOperation();
    schedule.addOperator(operator);
    //CompoundLikelihood
    List<Likelihood> likelihoods = new ArrayList<Likelihood>();
    likelihoods.add(coalescent);
    Likelihood prior = new CompoundLikelihood(0, likelihoods);
    prior.setId(CompoundLikelihoodParser.PRIOR);
    likelihoods.clear();
    likelihoods.add(treeLikelihood);
    Likelihood likelihood = new CompoundLikelihood(-1, likelihoods);
    likelihoods.clear();
    likelihoods.add(prior);
    likelihoods.add(likelihood);
    Likelihood posterior = new CompoundLikelihood(0, likelihoods);
    posterior.setId(CompoundLikelihoodParser.POSTERIOR);
    // Log
    ArrayLogFormatter formatter = new ArrayLogFormatter(false);
    MCLogger[] loggers = new MCLogger[2];
    loggers[0] = new MCLogger(formatter, 500, false);
    loggers[0].add(posterior);
    loggers[0].add(treeLikelihood);
    loggers[0].add(rootHeight);
    loggers[0].add(rateParameter);
    loggers[0].add(popSize);
    loggers[0].add(kappa);
    loggers[0].add(coalescent);
    loggers[1] = new MCLogger(new TabDelimitedFormatter(System.out), 10000, false);
    loggers[1].add(posterior);
    loggers[1].add(treeLikelihood);
    loggers[1].add(rootHeight);
    loggers[1].add(rateParameter);
    loggers[1].add(coalescent);
    // MCMC
    MCMC mcmc = new MCMC("mcmc1");
    MCMCOptions options = new MCMCOptions(1000000);
    mcmc.setShowOperatorAnalysis(true);
    mcmc.init(options, posterior, schedule, loggers);
    mcmc.run();
    // time
    System.out.println(mcmc.getTimer().toString());
    // Tracer
    List<Trace> traces = formatter.getTraces();
    ArrayTraceList traceList = new ArrayTraceList("RandomLocalClockTest", traces, 0);
    for (int i = 1; i < traces.size(); i++) {
        traceList.analyseTrace(i);
    }
    //        <expectation name="posterior" value="-3928.71"/>
    //        <expectation name="clock.rate" value="8.04835E-4"/>
    //        <expectation name="constant.popSize" value="37.3762"/>
    //        <expectation name="hky.kappa" value="18.2782"/>
    //        <expectation name="treeModel.rootHeight" value="69.0580"/>
    //        <expectation name="treeLikelihood" value="-3856.59"/>
    //        <expectation name="coalescent" value="-72.1285"/>
    TraceCorrelation likelihoodStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(CompoundLikelihoodParser.POSTERIOR));
    assertExpectation(CompoundLikelihoodParser.POSTERIOR, likelihoodStats, -3928.71);
    likelihoodStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(TreeLikelihoodParser.TREE_LIKELIHOOD));
    assertExpectation(TreeLikelihoodParser.TREE_LIKELIHOOD, likelihoodStats, -3856.59);
    TraceCorrelation treeHeightStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(TREE_HEIGHT));
    assertExpectation(TREE_HEIGHT, treeHeightStats, 69.0580);
    TraceCorrelation kappaStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(HKYParser.KAPPA));
    assertExpectation(HKYParser.KAPPA, kappaStats, 18.2782);
    TraceCorrelation rateStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(StrictClockBranchRates.RATE));
    assertExpectation(StrictClockBranchRates.RATE, rateStats, 8.04835E-4);
    TraceCorrelation popStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(ConstantPopulationModelParser.POPULATION_SIZE));
    assertExpectation(ConstantPopulationModelParser.POPULATION_SIZE, popStats, 37.3762);
    TraceCorrelation coalescentStats = traceList.getCorrelationStatistics(traceList.getTraceIndex("coalescent"));
    assertExpectation("coalescent", coalescentStats, -72.1285);
}
Also used : FrequencyModel(dr.oldevomodel.substmodel.FrequencyModel) CompoundLikelihood(dr.inference.model.CompoundLikelihood) Likelihood(dr.inference.model.Likelihood) TreeLikelihood(dr.oldevomodel.treelikelihood.TreeLikelihood) CoalescentLikelihood(dr.evomodel.coalescent.CoalescentLikelihood) TreeLikelihood(dr.oldevomodel.treelikelihood.TreeLikelihood) ExchangeOperator(dr.evomodel.operators.ExchangeOperator) ArrayList(java.util.ArrayList) MCMC(dr.inference.mcmc.MCMC) SubtreeSlideOperator(dr.evomodel.operators.SubtreeSlideOperator) StrictClockBranchRates(dr.evomodel.branchratemodel.StrictClockBranchRates) CoalescentLikelihood(dr.evomodel.coalescent.CoalescentLikelihood) MCMCOptions(dr.inference.mcmc.MCMCOptions) ArrayLogFormatter(dr.inference.loggers.ArrayLogFormatter) WilsonBalding(dr.evomodel.operators.WilsonBalding) TraceCorrelation(dr.inference.trace.TraceCorrelation) SitePatterns(dr.evolution.alignment.SitePatterns) ConstantPopulationModel(dr.evomodel.coalescent.ConstantPopulationModel) TaxonList(dr.evolution.util.TaxonList) CompoundLikelihood(dr.inference.model.CompoundLikelihood) TabDelimitedFormatter(dr.inference.loggers.TabDelimitedFormatter) Trace(dr.inference.trace.Trace) GammaSiteModel(dr.oldevomodel.sitemodel.GammaSiteModel) ArrayTraceList(dr.inference.trace.ArrayTraceList) HKY(dr.oldevomodel.substmodel.HKY) Parameter(dr.inference.model.Parameter) MCLogger(dr.inference.loggers.MCLogger)

Example 2 with MCMCOptions

use of dr.inference.mcmc.MCMCOptions 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);
}
Also used : TraceCorrelation(dr.inference.trace.TraceCorrelation) Likelihood(dr.inference.model.Likelihood) SpeciationLikelihood(dr.evomodel.speciation.SpeciationLikelihood) MCMC(dr.inference.mcmc.MCMC) SpeciationModel(dr.evomodel.speciation.SpeciationModel) TabDelimitedFormatter(dr.inference.loggers.TabDelimitedFormatter) SpeciationLikelihood(dr.evomodel.speciation.SpeciationLikelihood) Trace(dr.inference.trace.Trace) RandomLocalYuleModel(dr.evomodel.speciation.RandomLocalYuleModel) ArrayTraceList(dr.inference.trace.ArrayTraceList) MCMCOptions(dr.inference.mcmc.MCMCOptions) TreeLengthStatistic(dr.evomodel.tree.TreeLengthStatistic) TreeHeightStatistic(dr.evomodel.tree.TreeHeightStatistic) Parameter(dr.inference.model.Parameter) ArrayLogFormatter(dr.inference.loggers.ArrayLogFormatter) MCLogger(dr.inference.loggers.MCLogger)

Example 3 with MCMCOptions

use of dr.inference.mcmc.MCMCOptions 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);
}
Also used : TraceCorrelation(dr.inference.trace.TraceCorrelation) Likelihood(dr.inference.model.Likelihood) SpeciationLikelihood(dr.evomodel.speciation.SpeciationLikelihood) BirthDeathGernhard08Model(dr.evomodel.speciation.BirthDeathGernhard08Model) MCMC(dr.inference.mcmc.MCMC) SpeciationModel(dr.evomodel.speciation.SpeciationModel) TabDelimitedFormatter(dr.inference.loggers.TabDelimitedFormatter) SpeciationLikelihood(dr.evomodel.speciation.SpeciationLikelihood) Trace(dr.inference.trace.Trace) ArrayTraceList(dr.inference.trace.ArrayTraceList) MCMCOptions(dr.inference.mcmc.MCMCOptions) TreeLengthStatistic(dr.evomodel.tree.TreeLengthStatistic) TreeHeightStatistic(dr.evomodel.tree.TreeHeightStatistic) Parameter(dr.inference.model.Parameter) ArrayLogFormatter(dr.inference.loggers.ArrayLogFormatter) MCLogger(dr.inference.loggers.MCLogger)

Example 4 with MCMCOptions

use of dr.inference.mcmc.MCMCOptions 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);
    }
}
Also used : HashMap(java.util.HashMap) Likelihood(dr.inference.model.Likelihood) SpeciationLikelihood(dr.evomodel.speciation.SpeciationLikelihood) BirthDeathGernhard08Model(dr.evomodel.speciation.BirthDeathGernhard08Model) MCMC(dr.inference.mcmc.MCMC) SpeciationLikelihood(dr.evomodel.speciation.SpeciationLikelihood) TreeModel(dr.evomodel.tree.TreeModel) TreeLogger(dr.evomodel.tree.TreeLogger) MCMCOptions(dr.inference.mcmc.MCMCOptions) TreeLengthStatistic(dr.evomodel.tree.TreeLengthStatistic) FlexibleTree(dr.evolution.tree.FlexibleTree) Tree(dr.evolution.tree.Tree) FileReader(java.io.FileReader) HashSet(java.util.HashSet) NexusImporter(dr.evolution.io.NexusImporter) OperatorSchedule(dr.inference.operators.OperatorSchedule) SpeciationModel(dr.evomodel.speciation.SpeciationModel) TabDelimitedFormatter(dr.inference.loggers.TabDelimitedFormatter) FileOutputStream(java.io.FileOutputStream) TreeHeightStatistic(dr.evomodel.tree.TreeHeightStatistic) Parameter(dr.inference.model.Parameter) File(java.io.File) MCLogger(dr.inference.loggers.MCLogger)

Example 5 with MCMCOptions

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

the class LognormalPriorTest method testLognormalPrior.

public void testLognormalPrior() {
    //        ConstantPopulation constant = new ConstantPopulation(Units.Type.YEARS);
    //        constant.setN0(popSize); // popSize
    Parameter popSize = new Parameter.Default(6.0);
    popSize.setId(ConstantPopulationModelParser.POPULATION_SIZE);
    ConstantPopulationModel demo = new ConstantPopulationModel(popSize, Units.Type.YEARS);
    //Likelihood
    Likelihood dummyLikelihood = new DummyLikelihood(demo);
    // Operators
    OperatorSchedule schedule = new SimpleOperatorSchedule();
    MCMCOperator operator = new ScaleOperator(popSize, 0.75);
    operator.setWeight(1.0);
    schedule.addOperator(operator);
    // Log
    ArrayLogFormatter formatter = new ArrayLogFormatter(false);
    MCLogger[] loggers = new MCLogger[2];
    loggers[0] = new MCLogger(formatter, 1000, false);
    //        loggers[0].add(treeLikelihood);
    loggers[0].add(popSize);
    loggers[1] = new MCLogger(new TabDelimitedFormatter(System.out), 100000, false);
    //        loggers[1].add(treeLikelihood);
    loggers[1].add(popSize);
    // MCMC
    MCMC mcmc = new MCMC("mcmc1");
    MCMCOptions options = new MCMCOptions(1000000);
    // meanInRealSpace="false"
    DistributionLikelihood logNormalLikelihood = new DistributionLikelihood(new LogNormalDistribution(1.0, 1.0), 0);
    logNormalLikelihood.addData(popSize);
    List<Likelihood> likelihoods = new ArrayList<Likelihood>();
    likelihoods.add(logNormalLikelihood);
    Likelihood prior = new CompoundLikelihood(0, likelihoods);
    likelihoods.clear();
    likelihoods.add(dummyLikelihood);
    Likelihood likelihood = new CompoundLikelihood(-1, likelihoods);
    likelihoods.clear();
    likelihoods.add(prior);
    likelihoods.add(likelihood);
    Likelihood posterior = new CompoundLikelihood(0, likelihoods);
    mcmc.setShowOperatorAnalysis(true);
    mcmc.init(options, posterior, schedule, loggers);
    mcmc.run();
    // time
    System.out.println(mcmc.getTimer().toString());
    // Tracer
    List<Trace> traces = formatter.getTraces();
    ArrayTraceList traceList = new ArrayTraceList("LognormalPriorTest", traces, 0);
    for (int i = 1; i < traces.size(); i++) {
        traceList.analyseTrace(i);
    }
    //      <expectation name="param" value="4.48168907"/>
    TraceCorrelation popSizeStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(ConstantPopulationModelParser.POPULATION_SIZE));
    System.out.println("Expectation of Log-Normal(1,1) is e^(M+S^2/2) = e^(1.5) = " + Math.exp(1.5));
    assertExpectation(ConstantPopulationModelParser.POPULATION_SIZE, popSizeStats, Math.exp(1.5));
}
Also used : CompoundLikelihood(dr.inference.model.CompoundLikelihood) Likelihood(dr.inference.model.Likelihood) DistributionLikelihood(dr.inference.distribution.DistributionLikelihood) DummyLikelihood(dr.inference.model.DummyLikelihood) MCMC(dr.inference.mcmc.MCMC) ArrayList(java.util.ArrayList) LogNormalDistribution(dr.math.distributions.LogNormalDistribution) MCMCOptions(dr.inference.mcmc.MCMCOptions) DummyLikelihood(dr.inference.model.DummyLikelihood) ArrayLogFormatter(dr.inference.loggers.ArrayLogFormatter) TraceCorrelation(dr.inference.trace.TraceCorrelation) ConstantPopulationModel(dr.evomodel.coalescent.ConstantPopulationModel) OperatorSchedule(dr.inference.operators.OperatorSchedule) SimpleOperatorSchedule(dr.inference.operators.SimpleOperatorSchedule) CompoundLikelihood(dr.inference.model.CompoundLikelihood) TabDelimitedFormatter(dr.inference.loggers.TabDelimitedFormatter) Trace(dr.inference.trace.Trace) SimpleOperatorSchedule(dr.inference.operators.SimpleOperatorSchedule) ArrayTraceList(dr.inference.trace.ArrayTraceList) Parameter(dr.inference.model.Parameter) ScaleOperator(dr.inference.operators.ScaleOperator) DistributionLikelihood(dr.inference.distribution.DistributionLikelihood) MCMCOperator(dr.inference.operators.MCMCOperator) MCLogger(dr.inference.loggers.MCLogger)

Aggregations

MCMC (dr.inference.mcmc.MCMC)12 MCMCOptions (dr.inference.mcmc.MCMCOptions)12 MCLogger (dr.inference.loggers.MCLogger)11 ArrayLogFormatter (dr.inference.loggers.ArrayLogFormatter)10 TabDelimitedFormatter (dr.inference.loggers.TabDelimitedFormatter)10 ArrayTraceList (dr.inference.trace.ArrayTraceList)10 Trace (dr.inference.trace.Trace)10 TraceCorrelation (dr.inference.trace.TraceCorrelation)10 Likelihood (dr.inference.model.Likelihood)9 Parameter (dr.inference.model.Parameter)9 ArrayList (java.util.ArrayList)7 CompoundLikelihood (dr.inference.model.CompoundLikelihood)6 SitePatterns (dr.evolution.alignment.SitePatterns)5 ExchangeOperator (dr.evomodel.operators.ExchangeOperator)5 SubtreeSlideOperator (dr.evomodel.operators.SubtreeSlideOperator)5 WilsonBalding (dr.evomodel.operators.WilsonBalding)5 GammaSiteModel (dr.oldevomodel.sitemodel.GammaSiteModel)5 FrequencyModel (dr.oldevomodel.substmodel.FrequencyModel)5 TreeLikelihood (dr.oldevomodel.treelikelihood.TreeLikelihood)5 ConstantPopulationModel (dr.evomodel.coalescent.ConstantPopulationModel)4