Repository: StingraySoftware/stingray
Pull request #187 (2017-04-28 02:14:34)

Bounds and parameter links in Maximum likelihood fitting
At the moment, the only constraint on parameter evaluation that Stingray supports is the "freeze"; bounds can be accounted only through uniform priors (a little unintuitive) and no parameter linking is possible. However, we use Astropy models that, in principle, have all these functionalities. The problem is that, for fitting, we don't use Astropy but scipy.

Here I add some classes and methods to show one possible way to use the machinery that Astropy provides instead of Scipy's. 

For now, I only produced a Fitter à la Astropy that uses the correct PSDLogLikelihood and the Simplex minimization algorithm, plus a few methods and a test that show how to use it. The only key change is in `parameterestimation.py`. Elsewhere, I just duplicated a lot of code in order to change a few lines when needed without being destructive: if this PR is provisionally accepted, I will work towards harmonizing the code. 

What is there:
A fitting algorithm that correctly treats bounds and ties between parameters, for PSDLogLikelihood

To do:
- [x] Extend to other likelihoods
- [x] Errors on the fit parameters
- [x] Harmonize and merge the existing machinery for MCMC sampling

<!-- Reviewable:start -->
---
This change is [<img src="https://reviewable.io/review_button.svg" height="34" align="absmiddle" alt="Reviewable"/>](https://reviewable.io/reviews/stingraysoftware/stingray/187)
<!-- Reviewable:end -->



Patch:
diff --git a/stingray/modeling/parameterestimation.py b/stingray/modeling/parameterestimation.py
index f8e6271d..f48bc4fc 100644
--- a/stingray/modeling/parameterestimation.py
+++ b/stingray/modeling/parameterestimation.py
@@ -1,3 +1,4 @@
+from __future__ import print_function, division
 
 __all__ = ["OptimizationResults", "ParameterEstimation", "PSDParEst",
            "SamplingResults"]
@@ -38,8 +39,9 @@ try:
 except ImportError:
     comp_hessian = False
 
-from stingray.modeling.posterior import Posterior, PSDPosterior, LogLikelihood
-from astropy.modeling.fitting import _fitter_to_model_params
+from stingray.modeling.posterior import Posterior, LogLikelihood
+from astropy.modeling.fitting import _fitter_to_model_params, \
+    _model_to_fit_params, _validate_model, _convert_input
 
 
 class OptimizationResults(object):
@@ -132,11 +134,27 @@ class OptimizationResults(object):
         print("The best-fit model parameters plus errors are:")
 
         fixed = [lpost.model.fixed[n] for n in lpost.model.param_names]
-        parnames = [n for n, f in zip(lpost.model.param_names, fixed) \
+        tied = [lpost.model.tied[n] for n in lpost.model.param_names]
+        bounds = [lpost.model.bounds[n] for n in lpost.model.param_names]
+
+        parnames = [n for n, f in zip(lpost.model.param_names,
+                                      np.logical_or(fixed, tied)) \
                     if f is False]
 
-        for i, (x, y, p) in enumerate(zip(self.p_opt, self.err, parnames)):
-            print("%i) Parameter %s: %f.5 +/- %f.5f"%(i, p, x, y))
+        all_parnames = [n for n in lpost.model.param_names]
+        for i, par in enumerate(all_parnames):
+            print("{:3}) Parameter {:<20}: ".format(i, par), end="")
+
+            if par in parnames:
+                idx = parnames.index(par)
+                print("{:<20.5f} +/- {:<20.5f} ".format(self.p_opt[idx],
+                                                  self.err[idx]), end="")
+                print("[{:>10} {:>10}]".format(str(bounds[i][0]),
+                                               str(bounds[i][1])))
+            elif fixed[i]:
+                print("{:<20.5f} (Fixed) ".format(lpost.model.parameters[i]))
+            elif tied[i]:
+                print("{:<20.5f} (Tied) ".format(lpost.model.parameters[i]))
 
         print("\n")
 
@@ -233,7 +251,12 @@ class ParameterEstimation(object):
             raise TypeError("lpost must be a subclass of "
                             "Posterior or LogLikelihoood.")
 
-        if not len(t0) == lpost.npar:
+        newmod = lpost.model.copy()
+        newmod.parameters = t0
+        p0, _ = _model_to_fit_params(newmod)
+        # p0 will be shorter than t0, if there are any frozen/tied parameters
+        # this has to match with the npar attribute.
+        if not len(p0) == lpost.npar:
             raise ValueError("Parameter set t0 must be of right "
                              "length for model in lpost.")
 
@@ -249,19 +272,25 @@ class ParameterEstimation(object):
         # at least until scipy 0.11 is out
         funcval = 100.0
         i = 0
+
         while funcval == 100 or funcval == 200 or \
                 funcval == 0.0 or not np.isfinite(funcval):
 
             if i > 20:
                 raise Exception("Fitting unsuccessful!")
             # perturb parameters slightly
-            t0_p = np.random.multivariate_normal(t0, np.diag(np.abs(t0)/100.))
+            t0_p = np.random.multivariate_normal(p0, np.diag(np.abs(p0)/100.))
 
+            # print(lpost.model, dir(lpost.model), lpost.model.parameter_constraints, lpost.model.param_names)
+            params = [getattr(newmod,name) for name in newmod.param_names]
+            bounds = [p.bounds for p in params if not np.any([p.tied, p.fixed])]
+            # print(params, bounds)
             # if max_post is True, do the Maximum-A-Posteriori Fit
             if self.max_post:
                 opt = scipy.optimize.minimize(lpost, t0_p,
                                               method=self.fitmethod,
                                               args=args, tol=1.e-10,
+                                              bounds=bounds,
                                               **scipy_optimize_options)
 
             # if max_post is False, then do a Maximum Likelihood Fit
@@ -271,7 +300,8 @@ class ParameterEstimation(object):
                     opt = scipy.optimize.minimize(lpost.loglikelihood, t0_p,
                                                   method=self.fitmethod,
                                                   args=args, tol=1.e-10,
-                                                 **scipy_optimize_options)
+                                                  bounds=bounds,
+                                                  **scipy_optimize_options)
 
                 elif isinstance(lpost, LogLikelihood):
                     # Except this could be a `LogLikelihood object
@@ -280,10 +310,8 @@ class ParameterEstimation(object):
                     opt = scipy.optimize.minimize(lpost.evaluate, t0_p,
                                                   method=self.fitmethod,
                                                   args=args, tol=1.e-10,
+                                                  bounds=bounds,
                                                   **scipy_optimize_options)
-                else:
-                    raise TypeError("lpost must be a Posterior or LogLikelihood"
-                                    " object.")
 
             funcval = opt.fun
             i += 1
diff --git a/stingray/modeling/posterior.py b/stingray/modeling/posterior.py
index 01d5cc8a..81553bf7 100644
--- a/stingray/modeling/posterior.py
+++ b/stingray/modeling/posterior.py
@@ -3,6 +3,7 @@ from __future__ import division
 import abc
 
 import numpy as np
+import six
 from scipy.special import gamma as scipy_gamma
 from scipy.special import gammaln as scipy_gammaln
 from astropy.modeling.fitting import _fitter_to_model_params
@@ -139,8 +140,8 @@ def set_logprior(lpost, priors):
     return logprior
 
 
+@six.add_metaclass(abc.ABCMeta)
 class LogLikelihood(object):
-    __metaclass__ = abc.ABCMeta
 
     def __init__(self, x, y, model):
         """
@@ -301,7 +302,7 @@ class PSDLogLikelihood(LogLikelihood):
         self.m = m
         self.npar = 0
         for pname in self.model.param_names:
-            if not self.model.fixed[pname]:
+            if not self.model.fixed[pname] and not self.model.tied[pname]:
                 self.npar += 1
 
     def evaluate(self, pars, neg=False):
@@ -333,7 +334,6 @@ class PSDLogLikelihood(LogLikelihood):
         else:
             return loglike
 
-
 class Posterior(object):
 
     def __init__(self, x, y, model):
diff --git a/stingray/modeling/scripts.py b/stingray/modeling/scripts.py
index 1659dd39..bc66c16e 100644
--- a/stingray/modeling/scripts.py
+++ b/stingray/modeling/scripts.py
@@ -1,4 +1,4 @@
-
+from __future__ import division
 import numpy as np
 from astropy.modeling import models
 


Test patch:
diff --git a/stingray/modeling/tests/test_parameterestimation.py b/stingray/modeling/tests/test_parameterestimation.py
index 7717ac52..468666bb 100644
--- a/stingray/modeling/tests/test_parameterestimation.py
+++ b/stingray/modeling/tests/test_parameterestimation.py
@@ -1,6 +1,8 @@
+from __future__ import division, print_function
 import numpy as np
 import scipy.stats
 import os
+import re
 
 from astropy.tests.helper import pytest
 from astropy.modeling import models
@@ -572,6 +574,7 @@ class TestPSDParEst(object):
         cls.model.fwhm_0 = cls.fwhm_0
         cls.model.amplitude_0 = cls.amplitude_0
         cls.model.amplitude_1 = cls.amplitude_1
+        cls.model.x_0_0.fixed = True
 
         p = cls.model(freq)
 
@@ -614,6 +617,54 @@ class TestPSDParEst(object):
         cls.t0 = [cls.x_0_0, cls.fwhm_0, cls.amplitude_0, cls.amplitude_1]
         cls.neg = True
 
+    def test_fitting_with_ties_and_bounds(self, capsys):
+        double_f = lambda model : model.x_0_0 * 2
+        model = self.model.copy()
+        model =  self.model + models.Lorentz1D(amplitude=model.amplitude_0,
+                                   x_0 = model.x_0_0 * 2,
+                                   fwhm = model.fwhm_0)
+        model.x_0_0 = self.model.x_0_0
+        model.amplitude_0 = self.model.amplitude_0
+        model.amplitude_1 = self.model.amplitude_1
+        model.fwhm_0 = self.model.fwhm_0
+        model.x_0_2.tied = double_f
+        model.fwhm_0.bounds = [0, 10]
+        model.amplitude_0.fixed = True
+
+        p = model(self.ps.freq)
+
+        noise = np.random.exponential(size=len(p))
+        power = noise*p
+
+        ps = Powerspectrum()
+        ps.freq = self.ps.freq
+        ps.power = power
+        ps.m = self.ps.m
+        ps.df = self.ps.df
+        ps.norm = "leahy"
+
+        pe = PSDParEst(ps)
+        llike = PSDLogLikelihood(ps.freq, ps.power, model)
+
+        true_pars = [self.amplitude_0, self.x_0_0, self.fwhm_0,
+                     self.amplitude_1,
+                     model.amplitude_2.value, model.x_0_2.value,
+                     model.fwhm_2.value]
+        res = pe.fit(llike, true_pars)
+
+        res.print_summary(llike)
+        out, err = capsys.readouterr()
+        assert "100.00000            (Fixed)" in out
+        pattern = \
+            re.compile(r"5\) Parameter x_0_2\s+: [0-9]\.[0-9]{5}\s+\(Tied\)")
+        assert pattern.search(out)
+
+        compare_pars = [self.x_0_0, self.fwhm_0,
+                        self.amplitude_1,
+                        model.amplitude_2.value,
+                        model.fwhm_2.value]
+
+        assert np.all(np.isclose(compare_pars, res.p_opt, rtol=0.5))
 
     def test_par_est_initializes(self):
         pe = PSDParEst(self.ps)
@@ -664,7 +715,7 @@ class TestPSDParEst(object):
         assert isinstance(res, OptimizationResults), "res must be of type " \
                                                      "OptimizationResults"
 
-    def test_plotfits(self):
+    def test_plotfits_leahy(self):
         pe = PSDParEst(self.ps)
         t0 = [2.0, 1, 1, 1]
         lpost = PSDPosterior(self.ps, self.model, self.priors)
@@ -676,7 +727,7 @@ class TestPSDParEst(object):
         assert os.path.exists("test_ps_fit.png")
         os.unlink("test_ps_fit.png")
 
-    def test_plotfits_log(self):
+    def test_plotfits_log_leahy(self):
         pe = PSDParEst(self.ps)
         t0 = [2.0, 1, 1, 1]
         lpost = PSDPosterior(self.ps, self.model, self.priors)
@@ -688,7 +739,7 @@ class TestPSDParEst(object):
         assert os.path.exists("test_ps_fit.png")
         os.unlink("test_ps_fit.png")
 
-    def test_plotfits2(self):
+    def test_plotfits_rms(self):
         t0 = [2.0, 1, 1, 1]
         ps = Powerspectrum()
         ps.freq = self.ps.freq
@@ -706,7 +757,7 @@ class TestPSDParEst(object):
         assert os.path.exists("test_ps_fit.png")
         os.unlink("test_ps_fit.png")
 
-    def test_plotfits_log2(self):
+    def test_plotfits_log_rms(self):
         ps = Powerspectrum()
         ps.freq = self.ps.freq
         ps.power = self.ps.power
@@ -724,6 +775,42 @@ class TestPSDParEst(object):
         assert os.path.exists("test_ps_fit.png")
         os.unlink("test_ps_fit.png")
 
+    def test_plotfits_pow(self):
+        t0 = [2.0, 1, 1, 1]
+        ps = Powerspectrum()
+        ps.freq = self.ps.freq
+        ps.power = self.ps.power
+        ps.m = self.ps.m
+        ps.df = self.ps.df
+        ps.norm = "none"
+        pe = PSDParEst(ps)
+        lpost = PSDPosterior(self.ps, self.model, self.priors)
+
+        res = pe.fit(self.lpost, t0)
+
+        pe.plotfits(res, res2=res, save_plot=True)
+
+        assert os.path.exists("test_ps_fit.png")
+        os.unlink("test_ps_fit.png")
+
+    def test_plotfits_log_pow(self):
+        ps = Powerspectrum()
+        ps.freq = self.ps.freq
+        ps.power = self.ps.power
+        ps.m = self.ps.m
+        ps.df = self.ps.df
+        ps.norm = "none"
+        pe = PSDParEst(ps)
+        lpost = PSDPosterior(self.ps, self.model, self.priors)
+
+        t0 = [2.0, 1, 1, 1]
+        res = pe.fit(self.lpost, t0)
+
+        pe.plotfits(res, res2=res, save_plot=True, log=True)
+
+        assert os.path.exists("test_ps_fit.png")
+        os.unlink("test_ps_fit.png")
+
     def test_compute_lrt_fails_when_garbage_goes_in(self):
         pe = PSDParEst(self.ps)
         t0 = [2.0, 1, 1, 1]
@@ -766,3 +853,4 @@ class TestPSDParEst(object):
         os.unlink("test_corner.pdf")
         assert sample_res.acceptance > 0.25
         assert isinstance(sample_res, SamplingResults)
+
diff --git a/stingray/modeling/tests/test_posterior.py b/stingray/modeling/tests/test_posterior.py
index 8ce2830b..21572d9c 100644
--- a/stingray/modeling/tests/test_posterior.py
+++ b/stingray/modeling/tests/test_posterior.py
@@ -1,3 +1,4 @@
+from __future__ import division, print_function
 import numpy as np
 import scipy.stats
 import copy
@@ -11,9 +12,33 @@ from stingray.modeling import Posterior, PSDPosterior, PoissonPosterior, Gaussia
 from stingray.modeling import set_logprior
 from stingray.modeling.posterior import logmin
 from stingray.modeling.posterior import IncorrectParameterError
+from stingray.modeling.posterior import LogLikelihood
 
 np.random.seed(20150907)
 
+class TestMeta(object):
+    def test_use_loglikelihood_class_directly(self):
+        with pytest.raises(TypeError):
+            a = LogLikelihood(1, 2, models.Lorentz1D)
+
+    def test_inherit_loglikelihood_improperly(self):
+        class a(LogLikelihood):
+            def __init__(self, *args, **kwargs):
+                LogLikelihood.__init__(self, *args, **kwargs)
+
+        with pytest.raises(TypeError):
+            a(1, 2, models.Lorentz1D)
+
+    def test_inherit_loglikelihood_properly(self):
+        class a(LogLikelihood):
+            def __init__(self, *args, **kwargs):
+                LogLikelihood.__init__(self, *args, **kwargs)
+            def evaluate(self, parameters):
+                pass
+
+        a(1, 2, models.Lorentz1D)
+
+
 class TestSetPrior(object):
 
     @classmethod
diff --git a/stingray/modeling/tests/test_scripts.py b/stingray/modeling/tests/test_scripts.py
index 66e389d6..9904066e 100644
--- a/stingray/modeling/tests/test_scripts.py
+++ b/stingray/modeling/tests/test_scripts.py
@@ -1,19 +1,12 @@
+from __future__ import division, print_function
 import numpy as np
-import scipy.stats
-import copy
 
-from astropy.tests.helper import pytest
 from astropy.modeling import models
-from scipy.special import gammaln as scipy_gammaln
 
-from stingray import Lightcurve, Powerspectrum
-from stingray.modeling import Posterior, PSDPosterior, PoissonPosterior, GaussianPosterior
-from stingray.modeling import set_logprior
-from stingray.modeling.posterior import logmin
-from stingray.modeling.posterior import IncorrectParameterError
-
-from stingray.modeling.scripts import fit_lorentzians
+from stingray import Powerspectrum
 
+from stingray.modeling.scripts import fit_lorentzians#, fit_powerspectrum_bounds
+from stingray.modeling.scripts import fit_powerspectrum
 
 class TestFitLorentzians(object):
 
@@ -53,7 +46,8 @@ class TestFitLorentzians(object):
         cls.ps.df = cls.ps.freq[1] - cls.ps.freq[0]
         cls.ps.m = 1
 
-        cls.t0 = [200.0, 0.5, 0.1, 100.0, 2.0, 1.0, 50.0, 7.5, 0.5, 2.0]
+        cls.t0 = np.asarray([200.0, 0.5, 0.1, 100.0, 2.0, 1.0,
+                             50.0, 7.5, 0.5, 2.0])
 
         cls.parest, cls.res = fit_lorentzians(cls.ps, cls.nlor, cls.t0)
 
@@ -73,7 +67,43 @@ class TestFitLorentzians(object):
                      self.amplitude_2, self.x_0_2, self.fwhm_2,
                      self.whitenoise]
 
-        for t, p in zip(true_pars, self.res.p_opt):
-            print(str(t) + "\t" + str(p))
+        assert np.all(np.isclose(true_pars, self.res.p_opt, rtol=0.5))
+
+    def test_fitting_with_tied_pars(self):
+        double_f = lambda model : model.x_0_0 * 4
+        triple_f = lambda model : model.x_0_0 * 15
+        model = self.model.copy()
+        model.x_0_1.tied = double_f
+        model.x_0_2.tied = triple_f
+        model.amplitude_0 = self.t0[0]
+        # model.bounds = {}
+
+        parest, res = fit_powerspectrum(self.ps, model,
+                                np.random.normal(self.t0,
+                                                 self.t0 / 10))
+
+        true_pars = [self.amplitude_0,
+                     self.x_0_0, self.fwhm_0,
+                     self.amplitude_1, self.fwhm_1,
+                     self.amplitude_2, self.fwhm_2,
+                     self.whitenoise]
+
+        assert np.all(np.isclose(true_pars, res.p_opt, rtol=0.5))
+
+
+    def test_fitting_with_fixed_pars(self):
+        model = self.model.copy()
+        model.amplitude_0 = self.t0[0]
+        model.amplitude_0.fixed = True
+        # model.bounds = {}
+
+        parest, res = fit_powerspectrum(self.ps, model,
+                                  np.random.normal(self.t0,
+                                                   self.t0 / 10))
+
+        true_pars = [self.x_0_0, self.fwhm_0,
+                     self.amplitude_1, self.x_0_1, self.fwhm_1,
+                     self.amplitude_2, self.x_0_2, self.fwhm_2,
+                     self.whitenoise]
 
-        assert np.all(np.isclose(true_pars, self.res.p_opt, rtol=0.5))
\ No newline at end of file
+        assert np.all(np.isclose(true_pars, res.p_opt, rtol=0.5))