TSMapEstimator#

class gammapy.estimators.TSMapEstimator(model=None, kernel_width=None, downsampling_factor=None, n_sigma=1, n_sigma_ul=2, threshold=None, rtol=0.01, selection_optional=None, energy_edges=None, sum_over_energy_groups=True, n_jobs=None, parallel_backend=None, norm=None, max_niter=100)[source]#

Bases: Estimator, ParallelMixin

Compute test statistic map from a MapDataset using different optimization methods.

The map is computed fitting by a single parameter norm fit. The fit is simplified by finding roots of the derivative of the fit statistics using various root finding algorithms. The approach is described in Appendix A in Stewart (2009).

Parameters:

modelSkyModel

Source model kernel. If set to None, assume spatail model: point source model, PointSpatialModel. spectral model: PowerLawSpectral Model of index 2

kernel_widthAngle

Width of the kernel to use: the kernel will be truncated at this size

n_sigmaint

Number of sigma for flux error. Default is 1.

n_sigma_ulint

Number of sigma for flux upper limits. Default is 2.

downsampling_factorint

Sample down the input maps to speed up the computation. Only integer values that are a multiple of 2 are allowed. Note that the kernel is not sampled down, but must be provided with the downsampled bin size.

thresholdfloat, optional

If the test statistic value corresponding to the initial flux estimate is not above this threshold, the optimizing step is omitted to save computing time. Default is None.

rtolfloat

Relative precision of the flux estimate. Used as a stopping criterion for the norm fit. Default is 0.01.

selection_optionallist of str, optional

Which maps to compute besides TS, sqrt(TS), flux and symmetric error on flux. Available options are:

“all”: all the optional steps are executed

“errn-errp”: estimate asymmetric error on flux.

“ul”: estimate upper limits on flux.

“stat_scan”: estimate likelihood profile

Default is None so the optional steps are not executed.

energy_edgeslist of Quantity, optional

Edges of the target maps energy bins. The resulting bin edges won’t be exactly equal to the input ones, but rather the closest values to the energy axis edges of the parent dataset. Default is None: apply the estimator in each energy bin of the parent dataset. For further explanation see Estimators (DL4 to DL5, and DL6).

sum_over_energy_groupsbool

Whether to sum over the energy groups or fit the norm on the full energy cube.

normParameter or dict

Norm parameter used for the likelihood profile computation on a fixed norm range. Only used for “stat_scan” in selection_optional. Default is None and a new parameter is created automatically, with value=1, name=”norm”, scan_min=-100, scan_max=100, and values sampled such as we can probe a 0.1% relative error on the norm. If a dict is given the entries should be a subset of Parameter arguments.

n_jobsint

Number of processes used in parallel for the computation. Default is one, unless N_JOBS_DEFAULT was modified. The number of jobs limited to the number of physical CPUs.

parallel_backend{“multiprocessing”, “ray”}

Which backend to use for multiprocessing. Defaults to BACKEND_DEFAULT.

max_niterint

Maximal number of iterations used by the root finding algorithm. Default is 100.

Notes

Negative \(TS\) values are defined as following:

\[\begin{split}TS = \left \{ \begin{array}{ll} -TS \text{ if } F < 0 \\ TS \text{ else} \end{array} \right.\end{split}\]

Where \(F\) is the fitted flux norm.

References

[Stewart2009]

Examples

>>> import astropy.units as u
>>> from gammapy.estimators import TSMapEstimator
>>> from gammapy.datasets import MapDataset
>>> from gammapy.modeling.models import (SkyModel, PowerLawSpectralModel,PointSpatialModel)
>>> spatial_model = PointSpatialModel()
>>> spectral_model = PowerLawSpectralModel(amplitude="1e-22 cm-2 s-1 keV-1", index=2)
>>> model = SkyModel(spatial_model=spatial_model, spectral_model=spectral_model)
>>> dataset = MapDataset.read("$GAMMAPY_DATA/fermi-3fhl-gc/fermi-3fhl-gc.fits.gz")
>>> estimator = TSMapEstimator(
...            model, kernel_width="1 deg", energy_edges=[10, 100] * u.GeV, downsampling_factor=4
...        )
>>> maps = estimator.run(dataset)
>>> print(maps)
FluxMaps
--------

  geom                   : WcsGeom
  axes                   : ['lon', 'lat', 'energy']
  shape                  : (np.int64(400), np.int64(200), 1)
  quantities             : ['ts', 'norm', 'niter', 'norm_err', 'npred', 'npred_excess', 'stat', 'stat_null', 'success']
  ref. model             : pl
  n_sigma                : 1
  n_sigma_ul             : 2
  sqrt_ts_threshold_ul   : 2
  sed type init          : likelihood

Attributes Summary

`config_parameters`	Configuration parameters.
`n_jobs`	Number of jobs as an integer.
`parallel_backend`	Parallel backend as a string.
`selection_all`	Which quantities are computed.
`selection_optional`
`tag`

Methods Summary

`copy`()	Copy estimator.
`estimate_fit_input_maps`(dataset)	Estimate fit input maps.
`estimate_flux_default`(dataset[, kernel, ...])	Estimate default flux map using a given kernel.
`estimate_flux_map`(datasets)	Estimate flux and test statistic maps for single dataset.
`estimate_kernel`(dataset)	Get the convolution kernel for the input dataset.
`estimate_mask_default`(dataset)	Compute default mask where to estimate test statistic values.
`estimate_pad_width`(dataset[, kernel])	Estimate pad width of the dataset.
`run`(datasets)	Run test statistic map estimation.

Attributes Documentation

config_parameters#: Configuration parameters.

n_jobs#: Number of jobs as an integer.

parallel_backend#: Parallel backend as a string.

selection_all#: Which quantities are computed.

selection_optional#

tag = 'TSMapEstimator'#

Methods Documentation

copy()#: Copy estimator.

estimate_fit_input_maps(dataset)[source]#

Estimate fit input maps.

Parameters:

datasetMapDataset: Map dataset.

Returns:

mapsdict of Map: Maps dictionary.

estimate_flux_default(dataset, kernel=None, exposure=None)[source]#

Estimate default flux map using a given kernel.

Parameters:

datasetMapDataset: Input dataset.
kernelWcsNDMap: Source model kernel.
exposureWcsNDMap: Exposure map on reconstructed energy.

Returns:

fluxWcsNDMap: Approximate flux map.

estimate_flux_map(datasets)[source]#

Estimate flux and test statistic maps for single dataset.

Parameters:

datasetDatasets or MapDataset: Map dataset or Datasets (list of MapDataset with the same spatial geometry).

estimate_kernel(dataset)[source]#

Get the convolution kernel for the input dataset.

Convolves the model with the IRFs at the center of the dataset, or at the nearest position with non-zero exposure.

Parameters:

datasetMapDataset: Input dataset.

Returns:

kernelMap: Kernel map.

static estimate_mask_default(dataset)[source]#

Compute default mask where to estimate test statistic values.

Parameters:

datasetMapDataset: Input dataset.

Returns:

maskWcsNDMap: Mask map.

estimate_pad_width(dataset, kernel=None)[source]#

Estimate pad width of the dataset.

Parameters:

datasetMapDataset: Input MapDataset.
kernelWcsNDMap: Source model kernel.

Returns:

pad_widthtuple: Padding width.

run(datasets)[source]#

Run test statistic map estimation.

Requires a MapDataset with counts, exposure and background_model properly set to run.

Parameters:

datasetDatasets or MapDataset: Map dataset or Datasets (list of MapDataset with the same spatial geometry).

Returns:

mapsdict

Dictionary containing result maps. Keys are:

ts : delta(TS) map

sqrt_ts : sqrt(delta(TS)), or significance map

flux : flux map

flux_err : symmetric error map

flux_ul : upper limit map.

Notes

The progress bar can be displayed for this function.