Search in sources :

Example 1 with IllegalDimension

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;
}
Also used : MultivariateDistributionLikelihood(dr.inference.distribution.MultivariateDistributionLikelihood) MultivariateNormalGibbsOperator(dr.inference.operators.MultivariateNormalGibbsOperator) IllegalDimension(dr.math.matrixAlgebra.IllegalDimension)

Example 2 with IllegalDimension

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) {
    }
}
Also used : Matrix(dr.math.matrixAlgebra.Matrix) IllegalDimension(dr.math.matrixAlgebra.IllegalDimension) Vector(dr.math.matrixAlgebra.Vector)

Example 3 with 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();
}
Also used : IllegalDimension(dr.math.matrixAlgebra.IllegalDimension) SymmetricMatrix(dr.math.matrixAlgebra.SymmetricMatrix)

Example 4 with IllegalDimension

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;
}
Also used : Matrix(dr.math.matrixAlgebra.Matrix) IllegalDimension(dr.math.matrixAlgebra.IllegalDimension)

Example 5 with IllegalDimension

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;
}
Also used : Matrix(dr.math.matrixAlgebra.Matrix) IllegalDimension(dr.math.matrixAlgebra.IllegalDimension)

Aggregations

IllegalDimension (dr.math.matrixAlgebra.IllegalDimension)13 Matrix (dr.math.matrixAlgebra.Matrix)7 CholeskyDecomposition (dr.math.matrixAlgebra.CholeskyDecomposition)4 NodeRef (dr.evolution.tree.NodeRef)2 MultivariateNormalDistribution (dr.math.distributions.MultivariateNormalDistribution)2 SymmetricMatrix (dr.math.matrixAlgebra.SymmetricMatrix)2 Vector (dr.math.matrixAlgebra.Vector)2 MultivariateDistributionLikelihood (dr.inference.distribution.MultivariateDistributionLikelihood)1 MultivariateNormalGibbsOperator (dr.inference.operators.MultivariateNormalGibbsOperator)1 WishartSufficientStatistics (dr.math.distributions.WishartSufficientStatistics)1