Bug fix in MatrixNormal density with unit test for non-square matrices

alanisaac · Sep 6, 2014 · aeff418 · aeff418
1 parent 45de4a9
commit aeff418
Show file tree

Hide file tree

Showing 2 changed files with 32 additions and 2 deletions.
diff --git a/src/Numerics/Distributions/MatrixNormal.cs b/src/Numerics/Distributions/MatrixNormal.cs
@@ -207,8 +207,8 @@ public double Density(Matrix<double> x)
 
             return Math.Exp(-0.5*cholK.Solve(a.Transpose()*cholV.Solve(a)).Trace())
                    /Math.Pow(2.0*Constants.Pi, x.RowCount*x.ColumnCount/2.0)
-                   /Math.Pow(cholV.Determinant, x.RowCount/2.0)
-                   /Math.Pow(cholK.Determinant, x.ColumnCount/2.0);
+                   /Math.Pow(cholK.Determinant, x.RowCount/2.0)
+                   /Math.Pow(cholV.Determinant, x.ColumnCount/2.0);
         }
 
         /// <summary>

diff --git a/src/UnitTests/DistributionTests/Multivariate/MatrixNormalTests.cs b/src/UnitTests/DistributionTests/Multivariate/MatrixNormalTests.cs
@@ -250,6 +250,36 @@ public void ValidateDensity()
             AssertHelpers.AlmostEqualRelative(0.00015682927366491211, d.Density(x), 16);
         }
 
+        /// <summary>
+        /// Validate density with non-square matrices.
+        /// </summary>
+        [Test]
+        public void ValidateNonsquareDensity()
+        {
+            const int Rows = 2;
+            const int Cols = 1;
+            var m = Matrix<double>.Build.Dense(Rows, Cols);
+            m[0, 0] = 0.156065579983862;
+            m[1, 0] = -0.806288628097313;
+
+            var v = Matrix<double>.Build.Dense(Rows, Rows);
+            v[0, 0] = 0.674457817054746;
+            v[0, 1] = 0.878930403442185;
+            v[1, 0] = 0.878930403442185;
+            v[1, 1] = 1.76277498368061;
+
+            var k = Matrix<double>.Build.Dense(Cols, Cols);
+            k[0, 0] = 0.674457817054746;
+
+            var d = new MatrixNormal(m, v, k);
+
+            var x = Matrix<double>.Build.Dense(Rows, Cols);
+            x[0, 0] = 2;
+            x[1, 0] = 1.5;
+
+            AssertHelpers.AlmostEqualRelative(0.008613384131384546, d.Density(x), 12);
+        }
+
         /// <summary>
         /// Can sample.
         /// </summary>