Java MCMC類代碼示例

import dr.inference.mcmc.MCMC; //導入依賴的package包/類
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");
    }
 

import dr.inference.mcmc.MCMC; //導入依賴的package包/類
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();
        options.setChainLength(1000000);
        options.setUseCoercion(true); // autoOptimize = true
        options.setCoercionDelay(100);
        options.setTemperature(1.0);
        options.setFullEvaluationCount(2000);

        DistributionLikelihood logNormalLikelihood = new DistributionLikelihood(new LogNormalDistribution(1.0, 1.0), 0); // meanInRealSpace="false"
        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));
    }
 

import dr.inference.mcmc.MCMC; //導入依賴的package包/類
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);

        DistributionLikelihood logNormalLikelihood = new DistributionLikelihood(new LogNormalDistribution(1.0, 1.0), 0); // meanInRealSpace="false"
        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));
    }
 

import dr.inference.mcmc.MCMC; //導入依賴的package包/類
public MCMCMC(MCMC[] mcmcs, MCMCMCOptions mcmcmcOptions) {

        this.mcmcmcOptions = mcmcmcOptions;

        if (mcmcmcOptions.getChainTemperatures()[0] != 1.0) {
            throw new RuntimeException("The first chain in the array should be cold (temperature = 1.0)");
        }

        coldChain = 0;

        this.mcmcOptions = mcmcs[coldChain].getOptions();

        // Get all the loggers out of all the chains. We will only use the
        // loggers of the cold chain but we need to swap the formatters around
        // so that which every chain is cold always writes to the same destination.
        mcLoggers = new MCLogger[mcmcs.length][];
        for (int i = 0; i < mcmcs.length; i++) {
            Logger[] loggers = mcmcs[i].getLoggers();
            mcLoggers[i] = new MCLogger[loggers.length];
            for (int j = 0; j < loggers.length; j++) {
                mcLoggers[i][j] = (MCLogger) loggers[j];
            }
            if (mcLoggers[i] == null) {
                throw new RuntimeException("There are no loggers in the MCMC chains.");
            }
        }

        // Get all the operator schedules. The tuning values of these must be swapped
        // around as the temperatures are swapped.
        schedules = new OperatorSchedule[mcmcs.length];
        for (int i = 0; i < schedules.length; i++) {
            schedules[i] = mcmcs[i].getOperatorSchedule();
        }

        chains = new MarkovChain[mcmcs.length];

        chains[0] = mcmcs[0].getMarkovChain();
        for (int i = 1; i < chains.length; i++) {
            chains[i] = mcmcs[i].getMarkovChain();
            MCMCCriterion acceptor = ((MCMCCriterion) chains[i].getAcceptor());
            acceptor.setTemperature(mcmcmcOptions.getChainTemperatures()[i]);
        }

    }
 

import dr.inference.mcmc.MCMC; //導入依賴的package包/類
public String getParserName() {
    return MCMC;
}
 

import dr.inference.mcmc.MCMC; //導入依賴的package包/類
public Class getReturnType() {
    return MCMC.class;
}
 

import dr.inference.mcmc.MCMC; //導入依賴的package包/類
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);
    }
 

import dr.inference.mcmc.MCMC; //導入依賴的package包/類
private void yuleTester(TreeModel treeModel, OperatorSchedule schedule) {

        MCMC mcmc = new MCMC("mcmc1");
        MCMCOptions options = new MCMCOptions();
        options.setChainLength(1000000);
        options.setUseCoercion(true);
        options.setCoercionDelay(100);
        options.setTemperature(1.0);
        options.setFullEvaluationCount(2000);

        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);


    }
 

import dr.inference.mcmc.MCMC; //導入依賴的package包/類
private void randomLocalYuleTester(TreeModel treeModel, Parameter I, Parameter b, OperatorSchedule schedule) {

        MCMC mcmc = new MCMC("mcmc1");
        MCMCOptions options = new MCMCOptions();
        options.setChainLength(1000000);
        options.setUseCoercion(true);
        options.setCoercionDelay(100);
        options.setTemperature(1.0);
        options.setFullEvaluationCount(2000);

        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);
    }
 

import dr.inference.mcmc.MCMC; //導入依賴的package包/類
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);


    }