Likelihood Analysis from Python

Also see:

pyLikelihood Usage Notes

Installation

To have access to likelihood analysis modules from Python, set your INST_DIR environment variable to point to the directory where you have installed your ScienceTools distribution, and source or execute the appropriate python_setup file in the pyLikelihood package.  Following the example in the linux installation instructions, one would do:

noric07[jchiang] setenv INST_DIR ${HOME}/ST
noric07[jchiang] source ${INST_DIR}/pyLikelihood/v1r3p3/cmt/python_setup.csh

whereas under Windows, one might do:

C:\> set INST_DIR=%HOME%\ST
C:\> %INST_DIR%\pyLikelihood\v0r0p1\Visual\python_setup.bat

Note: For Windows, the python_setup.bat file lives in the Visual subdirectory, whereas for linux, the setup files, python_setup.[c]sh, live in the cmt subdirectory.

Run the python version that is in the distribution. For version v9r4p2 of the Science Tools is 2.5.1.

Package Description

Unbinned and binned analysis have been implemented in the UnbinnedAnalysis.py and BinnedAnalysis.py modules, respectively. Each module contains two classes: UnbinnedObs and UnbinnedAnalysis in the former, and BinnedObs and BinnedAnalysis in the latter.

The "Obs" classes differ in how they are constructed since the respective analyses require different sets of input data. Nonetheless, both implementations are intended to encapsulate the information associated with a specific observation. In this context, an observation is defined in terms of the extraction region in the data space of photon arrival time, measured energy, and direction. An observation also comprises ancillary information specific to the data selections such as the exposure information and the response functions to be used.

The signature of the UnbinnedObs constructor is: 

class UnbinnedObs(object):
    def __init__(self, eventFile=None, scFile=None, expMap=None,
                 expCube=None, irfs='HANDOFF', checkCuts=True):

and its parameters are:

  eventFile Event data file name(s), either a single file name or a tuple or list of file names. If given as a single file name, it can also be an ascii file of FITS file names.
  scFile Spacecraft data file name(s). This can also be an ascii file of FITS file names.
  expMap Exposure map (produced by gtexpmap). The geometry of this map matches that of the extraction region used to create the event files.
  expCube Live-time cube (produced by gtltcube). This is an off-axis histogram of live-times at each point in the sky, partitioned as nested HealPix.
  irfs Instrument response functions, e.g., HANDOFF
  checkCuts A flag used for debugging purposes. This should be set to True for a standard analysis.

The BinnedObs constructor signature is: 

class BinnedObs(object):
    def __init__(self, srcMaps, expCube, binnedExpMap=None, irfs='HANDOFF'):

and its parameters are:

  srcMaps This is a counts map file with axes in position and energy and source map extensions created using gtsrcmaps. The counts map may be created with gtcmap. If one provides a counts map without a complete set of source map extensions, the BinnedAnalysis class will examine the model definition file and will compute source maps in memory for all the sources that do not have source maps in the counts map file.
  expCube The live-time cube.
  binnedExpMap Exposure map appropriate for binned likelihood analysis. This is optional. If it is not specified, then Likelihood will compute an appropriate map and output it as binned_exposure.fits. These exposure maps are matched to the geometry of the counts map in srcMaps.
  irfs Instrument response functions, e.g., HANDOFF

Consult the unbinned and binned analysis tutorials for more detailed description of the inputs to these classes.

The "Analysis" classes have the same public interfaces and can therefore be used interchangably in scripts. The UnbinnedAnalysis constructor is defined as: 

class UnbinnedAnalysis(AnalysisBase):
    def __init__(self, observation, srcModel=None,  optimizer='Minuit'):

and its parameters are:

  observation An instance of the UnbinnedObs class.
  srcModel Name of the XML file containing the source model definition.
  optimizer Optimizer package to use: Minuit, NewMinuit Drmnfb, Drmngb, Lbfgs

The BinnedAnalysis constructor is identical except that it takes an instance of BinnedObs as the observation object.

Factoring into separate analysis and observation classes affords flexibility in mixing and matching observations and models in a single Python session or script while preserving computational resources, so that one can do something like this: 

>>> analysis1 = UnbinnedAnalysis(unbinnedObs, "model1.xml")
>>> analysis2 = UnbinnedAnalysis(unbinnedObs, "model2.xml")

Even though the two UnbinnedAnalysis instances access the same data in memory, distinct source models may be fit concurrently to those data without interfering with one another. This will be useful in comparing, via a likelihood ratio test, for example, how well one model compares to another.

Starting up from the Python command line

Here we perform an unbinned analysis. The steps that must be taken before proceeding, e.g., live-time, exposure and diffuse response calculations, are outlined in the unbinned analysis tutorial.

First, we import the data and analysis classes: 

>>> from UnbinnedAnalysis import *

then we create an UnbinnedObs object: 

>>> my_obs = 
UnbinnedObs(eventFiles, 'FT2.fits', expMap='expMap.fits', expCube='expCube.fits',
   irfs='HANDOFF')

Here, eventFiles is an ascii file containing the names of the event files: 

glitch [43] [cillis]: cat eventFiles 
eg_dif_filtered.fits
galdif_filtered.fits
ptsrcs_filtered.fits
glitch [44] [cillis]:

One could also have entered in a tuple or list if file names or use glob to generate the list. For these data, the following are equivalent to the above: 

>>> my_obs = UnbinnedObs(('eg_dif_filtered.fits', 'galdif_filtered.fits',
                          'ptsrcs_filtered.fits'),
                         'FT2.fits', expMap='expMap.fits', 
                         expCube ='expCube.fits', irfs='HANDOFF')

OR 

>>> my_obs = UnbinnedObs(glob.glob('*filtered.fits'), 'FT2.fits',
                         expMap='expMap.fits', expCube='expCube.fits', 
                         irfs='HANDOFF')

One may omit all of the arguments in the UnbinnedObs class constructor, in which case a small GUI dialog is launched that allows one to browse the file system, use wild cards for specifying groups of files, etc.. So, entering: 

>>> my_obs = UnbinnedObs()

launches this dialog:

Obs_dialog0.png


The buttons on the left will open file dialog boxes using the value shown in the text entry fields as a filter. One can include wild cards or comma-separated lists of files in the text entry field. Setting the entries as follows will create an UnbinnedObs object equivalent to the previous examples:

Obs_dialog1.png


Now we are ready to create an instance of UnbinnedAnalysis:

>>> analysis = UnbinnedAnalysis(my_obs, "srcModel.xml")

Here srcModel.xml is the xml file containing the model definition for Likelihood.

Just as with the UnbinnedObs class, a small GUI is provided that allows one to browse the file system in case one omits the name of the model XML file. For example: 

>>> analysis = UnbinnedAnalysis(my_obs)

launches:

Analysis_dialog.png

The UnbinnedObs and UnbinnedAnalysis classes both have been implemented to allow one to see easily what data these objects comprise: 

>>> print my_obs
Event file(s): ('eg_dif_filtered.fits', 'galdif_filtered.fits', 'ptsrcs_filtered.fits')
Spacecraft file(s): FT2.fits
Exposure map: expMap.fits
Exposure cube: expCube.fits
IRFs: HANDOFF
>>> 
>>> print analysis
Event file(s): ('eg_dif_filtered.fits', 'galdif_filtered.fits', 'ptsrcs_filtered.fits')
Spacecraft file(s): FT2.fits
Exposure map: expMap.fits
Exposure cube: expCube.fits
IRFs: HANDOFF
Source model file: srcModel.xml
Optimizer: Drmngb

>>>

Starting from a Python Script

It is convenient to prepare a script that creates the UnbinnedObs and UnbinnedAnalysis objects with the proper parameter values so that analysis of the data may be easily revisited in separate sessions:

glitch [61] [cillis]: cat my_analysis
import glob
from UnbinnedAnalysis import *

eventFiles = glob.glob('*filtered.fits')
obs = UnbinnedObs(eventFiles, 'FT2.fits', expMap='expMap.fits',
                  expCube='expCube.fits', irfs='HANDOFF')
like = UnbinnedAnalysis(obs, 'srcModel.xml')
glitch [62] [cillis]:

This just needs to be imported at the Python prompt: 

>>> from my_analysis import *

Viewing and Setting Parameters, Plotting, Fitting, etc.

Python provides various means of introspection. For example, the dir command shows an object's attributes:

>>> dir(like)

['Ts', '_Nobs', '__call__', '__class__', '__delattr__', '__dict__', '__doc__',
 '__getattribute__',  '__getitem__', '__hash__', '__init__', '__module__', '__new__',
 '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__setitem__', '__str__',
 '__weakref__', '_errors', '_importPlotter', '_inputs', '_isDiffuseOrNearby', '_normPar',
 '_npredValues', '_plotData', '_plotResiduals', '_plotSource', '_renorm', '_separation',
 '_srcCnts', '_srcDialog', 'addSource', 'covariance', 'deleteSource', 'disp', 'e_vals',
 'energies', 'fit', 'freeze', 'logLike', 'maxdist', 'model', 'nobs', 'observation',
 'oplot', 'optimizer', 'plot', 'resids', 'setPlotter', 'sourceNames', 'srcModel', 'state',
 'thaw', 'tol', 'writeCountsSpectra', 'writeXml']

>>>

The objects themselves may also offer introspection. Both UnbinnedAnalysis and BinnedAnalysis objects have a model attribute that allows one to view and manipulate the current state of the various fit parameters.

To view the current model parameters, do: 

>>> like.model
3C 273
   Spectrum: PowerLaw
0      Prefactor:  1.000e+01  0.000e+00  1.000e-03  1.000e+03 ( 1.000e-09)
1          Index: -2.100e+00  0.000e+00 -5.000e+00 -1.000e+00 ( 1.000e+00)
2          Scale:  1.000e+02  0.000e+00  3.000e+01  2.000e+03 ( 1.000e+00) fixed

3C 279
   Spectrum: PowerLaw
3      Prefactor:  1.000e+01  0.000e+00  1.000e-03  1.000e+03 ( 1.000e-09)
4          Index: -2.000e+00  0.000e+00 -5.000e+00 -1.000e+00 ( 1.000e+00)
5          Scale:  1.000e+02  0.000e+00  3.000e+01  2.000e+03 ( 1.000e+00) fixed

Extragalactic Diffuse
   Spectrum: PowerLaw
6      Prefactor:  1.600e+00  0.000e+00  1.000e-05  1.000e+02 ( 1.000e-07)
7          Index: -2.100e+00  0.000e+00 -3.500e+00 -1.000e+00 ( 1.000e+00) fixed
8          Scale:  1.000e+02  0.000e+00  5.000e+01  2.000e+02 ( 1.000e+00) fixed

GalProp Diffuse
   Spectrum: ConstantValue
9          Value:  1.000e+00  0.000e+00  0.000e+00  1.000e+01 ( 1.000e+00) fixed

>>>

The columns are:

  • parameter index
  • parameter name
  • parameter value
  • error estimate (zero before fitting)
  • lower bound
  • upper bound
  • parameter scaling (in parentheses)
  • fixed (or not)

By default, the plot() method uses pyROOT to create counts spectra plots for comparing the model fits to the data. The following creates two plot windows, one for the counts spectra and one for the fractional residuals of the fit: 

>>> like.plot()

 

One can set parameter values using the index: 

>>> like[6] = 1.595

optimizer::Parameter member functions are automatically dispatched to the underlying C++ class from Python so that one can set the "free" flag, set the bounds, and the scale factor 

  >>> like[7].setFree(1) >>> like[0]=1.01
  >>> like[6].setBounds(0.5,2)
  >>> like.model
  3C 273
  Spectrum: PowerLaw
  0      Prefactor:  1.010e+00  0.000e+00  1.000e-03  1.000e+03 ( 1.000e-09)
  1          Index: -2.100e+00  0.000e+00 -5.000e+00 -1.000e+00 ( 1.000e+00)
  2          Scale:  1.000e+02  0.000e+00  3.000e+01  2.000e+03 ( 1.000e+00) fixed
                          
  3C 279
  Spectrum: PowerLaw
  3      Prefactor:  1.000e+01  0.000e+00  1.000e-03  1.000e+03 ( 1.000e-09)
  4          Index: -2.000e+00  0.000e+00 -5.000e+00 -1.000e+00 ( 1.000e+00)
  5          Scale:  1.000e+02  0.000e+00  3.000e+01  2.000e+03 ( 1.000e+00) fixed
                          
  Extragalactic Diffuse
  Spectrum: PowerLaw
  6      Prefactor:  1.595e+00  0.000e+00  5.000e-01  2.000e+00 ( 1.000e-07)
  7          Index: -2.100e+00  0.000e+00 -3.500e+00 -1.000e+00 ( 1.000e+00)
  8          Scale:  1.000e+02  0.000e+00  5.000e+01  2.000e+02 ( 1.000e+00) fixed
                          
  GalProp Diffuse
  Spectrum: ConstantValue
  9          Value:  1.000e+00  0.000e+00  0.000e+00  1.000e+01 ( 1.000e+00) fixed

Now, perform a fit, suppressing Minuit's screen output; the resulting value -log(Likelihood) is returned:

>>> like.fit(verbosity=0)
9365.6790299170352

Overplot the new result in red, the default: 

>>> like.oplot()

The logLike attribute is a Likelihood::LogLike object, and so its member functions are exposed as Python methods: 

>>> like.logLike.writeXml("fitted_model.xml")

Owned by: Jim Chiang

Last updated by: Jim Chiang 09/29/2008