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