From 4933aa4920e55ff45794f87cc2c2706df2dc49ab Mon Sep 17 00:00:00 2001 From: Martin Roelfs Date: Wed, 6 Mar 2019 18:12:34 +0100 Subject: [PATCH 1/5] Added a test for bivariate normal likelihood fitting --- symfit/distributions.py | 21 +++++++++++++++++++++ tests/test_general.py | 35 ++++++++++++++++++++++++++++++++++- 2 files changed, 55 insertions(+), 1 deletion(-) diff --git a/symfit/distributions.py b/symfit/distributions.py index 1c6780ae..23364055 100644 --- a/symfit/distributions.py +++ b/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. diff --git a/tests/test_general.py b/tests/test_general.py index e5093d17..c2347b6a 100644 --- a/tests/test_general.py +++ b/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) + + pi = BivariateGaussian(x=x, mu_x=x0, sig_x=sig_x, y=y, mu_y=y0, + sig_y=sig_y, rho=rho) + + # Draw 10000 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(pi, 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. From fc877d6eac55e3f1cd090628804f7aacb54627ec Mon Sep 17 00:00:00 2001 From: Martin Roelfs Date: Wed, 6 Mar 2019 18:16:19 +0100 Subject: [PATCH 2/5] Added example file for bivariate normal distribution likelihood --- examples/bivariate_likelihood.py | 26 ++++++++++++++++++++++++++ tests/test_general.py | 6 +++--- 2 files changed, 29 insertions(+), 3 deletions(-) create mode 100644 examples/bivariate_likelihood.py diff --git a/examples/bivariate_likelihood.py b/examples/bivariate_likelihood.py new file mode 100644 index 00000000..bd330224 --- /dev/null +++ b/examples/bivariate_likelihood.py @@ -0,0 +1,26 @@ +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] +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() \ No newline at end of file diff --git a/tests/test_general.py b/tests/test_general.py index c2347b6a..7dec40e3 100644 --- a/tests/test_general.py +++ b/tests/test_general.py @@ -633,10 +633,10 @@ def test_likelihood_fitting_bivariate_gaussian(self): sig_y = Parameter('sig_y', value=0.05, max=1.0) rho = Parameter('rho', value=0.001, min=-1, max=1) - pi = BivariateGaussian(x=x, mu_x=x0, sig_x=sig_x, y=y, mu_y=y0, + pdf = BivariateGaussian(x=x, mu_x=x0, sig_x=sig_x, y=y, mu_y=y0, sig_y=sig_y, rho=rho) - # Draw 10000 samples from a bivariate distribution + # 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], @@ -644,7 +644,7 @@ def test_likelihood_fitting_bivariate_gaussian(self): np.random.seed(42) xdata, ydata = np.random.multivariate_normal(mean, cov, 100000).T - fit = Fit(pi, x=xdata, y=ydata, objective=LogLikelihood) + fit = Fit(pdf, x=xdata, y=ydata, objective=LogLikelihood) fit_result = fit.execute() self.assertAlmostEqual(fit_result.value(x0) / mean[0], 1, 2) From 03ca5e27c7c410304c8b5736c3819fafd09ec9f6 Mon Sep 17 00:00:00 2001 From: Martin Roelfs Date: Wed, 6 Mar 2019 18:24:23 +0100 Subject: [PATCH 3/5] Included the example in the docs --- docs/examples/ex_bivariate_likelihood.rst | 22 ++++++++++++++++++++++ examples/bivariate_likelihood.py | 3 ++- 2 files changed, 24 insertions(+), 1 deletion(-) create mode 100644 docs/examples/ex_bivariate_likelihood.rst diff --git a/docs/examples/ex_bivariate_likelihood.rst b/docs/examples/ex_bivariate_likelihood.rst new file mode 100644 index 00000000..f65b7b8f --- /dev/null +++ b/docs/examples/ex_bivariate_likelihood.rst @@ -0,0 +1,22 @@ +Example: Likelihood fitting a Bivariate Gaussian +================================================ + +In this example, we shall perform likelihood fitting to a `bivariate normal +distrubtion `_, +to demonstrate how ``symfit``'s API can easily be used to perform likelihood +fitting on multivariate problems. + +.. 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 \ No newline at end of file diff --git a/examples/bivariate_likelihood.py b/examples/bivariate_likelihood.py index bd330224..52e14515 100644 --- a/examples/bivariate_likelihood.py +++ b/examples/bivariate_likelihood.py @@ -23,4 +23,5 @@ xdata, ydata = np.random.multivariate_normal(mean, cov, 100000).T fit = Fit(pdf, x=xdata, y=ydata, objective=LogLikelihood) -fit_result = fit.execute() \ No newline at end of file +fit_result = fit.execute() +print(fit_result) \ No newline at end of file From c5a181ef3aec0ab0d046f42941a8c2a532d80b2d Mon Sep 17 00:00:00 2001 From: Martin Roelfs Date: Wed, 6 Mar 2019 18:30:10 +0100 Subject: [PATCH 4/5] Slightly updated the example --- docs/examples/ex_bivariate_likelihood.rst | 6 +++++- examples/bivariate_likelihood.py | 6 +++--- 2 files changed, 8 insertions(+), 4 deletions(-) diff --git a/docs/examples/ex_bivariate_likelihood.rst b/docs/examples/ex_bivariate_likelihood.rst index f65b7b8f..d62cad8b 100644 --- a/docs/examples/ex_bivariate_likelihood.rst +++ b/docs/examples/ex_bivariate_likelihood.rst @@ -2,10 +2,14 @@ Example: Likelihood fitting a Bivariate Gaussian ================================================ In this example, we shall perform likelihood fitting to a `bivariate normal -distrubtion `_, +distribution `_, 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 diff --git a/examples/bivariate_likelihood.py b/examples/bivariate_likelihood.py index 52e14515..55d18fdf 100644 --- a/examples/bivariate_likelihood.py +++ b/examples/bivariate_likelihood.py @@ -16,9 +16,9 @@ # 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]]) +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 From 5a72a468ceffc8ecadaabee67be929ef5dee5d5f Mon Sep 17 00:00:00 2001 From: Martin Roelfs Date: Fri, 8 Mar 2019 10:58:59 +0100 Subject: [PATCH 5/5] Added likelihood example to the toctree --- docs/examples/index.rst | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/docs/examples/index.rst b/docs/examples/index.rst index 4bd68384..e2f6df44 100644 --- a/docs/examples/index.rst +++ b/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 ------------------------