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

Example 16 with SymmetricMatrix

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

the class FactorTreeGibbsOperator method getMean.

double[] getMean(int column, double[][] precision) {
    Matrix variance = (new SymmetricMatrix(precision)).inverse();
    double[] midMean = new double[lfm.getLoadings().getColumnDimension()];
    double[] condMean = getTreeMean(column);
    double[][] condPrec = getTreePrec(column);
    for (int i = 0; i < midMean.length; i++) {
        //            for (int j = 0; j < midMean.length; j++) {
        midMean[i] += condPrec[i][i] * condMean[i];
    //            }
    }
    for (int i = 0; i < lfm.getLoadings().getRowDimension(); i++) {
        for (int j = 0; j < lfm.getLoadings().getColumnDimension(); j++) {
            if (missingIndicator == null || missingIndicator.getParameterValue(column * lfm.getScaledData().getRowDimension() + i) != 1)
                midMean[j] += lfm.getScaledData().getParameterValue(i, column) * errorPrec.getParameterValue(i, i) * lfm.getLoadings().getParameterValue(i, j) * pathParameter;
        }
    }
    double[] mean = new double[midMean.length];
    for (int i = 0; i < mean.length; i++) {
        for (int j = 0; j < mean.length; j++) {
            mean[i] += variance.component(i, j) * midMean[j];
        }
    }
    return mean;
}

Example 17 with SymmetricMatrix

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

the class FactorIndependenceOperator method doOperation.

@Override
public double doOperation() {
    setupParameters();
    getPrecision(precision);
    double[][] variance = (new SymmetricMatrix(precision)).inverse().toComponents();
    if (randomScan) {
        int i = MathUtils.nextInt(LFM.getFactors().getColumnDimension());
        randomDraw(i, variance);
    }
    for (int i = 0; i < LFM.getFactors().getColumnDimension(); i++) {
        randomDraw(i, variance);
    }
    LFM.getFactors().fireParameterChangedEvent();
    //To change body of implemented methods use File | Settings | File Templates.
    return 0;
}

Example 18 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 19 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()));
    }
}

Example 20 with SymmetricMatrix

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

the class LKJTransformConstrained method inverse.

// "Cholesky" transform from Stan manual
@Override
protected double[] inverse(double[] values) {
    WrappedMatrix.WrappedUpperTriangularMatrix L = fillDiagonal(super.inverse(values), dimVector);
    SymmetricMatrix R = L.transposedProduct();
    if (DEBUG) {
        System.err.println("Z: " + compoundCorrelationSymmetricMatrix(values, dimVector));
        System.err.println("R: " + R);
        try {
            if (!R.isPD()) {
                throw new RuntimeException("The LKJ transform should produce a Positive Definite matrix.");
            }
        } catch (IllegalDimension illegalDimension) {
            illegalDimension.printStackTrace();
        }
    }
    return extractUpperTriangular(R);
}