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

Example 16 with Matrix

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

the class CachedMatrixInverse method computeInverse.

private void computeInverse() {
    if (DEBUG) {
        System.err.println("CachedMatrixInverse.computeInverse()");
    }
    if (EMJL) {
        // TODO Avoid multiple copies
        DenseMatrix64F source = new DenseMatrix64F(base.getParameterAsMatrix());
        DenseMatrix64F destination = new DenseMatrix64F(getColumnDimension(), getColumnDimension());
        CommonOps.invert(source, destination);
        inverse = new WrappedMatrix.WrappedDenseMatrix(destination);
    } else {
        inverse = new WrappedMatrix.ArrayOfArray(new Matrix(base.getParameterAsMatrix()).inverse().toComponents());
    }
}

Example 17 with Matrix

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

the class MVOUCovarianceOperator method doOperation.

public double doOperation() {
    double[][] draw = WishartDistribution.nextWishart(priorDf, I);
    // double[][] good = varMatrix.getParameterAsMatrix();
    // double[][] saveOld = varMatrix.getParameterAsMatrix();
    // System.err.println("draw:\n"+new Matrix(draw));
    double[][] oldValue = varMatrix.getParameterAsMatrix();
    for (int i = 0; i < dim; i++) {
        Parameter column = varMatrix.getParameter(i);
        for (int j = 0; j < dim; j++) column.setParameterValue(j, mixingFactor * oldValue[j][i] + (1.0 - mixingFactor) * draw[j][i]);
    }
    // varMatrix.fireParameterChangedEvent();
    // calculate Hastings ratio
    // System.err.println("oldValue:\n"+new Matrix(oldValue).toString());
    // System.err.println("newValue:\n"+new Matrix(varMatrix.getParameterAsMatrix()).toString());
    Matrix forwardDrawMatrix = new Matrix(draw);
    for (int i = 0; i < dim; i++) {
        for (int j = 0; j < dim; j++) {
            // saveOld[i][j] *= - mixingFactor;
            // saveOld[i][j] += varMatrix.getParameterValue(i,j);
            // saveOld[i][j] /= 1.0 - mixingFactor;
            oldValue[i][j] -= mixingFactor * varMatrix.getParameterValue(i, j);
            oldValue[i][j] /= 1.0 - mixingFactor;
        }
    }
    // double[][] saveNew = varMatrix.getParameterAsMatrix();
    Matrix backwardDrawMatrix = new Matrix(oldValue);
    // System.err.println("forward:\n"+forwardDrawMatrix);
    // System.err.println("backward:\n"+backwardDrawMatrix);
    // System.err.println("calc start");
    // if( Math.abs(backwardDrawMatrix.component(0,0) + 0.251) < 0.001 ) {
    // System.err.println("found:\n"+backwardDrawMatrix);
    // 
    // System.err.println("original:\n"+new Matrix(good));
    // System.err.println("draw:\n"+new Matrix(draw));
    // System.err.println("proposed:\n"+new Matrix(varMatrix.getParameterAsMatrix()));
    // System.err.println("mixing = "+mixingFactor);
    // System.err.println("back[0][0] = "+backwardDrawMatrix.component(0,0));
    // System.err.println("saveOld[0][0] = "+saveOld[0][0]);
    // 
    // 
    // }
    double bProb = WishartDistribution.logPdf(backwardDrawMatrix, Iinv, priorDf, dim, // WishartDistribution.computeNormalizationConstant(Iinv,priorDf,dim));
    0);
    if (bProb == Double.NEGATIVE_INFINITY) {
        // not clear if this means a HR of -Inf or a RuntimeException
        return Double.NEGATIVE_INFINITY;
    }
    double fProb = WishartDistribution.logPdf(forwardDrawMatrix, Iinv, priorDf, dim, // WishartDistribution.computeNormalizationConstant(Iinv,priorDf,dim));
    0);
    return bProb - fProb;
}

Example 18 with Matrix

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

the class MultivariateNormalGibbsOperator method getMean.

private Vector getMean() throws IllegalDimension {
    Vector meanSum = new Vector(getMeanSum());
    Matrix workingPrecision = new Matrix(likelihoodPrecision.getParameterAsMatrix());
    Vector meanPart = workingPrecision.product(meanSum);
    meanPart = meanPart.add(priorPrecision.product(priorMean));
    Matrix varPart = getPrecision().inverse();
    Vector answer = varPart.product(meanPart);
    return answer;
}

Example 19 with Matrix

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

the class TraitValidationProvider method getReport.

@Override
public String getReport() {
    int iterations = 1000000;
    int nMissing = missingInds.length;
    double[] sums = new double[nMissing];
    double[][] sumSquares = new double[nMissing][nMissing];
    for (int iteration = 0; iteration < iterations; iteration++) {
        double[] extendedVals = extensionDelegate.getExtendedValues();
        for (int i = 0; i < nMissing; i++) {
            double vi = extendedVals[missingInds[i]];
            sums[i] += vi;
            sumSquares[i][i] += vi * vi;
            for (int j = 0; j < i; j++) {
                double vj = extendedVals[missingInds[j]];
                sumSquares[i][j] += vi * vj;
            }
        }
    }
    for (int i = 0; i < nMissing; i++) {
        sums[i] = sums[i] / iterations;
        sumSquares[i][i] = sumSquares[i][i] / iterations - sums[i] * sums[i];
        for (int j = 0; j < i; j++) {
            sumSquares[i][j] = sumSquares[i][j] / iterations - sums[i] * sums[j];
            sumSquares[j][i] = sumSquares[i][j];
        }
    }
    StringBuilder sb = new StringBuilder();
    sb.append("Mean:\t");
    sb.append(new Vector(sums));
    sb.append("\n");
    sb.append("Covariance:\t");
    sb.append(new Matrix(sumSquares));
    return sb.toString();
}

Example 20 with Matrix

use of dr.math.matrixAlgebra.Matrix 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;
}