commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From er...@apache.org
Subject svn commit: r1370984 - in /commons/proper/math/trunk/src: main/java/org/apache/commons/math3/distribution/MultivariateNormalDistribution.java test/java/org/apache/commons/math3/distribution/MultivariateNormalDistributionTest.java
Date Wed, 08 Aug 2012 21:52:22 GMT
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



Mime
View raw message