Examples with MultivariateNormalDistribution dr.math.distributions.MultivariateNormalDistribution used on opensource projects

Example 6 with MultivariateNormalDistribution

use of dr.math.distributions.MultivariateNormalDistribution in project beast-mcmc by beast-dev.

the class EllipticalSliceOperatorParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    final double weight = xo.getDoubleAttribute(MCMCOperator.WEIGHT);
    final Parameter variable = (Parameter) xo.getChild(Parameter.class);
    boolean drawByRowTemp = false;
    if (xo.hasAttribute(DRAW_BY_ROW)) {
        drawByRowTemp = xo.getBooleanAttribute(DRAW_BY_ROW);
    }
    final boolean drawByRow = drawByRowTemp;
    boolean signal = xo.getAttribute(SIGNAL_CONSTITUENT_PARAMETERS, true);
    if (!signal) {
        Parameter possiblyCompound = variable;
        if (variable instanceof MaskedParameter) {
            possiblyCompound = ((MaskedParameter) variable).getUnmaskedParameter();
        }
        if (!(possiblyCompound instanceof CompoundParameter)) {
            signal = true;
        }
    }
    double bracketAngle = xo.getAttribute(BRACKET_ANGLE, 0.0);
    boolean translationInvariant = xo.getAttribute(TRANSLATION_INVARIANT, false);
    boolean rotationInvariant = xo.getAttribute(ROTATION_INVARIANT, false);
    GaussianProcessRandomGenerator gaussianProcess = (GaussianProcessRandomGenerator) xo.getChild(GaussianProcessRandomGenerator.class);
    if (gaussianProcess == null) {
        final MultivariateDistributionLikelihood likelihood = (MultivariateDistributionLikelihood) xo.getChild(MultivariateDistributionLikelihood.class);
        if (!(likelihood.getDistribution() instanceof GaussianProcessRandomGenerator)) {
            throw new XMLParseException("Elliptical slice sampling only works for multivariate normally distributed random variables");
        }
        if (likelihood.getDistribution() instanceof MultivariateNormalDistribution)
            gaussianProcess = (MultivariateNormalDistribution) likelihood.getDistribution();
        if (likelihood.getDistribution() instanceof MultivariateNormalDistributionModel)
            gaussianProcess = (MultivariateNormalDistributionModel) likelihood.getDistribution();
    }
    EllipticalSliceOperator operator = new EllipticalSliceOperator(variable, gaussianProcess, drawByRow, signal, bracketAngle, translationInvariant, rotationInvariant);
    operator.setWeight(weight);
    return operator;
}

Example 7 with MultivariateNormalDistribution

use of dr.math.distributions.MultivariateNormalDistribution in project beast-mcmc by beast-dev.

the class FactorTreeGibbsOperator method doOperation.

@Override
public double doOperation() {
    if (randomScan) {
        int column = MathUtils.nextInt(factors.getColumnDimension());
        MultivariateNormalDistribution mvn = getMVN(column);
        double[] draw = (double[]) mvn.nextRandom();
        for (int i = 0; i < factors.getRowDimension(); i++) {
            factors.setParameterValue(i, column, draw[i]);
        }
    } else {
        for (int i = 0; i < factors.getColumnDimension(); i++) {
            MultivariateNormalDistribution mvn = getMVN(i);
            double[] draw = (double[]) mvn.nextRandom();
            for (int j = 0; j < factors.getRowDimension(); j++) {
                factors.setParameterValue(j, i, draw[j]);
            }
        }
    }
    return 0;
}

Example 8 with MultivariateNormalDistribution

use of dr.math.distributions.MultivariateNormalDistribution in project beast-mcmc by beast-dev.

the class FactorTreeGibbsOperator method doOperation.

@Override
public double doOperation() {
    if (randomScan) {
        int column = MathUtils.nextInt(factors.getColumnDimension());
        MultivariateNormalDistribution mvn = getMVN(column);
        double[] draw = (double[]) mvn.nextRandom();
        for (int i = 0; i < factors.getRowDimension(); i++) {
            factors.setParameterValue(i, column, draw[i]);
        }
    } else {
        for (int i = 0; i < factors.getColumnDimension(); i++) {
            MultivariateNormalDistribution mvn = getMVN(i);
            double[] draw = (double[]) mvn.nextRandom();
            for (int j = 0; j < factors.getRowDimension(); j++) {
                factors.setParameterValue(j, i, draw[j]);
            }
        }
    }
    return 0;
}

Example 9 with MultivariateNormalDistribution

use of dr.math.distributions.MultivariateNormalDistribution in project beast-mcmc by beast-dev.

the class MultivariateOUModel method calculateLogLikelihood.

public double calculateLogLikelihood() {
    double logLikelihood = 0;
    double[] previous = new double[K];
    double[] current = new double[K];
    double[] tmpHolder;
    double[][] G = gamma.getParameterAsMatrix();
    double[] theta = dependentParam.getParameterValues();
    double[] Xbeta = null;
    boolean hasEffects = getNumberOfFixedEffects() > 0;
    if (!conditionalPrecisionKnown)
        calculateConditionPrecision();
    try {
        if (new Matrix(G).determinant() < 0.01)
            return Double.NEGATIVE_INFINITY;
    } catch (IllegalDimension illegalDimension) {
        illegalDimension.printStackTrace();
    }
    int index = 0;
    if (!hasEffects) {
        for (int i = 0; i < K; i++) previous[i] = theta[index++];
    } else {
        Xbeta = getXBeta();
        for (int i = 0; i < K; i++) {
            previous[i] = theta[index] - Xbeta[index];
            index++;
        }
    }
    initialPrior = new MultivariateNormalDistribution(initialPriorMean, new Matrix(G).inverse().toComponents());
    logLikelihood += initialPrior.logPdf(previous);
    double save = logLikelihood;
    double save2 = 0;
    int oldMapTime = -1;
    double[][] conditionalPrecision = new double[K][K];
    for (int timeStep = 0; timeStep < numTimeSteps; timeStep++) {
        int thisMapTime = mapTime[timeStep];
        if (thisMapTime != oldMapTime) {
            System.arraycopy(Wt, Ksquared * thisMapTime, W, 0, Ksquared);
            for (int i = 0; i < K; i++) System.arraycopy(conditionPrecisionVector, Ksquared * thisMapTime + K * i, conditionalPrecision[i], 0, K);
            oldMapTime = thisMapTime;
        }
        double[] mean = new double[K];
        int u = 0;
        for (int i = 0; i < K; i++) {
            for (int j = 0; j < K; j++) mean[i] += W[u++] * previous[j];
        }
        if (!hasEffects) {
            for (int i = 0; i < K; i++) current[i] = theta[index++];
        } else {
            for (int i = 0; i < K; i++) {
                current[i] = theta[index] - Xbeta[index];
                index++;
            }
        }
        MultivariateNormalDistribution density = new MultivariateNormalDistribution(mean, conditionalPrecision);
        double partialLogLikelihood = density.logPdf(current);
        if (partialLogLikelihood > 10) {
            return Double.NEGATIVE_INFINITY;
        }
        logLikelihood += partialLogLikelihood;
        // move to next point
        tmpHolder = previous;
        previous = current;
        current = tmpHolder;
    }
    if (logLikelihood > 100) {
        System.err.println("got here end");
        System.err.println("save1 = " + save);
        System.err.println("save2 = " + save2);
        System.exit(-1);
    }
    likelihoodKnown = true;
    return logLikelihood;
}

Example 10 with MultivariateNormalDistribution

use of dr.math.distributions.MultivariateNormalDistribution in project beast-mcmc by beast-dev.

the class FullyConjugateMultivariateTraitLikelihood method getReport.

@Override
public String getReport() {
    StringBuilder sb = new StringBuilder();
    // sb.append(this.g)
    // System.err.println("Hello");
    sb.append("Tree:\n");
    sb.append(getId()).append("\t");
    sb.append(treeModel.toString());
    sb.append("\n\n");
    double[][] treeVariance = computeTreeVariance();
    double[][] traitPrecision = getDiffusionModel().getPrecisionmatrix();
    Matrix traitVariance = new Matrix(traitPrecision).inverse();
    double[][] jointVariance = KroneckerOperation.product(treeVariance, traitVariance.toComponents());
    sb.append("Tree variance:\n");
    sb.append(new Matrix(treeVariance));
    sb.append(matrixMin(treeVariance)).append("\t").append(matrixMax(treeVariance)).append("\t").append(matrixSum(treeVariance));
    sb.append("\n\n");
    sb.append("Trait variance:\n");
    sb.append(traitVariance);
    sb.append("\n\n");
    // sb.append("Joint variance:\n");
    // sb.append(new Matrix(jointVariance));
    // sb.append("\n\n");
    sb.append("Tree dim: ").append(treeVariance.length).append("\n");
    sb.append("data dim: ").append(jointVariance.length);
    sb.append("\n\n");
    double[] data = new double[jointVariance.length];
    System.arraycopy(meanCache, 0, data, 0, jointVariance.length);
    if (nodeToClampMap != null) {
        int offset = treeModel.getExternalNodeCount() * getDimTrait();
        for (Map.Entry<NodeRef, RestrictedPartials> clamps : nodeToClampMap.entrySet()) {
            double[] partials = clamps.getValue().getPartials();
            for (double partial : partials) {
                data[offset] = partial;
                ++offset;
            }
        }
    }
    sb.append("Data:\n");
    sb.append(new Vector(data)).append("\n");
    sb.append(data.length).append("\t").append(vectorMin(data)).append("\t").append(vectorMax(data)).append("\t").append(vectorSum(data));
    sb.append(treeModel.getNodeTaxon(treeModel.getExternalNode(0)).getId());
    sb.append("\n\n");
    MultivariateNormalDistribution mvn = new MultivariateNormalDistribution(new double[data.length], new Matrix(jointVariance).inverse().toComponents());
    double logDensity = mvn.logPdf(data);
    sb.append("logLikelihood: ").append(getLogLikelihood()).append(" == ").append(logDensity).append("\n\n");
    final WishartSufficientStatistics sufficientStatistics = getWishartStatistics();
    final double[] outerProducts = sufficientStatistics.getScaleMatrix();
    sb.append("Outer-products (DP):\n");
    sb.append(new Vector(outerProducts));
    sb.append(sufficientStatistics.getDf()).append("\n");
    Matrix treePrecision = new Matrix(treeVariance).inverse();
    final int n = data.length / traitPrecision.length;
    final int p = traitPrecision.length;
    double[][] tmp = new double[n][p];
    for (int i = 0; i < n; ++i) {
        System.arraycopy(data, i * p, tmp[i], 0, p);
    }
    Matrix y = new Matrix(tmp);
    Matrix S = null;
    try {
        // Using Matrix-Normal form
        S = y.transpose().product(treePrecision).product(y);
    } catch (IllegalDimension illegalDimension) {
        illegalDimension.printStackTrace();
    }
    sb.append("Outer-products (from tree variance:\n");
    sb.append(S);
    sb.append("\n\n");
    return sb.toString();
}