Repository: StingraySoftware/stingray
Pull request #246 (2017-08-15 21:46:54)

Added Laplacian log-likelihood
I added the Laplace log-likelihood that's used for Cospectra, but I'm not 100% sure it's correct.
The Laplace likelihood is simple. The log-likelihood should just be

$P(x_i | m_i) = -log(2*\sigma_i) - |x_i - m_i|/\sigma_i$ 

for a data point $x_i$ and a parametric model value at that position $m_i$, and a data uncertainty $\sigma_i$. 

Unless I've missed something.

@matteobachetti Are there any (astronomy) references I can look up?


Patch:
diff --git a/stingray/modeling/posterior.py b/stingray/modeling/posterior.py
index bc979fd0..ea0b1570 100644
--- a/stingray/modeling/posterior.py
+++ b/stingray/modeling/posterior.py
@@ -21,8 +21,8 @@ from stingray import Lightcurve, Powerspectrum
 #from stingray.modeling.parametricmodels import logmin
 
 __all__ = ["set_logprior", "Posterior", "PSDPosterior", "LogLikelihood",
-           "PSDLogLikelihood", "GaussianLogLikelihood",
-           "PoissonPosterior", "GaussianPosterior",
+           "PSDLogLikelihood", "GaussianLogLikelihood", "LaplaceLogLikelihood",
+           "PoissonPosterior", "GaussianPosterior", "LaplacePosterior",
            "PriorUndefinedError", "LikelihoodUndefinedError"]
 
 logmin = -10000000000000000.0
@@ -351,6 +351,59 @@ class PSDLogLikelihood(LogLikelihood):
         else:
             return loglike
 
+class LaplaceLogLikelihood(LogLikelihood):
+    def __init__(self, x, y, yerr, model):
+        """
+        A Gaussian likelihood.
+
+        Parameters
+        ----------
+        x : iterable
+            Array with independent variable
+
+        y : iterable
+            Array with dependent variable
+
+        model : an Astropy Model instance
+            The model to use in the likelihood.
+
+        yerr : iterable
+            Array with the uncertainties on `y`, in standard deviation
+
+        """
+
+        LogLikelihood.__init__(self, x, y, model)
+        self.yerr = yerr
+
+        self.npar = 0
+        for pname in self.model.param_names:
+            if not self.model.fixed[pname] and not self.model.tied[pname]:
+                self.npar += 1
+
+    def evaluate(self, pars, neg=False):
+
+        if np.size(pars) != self.npar:
+            raise IncorrectParameterError("Input parameters must" +
+                                          " match model parameters!")
+
+        _fitter_to_model_params(self.model, pars)
+
+        mean_model = self.model(self.x)
+
+        with warnings.catch_warnings(record=True) as out:
+
+                        loglike = np.sum(-np.log(2.*self.yerr) - \
+                                  (np.abs(self.y - mean_model)/self.yerr))
+
+        if not np.isfinite(loglike):
+            loglike = logmin
+
+        if neg:
+            return -loglike
+        else:
+            return loglike
+
+
 class Posterior(object):
 
     def __init__(self, x, y, model, **kwargs):
@@ -597,3 +650,38 @@ class GaussianPosterior(Posterior):
 
         if not priors is None:
             self.logprior = set_logprior(self, priors)
+
+class LaplacePosterior(Posterior):
+
+    def __init__(self, x, y, yerr, model, priors=None):
+        """
+        A general class for two-dimensional data following a Gaussian
+        sampling distribution.
+
+        Parameters
+        ----------
+        x: numpy.ndarray
+            independent variable
+
+        y: numpy.ndarray
+            dependent variable
+
+        yerr: numpy.ndarray
+            measurement uncertainties for y, in standard deviation
+
+        model: instance of any subclass of parameterclass.ParametricModel
+            The model for the power spectrum. Note that in order to define
+            the posterior properly, the ParametricModel subclass must be
+            instantiated with the hyperpars parameter set, or there won't
+            be a prior to be calculated! If all this object is used
+            for a maximum likelihood-style analysis, no prior is required.
+
+        """
+        self.loglikelihood = LaplaceLogLikelihood(x, y, yerr, model)
+
+        Posterior.__init__(self, x, y, model)
+
+        self.yerr = yerr
+
+        if not priors is None:
+            self.logprior = set_logprior(self, priors)


Test patch:
diff --git a/stingray/modeling/tests/test_posterior.py b/stingray/modeling/tests/test_posterior.py
index 336262d7..30c704d9 100644
--- a/stingray/modeling/tests/test_posterior.py
+++ b/stingray/modeling/tests/test_posterior.py
@@ -8,7 +8,9 @@ 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 Posterior, PSDPosterior, \
+    PoissonPosterior, GaussianPosterior, LaplacePosterior
+
 from stingray.modeling import set_logprior
 from stingray.modeling.posterior import logmin
 from stingray.modeling.posterior import IncorrectParameterError
@@ -494,6 +496,136 @@ class TestGaussianPosterior(object):
 
         assert np.isclose(lpost(t0), logmin, 1e-5)
 
+class TestLaplacePosterior(object):
+
+    @classmethod
+    def setup_class(cls):
+
+        nx = 1000000
+        cls.x = np.arange(nx)
+        cls.countrate = 10.0
+        cls.cerr = 2.0
+        cls.y = np.random.normal(cls.countrate, cls.cerr, size=cls.x.shape[0])
+        cls.yerr = np.ones_like(cls.y)*cls.cerr
+
+        cls.model = models.Const1D()
+
+        p_amplitude = lambda amplitude: \
+            scipy.stats.norm(loc=cls.countrate, scale=cls.cerr).pdf(amplitude)
+
+        cls.priors = {"amplitude":p_amplitude}
+
+    def test_logprior_fails_without_prior(self):
+        lpost = LaplacePosterior(self.x, self.y, self.yerr, self.model)
+
+        with pytest.raises(AttributeError):
+            lpost.logprior([10])
+
+    def test_making_posterior(self):
+        lpost = LaplacePosterior(self.x, self.y, self.yerr, self.model)
+        lpost.logprior = set_logprior(lpost, self.priors)
+
+        assert lpost.x.all() == self.x.all()
+        assert lpost.y.all() == self.y.all()
+
+    def test_correct_number_of_parameters(self):
+        lpost = LaplacePosterior(self.x, self.y, self.yerr, self.model)
+        lpost.logprior = set_logprior(lpost, self.priors)
+
+        with pytest.raises(IncorrectParameterError):
+            lpost([2,3])
+
+    def test_logprior(self):
+        t0 = [10.0]
+
+        lpost = LaplacePosterior(self.x, self.y, self.yerr, self.model)
+        lpost.logprior = set_logprior(lpost, self.priors)
+
+        lp_test = lpost.logprior(t0)
+        lp = np.log(scipy.stats.norm(self.countrate, self.cerr).pdf(t0))
+        assert lp == lp_test
+
+    def test_loglikelihood(self):
+        t0 = [10.0]
+        self.model.amplitude = t0[0]
+        mean_model = self.model(self.x)
+
+        loglike = np.sum(-np.log(2.0*self.yerr) -
+                         np.abs(self.y - mean_model)/self.yerr)
+
+        lpost = LaplacePosterior(self.x, self.y, self.yerr, self.model)
+        lpost.logprior = set_logprior(lpost, self.priors)
+
+        loglike_test = lpost.loglikelihood(t0, neg=False)
+
+        assert np.isclose(loglike, loglike_test)
+
+
+    def test_negative_loglikelihood(self):
+        t0 = [10.0]
+        self.model.amplitude = t0[0]
+        mean_model = self.model(self.x)
+
+        loglike = -np.sum(-np.log(2.0*self.yerr) -
+                         np.abs(self.y - mean_model)/self.yerr)
+
+        lpost = LaplacePosterior(self.x, self.y, self.yerr, self.model)
+        lpost.logprior = set_logprior(lpost, self.priors)
+
+        loglike_test = lpost.loglikelihood(t0, neg=True)
+
+        assert np.isclose(loglike, loglike_test)
+
+    def test_posterior(self):
+        t0 = [10.0]
+        self.model.amplitude = t0[0]
+        mean_model = self.model(self.x)
+
+        lpost = LaplacePosterior(self.x, self.y, self.yerr, self.model)
+        lpost.logprior = set_logprior(lpost, self.priors)
+
+        post_test = lpost(t0, neg=False)
+
+        loglike = np.sum(-np.log(2.0*self.yerr) -
+                         np.abs(self.y - mean_model)/self.yerr)
+
+        logprior = np.log(scipy.stats.norm(self.countrate, self.cerr).pdf(t0))
+
+        post = loglike + logprior
+
+        assert np.isclose(post_test, post, atol=1.e-10)
+
+    def test_negative_posterior(self):
+        t0 = [10.0]
+        self.model.amplitude = t0[0]
+        mean_model = self.model(self.x)
+
+        lpost = LaplacePosterior(self.x, self.y, self.yerr, self.model)
+        lpost.logprior = set_logprior(lpost, self.priors)
+
+        post_test = lpost(t0, neg=True)
+
+        loglike = np.sum(-np.log(2.0*self.yerr) -
+                         np.abs(self.y - mean_model)/self.yerr)
+
+        logprior = np.log(scipy.stats.norm(self.countrate, self.cerr).pdf(t0))
+
+        post = -loglike - logprior
+
+        assert np.isclose(post_test, post, atol=1.e-10)
+
+    def test_counts_are_nan(self):
+        y = np.nan * np.ones(self.x.shape[0])
+
+        t0 = [10.0]
+        self.model.amplitude = t0[0]
+
+        lpost = LaplacePosterior(self.x, y, self.yerr, self.model)
+        lpost.logprior = set_logprior(lpost, self.priors)
+
+        assert np.isclose(lpost(t0), logmin, 1e-5)
+
+
 
 
 class TestPerPosteriorAveragedPeriodogram(object):