Examples with SymmetricMatrix dr.math.matrixAlgebra.SymmetricMatrix used on opensource projects

Example 11 with SymmetricMatrix

use of dr.math.matrixAlgebra.SymmetricMatrix in project beast-mcmc by beast-dev.

the class LoadingsGibbsTruncatedOperator method getConditionalDistribution.

private NormalDistribution getConditionalDistribution(double[] meanArray, double[][] variance, int column, int row) {
    double[][] newVariance = new double[meanArray.length - 1][meanArray.length - 1];
    for (int i = 0; i < meanArray.length; i++) {
        for (int j = 0; j < meanArray.length; j++) {
            if (i < column && j < column) {
                newVariance[i][j] = variance[i][j];
            } else if (i < column && j > column) {
                newVariance[i][j - 1] = variance[i][j];
            } else if (i > column && j < column) {
                newVariance[i - 1][j] = variance[i][j];
            } else if (i > column && j > column) {
                newVariance[i - 1][j - 1] = variance[i][j];
            } else {
            }
        }
    }
    double[][] smallPrecision = (new SymmetricMatrix(newVariance)).inverse().toComponents();
    double[] meanStore1 = new double[meanArray.length - 1];
    double[] meanStore2 = new double[meanArray.length - 1];
    double[] precStore = new double[meanArray.length - 1];
    for (int i = 0; i < meanArray.length; i++) {
        if (i < column) {
            meanStore1[i] = LFM.getLoadings().getParameterValue(row, i) - meanArray[i];
        } else if (i > column) {
            meanStore1[i - 1] = LFM.getLoadings().getParameterValue(row, i) - meanArray[i];
        } else {
        }
    }
    for (int i = 0; i < meanArray.length - 1; i++) {
        for (int j = 0; j < meanArray.length - 1; j++) {
            meanStore2[i] += smallPrecision[i][j] * meanStore1[j];
        }
    }
    double mean = meanArray[column];
    for (int i = 0; i < meanArray.length - 1; i++) {
        if (i < column) {
            mean += meanStore2[i] * variance[i][column];
        } else {
            mean += meanStore2[i] * variance[i + 1][column];
        }
    }
    for (int i = 0; i < meanArray.length - 1; i++) {
        for (int j = 0; j < meanArray.length - 1; j++) {
            if (i < column)
                precStore[i] += smallPrecision[i][j] * variance[j][column];
            else
                precStore[i] += smallPrecision[i][j] * variance[j + 1][column];
        }
    }
    double varianceElement = variance[column][column];
    for (int i = 0; i < meanArray.length - 1; i++) {
        if (i < column)
            varianceElement -= precStore[i] * variance[i][column];
        else
            varianceElement -= precStore[i] * variance[i + 1][column];
    }
    return new NormalDistribution(mean, Math.sqrt(varianceElement));
}

Example 12 with SymmetricMatrix

use of dr.math.matrixAlgebra.SymmetricMatrix in project beast-mcmc by beast-dev.

the class LoadingsGibbsTruncatedOperator method drawI.

private void drawI(int i, int column) {
    double[] draws = null;
    precisionArray = new double[loadings.getColumnDimension()][loadings.getColumnDimension()];
    double[][] variance;
    meanMidArray = new double[loadings.getColumnDimension()];
    meanArray = new double[loadings.getColumnDimension()];
    double[][] cholesky = null;
    NormalDistribution conditioned;
    getPrecision(i, precisionArray);
    if (LFM.getLoadings().getParameterValue(i, column) != 0) {
        variance = (new SymmetricMatrix(precisionArray)).inverse().toComponents();
        //        try {
        //            cholesky = new CholeskyDecomposition(variance).getL();
        //        } catch (IllegalDimension illegalDimension) {
        //            illegalDimension.printStackTrace();
        //        }
        getMean(i, variance, meanMidArray, meanArray);
        if (LFM.getFactorDimension() != 1)
            conditioned = getConditionalDistribution(meanArray, variance, column, i);
        else
            conditioned = new NormalDistribution(meanArray[0], Math.sqrt(variance[0][0]));
    } else
        //TODO generify
        conditioned = new NormalDistribution(0, Math.sqrt(1 / priorPrecision));
    if (prior instanceof MomentDistributionModel) {
        if (MathUtils.nextDouble() < .5) {
            getTruncatedDraw(i, column, conditioned);
            getCutoffDraw(i, column, conditioned);
        } else {
            getCutoffDraw(i, column, conditioned);
            getTruncatedDraw(i, column, conditioned);
        }
    } else {
        loadings.setParameterValue(i, column, conditioned.quantile(MathUtils.nextDouble()));
    }
}