Repository: gammapy/gammapy
Pull request #5696 (2025-02-06 22:11:06)

Introducing FitStatistic
<!-- These are hidden commments. You can leave them or delete them. -->

<!-- Don't be shy! Github issues or pull requests are welcome any time. -->
<!-- But if you first want to chat, use the Gammapy Slack. -->
<!-- See https://gammapy.org/contact.html -->

**Description**
<!-- What is the motivation for this pull request? -->
<!-- Briefly: what changes are done here? -->
<!-- If this is to address a Github issue, mention its number and use the specific keywords to create a link -->
<!-- See docs.github.com/en/issues/tracking-your-work-with-issues/linking-a-pull-request-to-an-issue -->
<!-- Example: This PR is to resolve #42 -->

This pull request proposes another/simpler version of FitStatistic as discussed with @adonath @QRemy @AtreyeeS @bkhelifi .
All `stat_sum` and `stat_array` are now only on the base `Dataset` class.

**Dear reviewer**
<!-- Let the reviewer and Gammapy team know what you want: -->
<!-- * Is this ready for review? Or is it work in progress and you want some feedback? -->
<!-- * Do you want to go through review here? Or if someone should just finish this up and merge it? -->
<!-- Do you have any specific questions, e.g. about API or implementation? -->
<!-- Do you include a test executing the new code you're adding (to make sure it runs)? -->
<!-- Do you include some documentation? Is it needed? -->


codecov[bot]: ## [Codecov](https://app.codecov.io/gh/gammapy/gammapy/pull/5696?dropdown=coverage&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 `94.00000%` with `6 lines` in your changes missing coverage. Please review.
> Project coverage is 94.60%. Comparing base [(`e23a50a`)](https://app.codecov.io/gh/gammapy/gammapy/commit/e23a50a31d901cf828580c5f9444a06bf9a691e4?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy) to head [(`9003d96`)](https://app.codecov.io/gh/gammapy/gammapy/commit/9003d96b51b771b608f4c970d8358de16d77eac8?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy).
> Report is 5 commits behind head on main.

| [Files with missing lines](https://app.codecov.io/gh/gammapy/gammapy/pull/5696?dropdown=coverage&src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy) | Patch % | Lines |
|---|---|---|
| [gammapy/stats/fit\_statistics.py](https://app.codecov.io/gh/gammapy/gammapy/pull/5696?src=pr&el=tree&filepath=gammapy%2Fstats%2Ffit_statistics.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy#diff-Z2FtbWFweS9zdGF0cy9maXRfc3RhdGlzdGljcy5weQ==) | 93.33% | [6 Missing :warning: ](https://app.codecov.io/gh/gammapy/gammapy/pull/5696?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    #5696      +/-   ##
==========================================
- Coverage   94.60%   94.60%   -0.01%     
==========================================
  Files         238      238              
  Lines       37279    37319      +40     
==========================================
+ Hits        35269    35305      +36     
- Misses       2010     2014       +4     
```

</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/gammapy/gammapy/pull/5696?dropdown=coverage&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).


Patch:
diff --git a/gammapy/datasets/__init__.py b/gammapy/datasets/__init__.py
index 9d65a3279..27182a3f0 100644
--- a/gammapy/datasets/__init__.py
+++ b/gammapy/datasets/__init__.py
@@ -6,7 +6,6 @@ from .io import OGIPDatasetReader, OGIPDatasetWriter, FermipyDatasetsReader
 from .map import (
     MapDataset,
     MapDatasetOnOff,
-    MapDatasetWeighted,
     create_empty_map_dataset_from_irfs,
     create_map_dataset_from_observation,
     create_map_dataset_geoms,
@@ -20,7 +19,6 @@ DATASET_REGISTRY = Registry(
     [
         MapDataset,
         MapDatasetOnOff,
-        MapDatasetWeighted,
         SpectrumDataset,
         SpectrumDatasetOnOff,
         FluxPointsDataset,
@@ -41,7 +39,6 @@ __all__ = [
     "MapDataset",
     "MapDatasetEventSampler",
     "MapDatasetOnOff",
-    "MapDatasetWeighted",
     "ObservationEventSampler",
     "OGIPDatasetWriter",
     "OGIPDatasetReader",
diff --git a/gammapy/datasets/core.py b/gammapy/datasets/core.py
index b80a2607d..98253c8bf 100644
--- a/gammapy/datasets/core.py
+++ b/gammapy/datasets/core.py
@@ -10,6 +10,7 @@ from astropy.table import Table, vstack
 from gammapy.data import GTI
 from gammapy.modeling.models import DatasetModels, Models
 from gammapy.utils.scripts import make_name, make_path, read_yaml, to_yaml, write_yaml
+from gammapy.stats import FIT_STATISTICS_REGISTRY
 
 log = logging.getLogger(__name__)
 
@@ -36,6 +37,17 @@ class Dataset(abc.ABC):
         "diff/sqrt(model)": "(data - model) / sqrt(model)",
     }
 
+    @property
+    def stat_type(self):
+        """The Fit Statistic class used."""
+        return self._stat_type
+
+    @stat_type.setter
+    def stat_type(self, stat_type):
+        """Set the Fit Statistic."""
+        self._fit_statistic = FIT_STATISTICS_REGISTRY[stat_type]
+        self._stat_type = stat_type
+
     def _repr_html_(self):
         try:
             return self.to_html()
@@ -72,22 +84,16 @@ class Dataset(abc.ABC):
         prior_stat_sum = 0.0
         if self.models is not None:
             prior_stat_sum = self.models.parameters.prior_stat_sum()
-        return self._stat_sum_likelihood() + prior_stat_sum
+
+        return self._fit_statistic.stat_sum_dataset(self) + prior_stat_sum
 
     def _stat_sum_likelihood(self):
         """Total statistic given the current model parameters without the priors."""
-        stat = self.stat_array()
+        return self._fit_statistic.stat_sum_dataset(self)
 
-        if self.mask is not None:
-            if isinstance(self.mask, np.ndarray):
-                stat = stat[self.mask.astype(bool)]
-            else:
-                stat = stat[self.mask.data.astype(bool)]
-        return np.sum(stat, dtype=np.float64)
-
-    @abc.abstractmethod
     def stat_array(self):
         """Statistic array, one value per data point."""
+        return self._fit_statistic.stat_array_dataset(self)
 
     def copy(self, name=None):
         """A deep copy.
diff --git a/gammapy/datasets/flux_points.py b/gammapy/datasets/flux_points.py
index 8e22911c8..62cb133e6 100644
--- a/gammapy/datasets/flux_points.py
+++ b/gammapy/datasets/flux_points.py
@@ -1,7 +1,6 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 import logging
 import numpy as np
-from scipy.special import erfc
 from astropy import units as u
 from astropy.io import fits
 from astropy.table import Table
@@ -15,6 +14,7 @@ from gammapy.modeling.models import (
     SkyModel,
     TemplateSpatialModel,
 )
+from gammapy.stats import FIT_STATISTICS_REGISTRY
 from gammapy.utils.interpolation import interpolate_profile
 from gammapy.utils.scripts import make_name, make_path
 from .core import Dataset
@@ -141,12 +141,6 @@ class FluxPointsDataset(Dataset):
         self.models = models
         self.meta_table = meta_table
 
-        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
@@ -159,7 +153,7 @@ class FluxPointsDataset(Dataset):
 
     @property
     def available_stat_type(self):
-        return list(self._available_stat_type.keys())
+        return ["chi2", "distrib", "profile"]
 
     @property
     def stat_type(self):
@@ -183,6 +177,7 @@ class FluxPointsDataset(Dataset):
             self.stat_kwargs.setdefault("extrapolate", True)
             self._profile_interpolators = self._get_valid_profile_interpolators()
         self._stat_type = stat_type
+        self._fit_statistic = FIT_STATISTICS_REGISTRY[stat_type]
 
     @property
     def mask_valid(self):
@@ -377,30 +372,6 @@ class FluxPointsDataset(Dataset):
             flux += flux_model
         return flux
 
-    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.quantity
-        except AttributeError:
-            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]
@@ -420,46 +391,6 @@ class FluxPointsDataset(Dataset):
             )
         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]
-
-            mask_invalid = np.isnan(scale)
-            scale[mask_invalid] = self.data.dnde_err.data[mask_valid][mask_invalid]
-        except AttributeError:
-            scale = self.data.dnde_err.data[mask_valid]
-
-        stat[mask_valid] = ((value - loc) / scale) ** 2
-
-        mask_ul = self.data.is_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)
-        )
-
-        stat[np.isnan(stat.data)] = 0
-        return stat
-
     def residuals(self, method="diff"):
         """Compute flux point residuals.
 
diff --git a/gammapy/datasets/map.py b/gammapy/datasets/map.py
index 49211fbd3..3b8772421 100644
--- a/gammapy/datasets/map.py
+++ b/gammapy/datasets/map.py
@@ -15,11 +15,7 @@ from gammapy.modeling.models import DatasetModels, FoVBackgroundModel, Models
 from gammapy.stats import (
     CashCountsStatistic,
     WStatCountsStatistic,
-    cash,
-    cash_sum_cython,
-    weighted_cash_sum_cython,
     get_wstat_mu_bkg,
-    wstat,
 )
 from gammapy.utils.fits import HDULocation, LazyFitsData
 from gammapy.utils.random import get_random_state
@@ -33,7 +29,6 @@ from .utils import get_axes
 __all__ = [
     "MapDataset",
     "MapDatasetOnOff",
-    "MapDatasetWeighted",
     "create_empty_map_dataset_from_irfs",
     "create_map_dataset_geoms",
     "create_map_dataset_from_observation",
@@ -491,7 +486,6 @@ class MapDataset(Dataset):
     MapDatasetOnOff, SpectrumDataset, FluxPointsDataset.
     """
 
-    stat_type = "cash"
     tag = "MapDataset"
     counts = LazyFitsData(cache=True)
     exposure = LazyFitsData(cache=True)
@@ -528,6 +522,7 @@ class MapDataset(Dataset):
         meta_table=None,
         name=None,
         meta=None,
+        stat_type="cash",
     ):
         self._name = make_name(name)
         self._evaluators = {}
@@ -559,6 +554,8 @@ class MapDataset(Dataset):
         self.meta_table = meta_table
         self.meta = meta
 
+        self.stat_type = stat_type
+
     @property
     def _psf_kernel(self):
         """Precompute PSFkernel if there is only one spatial bin in the PSFmap"""
@@ -1206,10 +1203,6 @@ class MapDataset(Dataset):
         if self.meta and other.meta:
             self.meta.stack(other.meta)
 
-    def stat_array(self):
-        """Statistic function value per bin given the current model parameters."""
-        return cash(n_on=self.counts.data, mu_on=self.npred().data)
-
     def residuals(self, method="diff", **kwargs):
         """Compute residuals map.
 
@@ -1466,32 +1459,6 @@ class MapDataset(Dataset):
 
         return ax_spatial, ax_spectral
 
-    def stat_sum(self):
-        """Total statistic function value given the current model parameters and priors."""
-        prior_stat_sum = 0.0
-        if self.models is not None:
-            prior_stat_sum = self.models.parameters.prior_stat_sum()
-
-        counts, npred = self.counts.data.astype(float), self.npred().data
-
-        if self.mask is not None:
-            mask = ~(self.mask.data == False)  # noqa
-            counts = counts[mask]
-            npred = npred[mask]
-            if self.mask.data.dtype == bool or self.stat_type == "cash":
-                cash_sum = cash_sum_cython(counts, npred)
-            elif self.stat_type == "cash_weighted":
-                weight = self.mask.data[mask]
-                cash_sum = weighted_cash_sum_cython(counts, npred, weight)
-            else:
-                raise ValueError(
-                    f"'stat_type' must be a 'cash' or `cash_weighted`."
-                    f", got `{self.stat_type}` instead."
-                )
-        else:
-            cash_sum = cash_sum_cython(counts.ravel(), npred.ravel())
-        return cash_sum + prior_stat_sum
-
     def _to_asimov_dataset(self):
         """Create Asimov dataset from the current models."""
 
@@ -2387,11 +2354,6 @@ class MapDataset(Dataset):
         plot_mask(ax=axes[3], mask=self.mask_safe_image, hatches=["///"], colors="w")
 
 
-class MapDatasetWeighted(MapDataset):
-    stat_type = "cash_weighted"
-    tag = "MapDatasetWeighted"
-
-
 class MapDatasetOnOff(MapDataset):
     """Map dataset for on-off likelihood fitting.
 
@@ -2442,7 +2404,6 @@ class MapDatasetOnOff(MapDataset):
 
     """
 
-    stat_type = "wstat"
     tag = "MapDatasetOnOff"
 
     def __init__(
@@ -2461,6 +2422,7 @@ class MapDatasetOnOff(MapDataset):
         gti=None,
         meta_table=None,
         meta=None,
+        stat_type="wstat",
     ):
         self._name = make_name(name)
         self._evaluators = {}
@@ -2482,6 +2444,8 @@ class MapDatasetOnOff(MapDataset):
         else:
             self._meta = meta
 
+        self.stat_type = stat_type
+
     def __str__(self):
         str_ = super().__str__()
 
@@ -2586,17 +2550,6 @@ class MapDatasetOnOff(MapDataset):
             return None
         return self.alpha * self.counts_off
 
-    def stat_array(self):
-        """Statistic function value per bin given the current model parameters."""
-        mu_sig = self.npred_signal().data
-        on_stat_ = wstat(
-            n_on=self.counts.data,
-            n_off=self.counts_off.data,
-            alpha=list(self.alpha.data),
-            mu_sig=mu_sig,
-        )
-        return np.nan_to_num(on_stat_)
-
     @property
     def _counts_statistic(self):
         """Counts statistics of the dataset."""
@@ -2853,17 +2806,6 @@ class MapDatasetOnOff(MapDataset):
 
         super().stack(other, nan_to_num=nan_to_num)
 
-    def stat_sum(self):
-        """Total statistic function value given the current model parameters.
-
-        If the off counts are passed as None and the elements of the safe mask are False, zero will be returned.
-        Otherwise, the stat sum will be calculated and returned.
-        """
-        if self.counts_off is None and not np.any(self.mask_safe.data):
-            return 0
-        else:
-            return Dataset.stat_sum(self)
-
     def fake(self, npred_background, random_state="random-seed"):
         """Simulate fake counts (on and off) for the current model and reduced IRFs.
 
diff --git a/gammapy/datasets/spectrum.py b/gammapy/datasets/spectrum.py
index 06bc7638b..a4f218d5e 100644
--- a/gammapy/datasets/spectrum.py
+++ b/gammapy/datasets/spectrum.py
@@ -265,7 +265,6 @@ class SpectrumDataset(PlotMixin, MapDataset):
     For more information see :ref:`datasets`.
     """
 
-    stat_type = "cash"
     tag = "SpectrumDataset"
 
     def cutout(self, *args, **kwargs):
@@ -292,7 +291,6 @@ class SpectrumDatasetOnOff(PlotMixin, MapDatasetOnOff):
     For more information see :ref:`datasets`.
     """
 
-    stat_type = "wstat"
     tag = "SpectrumDatasetOnOff"
 
     def cutout(self, *args, **kwargs):
diff --git a/gammapy/makers/background/fov.py b/gammapy/makers/background/fov.py
index d18343afd..70ca1bd40 100644
--- a/gammapy/makers/background/fov.py
+++ b/gammapy/makers/background/fov.py
@@ -1,5 +1,6 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 """FoV background estimation."""
+
 import logging
 import numpy as np
 from gammapy.maps import Map, RegionGeom
@@ -132,7 +133,7 @@ class FoVBackgroundMaker(Maker):
         if self.exclusion_mask:
             mask = self.exclusion_mask.interp_to_geom(geom=geom)
         else:
-            mask = Map.from_geom(geom=geom, data=1, dtype=bool)
+            mask = Map.from_geom(geom=geom, data=True, dtype=bool)
         return mask
 
     def _make_masked_summed_counts(self, dataset):
diff --git a/gammapy/stats/__init__.py b/gammapy/stats/__init__.py
index a3a699122..a1d7b072d 100644
--- a/gammapy/stats/__init__.py
+++ b/gammapy/stats/__init__.py
@@ -2,7 +2,19 @@
 """Statistics."""
 
 from .counts_statistic import CashCountsStatistic, WStatCountsStatistic
-from .fit_statistics import cash, cstat, get_wstat_gof_terms, get_wstat_mu_bkg, wstat
+from .fit_statistics import (
+    cash,
+    cstat,
+    get_wstat_gof_terms,
+    get_wstat_mu_bkg,
+    wstat,
+    Chi2FitStatistic,
+    CashFitStatistic,
+    Chi2AsymmetricErrorFitStatistic,
+    ProfileFitStatistic,
+    WStatFitStatistic,
+    WeightedCashFitStatistic,
+)
 from .fit_statistics_cython import (
     weighted_cash_sum_cython,
     cash_sum_cython,
@@ -19,6 +31,15 @@ from .variability import (
     structure_function,
 )
 
+FIT_STATISTICS_REGISTRY = {
+    "cash": CashFitStatistic,
+    "wstat": WStatFitStatistic,
+    "chi2": Chi2FitStatistic,
+    "distrib": Chi2AsymmetricErrorFitStatistic,
+    "profile": ProfileFitStatistic,
+    "cash_weighted": WeightedCashFitStatistic,
+}
+
 __all__ = [
     "cash",
     "cash_sum_cython",
diff --git a/gammapy/stats/fit_statistics.py b/gammapy/stats/fit_statistics.py
index 962af8922..0dd310220 100644
--- a/gammapy/stats/fit_statistics.py
+++ b/gammapy/stats/fit_statistics.py
@@ -4,10 +4,28 @@
 see :ref:`fit-statistics`
 """
 
+from abc import ABC
 import numpy as np
-from gammapy.stats.fit_statistics_cython import TRUNCATION_VALUE
-
-__all__ = ["cash", "cstat", "wstat", "get_wstat_mu_bkg", "get_wstat_gof_terms"]
+from scipy.special import erfc
+from gammapy.maps import Map
+from gammapy.stats.fit_statistics_cython import (
+    TRUNCATION_VALUE,
+    cash_sum_cython,
+    weighted_cash_sum_cython,
+)
+
+
+__all__ = [
+    "cash",
+    "cstat",
+    "wstat",
+    "get_wstat_mu_bkg",
+    "get_wstat_gof_terms",
+    "CashFitStatistic",
+    "WStatFitStatistic",
+    "Chi2FitStatistic",
+    "Chi2AsymmetricErrorFitStatistic",
+]
 
 
 def cash(n_on, mu_on, truncation_value=TRUNCATION_VALUE):
@@ -228,3 +246,186 @@ def get_wstat_gof_terms(n_on, n_off):
     term += np.where(n_off == 0, 0, term2)
 
     return 2 * term
+
+
+class FitStatistic(ABC):
+    """Abstract base class for FitStatistic objects."""
+
+    @classmethod
+    def stat_sum_dataset(cls, dataset):
+        """Calculate -2 * sum log(L)."""
+        stat_array = cls.stat_array_dataset(dataset)
+        if dataset.mask is not None:
+            mask = dataset.mask.data if isinstance(dataset.mask, Map) else dataset.mask
+            stat_array = stat_array[mask]
+        return np.sum(stat_array)
+
+    @classmethod
+    def stat_array_dataset(cls, dataset):
+        """Calculate -2 * log(L)."""
+        raise NotImplementedError
+
+    @classmethod
+    def loglikelihood_dataset(cls, dataset):
+        """Calculate sum log(L)."""
+        return -0.5 * cls.stat_sum_dataset(dataset)
+
+
+class CashFitStatistic(FitStatistic):
+    """Cash statistic class for Poisson with known background."""
+
+    @classmethod
+    def stat_sum_dataset(cls, dataset):
+        mask = dataset.mask
+        counts, npred = dataset.counts.data, dataset.npred().data
+
+        if mask is not None:
+            mask = mask.data.astype("bool")
+            counts, npred = counts[mask], npred[mask]
+
+        counts = counts.astype(float)  # This might be done in the Dataset
+        return cash_sum_cython(counts.ravel(), npred.ravel())
+
+    @classmethod
+    def stat_array_dataset(cls, dataset):
+        counts, npred = dataset.counts.data, dataset.npred().data
+        return cash(n_on=counts, mu_on=npred)
+
+
+class WeightedCashFitStatistic(FitStatistic):
+    """Cash statistic class for Poisson with known background applying weights."""
+
+    @classmethod
+    def stat_sum_dataset(cls, dataset):
+        counts, npred = dataset.counts.data.astype(float), dataset.npred().data
+
+        if dataset.mask is not None:
+            mask = ~(dataset.mask.data == False)  # noqa
+            counts = counts[mask]
+            npred = npred[mask]
+
+            weights = dataset.mask.data[mask].astype("float")
+            return weighted_cash_sum_cython(counts, npred, weights)
+        else:
+            # No weights back to regular cash statistic
+            return cash_sum_cython(counts.ravel(), npred.ravel())
+
+    @classmethod
+    def stat_array_dataset(cls, dataset):
+        counts, npred = dataset.counts.data, dataset.npred().data
+        weights = 1.0
+        if dataset.mask is not None:
+            weights = dataset.mask.astype("float")
+        return cash(n_on=counts, mu_on=npred) * weights
+
+
+class WStatFitStatistic(FitStatistic):
+    """WStat fit statistic class for ON-OFF Poisson measurements."""
+
+    @classmethod
+    def stat_array_dataset(cls, dataset):
+        """Statistic function value per bin given the current model parameters."""
+        counts, counts_off, alpha = (
+            dataset.counts.data,
+            dataset.counts_off.data,
+            dataset.alpha.data,
+        )
+        npred_signal = dataset.npred_signal().data
+        on_stat_ = wstat(
+            n_on=counts,
+            n_off=counts_off,
+            alpha=alpha,
+            mu_sig=npred_signal,
+        )
+        return np.nan_to_num(on_stat_)
+
+    @classmethod
+    def stat_sum_dataset(cls, dataset):
+        """Statistic function value per bin given the current model parameters."""
+        if dataset.counts_off is None and not np.any(dataset.mask_safe.data):
+            return 0
+        else:
+            stat_array = cls.stat_array_dataset(dataset)
+            if dataset.mask is not None:
+                stat_array = stat_array[dataset.mask.data]
+            return np.sum(stat_array)
+
+
+class Chi2FitStatistic(FitStatistic):
+    """Chi2 fit statistic class for measurements with gaussian symmetric errors."""
+
+    @classmethod
+    def stat_array_dataset(cls, dataset):
+        """Statistic function value per bin given the current model."""
+        model = dataset.flux_pred()
+        data = dataset.data.dnde.quantity
+        try:
+            sigma = dataset.data.dnde_err.quantity
+        except AttributeError:
+            sigma = (dataset.data.dnde_errn + dataset.data.dnde_errp).quantity / 2
+        return ((data - model) / sigma).to_value("") ** 2
+
+
+class Chi2AsymmetricErrorFitStatistic(FitStatistic):
+    """Pseudo-Chi2 fit statistic class for measurements with gaussian asymmetric errors with upper limits.
+
+    Assumes that regular data follow asymmetric normal pdf and upper limits follow complementary error functions
+    """
+
+    @classmethod
+    def stat_array_dataset(cls, dataset):
+        """Estimate statistic from probability distributions,
+        assumes that flux points correspond to asymmetric gaussians
+        and upper limits complementary error functions.
+        """
+        model = np.zeros(dataset.data.dnde.data.shape) + dataset.flux_pred().to_value(
+            dataset.data.dnde.unit
+        )
+
+        stat = np.zeros(model.shape)
+
+        mask_valid = ~np.isnan(dataset.data.dnde.data)
+        loc = dataset.data.dnde.data[mask_valid]
+        value = model[mask_valid]
+        try:
+            mask_p = (model >= dataset.data.dnde.data)[mask_valid]
+            scale = np.zeros(mask_p.shape)
+            scale[mask_p] = dataset.data.dnde_errp.data[mask_valid][mask_p]
+            scale[~mask_p] = dataset.data.dnde_errn.data[mask_valid][~mask_p]
+
+            mask_invalid = np.isnan(scale)
+            scale[mask_invalid] = dataset.data.dnde_err.data[mask_valid][mask_invalid]
+        except AttributeError:
+            scale = dataset.data.dnde_err.data[mask_valid]
+
+        stat[mask_valid] = ((value - loc) / scale) ** 2
+
+        mask_ul = dataset.data.is_ul.data
+        value = model[mask_ul]
+        loc_ul = dataset.data.dnde_ul.data[mask_ul]
+        scale_ul = dataset.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)
+        )
+
+        stat[np.isnan(stat.data)] = 0
+        return stat
+
+
+class ProfileFitStatistic(FitStatistic):
+    """Pseudo-Chi2 fit statistic class for measurements with gaussian asymmetric errors with upper limits.
+
+    Assumes that regular data follow asymmetric normal pdf and upper limits follow complementary error functions
+    """
+
+    @classmethod
+    def stat_array_dataset(cls, dataset):
+        """Estimate statitistic from interpolation of the likelihood profile."""
+        model = np.zeros(dataset.data.dnde.data.shape) + (
+            dataset.flux_pred() / dataset.data.dnde_ref
+        ).to_value("")
+        stat = np.zeros(model.shape)
+        for idx in np.ndindex(dataset._profile_interpolators.shape):
+            stat[idx] = dataset._profile_interpolators[idx](model[idx])
+        return stat


Test patch:
diff --git a/gammapy/datasets/tests/test_map.py b/gammapy/datasets/tests/test_map.py
index e9c333f0b..69d34f64f 100644
--- a/gammapy/datasets/tests/test_map.py
+++ b/gammapy/datasets/tests/test_map.py
@@ -17,7 +17,6 @@ from gammapy.datasets import (
     Datasets,
     MapDataset,
     MapDatasetOnOff,
-    MapDatasetWeighted,
     create_empty_map_dataset_from_irfs,
     create_map_dataset_from_observation,
 )
@@ -230,7 +229,7 @@ def get_map_dataset(
     models = FoVBackgroundModel(dataset_name=name)
 
     if weighted:
-        return MapDatasetWeighted(
+        return MapDataset(
             models=models,
             exposure=exposure,
             background=background,
@@ -238,6 +237,7 @@ def get_map_dataset(
             edisp=edisp,
             mask_fit=mask_fit,
             name=name,
+            stat_type="cash_weighted",
             **kwargs,
         )
     else:
diff --git a/gammapy/stats/tests/test_fit_statistics.py b/gammapy/stats/tests/test_fit_statistics.py
index da9a9d5ab..a8a4d7877 100644
--- a/gammapy/stats/tests/test_fit_statistics.py
+++ b/gammapy/stats/tests/test_fit_statistics.py
@@ -3,6 +3,13 @@ import pytest
 import numpy as np
 from numpy.testing import assert_allclose
 from gammapy import stats
+from gammapy.stats.fit_statistics import (
+    CashFitStatistic,
+    WeightedCashFitStatistic,
+    WStatFitStatistic,
+    Chi2FitStatistic,
+    Chi2AsymmetricErrorFitStatistic,
+)
 
 
 @pytest.fixture
@@ -169,3 +176,213 @@ def test_wstat_corner_cases():
 
     actual = stats.get_wstat_mu_bkg(n_on=n_on, mu_sig=mu_sig, n_off=n_off, alpha=alpha)
     assert_allclose(actual, 0)
+
+
+class MockDataset:
+    @staticmethod
+    def create_region(size, axis_name="energy"):
+        from gammapy.maps import RegionGeom, MapAxis
+
+        axis = MapAxis.from_nodes(np.arange(size), name=axis_name, unit="TeV")
+        return RegionGeom.create(region=None, axes=[axis])
+
+
+class MockMapDataset(MockDataset):
+    """Mock dataset class"""
+
+    def __init__(self, counts, npred, mask=None):
+        from gammapy.maps import RegionNDMap
+
+        geom = self.create_region(len(counts))
+        self.counts = RegionNDMap(geom=geom, data=np.array(counts))
+        self.npred = lambda: RegionNDMap(
+            geom=geom, data=np.array(npred).astype("float")
+        )
+        self.mask = (
+            RegionNDMap(geom=geom, data=np.array(mask)) if mask is not None else None
+        )
+
+
+class MockMapDatasetOnOff(MockDataset):
+    """Mock dataset for ON-OFF Poisson measurements (WStat)."""
+
+    def __init__(self, counts, npred_signal, counts_off, alpha, mask=None):
+        from gammapy.maps import RegionNDMap
+
+        geom = self.create_region(len(counts))
+        self.counts = RegionNDMap(geom=geom, data=np.array(counts))
+        self.npred_signal = lambda: RegionNDMap(
+            geom=geom, data=np.array(npred_signal).astype("float")
+        )
+        self.mask = (
+            RegionNDMap(geom=geom, data=np.array(mask)) if mask is not None else None
+        )
+        self.counts_off = RegionNDMap(geom=geom, data=np.array(counts_off))
+        self.alpha = RegionNDMap(geom=geom, data=np.array(alpha).astype("float"))
+
+
+class MockFluxPointsDataset(MockDataset):
+    """Mock dataset for flux measurements (chi2, etc)."""
+
+    def __init__(
+        self,
+        norm,
+        norm_err,
+        norm_pred,
+        norm_errn=None,
+        norm_errp=None,
+        norm_ul=None,
+        is_ul=None,
+        model=None,
+        mask=None,
+    ):
+        from gammapy.maps import RegionNDMap
+        from gammapy.estimators import FluxPoints
+        from gammapy.modeling.models import ConstantSpectralModel
+
+        geom = self.create_region(len(norm))
+
+        data = {}
+        data["norm"] = RegionNDMap(geom=geom, data=np.array(norm))
+        data["norm_err"] = RegionNDMap(geom=geom, data=np.array(norm_err))
+        data["norm_errn"] = (
+            RegionNDMap(geom=geom, data=np.array(norm_errn)) if norm_errn else None
+        )
+        data["norm_errp"] = (
+            RegionNDMap(geom=geom, data=np.array(norm_errp)) if norm_errp else None
+        )
+        data["norm_ul"] = (
+            RegionNDMap(geom=geom, data=np.array(norm_ul)) if norm_ul else None
+        )
+        data["is_ul"] = RegionNDMap(geom=geom, data=np.array(is_ul)) if is_ul else None
+
+        model = model if model else ConstantSpectralModel()
+
+        self.data = FluxPoints(data, model)
+        new_geom = geom.copy(axes=geom.axes.rename_axes(["energy"], ["energy_true"]))
+        ref_flux = self.data.reference_model.evaluate_geom(new_geom)
+        self.flux_pred = lambda: np.array(norm_pred).reshape(ref_flux.shape) * ref_flux
+        self.mask = np.array(mask) if mask is not None else None
+
+
+@pytest.fixture
+def mock_map_dataset():
+    return MockMapDataset(counts=[1, 2, 3], npred=[1, 2, 3], mask=[True, False, True])
+
+
+def test_cash_fit_statistic_stat_sum_nomask(mock_map_dataset):
+    """Test the stat_sum_dataset method for CashFitStatistic."""
+    mock_map_dataset.mask = None
+    stat_sum = CashFitStatistic.stat_sum_dataset(mock_map_dataset)
+    assert_allclose(stat_sum, 2.63573737)
+
+
+def test_cash_fit_statistic_with_mask(mock_map_dataset):
+    """Test CashFitStatistic with a mask."""
+    stat_sum = CashFitStatistic.stat_sum_dataset(mock_map_dataset)
+    assert_allclose(stat_sum, 1.40832626799)
+
+
+def test_cash_fit_statistic_loglikelihood(mock_map_dataset):
+    """Test loglikelihood_dataset method."""
+    log_likelihood = CashFitStatistic.loglikelihood_dataset(mock_map_dataset)
+    assert_allclose(log_likelihood, 1.40832626799 * -0.5)
+
+
+def test_cash_fit_statistic_with_non_bool_mask(mock_map_dataset):
+    """Ensure stat_sum_dataset handles non-bool masks gracefully."""
+    mock_map_dataset.mask.data = np.array([1, 0, 2], dtype="float")
+
+    stat_sum = CashFitStatistic.stat_sum_dataset(mock_map_dataset)
+    assert_allclose(stat_sum, 1.40832626799)
+
+
+def test_weightedcash_fit_statistic(mock_map_dataset):
+    """Test WeightedCashFitStatistic."""
+    mock_map_dataset.mask.data = np.array([0.5, 0.5, 0.5], dtype="float")
+
+    stat_sum = WeightedCashFitStatistic.stat_sum_dataset(mock_map_dataset)
+    assert_allclose(stat_sum, 2.63573737 * 0.5)
+
+
+@pytest.fixture()
+def mock_map_dataset_onoff():
+    return MockMapDatasetOnOff(
+        counts=[3, 6, 9],
+        npred_signal=[1, 2, 3],
+        counts_off=[10, 10, 30],
+        alpha=[0.1, 0.2, 0.1],
+        mask=[True, False, True],
+    )
+
+
+def test_wstat_fit_statistic_stat_sum_nomask(mock_map_dataset_onoff):
+    """Test the stat_sum_dataset method for WStatFitStatistic."""
+    mock_map_dataset_onoff.mask = None
+    stat_sum = WStatFitStatistic.stat_sum_dataset(mock_map_dataset_onoff)
+    assert_allclose(stat_sum, 2.4088762929)
+
+
+def test_wstat_fit_statistic_with_mask(mock_map_dataset_onoff):
+    """Test WStatFitStatistic with a mask."""
+    stat_sum = WStatFitStatistic.stat_sum_dataset(mock_map_dataset_onoff)
+    assert_allclose(stat_sum, 1.63526008301)
+
+
+def test_wstat_fit_statistic_with_mask_false(mock_map_dataset_onoff):
+    """Test WStatFitStatistic with mask_safe handling."""
+    mask = [
+        False,
+    ] * 3
+    mock_map_dataset_onoff.mask.data = np.array(mask)
+    stat_sum = WStatFitStatistic.stat_sum_dataset(mock_map_dataset_onoff)
+    assert stat_sum == 0
+
+
+def test_wstat_fit_statistic_loglikelihood(mock_map_dataset_onoff):
+    """Test loglikelihood_dataset method."""
+    log_likelihood = WStatFitStatistic.loglikelihood_dataset(mock_map_dataset_onoff)
+    assert_allclose(log_likelihood, -0.817630041)
+
+
+@pytest.fixture()
+def mock_fp_dataset():
+    return MockFluxPointsDataset(
+        norm=[1.1, 0.9, 1.2, 0.8],
+        norm_err=[0.1, 0.1, 0.2, 0.2],
+        norm_pred=[1, 1, 1, 1],
+        norm_errn=[0.1, 0.1, 0.2, 0.2],
+        norm_errp=[0.1, 0.1, 0.1, 0.1],
+        norm_ul=[1.5, 1.5, 1.5, 1.5],
+        is_ul=[False, False, False, True],
+        mask=[True, True, False, True],
+    )
+
+
+def test_chi2_fit_statistic_stat_sum_nomask(mock_fp_dataset):
+    mock_fp_dataset.mask = None
+    stat_sum = Chi2FitStatistic.stat_sum_dataset(mock_fp_dataset)
+    assert_allclose(stat_sum, 4)
+
+
+def test_chi2_fit_statistic_with_mask(mock_fp_dataset):
+    stat_sum = Chi2FitStatistic.stat_sum_dataset(mock_fp_dataset)
+    assert_allclose(stat_sum, 3)
+
+
+def test_chi2_fit_statistic_with_mask_false(mock_fp_dataset):
+    mask = [False, False, False, False]
+    mock_fp_dataset.mask = np.array(mask)
+    stat_sum = Chi2FitStatistic.stat_sum_dataset(mock_fp_dataset)
+    assert stat_sum == 0
+
+
+def test_chi2_asym_fit_statistic_stat_sum_nomask(mock_fp_dataset):
+    mock_fp_dataset.mask = None
+    stat_sum = Chi2AsymmetricErrorFitStatistic.stat_sum_dataset(mock_fp_dataset)
+    assert_allclose(stat_sum, 5.798344)
+
+
+def test_chi2_asym_fit_statistic_with_mask(mock_fp_dataset):
+    stat_sum = Chi2AsymmetricErrorFitStatistic.stat_sum_dataset(mock_fp_dataset)
+    assert_allclose(stat_sum, 4.798344)