Merge 5a72a46 into 9aa0bd8

docs/examples/ex_bivariate_likelihood.rst

@@ -0,0 +1,26 @@
+Example: Likelihood fitting a Bivariate Gaussian
+================================================
+
+In this example, we shall perform likelihood fitting to a `bivariate normal
+distribution <http://mathworld.wolfram.com/BivariateNormalDistribution.html>`_,
+to demonstrate how ``symfit``'s API can easily be used to perform likelihood
+fitting on multivariate problems.
+
+In this example, we sample from a bivariate normal distribution with a
+significant correlation of :math:`\rho = 0.6` between :math:`x` and :math:`y`.
+We see that this is extracted from the data relatively straightforwardly.
+
+.. literalinclude:: ../../examples/bivariate_likelihood.py
+    :language: python
+
+This code prints::
+
+    Parameter Value        Standard Deviation
+	rho       6.026420e-01 2.013810e-03
+	sig_x     1.100898e-01 2.461684e-04
+	sig_y     2.303400e-01 5.150556e-04
+	x0        5.901317e-01 3.481346e-04
+	y0        8.014040e-01 7.283990e-04
+	Fitting status message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
+	Number of iterations:   35
+	Regression Coefficient: nan

docs/examples/index.rst

@@ -19,6 +19,17 @@ paper.
     ex_ODEModel
     ex_CallableNumericalModel
 
+Objective Examples
+------------------
+
+These are examples on how to use a different objective function from the
+standard least-squares objective, such as likelihood fitting.
+
+.. toctree::
+    :maxdepth: 1
+
+    ex_bivariate_likelihood
+
 Interactive Guess Module
 ------------------------
 

examples/bivariate_likelihood.py

@@ -0,0 +1,27 @@
+import numpy as np
+from symfit import Variable, Parameter, Fit
+from symfit.core.objectives import LogLikelihood
+from symfit.distributions import BivariateGaussian
+
+x = Variable('x')
+y = Variable('y')
+x0 = Parameter('x0', value=0.6, min=0.5, max=0.7)
+sig_x = Parameter('sig_x', value=0.1, max=1.0)
+y0 = Parameter('y0', value=0.7, min=0.6, max=0.9)
+sig_y = Parameter('sig_y', value=0.05, max=1.0)
+rho = Parameter('rho', value=0.001, min=-1, max=1)
+
+pdf = BivariateGaussian(x=x, mu_x=x0, sig_x=sig_x, y=y, mu_y=y0,
+                       sig_y=sig_y, rho=rho)
+
+# Draw 100000 samples from a bivariate distribution
+mean = [0.59, 0.8]
+corr = 0.6
+cov = np.array([[0.11 ** 2, 0.11 * 0.23 * corr],
+                [0.11 * 0.23 * corr, 0.23 ** 2]])
+np.random.seed(42)
+xdata, ydata = np.random.multivariate_normal(mean, cov, 100000).T
+
+fit = Fit(pdf, x=xdata, y=ydata, objective=LogLikelihood)
+fit_result = fit.execute()
+print(fit_result)

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) * (y - mu_y) / (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_general.py

@@ -14,7 +14,7 @@
 )
 from symfit.core.minimizers import MINPACK, LBFGSB, BoundedMinimizer, DifferentialEvolution
 from symfit.core.objectives import LogLikelihood
-from symfit.distributions import Gaussian, Exp
+from symfit.distributions import Gaussian, Exp, BivariateGaussian
 from tests.test_minimizers import subclasses
 
 if sys.version_info >= (3, 0):
@@ -620,6 +620,39 @@ def test_likelihood_fitting_gaussian(self):
         self.assertAlmostEqual(fit_result.stdev(mu) / mean_stdev, 1, 3)
         self.assertAlmostEqual(fit_result.value(sig) / np.std(xdata), 1, 6)
 
+    def test_likelihood_fitting_bivariate_gaussian(self):
+        """
+        Fit using the likelihood method.
+        """
+        # Make variables and parameters
+        x = Variable('x')
+        y = Variable('y')
+        x0 = Parameter('x0', value=0.6, min=0.5, max=0.7)
+        sig_x = Parameter('sig_x', value=0.1, max=1.0)
+        y0 = Parameter('y0', value=0.7, min=0.6, max=0.9)
+        sig_y = Parameter('sig_y', value=0.05, max=1.0)
+        rho = Parameter('rho', value=0.001, min=-1, max=1)
+
+        pdf = BivariateGaussian(x=x, mu_x=x0, sig_x=sig_x, y=y, mu_y=y0,
+                               sig_y=sig_y, rho=rho)
+
+        # Draw 100000 samples from a bivariate distribution
+        mean = [0.59, 0.8]
+        r = 0.6
+        cov = np.array([[0.11 ** 2, 0.11 * 0.23 * r],
+                        [0.11 * 0.23 * r, 0.23 ** 2]])
+        np.random.seed(42)
+        xdata, ydata = np.random.multivariate_normal(mean, cov, 100000).T
+
+        fit = Fit(pdf, x=xdata, y=ydata, objective=LogLikelihood)
+        fit_result = fit.execute()
+
+        self.assertAlmostEqual(fit_result.value(x0) / mean[0], 1, 2)
+        self.assertAlmostEqual(fit_result.value(y0) / mean[1], 1, 2)
+        self.assertAlmostEqual(fit_result.value(sig_x) / np.sqrt(cov[0, 0]), 1, 2)
+        self.assertAlmostEqual(fit_result.value(sig_y) / np.sqrt(cov[1, 1]), 1, 2)
+        self.assertAlmostEqual(fit_result.value(rho) / r, 1, 2)
+
     def test_evaluate_model(self):
         """
         Makes sure that models are callable and give the expected answer.