Author: erans
Date: Wed Aug 8 21:52:22 2012
New Revision: 1370984
URL: http://svn.apache.org/viewvc?rev=1370984&view=rev
Log:
MATH-815
Code update. Unit test.
Added:
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/MultivariateNormalDistributionTest.java
(with props)
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/MultivariateNormalDistribution.java
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/MultivariateNormalDistribution.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/MultivariateNormalDistribution.java?rev=1370984&r1=1370983&r2=1370984&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/MultivariateNormalDistribution.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/MultivariateNormalDistribution.java
Wed Aug 8 21:52:22 2012
@@ -1,7 +1,6 @@
package org.apache.commons.math3.distribution;
import org.apache.commons.math3.exception.DimensionMismatchException;
-import org.apache.commons.math3.exception.NotStrictlyPositiveException;
import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.EigenDecomposition;
import org.apache.commons.math3.linear.NonPositiveDefiniteMatrixException;
@@ -9,7 +8,6 @@ import org.apache.commons.math3.linear.R
import org.apache.commons.math3.linear.SingularMatrixException;
import org.apache.commons.math3.random.RandomGenerator;
import org.apache.commons.math3.random.Well19937c;
-import org.apache.commons.math3.stat.correlation.Covariance;
import org.apache.commons.math3.util.FastMath;
import org.apache.commons.math3.util.MathArrays;
@@ -143,7 +141,7 @@ public class MultivariateNormalDistribut
public RealMatrix getCovariances() {
return covarianceMatrix.copy();
}
-
+
/** {@inheritDoc} */
public double density(final double[] vals) throws
DimensionMismatchException {
final int dim = getDimensions();
@@ -151,11 +149,9 @@ public class MultivariateNormalDistribut
throw new DimensionMismatchException(vals.length, dim);
}
- final double kernel = getKernel(vals);
-
return FastMath.pow(2 * FastMath.PI, -dim / 2) *
FastMath.pow(covarianceMatrixDeterminant, -0.5) *
- FastMath.exp(kernel);
+ getExponentTerm(vals);
}
/**
@@ -193,19 +189,21 @@ public class MultivariateNormalDistribut
}
/**
- * Precomputes some of the multiplications used for determining densities.
+ * Computes the term used in the exponent (see definition of the
distribution).
*
* @param values Values at which to compute density.
* @return the multiplication factor of density calculations.
*/
- private double getKernel(final double[] values) {
- double k = 0;
- for (int col = 0; col < values.length; col++) {
- for (int v = 0; v < values.length; v++) {
- k += covarianceMatrixInverse.getEntry(v, col)
- * FastMath.pow(values[v] - means[v], 2);
- }
+ private double getExponentTerm(final double[] values) {
+ final double[] centered = new double[values.length];
+ for (int i = 0; i < centered.length; i++) {
+ centered[i] = values[i] - getMeans()[i];
+ }
+ final double[] preMultiplied =
covarianceMatrixInverse.preMultiply(centered);
+ double sum = 0;
+ for (int i = 0; i < preMultiplied.length; i++) {
+ sum += preMultiplied[i] * centered[i];
}
- return -0.5 * k;
+ return FastMath.exp(-0.5 * sum);
}
}
Added:
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/MultivariateNormalDistributionTest.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/MultivariateNormalDistributionTest.java?rev=1370984&view=auto
==============================================================================
---
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/MultivariateNormalDistributionTest.java
(added)
+++
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/MultivariateNormalDistributionTest.java
Wed Aug 8 21:52:22 2012
@@ -0,0 +1,134 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.math3.distribution;
+
+import org.apache.commons.math3.stat.correlation.Covariance;
+import org.apache.commons.math3.linear.RealMatrix;
+
+import org.junit.After;
+import org.junit.Assert;
+import org.junit.Before;
+import org.junit.Test;
+
+/**
+ * Test cases for {@link MultivariateNormalDistribution}.
+ */
+public class MultivariateNormalDistributionTest {
+ /**
+ * Test the ability of the distribution to report its mean value parameter.
+ */
+ @Test
+ public void testGetMean() {
+ final double[] mu = { -1.5, 2 };
+ final double[][] sigma = { { 2, -1.1 },
+ { -1.1, 2 } };
+ final MultivariateNormalDistribution d = new
MultivariateNormalDistribution(mu, sigma);
+
+ final double[] m = d.getMeans();
+ for (int i = 0; i < m.length; i++) {
+ Assert.assertEquals(mu[i], m[i], 0);
+ }
+ }
+
+ /**
+ * Test the ability of the distribution to report its covariance matrix
parameter.
+ */
+ @Test
+ public void testGetCovarianceMatrix() {
+ final double[] mu = { -1.5, 2 };
+ final double[][] sigma = { { 2, -1.1 },
+ { -1.1, 2 } };
+ final MultivariateNormalDistribution d = new
MultivariateNormalDistribution(mu, sigma);
+
+ final RealMatrix s = d.getCovariances();
+ final int dim = d.getDimensions();
+ for (int i = 0; i < dim; i++) {
+ for (int j = 0; j < dim; j++) {
+ Assert.assertEquals(sigma[i][j], s.getEntry(i, j), 0);
+ }
+ }
+ }
+
+ /**
+ * Test the accuracy of sampling from the distribution.
+ */
+ @Test
+ public void testSampling() {
+ final double[] mu = { -1.5, 2 };
+ final double[][] sigma = { { 2, -1.1 },
+ { -1.1, 2 } };
+ final MultivariateNormalDistribution d = new
MultivariateNormalDistribution(mu, sigma);
+ d.reseedRandomGenerator(50);
+
+ final int n = 30;
+
+ final double[][] samples = d.sample(n);
+ final int dim = d.getDimensions();
+ final double[] sampleMeans = new double[dim];
+
+ for (int i = 0; i < samples.length; i++) {
+ for (int j = 0; j < dim; j++) {
+ sampleMeans[j] += samples[i][j];
+ }
+ }
+
+ final double sampledMeanTolerance = 1e-1;
+ for (int j = 0; j < dim; j++) {
+ sampleMeans[j] /= samples.length;
+ Assert.assertEquals(mu[j], sampleMeans[j], sampledMeanTolerance);
+ }
+
+ final double sampledCovarianceTolerance = 2;
+ final double[][] sampleSigma = new
Covariance(samples).getCovarianceMatrix().getData();
+ for (int i = 0; i < dim; i++) {
+ for (int j = 0; j < dim; j++) {
+ Assert.assertEquals(sigma[i][j], sampleSigma[i][j],
sampledCovarianceTolerance);
+ }
+ }
+ }
+
+ /**
+ * Test the accuracy of the distribution when calculating densities.
+ */
+ @Test
+ public void testDensities() {
+ final double[] mu = { -1.5, 2 };
+ final double[][] sigma = { { 2, -1.1 },
+ { -1.1, 2 } };
+ final MultivariateNormalDistribution d = new
MultivariateNormalDistribution(mu, sigma);
+
+ final double[][] testValues = { { -1.5, 2 },
+ { 4, 4 },
+ { 1.5, -2 },
+ { 0, 0 } };
+ final double[] densities = new double[testValues.length];
+ for (int i = 0; i < densities.length; i++) {
+ densities[i] = d.density(testValues[i]);
+ }
+
+ // From dmvnorm function in R 2.15 CRAN package Mixtools v0.4.5
+ final double[] correctDensities = { 0.09528357207691344,
+ 5.80932710124009e-09,
+ 0.001387448895173267,
+ 0.03309922090210541 };
+
+ for (int i = 0; i < testValues.length; i++) {
+ Assert.assertEquals(correctDensities[i], densities[i], 1e-16);
+ }
+ }
+}
Propchange:
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/MultivariateNormalDistributionTest.java
------------------------------------------------------------------------------
svn:eol-style = native
