All Downloads are FREE. Search and download functionalities are using the official Maven repository.

gov.sandia.cognition.statistics.bayesian.conjugate.MultivariateGaussianMeanCovarianceBayesianEstimator Maven / Gradle / Ivy

There is a newer version: 4.0.1
Show newest version
/*
 * File:                MultivariateGaussianMeanCovarianceBayesianEstimator.java
 * Authors:             Kevin R. Dixon
 * Company:             Sandia National Laboratories
 * Project:             Cognitive Foundry
 * 
 * Copyright Mar 26, 2010, Sandia Corporation.
 * Under the terms of Contract DE-AC04-94AL85000, there is a non-exclusive
 * license for use of this work by or on behalf of the U.S. Government.
 * Export of this program may require a license from the United States
 * Government. See CopyrightHistory.txt for complete details.
 * 
 */

package gov.sandia.cognition.statistics.bayesian.conjugate;

import gov.sandia.cognition.annotation.PublicationReference;
import gov.sandia.cognition.annotation.PublicationReferences;
import gov.sandia.cognition.annotation.PublicationType;
import gov.sandia.cognition.collection.CollectionUtil;
import gov.sandia.cognition.math.MultivariateStatisticsUtil;
import gov.sandia.cognition.math.RingAccumulator;
import gov.sandia.cognition.math.matrix.Matrix;
import gov.sandia.cognition.math.matrix.MatrixFactory;
import gov.sandia.cognition.math.matrix.Vector;
import gov.sandia.cognition.statistics.bayesian.AbstractBayesianParameter;
import gov.sandia.cognition.statistics.bayesian.BayesianParameter;
import gov.sandia.cognition.statistics.distribution.MultivariateGaussian;
import gov.sandia.cognition.statistics.distribution.MultivariateStudentTDistribution;
import gov.sandia.cognition.statistics.distribution.NormalInverseWishartDistribution;
import gov.sandia.cognition.util.Pair;
import java.util.Arrays;

/**
 * Performs robust estimation of both the mean and covariance of a
 * MultivariateGaussian conditional distribution using the conjugate prior
 * Normal-Inverse-Wishart distribution.  The resulting predictive distribution
 * for future data is a multivariate Student's t-distribution.
 * @author Kevin R. Dixon
 * @since 3.0
 */
@PublicationReferences(
    references={
        @PublicationReference(
            author={
                "Andrew Gelman",
                "John B. Carlin",
                "Hal S. Stern",
                "Donald B. Rubin"
            },
            title="Bayesian Data Analysis, Second Edition",
            type=PublicationType.Book,
            year=2004,
            pages={87,88}
        )
        ,
        @PublicationReference(
            author="Wikipedia",
            title="Conjugate Prior",
            type=PublicationType.WebPage,
            year=2009,
            url="http://en.wikipedia.org/wiki/Conjugate_prior"
        )
    }
)
public class MultivariateGaussianMeanCovarianceBayesianEstimator 
    extends AbstractConjugatePriorBayesianEstimator
    implements ConjugatePriorBayesianEstimatorPredictor
{

    /** 
     * Creates a new instance of MultivariateGaussianMeanCovarianceBayesianEstimator 
     */
    public MultivariateGaussianMeanCovarianceBayesianEstimator()
    {
        this( new NormalInverseWishartDistribution() );
    }

    /**
     * Creates a new instance of MultivariateGaussianMeanCovarianceBayesianEstimator 
     * @param dimensionality
     * Dimensionality of the observations to consider.
     */
    public MultivariateGaussianMeanCovarianceBayesianEstimator(
        int dimensionality )
    {
        this( new NormalInverseWishartDistribution( dimensionality ) );
    }

    /** 
     * Creates a new instance of MultivariateGaussianMeanCovarianceBayesianEstimator 
     * @param belief 
     * Initial belief of the conditional parameters
     */
    public MultivariateGaussianMeanCovarianceBayesianEstimator(
        NormalInverseWishartDistribution belief )
    {
        this( new MultivariateGaussian( belief.getInputDimensionality() ), belief );
    }

    /**
     * Creates a new instance
     * @param prior
     * Default conjugate prior.
     * @param conditional
     * Conditional distribution of the conjugate prior.
     */
    public MultivariateGaussianMeanCovarianceBayesianEstimator(
        MultivariateGaussian conditional,
        NormalInverseWishartDistribution prior )
    {
        this( new MultivariateGaussianMeanCovarianceBayesianEstimator.Parameter(
            conditional, prior) );
    }
    
    /**
     * Creates a new instance
     * @param parameter
     * Bayesian hyperparameter relationship between the conditional
     * distribution and the conjugate prior distribution.
     */
    protected MultivariateGaussianMeanCovarianceBayesianEstimator(
        BayesianParameter parameter )
    {
        super( parameter );
    }

    public MultivariateGaussianMeanCovarianceBayesianEstimator.Parameter createParameter(
        MultivariateGaussian conditional,
        NormalInverseWishartDistribution prior)
    {
        return new MultivariateGaussianMeanCovarianceBayesianEstimator.Parameter( conditional, prior );
    }

    public void update(
        NormalInverseWishartDistribution target,
        Vector data)
    {
        this.update(target, Arrays.asList(data) );
    }

    @Override
    public void update(
        NormalInverseWishartDistribution prior,
        Iterable data)
    {
        final int n = CollectionUtil.size(data);

        Pair pair =
            MultivariateStatisticsUtil.computeMeanAndCovariance(data);
        Vector sampleMean = pair.getFirst();
        Matrix sampleCovariance = pair.getSecond();

        Vector lambda = prior.getGaussian().getMean();
        double nu = prior.getCovarianceDivisor();
        int alpha = prior.getInverseWishart().getDegreesOfFreedom();
        Matrix beta = prior.getInverseWishart().getInverseScale();

        int alphahat = alpha + n;
        double nuhat = nu+n;

        Vector lambdahat = lambda.scale(nu/n);
        lambdahat.plusEquals( sampleMean );
        lambdahat.scaleEquals( n/nuhat );

        Vector delta = sampleMean;
        delta.minusEquals(lambda);
        Matrix betahat = sampleCovariance;
        if( n > 1 )
        {
            betahat.scaleEquals(n);
        }
        betahat.plusEquals(beta);
        betahat.plusEquals( delta.outerProduct(delta.scale((n*nu)/nuhat)) );

        prior.getGaussian().setMean(lambdahat);
        prior.setCovarianceDivisor(nuhat);
        prior.getInverseWishart().setDegreesOfFreedom(alphahat);
        prior.getInverseWishart().setInverseScale(betahat);

    }

    public double computeEquivalentSampleSize(
        NormalInverseWishartDistribution belief)
    {
        return belief.getCovarianceDivisor();
    }

    public MultivariateStudentTDistribution createPredictiveDistribution(
        NormalInverseWishartDistribution posterior)
    {
        Vector mean = posterior.getGaussian().getMean();
        double dofs = posterior.getInverseWishart().getDegreesOfFreedom()
            - posterior.getInverseWishart().getInputDimensionality() + 1.0;

        Matrix covariance = posterior.getInverseWishart().getInverseScale().scale(
            (posterior.getCovarianceDivisor()+1.0) / (posterior.getCovarianceDivisor()*dofs) );
        Matrix precision = covariance.inverse();
        return new MultivariateStudentTDistribution( dofs, mean, precision );
    }

    /**
     * Parameter for this conjugate prior estimator.
     */
    public static class Parameter
        extends AbstractBayesianParameter
    {

        /**
         * Name of the parameter, {@value}.
         */
        public static final String NAME = "meanAndCovariance";

        /**
         * Creates a new instance
         * @param prior
         * Default conjugate prior.
         * @param conditional
         * Conditional distribution of the conjugate prior.
         */
        public Parameter(
            MultivariateGaussian conditional,
            NormalInverseWishartDistribution prior )
        {
            super( conditional, NAME, prior );
        }

        public void setValue(
            Matrix value)
        {
            final int dim = this.conditionalDistribution.getInputDimensionality();
            if( (value.getNumRows() != dim) ||
                (value.getNumColumns() != (dim+1) ) )
            {
                throw new IllegalArgumentException( "Expected (dim x dim+1) Matrix" );
            }
            Vector mean = value.getColumn(0);
            Matrix covariance = value.getSubMatrix(0, dim-1, 1,dim);
            this.conditionalDistribution.setMean(mean);
            this.conditionalDistribution.setCovariance(covariance);
        }

        public Matrix getValue()
        {
            final int dim = this.conditionalDistribution.getInputDimensionality();
            Matrix parameter = MatrixFactory.getDefault().createMatrix(dim, dim+1);
            parameter.setColumn(0, this.conditionalDistribution.getMean() );
            parameter.setSubMatrix(0,1, this.conditionalDistribution.getCovariance() );
            return parameter;
        }        

    }

}




© 2015 - 2024 Weber Informatics LLC | Privacy Policy