001/**
002 * Copyright (c) 2011, The University of Southampton and the individual contributors.
003 * All rights reserved.
004 *
005 * Redistribution and use in source and binary forms, with or without modification,
006 * are permitted provided that the following conditions are met:
007 *
008 *   *  Redistributions of source code must retain the above copyright notice,
009 *      this list of conditions and the following disclaimer.
010 *
011 *   *  Redistributions in binary form must reproduce the above copyright notice,
012 *      this list of conditions and the following disclaimer in the documentation
013 *      and/or other materials provided with the distribution.
014 *
015 *   *  Neither the name of the University of Southampton nor the names of its
016 *      contributors may be used to endorse or promote products derived from this
017 *      software without specific prior written permission.
018 *
019 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
020 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
021 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
022 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
023 * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
024 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
025 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
026 * ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
027 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
028 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
029 */
030package org.openimaj.math.statistics.distribution;
031
032import java.util.Arrays;
033import java.util.Random;
034
035import org.openimaj.math.matrix.MatrixUtils;
036
037import Jama.Matrix;
038import cern.jet.random.Normal;
039import cern.jet.random.engine.MersenneTwister;
040
041/**
042 * Implementation of a {@link MultivariateGaussian} with a diagonal covariance
043 * matrix.
044 *
045 * @author Jonathon Hare (jsh2@ecs.soton.ac.uk)
046 */
047public class DiagonalMultivariateGaussian extends AbstractMultivariateGaussian {
048        /**
049         * The diagonal of the covariance matrix
050         */
051        public double[] variance;
052
053        /**
054         * Construct the Gaussian with the provided center and covariance
055         *
056         * @param mean
057         *            centre of the Gaussian
058         * @param variance
059         *            variance of the Gaussian
060         */
061        public DiagonalMultivariateGaussian(Matrix mean, double[] variance) {
062                this.mean = mean;
063                this.variance = variance;
064        }
065
066        /**
067         * Construct the Gaussian with the zero mean and unit variance
068         *
069         * @param ndims
070         *            number of dimensions
071         */
072        public DiagonalMultivariateGaussian(int ndims) {
073                this.mean = new Matrix(1, ndims);
074                this.variance = new double[ndims];
075                Arrays.fill(variance, 1);
076        }
077
078        @Override
079        public Matrix getCovariance() {
080                return MatrixUtils.diag(variance);
081        }
082
083        @Override
084        public double getCovariance(int row, int col) {
085                if (row < 0 || row >= variance.length || col < 0 || col > variance.length)
086                        throw new IndexOutOfBoundsException();
087
088                if (row == col)
089                        return variance[row];
090                return 0;
091        }
092
093        @Override
094        public double estimateProbability(double[] sample) {
095                final int N = this.variance.length;
096                final double[] meanvector = mean.getArray()[0];
097
098                double det = variance[0];
099                for (int i = 1; i < N; i++)
100                        det *= variance[i];
101                final double pdf_const_factor = 1.0 / Math.sqrt((Math.pow((2 * Math.PI), N) * det));
102
103                double v = 0;
104                for (int i = 0; i < N; i++) {
105                        final double diff = sample[i] - meanvector[i];
106                        v += diff * diff / variance[i];
107                }
108
109                return pdf_const_factor * Math.exp(-0.5 * v);
110        }
111
112        @Override
113        public double estimateLogProbability(double[] sample) {
114                final int N = this.variance.length;
115                final double[] meanvector = mean.getArray()[0];
116
117                double log_sqrt_det = Math.log(Math.sqrt(variance[0]));
118                for (int i = 1; i < N; i++)
119                        log_sqrt_det += Math.log(Math.sqrt(variance[i]));
120                final double log_pdf_const_factor = -Math.log(Math.sqrt((Math.pow((2 * Math.PI), N)))) - log_sqrt_det;
121
122                double v = 0;
123                for (int i = 0; i < N; i++) {
124                        final double diff = sample[i] - meanvector[i];
125                        v += diff * diff / variance[i];
126                }
127
128                return log_pdf_const_factor + (-0.5 * v);
129        }
130
131        @Override
132        public double[] estimateLogProbability(double[][] samples) {
133                final int N = this.variance.length;
134                final double[] meanvector = mean.getArray()[0];
135
136                double log_sqrt_det = Math.log(Math.sqrt(variance[0]));
137                for (int i = 1; i < N; i++)
138                        log_sqrt_det += Math.log(Math.sqrt(variance[i]));
139                final double log_pdf_const_factor = -Math.log(Math.sqrt((Math.pow((2 * Math.PI), N)))) - log_sqrt_det;
140
141                final double[] lp = new double[samples.length];
142                for (int j = 0; j < samples.length; j++) {
143                        double v = 0;
144                        for (int i = 0; i < N; i++) {
145                                final double diff = samples[j][i] - meanvector[i];
146                                v += diff * diff / variance[i];
147                        }
148                        lp[j] = log_pdf_const_factor + (-0.5 * v);
149                }
150
151                return lp;
152        }
153
154        @Override
155        public double[][] sample(int nsamples, Random rng) {
156                if (nsamples == 0)
157                        return new double[0][0];
158
159                final Normal rng2 = new Normal(0, 1, new MersenneTwister());
160
161                final int N = mean.getColumnDimension();
162                final double[][] out = new double[nsamples][N];
163
164                final double[] meanv = mean.getArray()[0];
165                for (int i = 0; i < N; i++) {
166                        final double choli = Math.sqrt(this.variance[i]);
167
168                        for (int j = 0; j < nsamples; j++) {
169                                out[j][i] = choli * rng2.nextDouble() + meanv[i];
170                        }
171                }
172
173                return out;
174        }
175}