/
residmap.rst
145 lines (107 loc) · 5.12 KB
/
residmap.rst
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
.. _residmap:
Residual Map
============
:py:meth:`~fermipy.gtanalysis.GTAnalysis.residmap` calculates the
residual between smoothed data and model maps. Whereas a TS map is
only sensitive to positive deviations with respect to the model,
:py:meth:`~fermipy.gtanalysis.GTAnalysis.residmap` is sensitive to
both positive and negative residuals and therefore can be useful for
assessing the model goodness-of-fit. The significance of the
data/model residual at map position (*i*, *j*) is given by
.. math::
\sigma_{ij}^2 = 2 \mathrm{sgn}(\tilde{n}_{ij} - \tilde{m}_{ij})
\left(\ln L_{P}(\tilde{n}_{ij},\tilde{n}_{ij}) - \ln L_{P}(\tilde{n}_{ij},\tilde{m}_{ij})\right)
.. math::
\mathrm{with} \quad
\tilde{m}_{ij} = \sum_{k} (m_{k} \ast f_{k})_{ij} \quad \tilde{n}_{ij} = \sum_{k}(n_{k} \ast f_{k})_{ij}
\quad \ln L_{P}(n,m) = n\ln(m) - m
where *n*\ :sub:`k` and *m*\ :sub:`k` are the data and model maps at
energy plane *k* and *f*\ :sub:`k` is the convolution kernel. The
convolution kernel is proportional to the counts expectation at a
given pixel and normalized such that
.. math::
f_{ijk} = s_{ijk} \left(\sum_{ijk} s_{ijk}^{2}\right)^{-1}
where *s* is the expectation counts cube for a pure signal normalized to one.
Examples
--------
The spatial and spectral properties of the convolution kernel are
defined with the ``model`` dictionary argument. All source models are
supported as well as a gaussian kernel (defined by setting
*SpatialModel* to *Gaussian*).
.. code-block:: python
# Generate residual map for a Gaussian kernel with Index=2.0 and
# radius (R_68) of 0.3 degrees
model = {'Index' : 2.0,
'SpatialModel' : 'Gaussian', 'SpatialWidth' : 0.3 }
maps = gta.residmap('fit1',model=model)
# Generate residual map for a power-law point source with Index=2.0 for
# E > 3.16 GeV
model = {'Index' : 2.0, 'SpatialModel' : 'PointSource'}
maps = gta.residmap('fit1_emin35',model=model,loge_bounds=[3.5,None])
# Generate residual maps for a power-law point source with Index=1.5, 2.0, and 2.5
model={'SpatialModel' : 'PointSource'}
maps = []
for index in [1.5,2.0,2.5]:
model['Index'] = index
maps += [gta.residmap('fit1',model=model)]
:py:meth:`~fermipy.gtanalysis.GTAnalysis.residmap` returns a ``maps``
dictionary containing `~fermipy.skymap.Map` representations of the
residual significance and amplitude as well as the smoothed data and
model maps. The contents of the output dictionary are described in
the following table.
============= ====================== ======================================
Key Type Description
============= ====================== ======================================
sigma `~fermipy.skymap.Map` Residual significance in sigma.
excess `~fermipy.skymap.Map` Residual amplitude in counts.
data `~fermipy.skymap.Map` Smoothed counts map.
model `~fermipy.skymap.Map` Smoothed model map.
files dict File paths of the FITS image
files generated by this method.
src_dict dict Source dictionary with the
properties of the convolution kernel.
============= ====================== ======================================
The ``write_fits`` and ``write_npy`` options can used to write the
output to a FITS or numpy file. All output files are prepended with
the `prefix` argument.
Diagnostic plots can be generated by setting ``make_plots=True`` or by
passing the output dictionary to
`~fermipy.plotting.AnalysisPlotter.make_residmap_plots`:
.. code-block:: python
maps = gta.residmap('fit1',model=model, make_plots=True)
gta.plotter.make_residmap_plots(maps, roi=gta.roi)
This will generate the following plots:
* ``residmap_excess`` : Smoothed excess map (data-model).
* ``residmap_data`` : Smoothed data map.
* ``residmap_model`` : Smoothed model map.
* ``residmap_sigma`` : Map of residual significance. The color map is
truncated at -5 and 5 sigma with labeled isocontours at 2 sigma intervals
indicating values outside of this range.
* ``residmap_sigma_hist`` : Histogram of significance values for all
points in the map. Overplotted are distributions for the best-fit
Gaussian and a unit Gaussian.
.. |image_sigma| image:: residmap_gaussian_sigma.png
:width: 100%
.. |image_excess| image:: residmap_gaussian_excess.png
:width: 100%
.. |image_sigma_hist| image:: residmap_gaussian_sigma_hist.png
:width: 100%
.. csv-table::
:header: Residual Significance Map, Significance Histogram
:widths: 50, 50
|image_sigma|, |image_sigma_hist|
Configuration
-------------
The default configuration of the method is controlled with the
:ref:`config_residmap` section of the configuration file. The default
configuration can be overriden by passing the option as a *kwargs*
argument to the method.
.. csv-table:: *residmap* Options
:header: Option, Default, Description
:file: ../config/residmap.csv
:delim: tab
:widths: 10,10,80
Reference/API
-------------
.. automethod:: fermipy.gtanalysis.GTAnalysis.residmap
:noindex: