use of org.apache.commons.math.MaxIterationsExceededException in project beast-mcmc by beast-dev.
the class TwoStateOccupancyMarkovRewardsTest method run.
// public void testTwoStateSericolaRewards1() {
// final double rate = 0.0015;
//// final double prop = 0.5;
// final double prop = 0.66666;
//
// final double branchLength = 2000.0;
// final boolean print = false;
//
//// MarkovReward markovReward = createMarkovReward(rate, prop);
// MarkovReward markovReward = createSericolaMarkovReward(rate, prop);
//
// run(markovReward, rate, prop, branchLength, print, 1000);
// }
// public void testTwoStateSericolaRewards2() {
// final double rate = 0.0015;
// final double prop = 0.5;
//// final double prop = 0.66666;
// final double branchLength = 1000.0;
// final boolean print = false;
//
// MarkovReward markovReward = createMarkovReward(rate, prop);
//// MarkovReward markovReward = createSericolaMarkovReward(rate, prop);
//
// run(markovReward, rate, prop, branchLength, print, 1000);
// }
// public void testLatentStateBranchRateModel() throws FunctionEvaluationException, MaxIterationsExceededException {
//
// LatentStateBranchRateModel model = new LatentStateBranchRateModel(
// new Parameter.Default(0.001), new Parameter.Default(0.5));
//
// TrapezoidIntegrator integator = new TrapezoidIntegrator();
//
// final double branchLength = 2000;
// double integral = integator.integrate(new LatentStateDensityFunction(model, branchLength), 0.0, 1.0);
//
// System.out.println("testLatentStateBeanchRateModel");
// System.out.println("Integral = " + integral);
//
// assertEquals(integral, 1.0, tolerance);
// }
private void run(MarkovReward markovReward, double rate, double prop, double branchLength, boolean print, int length) {
DensityFunction densityFunction = new DensityFunction(markovReward, branchLength, rate, prop);
final double step = branchLength / length;
int i = 0;
double sum = 0.0;
double modeY = 0.0;
double modeX = 0.0;
for (double x = 0.0; x <= branchLength; x += step, ++i) {
double density = 0;
density = densityFunction.value(x);
if (x == 0.0) {
modeY = density;
} else {
if (density > modeY) {
modeY = density;
modeX = x;
}
}
if (x == 0.0 || x == branchLength) {
sum += density;
} else {
sum += 2.0 * density;
}
if (print) {
System.out.println(i + "\t" + String.format("%3.2f", x) + "\t" + String.format("%5.3e", density));
}
}
sum *= (branchLength / 2.0 / length);
// TODO Normalization is missing in LatentBranchRateModel
System.out.println("branchLength = " + branchLength);
System.out.println("rate = " + rate);
System.out.println("prop = " + prop);
System.out.println("Integral = " + sum);
System.out.println("Mode = " + String.format("%3.2e", modeY) + " at " + modeX);
assertEquals(sum, 1.0, tolerance);
TrapezoidIntegrator integrator = new TrapezoidIntegrator();
double integral = 0.0;
try {
integral = integrator.integrate(new UnitDensityFunction(markovReward, branchLength, rate, prop), 0.0, 1.0);
} catch (MaxIterationsExceededException e) {
e.printStackTrace();
} catch (FunctionEvaluationException e) {
e.printStackTrace();
}
System.out.println("unt int = " + integral);
assertEquals(integral, 1.0, tolerance);
System.out.println("\n");
}
use of org.apache.commons.math.MaxIterationsExceededException in project beast-mcmc by beast-dev.
the class GeneralizedIntegerGammaTest method testPdf.
public void testPdf() {
for (int i = 0; i < shape1s.length; ++i) {
final GeneralizedIntegerGammaDistribution gig = new GeneralizedIntegerGammaDistribution(shape1s[i], shape2s[i], rate1s[i], rate2s[i]);
TrapezoidIntegrator integrator = new TrapezoidIntegrator();
double m0 = 0.0;
double m1 = 0.0;
double m2 = 0.0;
try {
m0 = integrator.integrate(new UnivariateRealFunction() {
@Override
public double value(double x) throws FunctionEvaluationException {
final double pdf = gig.pdf(x);
return pdf;
}
}, 0.0, 1000.0);
m1 = integrator.integrate(new UnivariateRealFunction() {
@Override
public double value(double x) throws FunctionEvaluationException {
final double pdf = gig.pdf(x);
return x * pdf;
}
}, 0.0, 1000.0);
m2 = integrator.integrate(new UnivariateRealFunction() {
@Override
public double value(double x) throws FunctionEvaluationException {
final double pdf = gig.pdf(x);
return x * x * pdf;
}
}, 0.0, 1000.0);
} catch (MaxIterationsExceededException e) {
e.printStackTrace();
} catch (FunctionEvaluationException e) {
e.printStackTrace();
}
// Check normalization
assertEquals(1.0, m0, tolerance);
// Check mean
double mean = shape1s[i] / rate1s[i] + shape2s[i] / rate2s[i];
assertEquals(mean, m1, tolerance);
// Check variance
m2 -= m1 * m1;
double variance = shape1s[i] / rate1s[i] / rate1s[i] + shape2s[i] / rate2s[i] / rate2s[i];
// Harder to approximate
assertEquals(variance, m2, tolerance * 10);
}
}
Aggregations