Repository: gammapy/gammapy
Pull request #5289 (2024-06-03 16:29:41)

Add new methods for combine_flux_maps to use likelihood profile or its approximation
Compute combined  flux maps using likelihood profile information.
Requires #5285

codecov[bot]: ## [Codecov](https://app.codecov.io/gh/gammapy/gammapy/pull/5289?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.70588%` with `9 lines` in your changes missing coverage. Please review.
> Project coverage is 94.46%. Comparing base [(`f85fb7b`)](https://app.codecov.io/gh/gammapy/gammapy/commit/f85fb7b40717542f697d1707ce18d0aec211faaf?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy) to head [(`5d8b4a7`)](https://app.codecov.io/gh/gammapy/gammapy/commit/5d8b4a7b3aa4748c5ddeccd38b830209608ebcf7?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy).
> Report is 3 commits behind head on main.

| [Files](https://app.codecov.io/gh/gammapy/gammapy/pull/5289?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/estimators/map/ts.py](https://app.codecov.io/gh/gammapy/gammapy/pull/5289?src=pr&el=tree&filepath=gammapy%2Festimators%2Fmap%2Fts.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy#diff-Z2FtbWFweS9lc3RpbWF0b3JzL21hcC90cy5weQ==) | 88.23% | [8 Missing :warning: ](https://app.codecov.io/gh/gammapy/gammapy/pull/5289?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy) |
| [gammapy/estimators/utils.py](https://app.codecov.io/gh/gammapy/gammapy/pull/5289?src=pr&el=tree&filepath=gammapy%2Festimators%2Futils.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy#diff-Z2FtbWFweS9lc3RpbWF0b3JzL3V0aWxzLnB5) | 98.33% | [1 Missing :warning: ](https://app.codecov.io/gh/gammapy/gammapy/pull/5289?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    #5289      +/-   ##
==========================================
- Coverage   94.46%   94.46%   -0.01%     
==========================================
  Files         232      232              
  Lines       35636    35772     +136     
==========================================
+ Hits        33664    33792     +128     
- Misses       1972     1980       +8     
```



</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/gammapy/gammapy/pull/5289?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/estimators/map/core.py b/gammapy/estimators/map/core.py
index 77445a643..8f187db78 100644
--- a/gammapy/estimators/map/core.py
+++ b/gammapy/estimators/map/core.py
@@ -76,7 +76,6 @@ VALID_QUANTITIES = [
     "npred",
     "npred_excess",
     "stat",
-    "stat_scan",
     "stat_null",
     "stat_scan",
     "dnde_scan_values",
diff --git a/gammapy/estimators/utils.py b/gammapy/estimators/utils.py
index b5d4c0cb2..da76c2d40 100644
--- a/gammapy/estimators/utils.py
+++ b/gammapy/estimators/utils.py
@@ -1,7 +1,8 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 import numpy as np
 import scipy.ndimage
-from scipy import stats
+from scipy import special
+from scipy.interpolate import InterpolatedUnivariateSpline
 from astropy import units as u
 from astropy.coordinates import SkyCoord
 from astropy.table import Table
@@ -757,7 +758,9 @@ def get_combined_significance_maps(estimator, datasets):
     )
 
 
-def combine_flux_maps(maps, method="gaussian_errors", reference_model=None):
+def combine_flux_maps(
+    maps, method="gaussian_errors", reference_model=None, dnde_scan_axis=None
+):
     """Create a FluxMaps by combining a list of flux maps with the same geometry.
      This assumes the flux maps are independent measurements of the same true value.
      The GTI is stacked in the process.
@@ -770,83 +773,352 @@ def combine_flux_maps(maps, method="gaussian_errors", reference_model=None):
         * gaussian_errors :
             Under the gaussian error approximation the likelihood is given by the gaussian distibution.
             The product of gaussians is also a gaussian so can derive dnde, dnde_err, and ts.
+        * distrib :
+            Likelihood profile approximation assuming that probabilities distributions for
+            flux points correspond to asymmetric gaussians and for upper limits to complementary error functions.
+            Use available quantities among dnde, dnde_err, dnde_errp, dnde_errn, dnde_ul, and ts.
+        * profile :
+            Sum the likelihood profile maps.
+            The flux maps must contains the `stat_scan` maps.
+
+        Default is "gaussian_errors" which is the faster but least accurate solution,
+        "distrib"  will be more accurate if dnde_errp and dnde_errn are available,
+        "profile"  will be even more accurate if "stat_scan" is available.
+
     reference_model : `~gammapy.modeling.models.SkyModel`, optional
         Reference model to use for conversions.
-        If None, use the reference_model of the first FluxMaps in the list.
+        Default is None and is will use the reference_model of the first FluxMaps in the list.
+
+    dnde_scan_axis : `~gammapy.maps.MapAxis`
+        Map axis providing the dnde values used to compute the profile.
+        Default is None and it will be derived from the first FluxMaps in the list.
+        Used only if `method` is distrib or profile.
 
     Returns
     -------
     flux_maps : `~gammapy.estimators.FluxMaps`
         Joint flux map.
     """
-    if method != "gaussian_errors":
-        raise ValueError(
-            f'Invalid method : {method}, available methods are : "gaussian_errors"'
-        )
 
-    if isinstance(maps, FluxMaps):
-        geom = maps.dnde.sum_over_axes().geom
-        means = list(maps.dnde.copy().iter_by_image(keepdims=True))
-        sigmas = list(maps.dnde_err.copy().iter_by_image(keepdims=True))
-        mean = Map.from_geom(geom, data=means[0].data, unit=means[0].unit)
-        sigma = Map.from_geom(geom, data=sigmas[0].data, unit=sigmas[0].unit)
-        gti = maps.gti
-        meta = maps.meta
-        if reference_model is None:
-            reference_model = maps.reference_model
+    gtis = [map_.gti for map_ in maps if map_.gti is not None]
+    if np.any(gtis):
+        gti = gtis[0].copy()
+        for k in range(1, len(gtis)):
+            gti.stack(gtis[k])
     else:
+        gti = None
+    # TODO : change this once we have stackable metadata objets
+    metas = [map_.meta for map_ in maps if map_.meta is not None]
+    meta = {}
+    if np.any(metas):
+        for data in metas:
+            meta.update(data)
+    if reference_model is None:
+        reference_model = maps[0].reference_model
+
+    if method == "gaussian_errors":
+
         means = [map_.dnde.copy() for map_ in maps]
         sigmas = [map_.dnde_err.copy() for map_ in maps]
+        # compensate for the ts deviation from gaussian approximation expectation in each map
+        ts_diff = np.nansum(
+            [
+                map_.ts.data - (map_.dnde.data / map_.dnde_err.data) ** 2
+                for map_ in maps
+            ],
+            axis=0,
+        )
+
         mean = means[0]
         sigma = sigmas[0]
 
-        gtis = [map_.gti for map_ in maps if map_.gti is not None]
-        if np.any(gtis):
-            gti = gtis[0].copy()
-            for k in range(1, len(gtis)):
-                gti.stack(gtis[k])
-        else:
-            gti = None
+        for k in range(1, len(means)):
+
+            mean_k = means[k].quantity.to_value(mean.unit)
+            sigma_k = sigmas[k].quantity.to_value(sigma.unit)
+
+            mask_valid = np.isfinite(mean) & np.isfinite(sigma) & (sigma.data != 0)
+            mask_valid_k = np.isfinite(mean_k) & np.isfinite(sigma_k) & (sigma_k != 0)
+            mask = mask_valid & mask_valid_k
+            mask_k = ~mask_valid & mask_valid_k
+
+            mean.data[mask] = (
+                (mean.data * sigma_k**2 + mean_k * sigma.data**2)
+                / (sigma.data**2 + sigma_k**2)
+            )[mask]
+            sigma.data[mask] = (
+                sigma.data * sigma_k / np.sqrt(sigma.data**2 + sigma_k**2)
+            )[mask]
+
+            mean.data[mask_k] = mean_k[mask_k]
+            sigma.data[mask_k] = sigma_k[mask_k]
+
+        ts = mean * mean / sigma / sigma + ts_diff
+        ts.data[~np.isfinite(ts.data)] = np.nan
+
+        kwargs = dict(
+            sed_type="dnde", reference_model=reference_model, meta=meta, gti=gti
+        )
+        return FluxMaps.from_maps(dict(dnde=mean, dnde_err=sigma, ts=ts), **kwargs)
+
+    elif method in ["distrib", "profile"]:
+        if dnde_scan_axis is None:
+            dnde_scan_axis = _default_scan_map(maps[0]).geom.axes["dnde"]
+        for k, map_ in enumerate(maps):
+            if method == "profile":
+                map_stat_scan = interpolate_profile_map(map_, dnde_scan_axis)
+            else:
+                map_stat_scan = approximate_profile_map(map_, dnde_scan_axis)
+            map_stat_scan.data[np.isnan(map_stat_scan.data)] = 0.0
+            if k == 0:
+                stat_scan = map_stat_scan
+            else:
+                stat_scan += map_stat_scan
+
+        return get_flux_map_from_profile(
+            {"stat_scan": stat_scan},
+            reference_model=reference_model,
+            meta=meta,
+            gti=gti,
+        )
+
+    else:
+        raise ValueError(
+            f'Invalid method provided : {method}. Available methods are : "gaussian_errors", "distrib", "profile"'
+        )
+
+
+def _default_scan_map(flux_map, dnde_scan_axis=None):
+    if dnde_scan_axis is None:
+        dnde_scan_axis = MapAxis(
+            _generate_scan_values() * flux_map.dnde_ref.squeeze(),
+            interp="lin",
+            node_type="center",
+            name="dnde",
+            unit=flux_map.dnde_ref.unit,
+        )
+
+    geom = flux_map.dnde.geom
+    geom_scan = geom.to_image().to_cube([dnde_scan_axis] + list(geom.axes))
+    return Map.from_geom(geom_scan, data=np.nan, unit="")
+
+
+def interpolate_profile_map(flux_map, dnde_scan_axis=None):
+    """Interpolate sparse likelihood profile to regular grid
+
+    Parameters
+    ----------
+    flux_map : `~gammapy.estimators.FluxMaps`
+        Flux map.
+    dnde_scan_axis : `~gammapy.maps.MapAxis`
+        Map axis providing the dnde values used to compute the profile.
+        Default is None and it will be derived from the flux_map.
+
+    Returns
+    -------
+    scan_map: `~gammapy.estimators.Maps`
+        Likelihood profile map.
+
+    """
+
+    stat_scan = _default_scan_map(flux_map, dnde_scan_axis)
+    dnde_scan_axis = stat_scan.geom.axes["dnde"]
+
+    mask_valid = ~np.isnan(flux_map.dnde.data)
+    dnde_scan_values = flux_map.dnde_scan_values.quantity.to_value(dnde_scan_axis.unit)
+
+    for ij, il, ik in zip(*np.where(mask_valid)):
+        spline = InterpolatedUnivariateSpline(
+            dnde_scan_values[ij, :, il, ik],
+            flux_map.stat_scan.data[ij, :, il, ik],
+            k=1,
+            ext="raise",
+            check_finite=True,
+        )
+        stat_scan.data[ij, :, il, ik] = spline(dnde_scan_axis.center)
+    return stat_scan
+
+
+def approximate_profile_map(
+    flux_map, dnde_scan_axis=None, sqrt_ts_threshold_ul="ignore"
+):
+    """Likelihood profile approximation assuming that probabilities distributions for
+    flux points correspond to asymmetric gaussians and for upper limits to complementary error functions.
+    Use available quantities among dnde, dnde_err, dnde_errp, dnde_errn, dnde_ul and ts.
+
+    Parameters
+    ----------
+    flux_map : `~gammapy.estimators.FluxMaps`
+        Flux map.
+    dnde_scan_axis : `~gammapy.maps.MapAxis`
+        Map axis providing the dnde values used to compute the profile.
+        Default is None and it will be derived from the flux_map.
+    sqrt_ts_threshold_ul : int
+        Threshold value in sqrt(TS) for upper limits.
+        Default is `ignore` and no threshold is applied.
+        Setting to `None` will use the one of `flux_map`.
 
+    Returns
+    -------
+    scan_map: `~gammapy.estimators.Maps`
+        Likelihood profile map.
+
+    """
+
+    stat_approx = _default_scan_map(flux_map, dnde_scan_axis)
+    dnde_coord = stat_approx.geom.get_coord()["dnde"].value
+
+    if sqrt_ts_threshold_ul is None:
+        sqrt_ts_threshold_ul = flux_map.sqrt_ts_threshold_ul
+
+    mask_valid = ~np.isnan(flux_map.dnde.data)
+    ij, il, ik = np.where(mask_valid)
+    loc = flux_map.dnde.data[mask_valid][:, None]
+    value = dnde_coord[ij, :, il, ik]
+    try:
+        mask_p = dnde_coord >= flux_map.dnde.data
+        mask_p2d = mask_p[ij, :, il, ik]
+        new_axis = np.ones(mask_p2d.shape[1], dtype=bool)[None, :]
+        scale = np.zeros(mask_p2d.shape)
+        scale[mask_p2d] = (flux_map.dnde_errp.data[mask_valid][:, None] * new_axis)[
+            mask_p2d
+        ]
+        scale[~mask_p2d] = (flux_map.dnde_errn.data[mask_valid][:, None] * new_axis)[
+            ~mask_p2d
+        ]
+    except AttributeError:
+        scale = flux_map.dnde_err.data[mask_valid]
+        scale = scale[:, None]
+    stat_approx.data[ij, :, il, ik] = ((value - loc) / scale) ** 2
+
+    try:
+        invalid_value = 999
+        stat_min_p = (stat_approx.data + invalid_value * (~mask_p)).min(
+            axis=1, keepdims=True
+        )
+        stat_min_m = (stat_approx.data + invalid_value * mask_p).min(
+            axis=1, keepdims=True
+        )
+
+        mask_minp = mask_p & (stat_min_p > stat_min_m)
+        stat_approx.data[mask_minp] = (stat_approx.data + stat_min_m - stat_min_p)[
+            mask_minp
+        ]
+        mask_minn = ~mask_p & (stat_min_m >= stat_min_p)
+        stat_approx.data[mask_minn] = (stat_approx.data + stat_min_p - stat_min_m)[
+            mask_minn
+        ]
+    except NameError:
+        pass
+
+    if not sqrt_ts_threshold_ul == "ignore" and sqrt_ts_threshold_ul is not None:
+        mask_ul = (flux_map.sqrt_ts.data < sqrt_ts_threshold_ul) & ~np.isnan(
+            flux_map.dnde_ul.data
+        )
+        ij, il, ik = np.where(mask_ul)
+        value = dnde_coord[ij, :, il, ik]
+        loc_ul = flux_map.dnde_ul.data[mask_ul][:, None]
+        scale_ul = flux_map.dnde_ul.data[mask_ul][:, None]
+        stat_approx.data[ij, :, il, ik] = -2 * np.log(
+            (special.erfc((-loc_ul + value) / scale_ul) / 2)
+            / (special.erfc((-loc_ul + 0) / scale_ul) / 2)
+        )
+
+    stat_approx.data[np.isnan(stat_approx.data)] = np.inf
+    stat_approx.data += -flux_map.ts.data - stat_approx.data.min(axis=1)
+    return stat_approx
+
+
+def get_flux_map_from_profile(
+    flux_map, n_sigma=1, n_sigma_ul=2, reference_model=None, meta=None, gti=None
+):
+    """Create a new flux map using the likehood profile (stat_scan)
+    to get ts, dnde, dnde_err, dnde_errp, dnde_errn, and dnde_ul.
+
+    Parameters
+    ----------
+    flux_maps : `~gammapy.estimators.FluxMaps` or dict of `~gammapy.maps.WcsNDMap`
+        Flux map or dict containing  a `stat_scan` entry
+    n_sigma : int
+        Number of sigma for flux error. Default is 1.
+    n_sigma_ul : int
+        Number of sigma for flux upper limits. Default is 2.
+    reference_model : `~gammapy.modeling.models.SkyModel`, optional
+        The reference model to use for conversions. If None, a model consisting
+        of a point source with a power law spectrum of index 2 is assumed.
+        Default is None and and the one of `flux_map` will be used if available
+    meta : dict, optional
+        Dict of metadata.
+        Default is None and and the oneof `flux_map` will be used if available
+    gti : `~gammapy.data.GTI`, optional
+        Maps GTI information.
+        Default is None and and the one of `flux_map` will be used if available
+
+    Returns
+    -------
+    flux_maps : `~gammapy.estimators.FluxMaps`
+        Flux map.
+
+    """
+
+    if isinstance(flux_map, dict):
+        output_maps = flux_map
+    else:
         if reference_model is None:
-            reference_model = maps[0].reference_model
-
-        # TODO : change this once we have stackable metadata objets
-        metas = [map_.meta for map_ in maps if map_.meta is not None]
-        meta = {}
-        if np.any(metas):
-            for data in metas:
-                meta.update(data)
-
-    for k in range(1, len(means)):
-
-        mean_k = means[k].quantity.to_value(mean.unit)
-        sigma_k = sigmas[k].quantity.to_value(sigma.unit)
-
-        mask_valid = np.isfinite(mean) & np.isfinite(sigma) & (sigma.data != 0)
-        mask_valid_k = np.isfinite(mean_k) & np.isfinite(sigma_k) & (sigma_k != 0)
-        mask = mask_valid & mask_valid_k
-        mask_k = ~mask_valid & mask_valid_k
-
-        mean.data[mask] = (
-            (mean.data * sigma_k**2 + mean_k * sigma.data**2)
-            / (sigma.data**2 + sigma_k**2)
-        )[mask]
-        sigma.data[mask] = (
-            sigma.data * sigma_k / np.sqrt(sigma.data**2 + sigma_k**2)
-        )[mask]
-
-        mean.data[mask_k] = mean_k[mask_k]
-        sigma.data[mask_k] = sigma_k[mask_k]
-
-    ts = -2 * np.log(
-        stats.norm.pdf(0, loc=mean, scale=sigma)
-        / stats.norm.pdf(mean, loc=mean, scale=sigma)
+            reference_model = flux_map.reference_model
+        if gti is None:
+            gti = flux_map.gti
+        if meta is None:
+            meta = flux_map.meta
+
+        output_maps = dict(
+            stat_scan=flux_map.stat_scan, dnde_scan_values=flux_map.dnde_scan_values
+        )
+
+    if getattr(flux_map, "dnde_scan_values", False):
+        dnde_coord = flux_map["dnde_scan_values"].quantity
+    else:
+        dnde_coord = flux_map["stat_scan"].geom.get_coord()["dnde"]
+
+    geom = (
+        flux_map["stat_scan"]
+        .geom.to_image()
+        .to_cube([flux_map["stat_scan"].geom.axes["energy"]])
     )
-    ts = Map.from_geom(mean.geom, data=ts)
 
+    ts = -flux_map["stat_scan"].data.min(axis=1) * u.Unit("")
+    ind = flux_map["stat_scan"].data.argmin(axis=1)
+    ij, ik, il = np.indices(ind.shape)
+    dnde = dnde_coord[ij, ind, ik, il]
+
+    maskp = dnde_coord > dnde
+    stat_diff = flux_map["stat_scan"].data - flux_map["stat_scan"].data.min(axis=1)
+
+    invalid_value = 999
+    ind = np.abs(stat_diff + invalid_value * maskp - n_sigma**2).argmin(axis=1)
+    dnde_errn = dnde - dnde_coord[ij, ind, ik, il]
+
+    ind = np.abs(stat_diff + invalid_value * (~maskp) - n_sigma**2).argmin(axis=1)
+    dnde_errp = dnde_coord[ij, ind, ik, il] - dnde
+
+    ind = np.abs(stat_diff + invalid_value * (~maskp) - n_sigma_ul**2).argmin(axis=1)
+    dnde_ul = dnde_coord[ij, ind, ik, il]
+
+    dnde_err = (dnde_errn + dnde_errp) / 2
+
+    maps = dict(
+        ts=ts,
+        dnde=dnde,
+        dnde_err=dnde_err,
+        dnde_errn=dnde_errn,
+        dnde_errp=dnde_errp,
+        dnde_ul=dnde_ul,
+    )
+    for key in maps.keys():
+        maps[key] = Map.from_geom(geom, data=maps[key].value, unit=maps[key].unit)
     kwargs = dict(sed_type="dnde", gti=gti, reference_model=reference_model, meta=meta)
-    return FluxMaps.from_maps(dict(dnde=mean, dnde_err=sigma, ts=ts), **kwargs)
+    output_maps.update(maps)
+    return FluxMaps.from_maps(output_maps, **kwargs)
 
 
 def _generate_scan_values(power_min=-6, power_max=2, relative_error=1e-2):


Test patch:
diff --git a/gammapy/estimators/map/tests/test_core.py b/gammapy/estimators/map/tests/test_core.py
index b02ba2a6a..ee156c14f 100644
--- a/gammapy/estimators/map/tests/test_core.py
+++ b/gammapy/estimators/map/tests/test_core.py
@@ -217,8 +217,14 @@ def test_combine_flux_maps(map_flux_estimate, wcs_flux_map, reference_model):
     stop = u.Quantity([1.5, 2.5], "min")
     gti = GTI.create(start, stop)
 
+    data = map_flux_estimate.copy()
+    ts = data["norm"].copy()
+    data["norm_err"].data[0, 0, 0] = 0.0
+    ts.data = (data["norm"].data / data["norm_err"].data) ** 2
+    data["ts"] = ts
+
     model = SkyModel(PowerLawSpectralModel(amplitude="1e-10 cm-2s-1TeV-1", index=2))
-    fe = FluxMaps(data=map_flux_estimate, reference_model=model, gti=gti)
+    fe = FluxMaps(data=data, reference_model=model, gti=gti)
     iE = 0
     energy = fe.geom.axes[0].center[iE]
 
@@ -245,7 +251,8 @@ def test_combine_flux_maps(map_flux_estimate, wcs_flux_map, reference_model):
     )
     assert fe_new.meta == fe.meta
 
-    fe = FluxMaps(wcs_flux_map, reference_model)
+    fe = FluxMaps(data, reference_model)
+
     fe_new = combine_flux_maps([fe, fe], reference_model=reference_model)
     assert_allclose(fe_new.dnde.quantity, fe.dnde.quantity)
     assert_allclose(
@@ -253,14 +260,16 @@ def test_combine_flux_maps(map_flux_estimate, wcs_flux_map, reference_model):
     )
     assert_allclose(fe_new.norm.quantity, fe.norm.quantity)
 
-    ts = -2 * np.log(
-        stats.norm.pdf(0, loc=fe.dnde.data, scale=fe.dnde_err.data)
-        / stats.norm.pdf(fe.dnde.data, loc=fe.dnde.data, scale=fe.dnde_err.data)
+    ts = (
+        -2
+        * np.log(
+            stats.norm.pdf(0, loc=fe.dnde.data, scale=fe.dnde_err.data)
+            / stats.norm.pdf(fe.dnde.data, loc=fe.dnde.data, scale=fe.dnde_err.data)
+        )
+        + (fe.ts - (fe.norm.data / fe.norm_err.data) ** 2) * 2
     )
     assert_allclose(fe_new.ts, ts * 2)
 
-    fe_new = combine_flux_maps(fe)
-
 
 def test_flux_map_properties(wcs_flux_map, reference_model):
     fluxmap = FluxMaps(wcs_flux_map, reference_model)
diff --git a/gammapy/estimators/map/tests/test_ts.py b/gammapy/estimators/map/tests/test_ts.py
index 7a531abb4..557e882a8 100644
--- a/gammapy/estimators/map/tests/test_ts.py
+++ b/gammapy/estimators/map/tests/test_ts.py
@@ -1,4 +1,5 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
+from copy import deepcopy
 import pytest
 import numpy as np
 from numpy.testing import assert_allclose
@@ -6,7 +7,12 @@ import astropy.units as u
 from astropy.coordinates import Angle, SkyCoord
 from gammapy.datasets import Datasets, MapDataset, MapDatasetOnOff
 from gammapy.estimators import TSMapEstimator
-from gammapy.estimators.utils import get_combined_significance_maps
+from gammapy.estimators.utils import (
+    approximate_profile_map,
+    combine_flux_maps,
+    get_combined_significance_maps,
+    get_flux_map_from_profile,
+)
 from gammapy.irf import EDispKernelMap, PSFMap
 from gammapy.maps import Map, MapAxis, WcsGeom
 from gammapy.modeling.models import (
@@ -304,14 +310,54 @@ def test_ts_map_stat_scan(fake_dataset):
     ts = np.abs(maps["stat_scan"].data.min(axis=1))
     assert_allclose(ts[success], maps_ref.ts.data[success], rtol=1e-3)
 
-    ind_best = maps.stat_scan.data.argmin(axis=1)
-    ij, ik, il = np.indices(ind_best.shape)
-
     dnde_ref = maps.dnde_ref.squeeze()
     assert maps.dnde_scan_values.unit == dnde_ref.unit
+
+    ind_best = maps.stat_scan.data.argmin(axis=1)
+    ij, ik, il = np.indices(ind_best.shape)
     norm = maps.dnde_scan_values.data[ij, ind_best, ik, il] / dnde_ref.value
     assert_allclose(norm[success], maps_ref.norm.data[success], rtol=1e-5)
 
+    stat_scan_aprrox = approximate_profile_map(maps, sqrt_ts_threshold_ul="ignore")
+    ts_aprrox = np.abs(stat_scan_aprrox.data.min(axis=1))
+    assert_allclose(ts_aprrox[success], maps_ref.ts.data[success], rtol=1e-3)
+
+    stat_scan_aprrox = approximate_profile_map(maps, sqrt_ts_threshold_ul=None)
+    ts_aprrox = np.abs(stat_scan_aprrox.data.min(axis=1))
+    assert_allclose(ts_aprrox[success], maps_ref.ts.data[success], rtol=1e-3)
+
+    maps_from_scan = get_flux_map_from_profile(maps)
+    assert_allclose(
+        maps_from_scan.ts.data[success], maps_ref.ts.data[success], rtol=1e-3
+    )
+
+    combined_map = combine_flux_maps([maps, maps], method="profile")
+    assert_allclose(combined_map.ts.data, 2 * ts, rtol=1e-4)
+    assert_allclose(combined_map.norm.data[success], norm[success], rtol=5e-2)
+
+    maps1 = deepcopy(maps)
+    combined_map = combine_flux_maps([maps, maps1], method="gaussian_errors")
+    assert_allclose(combined_map.ts.data[success], 2 * ts[success], rtol=1e-4)
+    assert_allclose(combined_map.norm.data[success], norm[success], rtol=5e-2)
+
+    combined_map = combine_flux_maps([maps, maps1], method="distrib")
+    assert_allclose(combined_map.ts.data, 2 * ts, rtol=1e-4)
+    assert_allclose(combined_map.norm.data[success], norm[success], rtol=5e-2)
+
+    combined_map = combine_flux_maps([maps, maps1], method="profile")
+    assert_allclose(combined_map.ts.data, 2 * ts, rtol=1e-4)
+    assert_allclose(combined_map.norm.data[success], norm[success], rtol=5e-2)
+
+    maps1._reference_model.parameters["amplitude"].value = 1
+
+    combined_map = combine_flux_maps([maps, maps1], method="distrib")
+    assert_allclose(combined_map.ts.data, 2 * ts, rtol=1e-4)
+    assert_allclose(combined_map.norm.data[success], norm[success], rtol=5e-2)
+
+    combined_map = combine_flux_maps([maps, maps1], method="profile")
+    assert_allclose(combined_map.ts.data, 2 * ts, rtol=1e-4)
+    assert_allclose(combined_map.norm.data[success], norm[success], rtol=5e-2)
+
     assert_allclose(maps.norm.data[success], maps_ref.norm.data[success], rtol=1e-4)
     assert_allclose(
         maps.norm_errn.data[success], maps_ref.norm_errn.data[success], rtol=5e-2
@@ -322,6 +368,8 @@ def test_ts_map_stat_scan(fake_dataset):
     assert_allclose(
         maps.norm_ul.data[success], maps_ref.norm_ul.data[success], rtol=5e-2
     )
+    with pytest.raises(ValueError):
+        combine_flux_maps([maps, maps1], method="test")
 
 
 def test_ts_map_with_model(fake_dataset):