Examples with MCMC dr.inference.mcmc.MCMC used on opensource projects

Example 1 with MCMC

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

the class Tutorial1 method main.

public static void main(String[] arg) throws IOException, TraceException {
    // constructing random variable representing mean of normal distribution
    Variable.D mean = new Variable.D("mean", 1.0);
    // give mean a uniform prior [-1000, 1000]
    mean.addBounds(new Parameter.DefaultBounds(1000, -1000, 1));
    // constructing random variable representing stdev of normal distribution
    Variable.D stdev = new Variable.D("stdev", 1.0);
    // give stdev a uniform prior [0, 1000]
    stdev.addBounds(new Parameter.DefaultBounds(1000, 0, 1));
    // construct normal distribution model
    NormalDistributionModel normal = new NormalDistributionModel(mean, stdev);
    // construct a likelihood for normal distribution
    DistributionLikelihood likelihood = new DistributionLikelihood(normal);
    // construct data
    Attribute.Default<double[]> d = new Attribute.Default<double[]>("x", new double[] { 1, 2, 3, 4, 5, 6, 7, 8, 9 });
    // add data (representing 9 independent observations) to likelihood
    likelihood.addData(d);
    // construct two "operators" to be used as the proposal distribution
    MCMCOperator meanMove = new ScaleOperator(mean, 0.75);
    MCMCOperator stdevMove = new ScaleOperator(stdev, 0.75);
    // construct a logger to log progress of MCMC run to stdout (screen)
    MCLogger logger1 = new MCLogger(100);
    logger1.add(mean);
    logger1.add(stdev);
    // construct a logger to log to a log file for later analysis
    MCLogger logger2 = new MCLogger("tutorial1.log", 100, false, 0);
    logger2.add(mean);
    logger2.add(stdev);
    // construct MCMC object
    MCMC mcmc = new MCMC("tutorial1:normal");
    // initialize MCMC with chain length, likelihood, operators and loggers
    mcmc.init(100000, likelihood, new MCMCOperator[] { meanMove, stdevMove }, new Logger[] { logger1, logger2 });
    // run the mcmc
    mcmc.chain();
    // perform post-analysis
    TraceAnalysis.report("tutorial1.log");
}

Example 2 with MCMC

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

the class MCMCParser method parseMCMC.

/**
 * Parse the MCMC object.
 * @param xo the XML object
 * @return the MXMX object
 * @throws XMLParseException an exception of there is an XML parse error
 */
private MCMC parseMCMC(String id, XMLObject xo) throws XMLParseException {
    MCMC mcmc = new MCMC(id);
    long chainLength = xo.getLongIntegerAttribute(CHAIN_LENGTH);
    boolean useAdaptation = xo.getAttribute(ADAPTATION, true) || xo.getAttribute(AUTO_OPTIMIZE, true);
    if (System.getProperty("mcmc.use_adaptation") != null) {
        useAdaptation = Boolean.parseBoolean(System.getProperty("mcmc.use_adaptation"));
    }
    long adaptationDelay = chainLength / 100;
    adaptationDelay = xo.getAttribute(ADAPTATION_DELAY, xo.getAttribute(AUTO_OPTIMIZE_DELAY, xo.getAttribute(PRE_BURNIN, adaptationDelay)));
    double adaptationTarget = 0.234;
    if (System.getProperty("mcmc.adaptation_target") != null) {
        adaptationTarget = Double.parseDouble(System.getProperty("mcmc.adaptation_target"));
    }
    boolean useSmoothAcceptanceRatio = xo.getAttribute(SMOOTHED_ACCEPTANCE_RATIO, false);
    double temperature = xo.getAttribute(TEMPERATURE, 1.0);
    long fullEvaluationCount = xo.getAttribute(FULL_EVALUATION, DEFAULT_FULL_EVALUATION_COUNT);
    if (System.getProperty("mcmc.evaluation.count") != null) {
        fullEvaluationCount = Long.parseLong(System.getProperty("mcmc.evaluation.count"));
    }
    double evaluationTestThreshold = MarkovChain.EVALUATION_TEST_THRESHOLD;
    evaluationTestThreshold = xo.getAttribute(EVALUATION_THRESHOLD, evaluationTestThreshold);
    if (System.getProperty("mcmc.evaluation.threshold") != null) {
        evaluationTestThreshold = Double.parseDouble(System.getProperty("mcmc.evaluation.threshold"));
    }
    int minOperatorCountForFullEvaluation = xo.getAttribute(MIN_OPS_EVALUATIONS, 1);
    MCMCOptions options = new MCMCOptions(chainLength, fullEvaluationCount, minOperatorCountForFullEvaluation, evaluationTestThreshold, useAdaptation, adaptationDelay, adaptationTarget, useSmoothAcceptanceRatio, temperature);
    OperatorSchedule opsched = (OperatorSchedule) xo.getChild(OperatorSchedule.class);
    Likelihood likelihood = (Likelihood) xo.getChild(Likelihood.class);
    likelihood.setUsed();
    if (Boolean.valueOf(System.getProperty("show_warnings", "false"))) {
        // check that all models, parameters and likelihoods are being used
        for (Likelihood l : Likelihood.FULL_LIKELIHOOD_SET) {
            if (!l.isUsed()) {
                java.util.logging.Logger.getLogger("dr.inference").warning("Likelihood, " + l.getId() + ", of class " + l.getClass().getName() + " is not being handled by the MCMC.");
            }
        }
        for (Model m : Model.FULL_MODEL_SET) {
            if (!m.isUsed()) {
                java.util.logging.Logger.getLogger("dr.inference").warning("Model, " + m.getId() + ", of class " + m.getClass().getName() + " is not being handled by the MCMC.");
            }
        }
        for (Parameter p : Parameter.FULL_PARAMETER_SET) {
            if (!p.isUsed()) {
                java.util.logging.Logger.getLogger("dr.inference").warning("Parameter, " + p.getId() + ", of class " + p.getClass().getName() + " is not being handled by the MCMC.");
            }
        }
    }
    ArrayList<Logger> loggers = new ArrayList<Logger>();
    for (int i = 0; i < xo.getChildCount(); i++) {
        Object child = xo.getChild(i);
        if (child instanceof Logger) {
            loggers.add((Logger) child);
        }
    }
    mcmc.setShowOperatorAnalysis(true);
    if (xo.hasAttribute(OPERATOR_ANALYSIS)) {
        mcmc.setOperatorAnalysisFile(XMLParser.getLogFile(xo, OPERATOR_ANALYSIS));
    }
    Logger[] loggerArray = new Logger[loggers.size()];
    loggers.toArray(loggerArray);
    java.util.logging.Logger.getLogger("dr.inference").info("\nCreating the MCMC chain:" + "\n  chain length = " + options.getChainLength() + "\n  operator adaption = " + options.useAdaptation() + (options.useAdaptation() ? "\n  adaptation delayed for " + options.getAdaptationDelay() + " steps" : "") + (options.getFullEvaluationCount() == 0 ? "\n  full evaluation test off" : ""));
    mcmc.init(options, likelihood, opsched, loggerArray);
    MarkovChain mc = mcmc.getMarkovChain();
    double initialScore = mc.getCurrentScore();
    if (initialScore == Double.NEGATIVE_INFINITY) {
        String message = "The initial posterior is zero";
        if (likelihood instanceof CompoundLikelihood) {
            message += ": " + ((CompoundLikelihood) likelihood).getDiagnosis(2);
        } else {
            message += "!";
        }
        throw new IllegalArgumentException(message);
    }
    if (!xo.getAttribute(SPAWN, true))
        mcmc.setSpawnable(false);
    return mcmc;
}

Example 3 with MCMC

use of dr.inference.mcmc.MCMC 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 4 with MCMC

use of dr.inference.mcmc.MCMC 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 5 with MCMC

use of dr.inference.mcmc.MCMC 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);
    TreeIntervals intervalList = new TreeIntervals(treeModel, null, null);
    CoalescentLikelihood coalescent = new CoalescentLikelihood(intervalList, 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 = ((DefaultTreeModel) treeModel).createNodeHeightsParameter(true, true, false);
    operator = new UpDownOperator(new Scalable[] { new Scalable.Default(rateParameter) }, new Scalable[] { new Scalable.Default(allInternalHeights) }, 0.75, 3.0, AdaptationMode.ADAPTATION_ON);
    schedule.addOperator(operator);
    operator = new ScaleOperator(popSize, 0.75);
    operator.setWeight(3.0);
    schedule.addOperator(operator);
    Parameter rootHeight = ((DefaultTreeModel) treeModel).getRootHeightParameter();
    rootHeight.setId(TREE_HEIGHT);
    operator = new ScaleOperator(rootHeight, 0.75);
    operator.setWeight(3.0);
    schedule.addOperator(operator);
    Parameter internalHeights = ((DefaultTreeModel) treeModel).createNodeHeightsParameter(false, true, false);
    operator = new UniformOperator(internalHeights, 30.0);
    schedule.addOperator(operator);
    operator = new SubtreeSlideOperator(((DefaultTreeModel) treeModel), 15.0, 1.0, true, false, false, false, AdaptationMode.ADAPTATION_ON, AdaptableMCMCOperator.DEFAULT_ADAPTATION_TARGET);
    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);
}