Improved TestFitResults.test_fitting_2 by adding a bivariate normal d…

…istribution.

symfit/distributions.py

@@ -18,6 +18,27 @@ def Gaussian(x, mu, sig):
     """
     return sympy.exp(-(x - mu)**2/(2*sig**2))/sympy.sqrt(2*sympy.pi*sig**2)
 
+def BivariateGaussian(x, y, mu_x, mu_y, sig_x, sig_y, rho):
+    """
+    Bivariate Gaussian pdf.
+
+    :param x: :class:`symfit.core.argument.Variable`
+    :param y: :class:`symfit.core.argument.Variable`
+    :param mu_x: :class:`symfit.core.argument.Parameter` for the mean of `x`
+    :param mu_y: :class:`symfit.core.argument.Parameter` for the mean of `y`
+    :param sig_x: :class:`symfit.core.argument.Parameter` for the standard
+        deviation of `x`
+    :param sig_y: :class:`symfit.core.argument.Parameter` for the standard
+        deviation of `y`
+    :param rho: :class:`symfit.core.argument.Parameter` for the correlation
+        between `x` and `y`.
+    :return: sympy expression for a Bivariate Gaussian pdf.
+    """
+    exponent = - 1 / (2 * (1 - rho**2))
+    exponent *= (x - mu_x)**2 / sig_x**2 + (y - mu_y)**2 / sig_y**2 \
+                - 2 * rho * (x - mu_x) * (x - mu_x) / (sig_x * sig_y)
+    return sympy.exp(exponent) / (2 * sympy.pi * sig_x * sig_y * sympy.sqrt(1 - rho**2))
+
 def Exp(x, l):
     """
     Exponential Distribution pdf.

tests/test_fit_result.py

@@ -6,8 +6,8 @@
 from collections import OrderedDict
 
 import numpy as np
-from symfit import Variable, Parameter, Fit, FitResults
-from symfit.distributions import Gaussian
+from symfit import Variable, Parameter, Fit, FitResults, GradientModel
+from symfit.distributions import Gaussian, BivariateGaussian
 from symfit.core.minimizers import MINPACK
 
 class TestFitResults(unittest.TestCase):
@@ -66,55 +66,68 @@ def test_fitting(self):
         self.assertRaises(AttributeError, getattr, *[fit_result.params, 'a__stdev'])
 
     def test_fitting_2(self):
-        np.random.seed(4242)
-        mean = (0.33, 0.28)  # x, y mean 0.3, 0.3
+        np.random.seed(43)
+        mean = (0.62, 0.71)  # x, y mean 0.7, 0.7
         cov = [
-            [0.01**2, 0.39],
-            [0.41, 0.0101**2]
+            [0.102**2, 0],
+            [0, 0.07**2]
+        ]
+        data_1 = np.random.multivariate_normal(mean, cov, 10**5)
+        mean = (0.33, 0.28)  # x, y mean 0.3, 0.3
+        cov = [  # rho = 0.25
+            [0.05 ** 2, 0.25 * 0.05 * 0.101],
+            [0.25 * 0.05 * 0.101, 0.101 ** 2]
         ]
-        data = np.random.multivariate_normal(mean, cov, 1000000)
-        mean = (0.68, 0.71)  # x, y mean 0.6, 0.4
-        cov = [[0.0102**2, 0.00], [1e-9, 0.010**2]]
-        data_2 = np.random.multivariate_normal(mean, cov, 1000000)
-        data = np.vstack((data, data_2))
+        data_2 = np.random.multivariate_normal(mean, cov, 10**5)
+        data = np.vstack((data_1, data_2))
 
         # Insert them as y,x here as np fucks up cartesian conventions.
-        ydata, xedges, yedges = np.histogram2d(data[:,1], data[:,0], bins=200, range=[[0.0, 1.0], [0.0, 1.0]])
+        ydata, xedges, yedges = np.histogram2d(data[:, 1], data[:, 0], bins=200,
+                                               range=[[0.0, 1.0], [0.0, 1.0]],
+                                               density=True)
         xcentres = (xedges[:-1] + xedges[1:]) / 2
         ycentres = (yedges[:-1] + yedges[1:]) / 2
 
         # Make a valid grid to match ydata
         xx, yy = np.meshgrid(xcentres, ycentres, sparse=False)
-        # xdata = np.dstack((xx, yy)).T
 
         x = Variable('x')
         y = Variable('y')
 
-        x0_1 = Parameter('x0_1', value=0.7, min=0.6, max=0.8)
-        sig_x_1 = Parameter('sig_x_1', value=0.01, min=0.0, max=0.2)
+        x0_1 = Parameter('x0_1', value=0.6, min=0.5, max=0.7)
+        sig_x_1 = Parameter('sig_x_1', value=0.1, min=0.0, max=0.2)
         y0_1 = Parameter('y0_1', value=0.7, min=0.6, max=0.8)
-        sig_y_1 = Parameter('sig_y_1', value=0.01, min=0.0, max=0.2)
-        A_1 = Parameter('A_1')
-        g_1 = A_1 * Gaussian(x, x0_1, sig_x_1) * Gaussian(y, y0_1, sig_y_1)
+        sig_y_1 = Parameter('sig_y_1', value=0.05, min=0.0, max=0.2)
+        rho_1 = Parameter('rho_1', value=0.0, min=-0.5, max=0.5)
+        A_1 = Parameter('A_1', value=0.5, min=0.3, max=0.7)
+        g_1 = A_1 * BivariateGaussian(x=x, y=y, mu_x=x0_1, mu_y=y0_1,
+                                      sig_x=sig_x_1, sig_y=sig_y_1, rho=rho_1)
 
         x0_2 = Parameter('x0_2', value=0.3, min=0.2, max=0.4)
-        sig_x_2 = Parameter('sig_x_2', value=0.01, min=0.0, max=0.2)
+        sig_x_2 = Parameter('sig_x_2', value=0.05, min=0.0, max=0.2)
         y0_2 = Parameter('y0_2', value=0.3, min=0.2, max=0.4)
-        sig_y_2 = Parameter('sig_y_2', value=0.01, min=0.0, max=0.2)
-        A_2 = Parameter('A_2')
-        g_2 = A_2 * Gaussian(x, x0_2, sig_x_2) * Gaussian(y, y0_2, sig_y_2)
+        sig_y_2 = Parameter('sig_y_2', value=0.1, min=0.0, max=0.2)
+        rho_2 = Parameter('rho_2', value=0.26, min=0.0, max=0.8)
+        A_2 = Parameter('A_2', value=0.5, min=0.3, max=0.7)
+        g_2 = A_2 * BivariateGaussian(x=x, y=y, mu_x=x0_2, mu_y=y0_2,
+                                      sig_x=sig_x_2, sig_y=sig_y_2, rho=rho_2)
 
         model = g_1 + g_2
         fit = Fit(model, xx, yy, ydata)
         fit_result = fit.execute()
 
+        self.assertGreater(fit_result.r_squared, 0.95)
         for param in fit.model.params:
-            self.assertAlmostEqual(fit_result.stdev(param)**2, fit_result.variance(param))
+            try:
+                self.assertAlmostEqual(fit_result.stdev(param)**2 / fit_result.variance(param), 1.0)
+            except AssertionError:
+                self.assertLessEqual(fit_result.variance(param), 0.0)
+                self.assertTrue(np.isnan(fit_result.stdev(param)))
 
         # Covariance matrix should be symmetric
         for param_1 in fit.model.params:
             for param_2 in fit.model.params:
-                self.assertAlmostEqual(fit_result.covariance(param_1, param_2) / fit_result.covariance(param_2, param_1), 1.0, 6)
+                self.assertAlmostEqual(fit_result.covariance(param_1, param_2) / fit_result.covariance(param_2, param_1), 1.0, 3)
 
     def test_minimizer_included(self):
         """"The minimizer used should be included in the results."""