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}