Search in sources :

Example 1 with LogNormalDistribution

use of dr.math.distributions.LogNormalDistribution 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)

Example 2 with LogNormalDistribution

use of dr.math.distributions.LogNormalDistribution 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());
}
Also used : BooleanLikelihood(dr.inference.model.BooleanLikelihood) CompoundLikelihood(dr.inference.model.CompoundLikelihood) Likelihood(dr.inference.model.Likelihood) DistributionLikelihood(dr.inference.distribution.DistributionLikelihood) SpeciationLikelihood(dr.evomodel.speciation.SpeciationLikelihood) BooleanLikelihood(dr.inference.model.BooleanLikelihood) BirthDeathGernhard08Model(dr.evomodel.speciation.BirthDeathGernhard08Model) MCMC(dr.inference.mcmc.MCMC) ArrayList(java.util.ArrayList) LogNormalDistribution(dr.math.distributions.LogNormalDistribution) SpeciationLikelihood(dr.evomodel.speciation.SpeciationLikelihood) Taxa(dr.evolution.util.Taxa) MCMCOptions(dr.inference.mcmc.MCMCOptions) ArrayLogFormatter(dr.inference.loggers.ArrayLogFormatter) TraceCorrelation(dr.inference.trace.TraceCorrelation) Taxon(dr.evolution.util.Taxon) CompoundLikelihood(dr.inference.model.CompoundLikelihood) SpeciationModel(dr.evomodel.speciation.SpeciationModel) Trace(dr.inference.trace.Trace) ArrayTraceList(dr.inference.trace.ArrayTraceList) DistributionLikelihood(dr.inference.distribution.DistributionLikelihood) MCLogger(dr.inference.loggers.MCLogger) NumberFormatter(dr.util.NumberFormatter) NumberFormat(java.text.NumberFormat)

Example 3 with LogNormalDistribution

use of dr.math.distributions.LogNormalDistribution in project beast-mcmc by beast-dev.

the class LogNormalDistributionTest method testCDFAndQuantile2.

public void testCDFAndQuantile2() {
    final LogNormalDistribution f = new LogNormalDistribution(1, 1);
    for (double i = 0.01; i < 0.95; i += 0.01) {
        final double y = i;
        BisectionZeroFinder zeroFinder = new BisectionZeroFinder(new OneVariableFunction() {

            public double value(double x) {
                return f.cdf(x) - y;
            }
        }, 0.01, 100);
        zeroFinder.evaluate();
        assertEquals(f.quantile(i), zeroFinder.getResult(), 1e-6);
    }
}
Also used : OneVariableFunction(dr.math.interfaces.OneVariableFunction) BisectionZeroFinder(dr.math.iterations.BisectionZeroFinder) LogNormalDistribution(dr.math.distributions.LogNormalDistribution)

Aggregations

LogNormalDistribution (dr.math.distributions.LogNormalDistribution)3 DistributionLikelihood (dr.inference.distribution.DistributionLikelihood)2 ArrayLogFormatter (dr.inference.loggers.ArrayLogFormatter)2 MCLogger (dr.inference.loggers.MCLogger)2 MCMC (dr.inference.mcmc.MCMC)2 MCMCOptions (dr.inference.mcmc.MCMCOptions)2 CompoundLikelihood (dr.inference.model.CompoundLikelihood)2 Likelihood (dr.inference.model.Likelihood)2 ArrayTraceList (dr.inference.trace.ArrayTraceList)2 Trace (dr.inference.trace.Trace)2 TraceCorrelation (dr.inference.trace.TraceCorrelation)2 ArrayList (java.util.ArrayList)2 Taxa (dr.evolution.util.Taxa)1 Taxon (dr.evolution.util.Taxon)1 ConstantPopulationModel (dr.evomodel.coalescent.ConstantPopulationModel)1 BirthDeathGernhard08Model (dr.evomodel.speciation.BirthDeathGernhard08Model)1 SpeciationLikelihood (dr.evomodel.speciation.SpeciationLikelihood)1 SpeciationModel (dr.evomodel.speciation.SpeciationModel)1 TabDelimitedFormatter (dr.inference.loggers.TabDelimitedFormatter)1 BooleanLikelihood (dr.inference.model.BooleanLikelihood)1