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

Example 1 with MultivariateNormalDistribution

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

the class FullyConjugateMultivariateTraitLikelihood method checkLogLikelihood.

protected void checkLogLikelihood(double logLikelihood, double logRemainders, double[] conditionalRootMean, double conditionalRootPrecision, double[][] traitPrecision) {
    // System.err.println("root cMean    : " + new Vector(conditionalRootMean));
    // System.err.println("root cPrec    : " + conditionalRootPrecision);
    // System.err.println("diffusion prec: " + new Matrix(traitPrecision));
    // 
    // System.err.println("prior mean    : " + new Vector(rootPriorMean));
    // System.err.println("prior prec    : " + rootPriorSampleSize);
    double upperPrecision = conditionalRootPrecision * rootPriorSampleSize / (conditionalRootPrecision + rootPriorSampleSize);
    // System.err.println("root cPrec    : " + upperPrecision);
    double[][] newPrec = new double[traitPrecision.length][traitPrecision.length];
    for (int i = 0; i < traitPrecision.length; ++i) {
        for (int j = 0; j < traitPrecision.length; ++j) {
            newPrec[i][j] = traitPrecision[i][j] * upperPrecision;
        }
    }
    MultivariateNormalDistribution mvn = new MultivariateNormalDistribution(rootPriorMean, newPrec);
    double logPdf = mvn.logPdf(conditionalRootMean);
    if (Math.abs(logLikelihood - logRemainders - logPdf) > 1E-3) {
        System.err.println("Got here subclass: " + logLikelihood);
        System.err.println("logValue         : " + (logRemainders + logPdf));
        System.err.println("logRemainder = " + logRemainders);
        System.err.println("");
    }
// System.err.println("logRemainders    : " + logRemainders);
// System.err.println("logPDF           : " + logPdf);
// System.exit(-1);
}

Example 2 with MultivariateNormalDistribution

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

the class ContinuousDataLikelihoodDelegate method getReport.

@Override
public String getReport() {
    StringBuilder sb = new StringBuilder();
    sb.append("Tree:\n");
    sb.append(callbackLikelihood.getId()).append("\t");
    sb.append(cdi.getReport());
    final Tree tree = callbackLikelihood.getTree();
    sb.append(tree.toString());
    sb.append("\n\n");
    final double normalization = rateTransformation.getNormalization();
    final double priorSampleSize = rootProcessDelegate.getPseudoObservations();
    double[][] treeStructure = MultivariateTraitDebugUtilities.getTreeVariance(tree, callbackLikelihood.getBranchRateModel(), 1.0, Double.POSITIVE_INFINITY);
    sb.append("Tree structure:\n");
    sb.append(new Matrix(treeStructure));
    sb.append("\n\n");
    double[][] treeSharedLengths = MultivariateTraitDebugUtilities.getTreeVariance(tree, callbackLikelihood.getBranchRateModel(), rateTransformation.getNormalization(), Double.POSITIVE_INFINITY);
    double[][] treeVariance = getTreeVariance();
    double[][] traitPrecision = getDiffusionModel().getPrecisionmatrix();
    Matrix traitVariance = new Matrix(traitPrecision).inverse();
    double[][] jointVariance = diffusionProcessDelegate.getJointVariance(priorSampleSize, treeVariance, treeSharedLengths, traitVariance.toComponents());
    if (dataModel instanceof RepeatedMeasuresTraitDataModel) {
        for (int tip = 0; tip < tipCount; ++tip) {
            double[] partial = dataModel.getTipPartial(tip, false);
            WrappedMatrix tipVariance = new WrappedMatrix.Raw(partial, dimTrait + dimTrait * dimTrait, dimTrait, dimTrait);
            for (int row = 0; row < dimTrait; ++row) {
                for (int col = 0; col < dimTrait; ++col) {
                    jointVariance[tip * dimTrait + row][tip * dimTrait + col] += tipVariance.get(row, col);
                }
            }
        }
    }
    Matrix treeV = new Matrix(treeVariance);
    Matrix treeP = treeV.inverse();
    sb.append("Tree variance:\n");
    sb.append(treeV);
    sb.append("Tree precision:\n");
    sb.append(treeP);
    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");
    double[] priorMean = rootPrior.getMean();
    sb.append("prior mean: ").append(new dr.math.matrixAlgebra.Vector(priorMean));
    sb.append("\n\n");
    sb.append("Joint variance:\n");
    sb.append(new Matrix(jointVariance));
    sb.append("\n\n");
    sb.append("Joint precision:\n");
    sb.append(new Matrix(getTreeTraitPrecision()));
    sb.append("\n\n");
    double[][] treeDrift = MultivariateTraitDebugUtilities.getTreeDrift(tree, priorMean, cdi, diffusionProcessDelegate);
    if (diffusionProcessDelegate.hasDrift()) {
        sb.append("Tree drift (including root mean):\n");
        sb.append(new Matrix(treeDrift));
        sb.append("\n\n");
    }
    double[] drift = KroneckerOperation.vectorize(treeDrift);
    final int datumLength = tipCount * dimProcess;
    sb.append("Tree dim : ").append(treeStructure.length).append("\n");
    sb.append("dimTrait : ").append(dimTrait).append("\n");
    sb.append("numTraits: ").append(numTraits).append("\n");
    sb.append("jVar dim : ").append(jointVariance.length).append("\n");
    sb.append("datum dim: ").append(datumLength);
    sb.append("\n\n");
    double[] data = dataModel.getParameter().getParameterValues();
    if (dataModel instanceof ContinuousTraitDataModel) {
        for (int tip = 0; tip < tipCount; ++tip) {
            double[] tipData = ((ContinuousTraitDataModel) dataModel).getTipObservation(tip, precisionType);
            System.arraycopy(tipData, 0, data, tip * numTraits * dimProcess, numTraits * dimProcess);
        }
    }
    sb.append("data: ").append(new dr.math.matrixAlgebra.Vector(data));
    sb.append("\n\n");
    double[][] graphStructure = MultivariateTraitDebugUtilities.getGraphVariance(tree, callbackLikelihood.getBranchRateModel(), normalization, priorSampleSize);
    double[][] jointGraphVariance = KroneckerOperation.product(graphStructure, traitVariance.toComponents());
    sb.append("graph structure:\n");
    sb.append(new Matrix(graphStructure));
    sb.append("\n\n");
    for (int trait = 0; trait < numTraits; ++trait) {
        sb.append("Trait #").append(trait).append("\n");
        double[] rawDatum = new double[datumLength];
        List<Integer> missing = new ArrayList<Integer>();
        int index = 0;
        for (int tip = 0; tip < tipCount; ++tip) {
            for (int dim = 0; dim < dimProcess; ++dim) {
                double d = data[tip * dimProcess * numTraits + trait * dimProcess + dim];
                rawDatum[index] = d;
                if (Double.isNaN(d)) {
                    missing.add(index);
                }
                ++index;
            }
        }
        double[][] varianceDatum = jointVariance;
        double[] datum = rawDatum;
        double[] driftDatum = drift;
        int[] notMissingIndices;
        notMissingIndices = new int[datumLength - missing.size()];
        int offsetNotMissing = 0;
        for (int i = 0; i < datumLength; ++i) {
            if (!missing.contains(i)) {
                notMissingIndices[offsetNotMissing] = i;
                ++offsetNotMissing;
            }
        }
        datum = Matrix.gatherEntries(rawDatum, notMissingIndices);
        varianceDatum = Matrix.gatherRowsAndColumns(jointVariance, notMissingIndices, notMissingIndices);
        driftDatum = Matrix.gatherEntries(drift, notMissingIndices);
        sb.append("datum : ").append(new dr.math.matrixAlgebra.Vector(datum)).append("\n");
        sb.append("drift : ").append(new dr.math.matrixAlgebra.Vector(driftDatum)).append("\n");
        sb.append("variance:\n");
        sb.append(new Matrix(varianceDatum));
        MultivariateNormalDistribution mvn = new MultivariateNormalDistribution(driftDatum, new Matrix(varianceDatum).inverse().toComponents());
        double logDensity = mvn.logPdf(datum);
        sb.append("\n\n");
        sb.append("logDatumLikelihood: ").append(logDensity).append("\n\n");
        // Compute joint for internal nodes
        int[] cNotMissingJoint = new int[dimProcess * tipCount];
        int[] cMissingJoint = new int[dimProcess * (tipCount - 1)];
        // External nodes
        for (int tipTrait = 0; tipTrait < dimProcess * tipCount; ++tipTrait) {
            cNotMissingJoint[tipTrait] = tipTrait;
        }
        // Internal nodes
        for (int tipTrait = dimProcess * tipCount; tipTrait < dimProcess * (2 * tipCount - 1); ++tipTrait) {
            cMissingJoint[tipTrait - dimProcess * tipCount] = tipTrait;
        }
        double[] rawDatumJoint = new double[dimProcess * (2 * tipCount - 1)];
        System.arraycopy(rawDatum, 0, rawDatumJoint, 0, rawDatum.length);
        double[][] driftJointMatrix = MultivariateTraitDebugUtilities.getGraphDrift(tree, cdi, diffusionProcessDelegate);
        double[] driftJoint = KroneckerOperation.vectorize(driftJointMatrix);
        for (int idx = 0; idx < driftJoint.length / dimProcess; ++idx) {
            for (int dim = 0; dim < dimProcess; ++dim) {
                driftJoint[idx * dimProcess + dim] += priorMean[dim];
            }
        }
        ConditionalVarianceAndTransform cVarianceJoint = new ConditionalVarianceAndTransform(new Matrix(jointGraphVariance), cMissingJoint, cNotMissingJoint);
        double[] cMeanJoint = cVarianceJoint.getConditionalMean(rawDatumJoint, 0, driftJoint, 0);
        sb.append("cDriftJoint: ").append(new dr.math.matrixAlgebra.Vector(driftJoint)).append("\n\n");
        sb.append("cMeanInternalJoint: ").append(new dr.math.matrixAlgebra.Vector(cMeanJoint)).append("\n\n");
        // Compute full conditional distributions
        sb.append("Full conditional distributions:\n");
        int offsetNotMissing2 = 0;
        for (int tip = 0; tip < tipCount; ++tip) {
            int offset = tip * dimProcess;
            int dimTip = 0;
            for (int cTrait = 0; cTrait < dimProcess; cTrait++) {
                if ((offsetNotMissing2 + cTrait < notMissingIndices.length) && notMissingIndices[offsetNotMissing2 + cTrait] < offset + dimProcess) {
                    dimTip++;
                }
            }
            int[] cMissing = new int[dimProcess];
            int[] cNotMissing = new int[notMissingIndices.length - dimTip];
            for (int cTrait = 0; cTrait < dimProcess; ++cTrait) {
                cMissing[cTrait] = offset + cTrait;
            }
            for (int m = 0; m < offsetNotMissing2; ++m) {
                cNotMissing[m] = notMissingIndices[m];
            }
            offsetNotMissing2 += dimTip;
            for (int m = offsetNotMissing2; m < notMissingIndices.length; ++m) {
                cNotMissing[m - dimTip] = notMissingIndices[m];
            }
            ConditionalVarianceAndTransform cVariance = new ConditionalVarianceAndTransform(new Matrix(jointVariance), cMissing, cNotMissing);
            double[] cMean = cVariance.getConditionalMean(rawDatum, 0, drift, 0);
            Matrix cVar = cVariance.getConditionalVariance();
            sb.append("cMean #").append(tip).append(" ").append(new dr.math.matrixAlgebra.Vector(cMean)).append(" cVar [").append(cVar).append("]\n");
        }
    }
    return sb.toString();
}

Example 3 with MultivariateNormalDistribution

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

the class IntegratedFactorAnalysisLikelihood method getReport.

@Override
public String getReport() {
    StringBuilder sb = new StringBuilder();
    double logComponents = 0;
    if (delegate != null) {
        double logInc = delegate.getCallbackLikelihood().getLogLikelihood();
        final Tree tree = delegate.getCallbackLikelihood().getTree();
        final BranchRates branchRates = delegate.getCallbackLikelihood().getBranchRateModel();
        sb.append(tree.toString());
        sb.append("\n\n");
        final double normalization = delegate.getRateTransformation().getNormalization();
        final double priorSampleSize = delegate.getRootProcessDelegate().getPseudoObservations();
        double[][] treeStructure = MultivariateTraitDebugUtilities.getTreeVariance(tree, branchRates, 1.0, Double.POSITIVE_INFINITY);
        sb.append("Tree structure:\n");
        sb.append(new Matrix(treeStructure));
        sb.append("\n\n");
        double[][] treeSharedLengths = MultivariateTraitDebugUtilities.getTreeVariance(tree, branchRates, normalization, Double.POSITIVE_INFINITY);
        double[][] treeVariance = MultivariateTraitDebugUtilities.getTreeVariance(tree, branchRates, normalization, priorSampleSize);
        Matrix treeV = new Matrix(treeVariance);
        Matrix treeP = treeV.inverse();
        sb.append("Tree variance:\n");
        sb.append(treeV);
        sb.append("Tree precision:\n");
        sb.append(treeP);
        sb.append("\n\n");
        Matrix Lt = new Matrix(loadingsTransposed.getParameterAsMatrix());
        sb.append("Loadings:\n");
        sb.append(Lt.transpose());
        sb.append("\n\n");
        double[][] diffusionPrecision = delegate.getDiffusionModel().getPrecisionmatrix();
        Matrix diffusionVariance = new Matrix(diffusionPrecision).inverse();
        Matrix loadingsVariance = null;
        try {
            loadingsVariance = Lt.product(diffusionVariance.product(Lt.transpose()));
        } catch (IllegalDimension illegalDimension) {
            illegalDimension.printStackTrace();
        }
        sb.append("Loadings variance:\n");
        sb.append(loadingsVariance);
        sb.append("\n\n");
        assert (loadingsVariance != null);
        Matrix loadingsFactorsVariance = MultivariateTraitDebugUtilities.getJointVarianceFactor(priorSampleSize, treeVariance, treeSharedLengths, loadingsVariance.toComponents(), diffusionVariance.toComponents(), delegate.getDiffusionProcessDelegate(), Lt);
        Matrix gamma = buildDiagonalMatrix(traitPrecision.getParameterValues());
        sb.append("Trait precision:\n");
        sb.append(gamma);
        sb.append("\n\n");
        Matrix gammaVariance = gamma.inverse();
        double[] tmp = new double[tree.getExternalNodeCount()];
        Arrays.fill(tmp, 1.0);
        Matrix identity = buildDiagonalMatrix(tmp);
        Matrix errorVariance = new Matrix(KroneckerOperation.product(identity.toComponents(), gammaVariance.toComponents()));
        sb.append("Loadings-factors variance:\n");
        sb.append(loadingsFactorsVariance);
        sb.append("\n\n");
        sb.append("Error variance\n");
        sb.append(errorVariance);
        sb.append("\n\n");
        Matrix totalVariance = null;
        try {
            totalVariance = loadingsFactorsVariance.add(errorVariance);
        } catch (IllegalDimension illegalDimension) {
            illegalDimension.printStackTrace();
        }
        double[] allData = getParameter().getParameterValues();
        List<Integer> notMissing = new ArrayList<>();
        for (int taxon = 0; taxon < numTaxa; ++taxon) {
            double[] observed = observedIndicators[taxon];
            for (int trait = 0; trait < dimTrait; ++trait) {
                if (observed[trait] == 0.0) {
                    System.err.println("Missing taxon " + taxon + " trait " + trait);
                } else {
                    notMissing.add(taxon * dimTrait + trait);
                }
            }
        }
        double[] priorMean = delegate.getRootPrior().getMean();
        Matrix treeDrift = new Matrix(MultivariateTraitDebugUtilities.getTreeDrift(tree, priorMean, delegate.getIntegrator(), delegate.getDiffusionProcessDelegate()));
        if (delegate.getDiffusionProcessDelegate().hasDrift()) {
            sb.append("Tree drift (including root mean):\n");
            sb.append(new Matrix(treeDrift.toComponents()));
            sb.append("\n\n");
        }
        try {
            loadingsFactorsVariance = treeDrift.product(Lt.transpose());
        } catch (IllegalDimension illegalDimension) {
            illegalDimension.printStackTrace();
        }
        double[] drift = KroneckerOperation.vectorize(loadingsFactorsVariance.toComponents());
        int[] notMissingIndices = new int[notMissing.size()];
        double[] data = new double[notMissing.size()];
        for (int i = 0; i < notMissing.size(); ++i) {
            notMissingIndices[i] = notMissing.get(i);
            data[i] = allData[notMissing.get(i)];
        }
        if (totalVariance != null) {
            totalVariance = new Matrix(Matrix.gatherRowsAndColumns(totalVariance.toComponents(), notMissingIndices, notMissingIndices));
        }
        Matrix totalPrecision = null;
        if (totalVariance != null) {
            totalPrecision = totalVariance.inverse();
        }
        drift = Matrix.gatherEntries(drift, notMissingIndices);
        sb.append("Total variance:\n");
        sb.append(totalVariance);
        sb.append("\n\n");
        sb.append("Total precision:\n");
        sb.append(totalPrecision);
        sb.append("\n\n");
        sb.append("Data:\n");
        sb.append(new Vector(data));
        sb.append("\n\n");
        sb.append("Expectations:\n");
        sb.append(new Vector(drift));
        sb.append("\n\n");
        MultivariateNormalDistribution mvn = null;
        if (totalPrecision != null) {
            mvn = new MultivariateNormalDistribution(drift, totalPrecision.toComponents());
        }
        double logDensity = 0;
        if (mvn != null) {
            logDensity = mvn.logPdf(data);
        }
        sb.append("logMultiVariateNormalDensity = ").append(logDensity).append("\n\n");
        sb.append("traitDataLikelihood = ").append(logInc).append("\n");
        logComponents += logInc;
    }
    sb.append("logLikelihood = ").append(getLogLikelihood()).append("\n");
    if (logComponents != 0.0) {
        sb.append("total likelihood = ").append((getLogLikelihood() + logComponents)).append("\n");
    }
    return sb.toString();
}

Example 4 with MultivariateNormalDistribution

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

the class NumericalSpaceTimeProbs2D method main.

public static void main(String[] args) throws FileNotFoundException {
    Rectangle2D bounds = new Rectangle2D.Double(0, 0, 1, 1);
    MultivariateNormalDistribution D = new MultivariateNormalDistribution(new double[] { 0.0 }, new double[][] { { 1, 0 }, { 0, 1 } });
    NumericalSpaceTimeProbs2D nstp2D = new NumericalSpaceTimeProbs2D(50, 50, 50, 1, 0.02, bounds, D, SpaceTimeRejector.Utils.createSimpleBounds2D(bounds));
    long startTime = System.currentTimeMillis();
    nstp2D.populate(0, 0, 1000, true);
    long stopTime = System.currentTimeMillis();
    System.out.println("Time taken = " + (stopTime - startTime) / 1000 + " seconds");
    for (int i = 0; i < 10; i++) {
        Point2D start = new Point2D.Double(Math.random(), Math.random());
        Point2D end = new Point2D.Double(Math.random(), Math.random());
        double time = Math.random();
        double p = nstp2D.getProb(start, end, time);
        System.out.println("Pr(" + end.getX() + ", " + end.getY() + " | " + start.getX() + ", " + start.getY() + ", t=" + time + ") = " + p);
    }
}

Example 5 with MultivariateNormalDistribution

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

the class EllipticalSliceOperator method main.

public static void main(String[] arg) {
    // Define normal model
    // Starting values
    Parameter thetaParameter = new Parameter.Default(new double[] { 1.0, 0.0 });
    MaskedParameter meanParameter = new MaskedParameter(thetaParameter, new Parameter.Default(new double[] { 1.0, 0.0 }), true);
    TransformedParameter precParameter = new TransformedParameter(new MaskedParameter(thetaParameter, new Parameter.Default(new double[] { 0.0, 1.0 }), true), new Transform.LogTransform(), true);
    // System.err.println(thetaParameter);
    // System.err.println(meanParameter);
    // System.err.println(precParameter);
    ParametricDistributionModel densityModel = new NormalDistributionModel(meanParameter, precParameter, true);
    DistributionLikelihood likelihood = new DistributionLikelihood(densityModel);
    // Define prior
    MultivariateNormalDistribution priorDistribution = new MultivariateNormalDistribution(new double[] { 0.0, 0.0 }, new double[][] { { 0.001, 0.0 }, { 0.0, 0.001 } });
    MultivariateDistributionLikelihood prior = new MultivariateDistributionLikelihood(priorDistribution);
    prior.addData(thetaParameter);
    // Define data
    // likelihood.addData(new Attribute.Default<double[]>("Data", new double[] {0.0, 2.0, 4.0}));
    likelihood.addData(new Attribute.Default<double[]>("Data", new double[] { 1, 2, 3, 4, 5, 6, 7, 8, 9 }));
    List<Likelihood> list = new ArrayList<Likelihood>();
    list.add(likelihood);
    list.add(prior);
    CompoundLikelihood posterior = new CompoundLikelihood(0, list);
    EllipticalSliceOperator sliceSampler = new EllipticalSliceOperator(thetaParameter, priorDistribution, false, true);
    final int dim = thetaParameter.getDimension();
    final int length = 100000;
    double[] mean = new double[dim];
    double[] variance = new double[dim];
    Parameter[] log = new Parameter[dim];
    log[0] = meanParameter;
    log[1] = precParameter;
    for (int i = 0; i < length; i++) {
        sliceSampler.doOperation(posterior);
        for (int j = 0; j < dim; ++j) {
            double x = log[j].getValue(0);
            mean[j] += x;
            variance[j] += x * x;
        }
    }
    for (int j = 0; j < dim; ++j) {
        mean[j] /= length;
        variance[j] /= length;
        variance[j] -= mean[j] * mean[j];
    }
    System.out.println("E(x)\tStErr(x)");
    for (int j = 0; j < dim; ++j) {
        System.out.println(mean[j] + " " + Math.sqrt(variance[j]));
    }
}