Search in sources :

Example 11 with MCMCOptions

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

the class MCMCTest method testMCMC.

public void testMCMC() {
    // Sub model
    //new double[]{0.25, 0.25, 0.25, 0.25});
    Parameter freqs = new Parameter.Default(alignment.getStateFrequencies());
    Parameter kappa = new Parameter.Default(HKYParser.KAPPA, 1.0, 1.0E-8, Double.POSITIVE_INFINITY);
    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, null, null, false, false, true, false, false);
    treeLikelihood.setId(TreeLikelihoodParser.TREE_LIKELIHOOD);
    // Operators
    OperatorSchedule schedule = new SimpleOperatorSchedule();
    MCMCOperator operator = new ScaleOperator(kappa, 0.5);
    operator.setWeight(1.0);
    schedule.addOperator(operator);
    //        Parameter rootParameter = treeModel.createNodeHeightsParameter(true, false, false);
    //        ScaleOperator scaleOperator = new ScaleOperator(rootParameter, 0.75, CoercionMode.COERCION_ON, 1.0);
    Parameter rootHeight = treeModel.getRootHeightParameter();
    rootHeight.setId(TREE_HEIGHT);
    operator = new ScaleOperator(rootHeight, 0.5);
    operator.setWeight(1.0);
    schedule.addOperator(operator);
    Parameter internalHeights = treeModel.createNodeHeightsParameter(false, true, false);
    operator = new UniformOperator(internalHeights, 10.0);
    schedule.addOperator(operator);
    operator = new SubtreeSlideOperator(treeModel, 1, 1, true, false, false, false, CoercionMode.COERCION_ON);
    schedule.addOperator(operator);
    operator = new ExchangeOperator(ExchangeOperator.NARROW, treeModel, 1.0);
    //        operator.doOperation();
    schedule.addOperator(operator);
    operator = new ExchangeOperator(ExchangeOperator.WIDE, treeModel, 1.0);
    //        operator.doOperation();
    schedule.addOperator(operator);
    operator = new WilsonBalding(treeModel, 1.0);
    //        operator.doOperation();
    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(rootHeight);
    loggers[0].add(kappa);
    loggers[1] = new MCLogger(new TabDelimitedFormatter(System.out), 100000, false);
    loggers[1].add(treeLikelihood);
    loggers[1].add(rootHeight);
    loggers[1].add(kappa);
    // MCMC
    MCMC mcmc = new MCMC("mcmc1");
    MCMCOptions options = new MCMCOptions(10000000);
    mcmc.setShowOperatorAnalysis(true);
    mcmc.init(options, treeLikelihood, schedule, loggers);
    mcmc.run();
    // time
    System.out.println(mcmc.getTimer().toString());
    // Tracer
    List<Trace> traces = formatter.getTraces();
    ArrayTraceList traceList = new ArrayTraceList("MCMCTest", traces, 0);
    for (int i = 1; i < traces.size(); i++) {
        traceList.analyseTrace(i);
    }
    //      <expectation name="likelihood" value="-1815.75"/>
    //		<expectation name="treeModel.rootHeight" value="6.42048E-2"/>
    //		<expectation name="hky.kappa" value="32.8941"/>
    TraceCorrelation likelihoodStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(TreeLikelihoodParser.TREE_LIKELIHOOD));
    assertExpectation(TreeLikelihoodParser.TREE_LIKELIHOOD, likelihoodStats, -1815.75);
    TraceCorrelation treeHeightStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(TREE_HEIGHT));
    assertExpectation(TREE_HEIGHT, treeHeightStats, 6.42048E-2);
    TraceCorrelation kappaStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(HKYParser.KAPPA));
    assertExpectation(HKYParser.KAPPA, kappaStats, 32.8941);
}
Also used : FrequencyModel(dr.oldevomodel.substmodel.FrequencyModel) TreeLikelihood(dr.oldevomodel.treelikelihood.TreeLikelihood) ExchangeOperator(dr.evomodel.operators.ExchangeOperator) MCMC(dr.inference.mcmc.MCMC) SubtreeSlideOperator(dr.evomodel.operators.SubtreeSlideOperator) 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) 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 12 with MCMCOptions

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

the class PMDTestProblem method testPMD.

public void testPMD() throws Exception {
    Parameter popSize = new Parameter.Default(ConstantPopulationModelParser.POPULATION_SIZE, 496432.69917113904, 0, Double.POSITIVE_INFINITY);
    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, 4.0E-7, 0, 100.0);
    StrictClockBranchRates branchRateModel = new StrictClockBranchRates(rateParameter);
    // Sub model
    Parameter freqs = new Parameter.Default(new double[] { 0.25, 0.25, 0.25, 0.25 });
    Parameter kappa = new Parameter.Default(HKYParser.KAPPA, 1.0, 1.0E-8, Double.POSITIVE_INFINITY);
    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);
    // SequenceErrorModel
    Parameter ageRelatedRateParameter = new Parameter.Default(SequenceErrorModelParser.AGE_RELATED_RATE, 4.0E-7, 0, 100.0);
    TipStatesModel aDNADamageModel = new SequenceErrorModel(null, null, SequenceErrorModel.ErrorType.TRANSITIONS_ONLY, null, ageRelatedRateParameter, null);
    //treeLikelihood
    SitePatterns patterns = new SitePatterns(alignment, null, 0, -1, 1, true);
    TreeLikelihood treeLikelihood = new TreeLikelihood(patterns, treeModel, siteModel, branchRateModel, aDNADamageModel, 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);
    operator = new ScaleOperator(ageRelatedRateParameter, 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, 49643.2699171139, 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);
    // ??? correct?
    operator = new DeltaExchangeOperator(freqs, new int[] { 1, 1, 1, 1 }, 0.01, 1.0, false, CoercionMode.COERCION_ON);
    schedule.addOperator(operator);
    //CompoundLikelihood
    OneOnXPrior likelihood1 = new OneOnXPrior();
    likelihood1.addData(popSize);
    OneOnXPrior likelihood2 = new OneOnXPrior();
    likelihood2.addData(kappa);
    List<Likelihood> likelihoods = new ArrayList<Likelihood>();
    likelihoods.add(likelihood1);
    likelihoods.add(likelihood2);
    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, 1000, false);
    loggers[0].add(posterior);
    loggers[0].add(treeLikelihood);
    loggers[0].add(rootHeight);
    loggers[0].add(rateParameter);
    loggers[0].add(ageRelatedRateParameter);
    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);
    // 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("PMDTest", traces, 0);
    for (int i = 1; i < traces.size(); i++) {
        traceList.analyseTrace(i);
    }
    //        <expectation name="clock.rate" value="1.5E-7"/>
    //        <expectation name="errorModel.ageRate" value="0.7E-7"/>
    //        <expectation name="hky.kappa" value="10"/>
    TraceCorrelation kappaStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(HKYParser.KAPPA));
    assertExpectation(HKYParser.KAPPA, kappaStats, 10);
    TraceCorrelation rateStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(StrictClockBranchRates.RATE));
    assertExpectation(StrictClockBranchRates.RATE, rateStats, 1.5E-7);
    TraceCorrelation ageRateStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(SequenceErrorModelParser.AGE_RELATED_RATE));
    assertExpectation(SequenceErrorModelParser.AGE_RELATED_RATE, ageRateStats, 0.7E-7);
}
Also used : FrequencyModel(dr.oldevomodel.substmodel.FrequencyModel) OneOnXPrior(dr.inference.model.OneOnXPrior) 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) SequenceErrorModel(dr.evomodel.tipstatesmodel.SequenceErrorModel) 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) TipStatesModel(dr.evomodel.tipstatesmodel.TipStatesModel) 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)

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