Repository: gammapy/gammapy
Pull request #5135 (2024-02-29 16:24:19)

Add different options to compute stat_array on FluxPointsDatasets
Add different options to compute stat_array on FluxPointsDatasets :
- chi2 : etimate from chi2 statistics.
- profile : estimate from interpolation of the likelihood profile.
- distrib : estimate from probability distributions,
             assumes that flux points correspond to asymmetric gaussians and upper limits complemantary error functions.

 Default is `chi2`, in that case upper limits are ignored and the mean of asymetrics error is used.
The `distrib` case provides an approximation if the `profile` is not available
 which allows to take into accounts upper limits and asymetrics errors.

Updated /tutorials/analysis-1d/spectral_analysis.py to present the different cases:
![image](https://github.com/gammapy/gammapy/assets/25845654/3d997a81-d4ca-40f4-a708-af6d66bb8cd1)
Without ul the results are the same:
![image](https://github.com/gammapy/gammapy/assets/25845654/f3ec38f1-d893-4605-8baf-9a4e1b1c7533)


codecov[bot]: ## [Codecov](https://app.codecov.io/gh/gammapy/gammapy/pull/5135?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy) Report
Attention: Patch coverage is `43.37349%` with `47 lines` in your changes are missing coverage. Please review.
> Project coverage is 75.15%. Comparing base [(`f462141`)](https://app.codecov.io/gh/gammapy/gammapy/commit/f4621417c1dc0edb67110a94a1fe99dbbda45fe8?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy) to head [(`522482a`)](https://app.codecov.io/gh/gammapy/gammapy/pull/5135?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy).
> Report is 16 commits behind head on main.

| [Files](https://app.codecov.io/gh/gammapy/gammapy/pull/5135?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy) | Patch % | Lines |
|---|---|---|
| [gammapy/datasets/flux\_points.py](https://app.codecov.io/gh/gammapy/gammapy/pull/5135?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy#diff-Z2FtbWFweS9kYXRhc2V0cy9mbHV4X3BvaW50cy5weQ==) | 41.55% | [45 Missing :warning: ](https://app.codecov.io/gh/gammapy/gammapy/pull/5135?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy) |
| [gammapy/utils/interpolation.py](https://app.codecov.io/gh/gammapy/gammapy/pull/5135?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy#diff-Z2FtbWFweS91dGlscy9pbnRlcnBvbGF0aW9uLnB5) | 66.66% | [2 Missing :warning: ](https://app.codecov.io/gh/gammapy/gammapy/pull/5135?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy) |

<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff             @@
##             main    #5135      +/-   ##
==========================================
- Coverage   75.23%   75.15%   -0.08%     
==========================================
  Files         234      234              
  Lines       34845    34920      +75     
==========================================
+ Hits        26216    26245      +29     
- Misses       8629     8675      +46     
```



</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/gammapy/gammapy/pull/5135?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy).

QRemy: Thanks @AtreyeeS I implemented most of the comments.

> I don't understand the implementation of the `distrib` option. Maybe you can add some information in `gammapy.stat` docs page?

Assuming gaussian errors the likelihood is given by the probability density function of the normal distribution.
For the upper limit case it is necessary to marginalize over the unknown measurement so we integrate the normal distribution up to the upper limit value which gives the complementary error function.

![Screenshot 2024-03-06 at 17 03 23](https://github.com/gammapy/gammapy/assets/25845654/aefae4dd-8587-4b7e-aeae-251624d09945)
see eq. C7 from https://iopscience.iop.org/article/10.1088/0004-637X/773/2/168/pdf
(in eq C6 they are doing something more complicated as they also marginalised over the unknown true value)
adonath: Thanks @QRemy, the statistical approach to handle the ULs looks reasonable to me. I was wonderign whether we can maybe find a more descriptive name. I think even "chi2-marginalize-ul" would be better. We should definitely give the paper reference in addition.  

Patch:
diff --git a/examples/tutorials/analysis-1d/spectral_analysis.py b/examples/tutorials/analysis-1d/spectral_analysis.py
index 0022b2677..5ef9f5a64 100644
--- a/examples/tutorials/analysis-1d/spectral_analysis.py
+++ b/examples/tutorials/analysis-1d/spectral_analysis.py
@@ -380,23 +380,13 @@ plt.show()
 #
 # To round up our analysis we can compute flux points by fitting the norm
 # of the global model in energy bands.
-# We can utilise the `~gammapy.estimators.utils.resample_energy_edges`
-# for defining the energy bins in which the minimum number of `sqrt_ts` is 2.
-# To do so we first stack the individual datasets, only for obtaining the energies:
-#
-
-dataset_stacked = Datasets(datasets).stack_reduce()
-energy_edges = resample_energy_edges(dataset_stacked, conditions={"sqrt_ts_min": 2})
-
-
-######################################################################
-# Now we create an instance of the
+# We create an instance of the
 # `~gammapy.estimators.FluxPointsEstimator`, by passing the dataset and
 # the energy binning:
 #
 
 fpe = FluxPointsEstimator(
-    energy_edges=energy_edges, source="crab", selection_optional="all"
+    energy_edges=energy_axis.edges, source="crab", selection_optional="all"
 )
 flux_points = fpe.run(datasets=datasets)
 
@@ -416,6 +406,7 @@ display(flux_points.to_table(sed_type="dnde", formatted=True))
 fig, ax = plt.subplots()
 flux_points.plot(ax=ax, sed_type="e2dnde", color="darkorange")
 flux_points.plot_ts_profiles(ax=ax, sed_type="e2dnde")
+ax.set_xlim(0.6, 40)
 plt.show()
 
 ######################################################################
@@ -430,8 +421,11 @@ plt.show()
 # quickly made like this:
 #
 
-flux_points_dataset = FluxPointsDataset(data=flux_points, models=model_best_joint)
-flux_points_dataset.plot_fit()
+flux_points_dataset = FluxPointsDataset(
+    data=flux_points, models=model_best_joint.copy()
+)
+ax, _ = flux_points_dataset.plot_fit()
+ax.set_xlim(0.6, 40)
 plt.show()
 
 
@@ -503,6 +497,82 @@ plt.show()
 
 # sphinx_gallery_thumbnail_number = 5
 
+######################################################################
+# A note on statistics
+# --------------------
+#
+# Different statistic are available for the FluxPointDataset :
+# * chi2 : estimate from chi2 statistics.
+# * profile : estimate from interpolation of the likelihood profile.
+# * distrib : estimate from probability distributions,
+#             assuming that flux points correspond to asymmetric gaussians
+#             and upper limits complementary error functions.
+# Default is `chi2`, in that case upper limits are ignored and the mean of asymetrics error is used.
+# So it is recommended to use `profile` if `stat_scan` is available on flux points.
+# The `distrib` case provides an approximation if the `profile` is not available
+# which allows to take into accounts upper limit and asymetrics error.
+#
+# In the example below we can see that the `profile` case matches exactly the result
+# from the joint analysis of the ON/OFF datasets using `wstat` (as labelled).
+
+
+def plot_stat(fp_dataset):
+    fig, ax = plt.subplots()
+
+    plot_kwargs = {
+        "energy_bounds": [0.1, 30] * u.TeV,
+        "sed_type": "e2dnde",
+        "ax": ax,
+    }
+
+    fp_dataset.data.plot(energy_power=2, ax=ax)
+    model_best_joint.spectral_model.plot(
+        color="b", lw=0.5, **plot_kwargs, label="wstat"
+    )
+
+    stat_types = ["chi2", "profile", "distrib"]
+    colors = ["red", "g", "c"]
+    lss = ["--", ":", "--"]
+
+    for ks, stat in enumerate(stat_types):
+
+        fp_dataset.stat_type = stat
+
+        fit = Fit()
+        fit.run([fp_dataset])
+
+        fp_dataset.models[0].spectral_model.plot(
+            color=colors[ks], ls=lss[ks], **plot_kwargs, label=stat
+        )
+        fp_dataset.models[0].spectral_model.plot_error(
+            facecolor=colors[ks], **plot_kwargs
+        )
+        plt.legend()
+
+
+plot_stat(flux_points_dataset)
+
+######################################################################
+
+# In order to avoid discrepancies due to the treatment of upper limits
+# we can utilise the `~gammapy.estimators.utils.resample_energy_edges`
+# for defining energy bins in which the minimum number of `sqrt_ts` is 2.
+# In that case all the statistics definitions give equivalent results.
+#
+
+energy_edges = resample_energy_edges(dataset_stacked, conditions={"sqrt_ts_min": 2})
+
+fpe_no_ul = FluxPointsEstimator(
+    energy_edges=energy_edges, source="crab", selection_optional="all"
+)
+flux_points_no_ul = fpe_no_ul.run(datasets=datasets)
+flux_points_dataset_no_ul = FluxPointsDataset(
+    data=flux_points_no_ul,
+    models=model_best_joint.copy(),
+)
+
+plot_stat(flux_points_dataset_no_ul)
+
 ######################################################################
 # Exercises
 # ---------
diff --git a/gammapy/datasets/flux_points.py b/gammapy/datasets/flux_points.py
index 121777302..2ff79b2f2 100644
--- a/gammapy/datasets/flux_points.py
+++ b/gammapy/datasets/flux_points.py
@@ -1,6 +1,8 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 import logging
 import numpy as np
+from scipy.special import erfc
+from scipy.stats import norm
 from astropy import units as u
 from astropy.io import fits
 from astropy.table import Table
@@ -14,6 +16,7 @@ from gammapy.modeling.models import (
     SkyModel,
     TemplateSpatialModel,
 )
+from gammapy.utils.interpolation import interpolate_profile
 from gammapy.utils.scripts import make_name, make_path
 from .core import Dataset
 
@@ -55,6 +58,25 @@ class FluxPointsDataset(Dataset):
     meta_table : `~astropy.table.Table`
         Table listing information on observations used to create the dataset.
         One line per observation for stacked datasets.
+    stat_type : str
+        Method used to compute the statistics:
+        * chi2 : estimate from chi2 statistics.
+        * profile : estimate from interpolation of the likelihood profile.
+        * distrib : Assuming gaussian errors the likelihood is given by the
+                    probability density function of the normal distribution.
+                    For the upper limit case it is necessary to marginalize over the unknown measurement,
+                    So we integrate the normal distribution up to the upper limit value
+                    which gives the complementary error function.
+                    See eq. C7 of Mohanty et al (2013) :
+                    https://iopscience.iop.org/article/10.1088/0004-637X/773/2/168/pdf
+
+        Default is `chi2`, in that case upper limits are ignored and the mean of asymetrics error is used.
+        However it is recommended to use `profile` if `stat_scan` is available on flux points.
+        The `distrib` case provides an approximation if the profile is not available.
+    stat_kwargs : dict
+        Extra arguments specifying the interpolation scheme of the likelihood profile.
+        Used only if `stat_type=="profile"`. In that case the default is :
+        `stat_kwargs={"interp_scale":"sqrt", "extrapolate":True}
 
     Examples
     --------
@@ -102,7 +124,6 @@ class FluxPointsDataset(Dataset):
     ``gammapy download datasets --tests --out $GAMMAPY_DATA``
     """
 
-    stat_type = "chi2"
     tag = "FluxPointsDataset"
 
     def __init__(
@@ -113,6 +134,8 @@ class FluxPointsDataset(Dataset):
         mask_safe=None,
         name=None,
         meta_table=None,
+        stat_type="chi2",
+        stat_kwargs=None,
     ):
         if not data.geom.has_energy_axis:
             raise ValueError("FluxPointsDataset needs an energy axis")
@@ -122,11 +145,62 @@ class FluxPointsDataset(Dataset):
         self.models = models
         self.meta_table = meta_table
 
-        if mask_safe is None:
-            mask_safe = (~data.is_ul).data
+        self._available_stat_type = dict(
+            chi2=self._stat_array_chi2,
+            profile=self._stat_array_profile,
+            distrib=self._stat_array_distrib,
+        )
+
+        if stat_kwargs is None:
+            stat_kwargs = dict()
+        self.stat_kwargs = stat_kwargs
+
+        self.stat_type = stat_type
 
+        if mask_safe is None:
+            mask_safe = np.ones(self.data.dnde.data.shape, dtype=bool)
         self.mask_safe = mask_safe
 
+    @property
+    def available_stat_type(self):
+        return list(self._available_stat_type.keys())
+
+    @property
+    def stat_type(self):
+        return self._stat_type
+
+    @stat_type.setter
+    def stat_type(self, stat_type):
+
+        if stat_type not in self.available_stat_type:
+            raise ValueError(
+                f"Invalid stat_type: possible options are {self.available_stat_type}"
+            )
+
+        if stat_type == "chi2":
+            self._mask_valid = (~self.data.is_ul).data & np.isfinite(self.data.dnde)
+        elif stat_type == "distrib":
+            self._mask_valid = (
+                self.data.is_ul.data & np.isfinite(self.data.dnde_ul)
+            ) | np.isfinite(self.data.dnde)
+        elif stat_type == "profile":
+            self.stat_kwargs.setdefault("interp_scale", "sqrt")
+            self.stat_kwargs.setdefault("extrapolate", True)
+            self._profile_interpolators = self._get_valid_profile_interpolators()
+        self._stat_type = stat_type
+
+    @property
+    def mask_valid(self):
+        return self._mask_valid
+
+    @property
+    def mask_safe(self):
+        return self._mask_safe & self.mask_valid
+
+    @mask_safe.setter
+    def mask_safe(self, mask_safe):
+        self._mask_safe = mask_safe
+
     @property
     def name(self):
         return self._name
@@ -310,13 +384,81 @@ class FluxPointsDataset(Dataset):
 
     def stat_array(self):
         """Fit statistic array."""
+        return self._available_stat_type[self.stat_type]()
+
+    def _stat_array_chi2(self):
+        """Chi2 statistics."""
         model = self.flux_pred()
         data = self.data.dnde.quantity
         try:
-            sigma = self.data.dnde_err
+            sigma = self.data.dnde_err.quantity
         except AttributeError:
-            sigma = (self.data.dnde_errn + self.data.dnde_errp) / 2
-        return ((data - model) / sigma.quantity).to_value("") ** 2
+            sigma = (self.data.dnde_errn + self.data.dnde_errp).quantity / 2
+        return ((data - model) / sigma).to_value("") ** 2
+
+    def _stat_array_profile(self):
+        """Estimate statitistic from interpolation of the likelihood profile."""
+        model = np.zeros(self.data.dnde.data.shape) + (
+            self.flux_pred() / self.data.dnde_ref
+        ).to_value("")
+        stat = np.zeros(model.shape)
+        for idx in np.ndindex(self._profile_interpolators.shape):
+            stat[idx] = self._profile_interpolators[idx](model[idx])
+        return stat
+
+    def _get_valid_profile_interpolators(self):
+        value_scan = self.data.stat_scan.geom.axes["norm"].center
+        shape_axes = self.data.stat_scan.geom._shape[slice(3, None)][::-1]
+        interpolators = np.empty(shape_axes, dtype=object)
+        self._mask_valid = np.ones(self.data.dnde.data.shape, dtype=bool)
+        for idx in np.ndindex(shape_axes):
+            stat_scan = np.abs(
+                self.data.stat_scan.data[idx].squeeze()
+                - self.data.stat.data[idx].squeeze()
+            )
+            self._mask_valid[idx] = np.all(np.isfinite(stat_scan))
+            interpolators[idx] = interpolate_profile(
+                value_scan,
+                stat_scan,
+                interp_scale=self.stat_kwargs["interp_scale"],
+                extrapolate=self.stat_kwargs["extrapolate"],
+            )
+        return interpolators
+
+    def _stat_array_distrib(self):
+        """Estimate statistic from probability distributions,
+        assumes that flux points correspond to asymmetric gaussians
+        and upper limits complementary error functions.
+        """
+        model = np.zeros(self.data.dnde.data.shape) + self.flux_pred().to_value(
+            self.data.dnde.unit
+        )
+
+        stat = np.zeros(model.shape)
+
+        mask_valid = ~np.isnan(self.data.dnde.data)
+        loc = self.data.dnde.data[mask_valid]
+        value = model[mask_valid]
+        try:
+            mask_p = (model >= self.data.dnde.data)[mask_valid]
+            scale = np.zeros(mask_p.shape)
+            scale[mask_p] = self.data.dnde_errp.data[mask_valid][mask_p]
+            scale[~mask_p] = self.data.dnde_errn.data[mask_valid][~mask_p]
+        except AttributeError:
+            scale = self.data.dnde_err.data[mask_valid]
+        stat[mask_valid] = -2 * np.log(
+            norm.pdf(value, loc=loc, scale=scale) / norm.pdf(loc, loc=loc, scale=scale)
+        )
+
+        mask_ul = self.data.is_ul.data & ~np.isnan(self.data.dnde_ul.data)
+        value = model[mask_ul]
+        loc_ul = self.data.dnde_ul.data[mask_ul]
+        scale_ul = self.data.dnde_ul.data[mask_ul]
+        stat[mask_ul] = 2 * np.log(
+            (erfc((loc_ul - value) / scale_ul) / 2)
+            / (erfc((loc_ul - 0) / scale_ul) / 2)
+        )
+        return stat
 
     def residuals(self, method="diff"):
         """Compute flux point residuals.
diff --git a/gammapy/utils/interpolation.py b/gammapy/utils/interpolation.py
index 3d8968079..7c48c52fb 100644
--- a/gammapy/utils/interpolation.py
+++ b/gammapy/utils/interpolation.py
@@ -240,7 +240,7 @@ class LinearScale(InterpolationScale):
         return values
 
 
-def interpolate_profile(x, y, interp_scale="sqrt"):
+def interpolate_profile(x, y, interp_scale="sqrt", extrapolate=False):
     """Helper function to interpolate one-dimensional profiles.
 
     Parameters
@@ -254,6 +254,9 @@ def interpolate_profile(x, y, interp_scale="sqrt"):
         of parabolic shape, a "sqrt" scaling is recommended. In other cases or
         for fine sampled profiles a "lin" can also be used.
         Default is "sqrt".
+    extrapolate : bool
+        Extrapolate or not if the evaluation value is outside the range of x values.
+        Default is False.
 
     Returns
     -------
@@ -261,4 +264,8 @@ def interpolate_profile(x, y, interp_scale="sqrt"):
         Interpolator.
     """
     method_dict = {"sqrt": "quadratic", "lin": "linear"}
-    return scipy.interpolate.interp1d(x, y, kind=method_dict[interp_scale])
+    kwargs = dict(kind=method_dict[interp_scale])
+    if extrapolate:
+        kwargs["bounds_error"] = False
+        kwargs["fill_value"] = "extrapolate"
+    return scipy.interpolate.interp1d(x, y, **kwargs)


Test patch:
diff --git a/gammapy/datasets/tests/test_flux_points.py b/gammapy/datasets/tests/test_flux_points.py
index ad1333bf9..9d163b196 100644
--- a/gammapy/datasets/tests/test_flux_points.py
+++ b/gammapy/datasets/tests/test_flux_points.py
@@ -102,6 +102,19 @@ def test_flux_point_dataset_flux_pred(dataset):
     assert_allclose(dataset.flux_pred()[0].value, 0.000472, rtol=1e-3)
 
 
+@requires_data()
+def test_flux_point_dataset_stat(dataset):
+    dataset.stat_type = "chi2"
+    fit = Fit()
+    fit.run([dataset])
+    assert_allclose(dataset.stat_sum(), 25.205933, rtol=1e-3)
+
+    dataset.stat_type = "distrib"
+    fit = Fit()
+    fit.run([dataset])
+    assert_allclose(dataset.stat_sum(), 36.153428, rtol=1e-3)
+
+
 def test_flux_point_dataset_with_time_axis(tmp_path):
     meta = dict(TIMESYS="utc", SED_TYPE="flux")
 
@@ -158,6 +171,9 @@ def test_flux_point_dataset_with_time_axis(tmp_path):
     with pytest.raises(ValueError):
         flux_points_dataset.plot_residuals()
 
+    flux_points_dataset.stat_type = "distrib"
+    assert_allclose(flux_points_dataset.stat_sum(), 193.8093, rtol=1e-3)
+
 
 @requires_data()
 class TestFluxPointFit:
diff --git a/gammapy/estimators/points/tests/test_lightcurve.py b/gammapy/estimators/points/tests/test_lightcurve.py
index 716ce5f7b..f97ed28d2 100644
--- a/gammapy/estimators/points/tests/test_lightcurve.py
+++ b/gammapy/estimators/points/tests/test_lightcurve.py
@@ -7,7 +7,7 @@ from astropy.table import Column, Table
 from astropy.time import Time
 from astropy.timeseries import BinnedTimeSeries, BoxLeastSquares
 from gammapy.data import GTI
-from gammapy.datasets import Datasets
+from gammapy.datasets import Datasets, FluxPointsDataset
 from gammapy.datasets.actors import DatasetsActor
 from gammapy.estimators import FluxPoints, LightCurveEstimator
 from gammapy.estimators.points.lightcurve import parallel as parallel
@@ -738,6 +738,43 @@ def test_recompute_ul():
         lightcurve.recompute_ul(n_sigma_ul=4)
 
 
+@requires_data()
+def test_dataset_stat_type():
+    datasets = get_spectrum_datasets()
+    selection = ["all"]
+    estimator = LightCurveEstimator(
+        energy_edges=[1, 3, 30] * u.TeV, selection_optional=selection, n_sigma_ul=2
+    )
+    lightcurve = estimator.run(datasets)
+    assert_allclose(
+        lightcurve.dnde_ul.data[0], [[[3.260703e-13]], [[1.159354e-14]]], rtol=1e-3
+    )
+
+    lightcurve_dataset = FluxPointsDataset(
+        data=lightcurve, models=[lightcurve.reference_model]
+    )
+    lightcurve_dataset.stat_type = "profile"
+    assert_allclose(lightcurve_dataset.stat_sum(), 5.741539, rtol=1e-5)
+    lightcurve_dataset.stat_type = "distrib"
+    assert_allclose(lightcurve_dataset.stat_sum(), 5.824218, rtol=1e-5)
+    lightcurve_dataset.stat_type = "chi2"
+    assert_allclose(lightcurve_dataset.stat_sum(), 6.014128, rtol=1e-5)
+
+    # test if scan is not present
+    selection = []
+    estimator = LightCurveEstimator(
+        energy_edges=[1, 30] * u.TeV, selection_optional=selection
+    )
+    lightcurve = estimator.run(datasets)
+
+    lightcurve_dataset = FluxPointsDataset(
+        data=lightcurve, models=[lightcurve.reference_model]
+    )
+    with pytest.raises(AttributeError):
+        lightcurve_dataset.stat_type = "profile"
+        lightcurve_dataset.stat_sum()
+
+
 @requires_data()
 def test_lightcurve_parallel_multiprocessing():
     datasets = get_spectrum_datasets()
diff --git a/gammapy/estimators/points/tests/test_sed.py b/gammapy/estimators/points/tests/test_sed.py
index 4a8b23d10..056669323 100644
--- a/gammapy/estimators/points/tests/test_sed.py
+++ b/gammapy/estimators/points/tests/test_sed.py
@@ -309,6 +309,15 @@ def test_run_pwl(fpe_pwl, tmpdir):
     fp_new = FluxPoints.read(tmpdir / "test.fits")
     assert fp_new.meta["sed_type_init"] == "likelihood"
 
+    # test datasets stat
+    fp_dataset = FluxPointsDataset(data=fp, models=fp.reference_model)
+    fp_dataset.stat_type = "chi2"
+    assert_allclose(fp_dataset.stat_sum(), 3.82374, rtol=1e-4)
+    fp_dataset.stat_type = "profile"
+    assert_allclose(fp_dataset.stat_sum(), 3.790053, rtol=1e-4)
+    fp_dataset.stat_type = "distrib"
+    assert_allclose(fp_dataset.stat_sum(), 3.783325, rtol=1e-4)
+
 
 def test_run_ecpl(fpe_ecpl, tmpdir):
     datasets, fpe = fpe_ecpl