use of dr.math.matrixAlgebra.IllegalDimension in project beast-mcmc by beast-dev.
the class MultivariateNormalGibbsOperatorParser method parseXMLObject.
@Override
public Object parseXMLObject(XMLObject xo) throws XMLParseException {
MultivariateDistributionLikelihood prior = (MultivariateDistributionLikelihood) xo.getChild(PRIOR).getChild(MultivariateDistributionLikelihood.class);
MultivariateDistributionLikelihood likelihood = (MultivariateDistributionLikelihood) xo.getChild(LIKELIHOOD).getChild(MultivariateDistributionLikelihood.class);
// CompoundParameter data = (CompoundParameter) xo.getChild(CompoundParameter.class);
String weightTemp = (String) xo.getAttribute(WEIGHT);
Double weight = Double.parseDouble(weightTemp);
try {
//To change body of implemented methods use File | Settings | File Templates.
return new MultivariateNormalGibbsOperator(likelihood, prior, weight);
} catch (IllegalDimension illegalDimension) {
//To change body of catch statement use File | Settings | File Templates.
illegalDimension.printStackTrace();
}
return null;
}
use of dr.math.matrixAlgebra.IllegalDimension in project beast-mcmc by beast-dev.
the class BinaryCovarionModelTest method testEquilibriumDistribution.
/**
* Tests that pi*Q = 0
*/
public void testEquilibriumDistribution() {
alpha.setParameterValue(0, 0.1);
switchingRate.setParameterValue(0, 1.0);
model.setupMatrix();
double[] pi = model.getFrequencyModel().getFrequencies();
try {
Matrix m = new Matrix(model.getQ());
Vector p = new Vector(pi);
Vector y = m.product(p);
assertEquals(0.0, y.norm(), 1e-14);
} catch (IllegalDimension illegalDimension) {
}
}
use of dr.math.matrixAlgebra.IllegalDimension in project beast-mcmc by beast-dev.
the class PrecisionMatrixGibbsOperator method getOperationScaleMatrixAndSetObservationCount.
public double[][] getOperationScaleMatrixAndSetObservationCount() {
// calculate sum-of-the-weighted-squares matrix over tree
double[][] S = new double[dim][dim];
SymmetricMatrix S2;
SymmetricMatrix inverseS2 = null;
// Need to reset, as incrementOuterProduct can be recursive
numberObservations = 0;
if (isSampledTraitLikelihood) {
incrementOuterProduct(S, treeModel.getRoot());
} else {
// IntegratedTraitLikelihood
if (traitModel != null) {
// is a tree
// TODO deprecate usage
incrementOuterProduct(S, (ConjugateWishartStatisticsProvider) traitModel);
} else if (conjugateWishartProvider != null) {
incrementOuterProduct(S, conjugateWishartProvider);
} else {
// is a normal-normal-wishart model
incrementOuterProduct(S, multivariateLikelihood);
}
}
try {
S2 = new SymmetricMatrix(S);
if (pathWeight != 1.0) {
S2 = (SymmetricMatrix) S2.product(pathWeight);
}
if (priorInverseScaleMatrix != null)
S2 = priorInverseScaleMatrix.add(S2);
inverseS2 = (SymmetricMatrix) S2.inverse();
// if (S[0][0] < 0.0) {
// System.err.println("ERROR A");
// System.err.println(new Matrix(S));
// }
//
// if (S2.component(0, 0) < 0.0) {
// System.err.println("ERROR B");
// System.err.println(S2);
// }
//
// if (inverseS2.component(0, 0) < 0.0) {
// System.err.println("ERROR C");
// System.err.println("S:\n" + new Matrix(S));
// System.err.println("S2:\n" + S2);
// System.err.println(inverseS2);
// }
} catch (IllegalDimension illegalDimension) {
illegalDimension.printStackTrace();
}
assert inverseS2 != null;
return inverseS2.toComponents();
}
use of dr.math.matrixAlgebra.IllegalDimension in project beast-mcmc by beast-dev.
the class WishartDistribution method logPdf.
public static double logPdf(Matrix W, Matrix Sinv, double df, int dim, double logNormalizationConstant) {
double logDensity = 0;
try {
// if (!W.isPD()) { // TODO isPD() does not appear to work
// return Double.NEGATIVE_INFINITY;
// }
// Returns NaN is W is not positive-definite.
logDensity = W.logDeterminant();
if (Double.isInfinite(logDensity) || Double.isNaN(logDensity)) {
return Double.NEGATIVE_INFINITY;
}
logDensity *= 0.5;
logDensity *= df - dim - 1;
// the whole matrix
if (Sinv != null) {
Matrix product = Sinv.product(W);
for (int i = 0; i < dim; i++) logDensity -= 0.5 * product.component(i, i);
}
} catch (IllegalDimension illegalDimension) {
illegalDimension.printStackTrace();
}
logDensity += logNormalizationConstant;
return logDensity;
}
use of dr.math.matrixAlgebra.IllegalDimension in project beast-mcmc by beast-dev.
the class InverseWishartDistribution method logPdf.
public double logPdf(double[] x) {
Matrix W = new Matrix(x, dim, dim);
double logDensity = 0;
try {
logDensity = Math.log(W.determinant());
logDensity *= -0.5;
logDensity *= df + dim + 1;
Matrix product = S.product(W.inverse());
for (int i = 0; i < dim; i++) logDensity -= 0.5 * product.component(i, i);
} catch (IllegalDimension illegalDimension) {
illegalDimension.printStackTrace();
}
logDensity += logNormalizationConstant;
return logDensity;
}
Aggregations