#
# Copyright (C) 2009, 2015, 2016, 2017, 2018, 2019, 2020, 2022
# Smithsonian Astrophysical Observatory
#
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
#
import warnings
from configparser import ConfigParser
import numpy
from sherpa.utils import NoNewAttributesAfterInit, igamc
from sherpa.utils.err import FitErr, StatErr
from sherpa.data import DataSimulFit
from sherpa.models import SimulFitModel
from . import _statfcts
from sherpa import get_config
__all__ = ('Stat', 'Cash', 'CStat', 'LeastSq',
'Chi2Gehrels', 'Chi2ConstVar', 'Chi2DataVar', 'Chi2ModVar',
'Chi2XspecVar', 'Chi2',
'UserStat', 'WStat')
config = ConfigParser()
config.read(get_config())
# truncation_flag indicates whether or not model truncation
# should be performed. If true, use the truncation_value from
# the config file.
truncation_flag = config.get('statistics', 'truncate',
fallback='True').upper()
truncation_value = float(config.get('statistics', 'trunc_value',
fallback=1.0e-25))
if (bool(truncation_flag) is False or truncation_flag == "FALSE" or
truncation_flag == "NONE" or truncation_flag == "0"):
truncation_value = 1.0e-25
[docs]class Stat(NoNewAttributesAfterInit):
"""The base class for calculating a statistic given data and model."""
# Used by calc_stat
#
_calc = None
# This should be overridden by derived classes and set to True if the rstat and qvalue
# figures can be calculated for that class of statistics.
_can_calculate_rstat = None
def __init__(self, name):
self.name = name
NoNewAttributesAfterInit.__init__(self)
def __repr__(self):
if self.__doc__ is not None:
return self.__doc__
return ("<%s statistic instance '%s'>" %
(type(self).__name__, self.name))
@staticmethod
def _bundle_inputs(data, model):
"""Convert input into SimulFit instances.
Convert the inputs into `sherpa.data.DataSimulFit` and
`sherpa.models.model.SimulFitModel`
instances.
Parameters
----------
data : `sherpa.data.Data` or `sherpa.data.DataSimulFit`
The data set, or sets, to use.
model : `sherpa.models.model.Model` or `sherpa.models.model.SimulFitModel`
The model expression, or expressions. If a
`~sherpa.models.model.SimulFitModel`
is given then it must match the number of data sets in the
data parameter.
Returns
-------
data : `sherpa.data.DataSimulFit`
If the input was a `~sherpa.data.DataSimulFit` object
then this is just the input value.
model : `sherpa.models.model.SimulFitModel`
If the input was a `~sherpa.models.model.SimulFitModel`
object then this is just the input value.
"""
if not isinstance(data, DataSimulFit):
data = DataSimulFit('simulfit data', (data,))
if not isinstance(model, SimulFitModel):
model = SimulFitModel('simulfit model', (model,))
return data, model
@staticmethod
def _check_has_bins(data):
"""Raise an error if there are no noticed bins in the dataset.
Parameters
----------
data : `sherpa.data.DataSimulFit`
The data sets to use.
Raises
------
FitErr
Notes
-----
It is unclear whether this should error out if any particular
dataset has no noticed bins, or only if all the datasets have
no noticed bins (the latter approach is taken).
"""
for dset in data.datasets:
# Assume that the error column does not need to be
# calculated for this check, so staterrfunc can be
# None.
dep, _, _ = dset.to_fit(staterrfunc=None)
if numpy.iterable(dep) and len(dep) > 0:
return
raise FitErr('nobins')
@staticmethod
def _check_sizes_match(data, model):
"""Raise an error if number of datasets and models do not match.
Parameters
----------
data : `sherpa.data.DataSimulFit`
The data sets to use.
model : `sherpa.models.model.SimulFitModel`
The model expressions for each data set. It must match
the data parameter (the models are in the same order
as the data objects).
Raises
------
StatErr
"""
ndata = len(data.datasets)
nmdl = len(model.parts)
if ndata != nmdl:
raise StatErr('mismatch',
'number of data sets', ndata,
'model expressions', nmdl)
def _validate_inputs(self, data, model):
"""Ensure that the inputs are correct for the statistic.
The default behavior is to check that the data contains
at least one bin and that the number of datasets matches
the number of models. It also converts single values to
simultaneous objects, if necessary.
Parameters
----------
data : `sherpa.data.Data` or `sherpa.data.DataSimulFit`
The data set, or sets, to use.
model : `sherpa.models.model.Model` or `sherpa.models.model.SimulFitModel`
The model expressions for each data set. It must match
the data parameter (the models are in the same order
as the data objects).
Returns
-------
data : `sherpa.data.DataSimulFit`
If the input was a `~sherpa.data.DataSimulFit` object
then this is just the input value.
model : `sherpa.models.model.SimulFitModel`
If the input was a `~sherpa.models.model.SimulFitModel`
object then this is just the input value.
"""
if self._calc is None:
# This is a programmer error rather than a user error,
# so use NotImplementedError rather than
# StatErr('nostat', self.name, '_calc')
#
raise NotImplementedError("_calc method has not been set")
data, model = self._bundle_inputs(data, model)
self._check_has_bins(data)
self._check_sizes_match(data, model)
return data, model
def _get_fit_model_data(self, data, model):
data, model = self._validate_inputs(data, model)
fitdata = data.to_fit(staterrfunc=self.calc_staterror)
modeldata = data.eval_model_to_fit(model)
return fitdata, modeldata
# TODO:
# - should this accept sherpa.data.Data input instead of
# "raw" data (i.e. to match calc_stat)
# - should this be moved out of the base Stat class since
# isn't relevant for likelihood statistics?
#
[docs] def calc_staterror(self, data):
"""Return the statistic error values for the data.
Parameters
----------
data : scalar or 1D array of numbers
The data values.
Returns
-------
staterror : scalar or array of numbers
The errors for the input data values (matches the data
argument).
"""
raise NotImplementedError
# TODO: add *args, **kwargs?
[docs] def calc_stat(self, data, model):
"""Return the statistic value for the data and model.
Parameters
----------
data : `sherpa.data.Data` or `sherpa.data.DataSimulFit`
The data set, or sets, to use.
model : `sherpa.models.model.Model` or `sherpa.models.model.SimulFitModel`
The model expression, or expressions. If a
`sherpa.models.model.SimulFitModel`
is given then it must match the number of data sets in the
data parameter.
Returns
-------
statval : number
The value of the statistic.
fvec : array of numbers
The per-bin "statistic" value.
"""
raise NotImplementedError
[docs] def goodness_of_fit(self, statval, dof):
"""Return the reduced statistic and q value.
The reduced statisitc is conceptually simple, as it is just
statistic / degrees-of-freedom, but it is not meaningful for
all statistics, and it is only valid if there are any degrees
of freedom.
Parameters
----------
statval : float
The statistic value. It is assumed to be finite.
dof : int
The number of degrees of freedom, which may be 0 or negative.
Returns
-------
rstat : float or NaN or None
The reduced statistic. If the statistic does not support
a goodness of fit then the return value is `None`. If it does
then NaN is returned if either the number of degrees of freedom is
0 (or less), or the statistic value is less than 0.
qval : float or NaN or None
The q value. If the statistic does not support
a goodness of fit then the return values are `None`. If it does
then NaN is returned if either the number of degrees of freedom is
0 (or less), or the statistic value is less than 0.
"""
if not self._can_calculate_rstat:
return None, None
if dof > 0 and statval >= 0.0:
qval = igamc(dof / 2.0, statval / 2.0)
rstat = statval / dof
return rstat, qval
return numpy.nan, numpy.nan
[docs]class Likelihood(Stat):
"""Likelihood functions"""
def __init__(self, name='likelihood'):
Stat.__init__(self, name)
[docs] @staticmethod
def calc_staterror(data):
# Likelihood stats do not have 'errors' associated with them.
# return 1 to avoid dividing by 0 by some optimization methods.
return numpy.ones_like(data)
def _check_background_subtraction(self, data):
"""Raise an error if any dataset has been background subtracted.
Parameters
----------
data : a DataSimulFit instance
The data sets to use.
Raises
------
FitErr
"""
for dobj in data.datasets:
if getattr(dobj, 'subtracted', False):
# TODO:
# Historically this has been a FitErr, but it
# would make more sense to be a StatErr. This could
# break people's code, so hold off for now
# (October 2016).
#
raise FitErr('statnotforbackgsub', self.name)
def _validate_inputs(self, data, model):
data, model = Stat._validate_inputs(self, data, model)
self._check_background_subtraction(data)
return data, model
[docs] def calc_stat(self, data, model):
fitdata, modeldata = self._get_fit_model_data(data, model)
return self._calc(fitdata[0], modeldata, None,
truncation_value)
# DOC-TODO: where is the truncate/trunc_value stored for objects
# AHA: it appears to be taken straight from the config
# file rather than associated with the Stat class
# DOC-TODO: where to talk about the .sherpa.rc config file?
[docs]class Cash(Likelihood):
"""Poisson Log-likelihood function.
Counts are sampled from the Poisson distribution, and so the best
way to assess the quality of model fits is to use the product of
individual Poisson probabilities computed in each bin i, or the
likelihood L:
L = (product)_i [ M(i)^D(i)/D(i)! ] * exp[-M(i)]
where M(i) = S(i) + B(i) is the sum of source and background model
amplitudes, and D(i) is the number of observed counts, in bin i.
The Cash statistic [1]_ is derived by (1) taking the logarithm of
the likelihood function, (2) changing its sign, (3) dropping the
factorial term (which remains constant during fits to the same
dataset), and (4) multiplying by two:
C = 2 * (sum)_i [ M(i) - D(i) log M(i) ]
The factor of two exists so that the change in cash statistic from
one model fit to the next, (Delta)C, is distributed approximately
as (Delta)chi-square when the number of counts in each bin is
high. One can then in principle use (Delta)C instead of
(Delta)chi-square in certain model comparison tests. However,
unlike chi-square, the cash statistic may be used regardless of
the number of counts in each bin.
The magnitude of the Cash statistic depends upon the number of
bins included in the fit and the values of the data
themselves. Hence one cannot analytically assign a
goodness-of-fit measure to a given value of the Cash statistic.
Such a measure can, in principle, be computed by performing
Monte Carlo simulations. One would repeatedly sample new
datasets from the best-fit model, fit them, and note where the
observed Cash statistic lies within the derived distribution
of Cash statistics. Alternatively, the `cstat` statistic can
be used.
Notes
-----
The background should not be subtracted from the data when this
statistic is used. It should be modeled simultaneously with the
source.
The Cash statistic function evaluates the logarithm of each data
point. If the number of counts is zero or negative, it's not
possible to take the log of that number. The behavior in this case
is controlled by the `truncate` and `trunc_value` settings in the
.sherpa.rc file:
- if `truncate` is `True` (the default value), then
`log(trunc_value)` is used whenever the data value is <= 0. The
default is `trunc_value=1.0e-25`.
- when `truncate` is `False` an error is raised.
References
----------
.. [1] "Parameter estimation in astronomy through application of
the likelihood ratio", Cash, W. 1979, ApJ 228, 939
http://adsabs.harvard.edu/abs/1979ApJ...228..939C
"""
_calc = _statfcts.calc_cash_stat
def __init__(self, name='cash'):
Likelihood.__init__(self, name)
[docs]class CStat(Likelihood):
"""Poisson Log-likelihood function (XSPEC style).
This is equivalent to the XSPEC implementation of the
Cash statistic [1]_ except that it requires a model to be fit
to the background. To handle the background in the same manner
as XSPEC, use the WStat statistic.
Counts are sampled from the Poisson distribution, and so the best
way to assess the quality of model fits is to use the product of
individual Poisson probabilities computed in each bin i, or the
likelihood L:
L = (product)_i [ M(i)^D(i)/D(i)! ] * exp[-M(i)]
where M(i) = S(i) + B(i) is the sum of source and background model
amplitudes, and D(i) is the number of observed counts, in bin i.
The cstat statistic is derived by (1) taking the logarithm of the
likelihood function, (2) changing its sign, (3) dropping the
factorial term (which remains constant during fits to the same
dataset), (4) adding an extra data-dependent term (this is what
makes it different to `Cash`, and (5) multiplying by two:
C = 2 * (sum)_i [ M(i) - D(i) + D(i)*[log D(i) - log M(i)] ]
The factor of two exists so that the change in the cstat statistic
from one model fit to the next, (Delta)C, is distributed
approximately as (Delta)chi-square when the number of counts in
each bin is high. One can then in principle use (Delta)C instead
of (Delta)chi-square in certain model comparison tests. However,
unlike chi-square, the cstat statistic may be used regardless of
the number of counts in each bin.
The inclusion of the data term in the expression means that,
unlike the Cash statistic, one can assign an approximate
goodness-of-fit measure to a given value of the cstat statistic,
i.e. the observed statistic, divided by the number of degrees of
freedom, should be of order 1 for good fits.
Notes
-----
The background should not be subtracted from the data when this
statistic is used. It should be modeled simultaneously with the
source.
The cstat statistic function evaluates the logarithm of each data
point. If the number of counts is zero or negative, it's not
possible to take the log of that number. The behavior in this case
is controlled by the `truncate` and `trunc_value` settings in the
.sherpa.rc file:
- if `truncate` is `True` (the default value), then
`log(trunc_value)` is used whenever the data value is <= 0. The
default is `trunc_value=1.0e-25`.
- when `truncate` is `False` an error is raised.
References
----------
.. [1] The description of the Cash statistic (`cstat`) in
https://heasarc.gsfc.nasa.gov/xanadu/xspec/manual/XSappendixStatistics.html
"""
_calc = _statfcts.calc_cstat_stat
_can_calculate_rstat = True
def __init__(self, name='cstat'):
Likelihood.__init__(self, name)
[docs]class Chi2(Stat):
"""A Gaussian Log-likelihood function.
It is assumed that the counts are sampled from the Gaussian
(Normal) distribution and so the best way to assess the quality of
model fit is to use the product of individual Gaussian probabilities
computed in each bin i, or the likelihood:
L = (prod)_i 1/(sigma^2 sqrt(2 pi)) exp[(N(i) - M(i))^2/2 sigma(i)^2]
where M(i) = S(i) + B(i) is the sum of source and background model
amplitudes, and N(i) is the total number of observed counts in bin i.
The chi-square statistic is:
chi^2 = (sum)_i [ [ N(i,S) - B(i,x,pB) - S(i,x,pS) ]^2 / sigma(i)^2 ]
where N(i,S) is the total number of observed counts in bin i of
the on-source region; B(i,x,pB) is the number of predicted
background model counts in bin i of the on-source region (zero for
background-subtracted data), rescaled from bin i of the off-source
region, and computed as a function of the model argument x(i)
(e.g., energy or time) and set of background model parameter
values pB; S(i,x,pS) is the number of predicted source model
counts in bin i, as a function of the model argument x(i) and set
of source model parameter values pS; and sigma(i) is the error in
bin i.
Note that there are several weightings of this statistics depending
on calculation of sigma(i). N(i,S) contains the background counts and
in a case of background subtraction the number of contributing
background counts needs to be estimated from the background, so an
off-source region. In such case, N(i,B) is the total number of observed
counts in bin i of the off-source region; A(B) is the off-source "area",
which could be the size of the region from which the background is extracted, or
the length of a background time segment, or a product of the two,
etc.; and A(S) is the on-source "area". These terms may be defined
for a particular type of data: for example, PHA data sets A(B) to
`BACKSCAL * EXPOSURE` from the background data set and A(S) to
`BACKSCAL * EXPOSURE` from the source data set.
There are different ways of defining the sigma(i) terms,
supported by the sub-classes.
Notes
-----
It is assumed that there is a one-to-one mapping between a given
background region bin and a given source region bin. For
instance, in the analysis of PHA data, it is assumed that the
input background counts spectrum is binned in exactly the same way
as the input source counts spectrum, and any filter applied to the
source spectrum automatically applied to the background spectrum.
This means that the user cannot, for example, specify arbitrary
background and source regions in two dimensions and get correct
results. This limitation *only* applies to backgrounds included
included as part of the data set - e.g. as with PHA files - and
can be avoided by treating the background as a separate data set.
"""
_calc = _statfcts.calc_chi2_stat
_can_calculate_rstat = True
def __init__(self, name='chi2'):
Stat.__init__(self, name)
[docs] @staticmethod
def calc_staterror(data):
raise StatErr('chi2noerr')
[docs] def calc_stat(self, data, model):
fitdata, modeldata = self._get_fit_model_data(data, model)
return self._calc(fitdata[0], modeldata,
fitdata[1], fitdata[2],
None, # TODO: weights
truncation_value)
[docs] def calc_chisqr(self, data, model):
"""Return the chi-square value for each bin.
Parameters
----------
data : `sherpa.data.Data` or `sherpa.data.DataSimulFit`
The data set, or sets, to use.
model : `sherpa.models.model.Model` or `sherpa.models.model.SimulFitModel`
The model expression, or expressions. If a
`sherpa.models.model.SimulFitModel`
is given then it must match the number of data sets in the
data parameter.
Returns
-------
chisqr : array of numbers
The per-bin chi-square values.
"""
_, fvec = self.calc_stat(data, model)
return fvec * fvec
[docs]class LeastSq(Chi2):
"""Least Squared Statistic.
The least-square statistic is equivalent to a chi-square
statistic where the error on each point - sigma(i) - is 1.
"""
_calc = _statfcts.calc_lsq_stat
_can_calculate_rstat = False
def __init__(self, name='leastsq'):
Stat.__init__(self, name)
[docs] @staticmethod
def calc_staterror(data):
return numpy.ones_like(data)
[docs]class Chi2Gehrels(Chi2):
"""Chi Squared with Gehrels variance.
The variance is estimated from the number of counts in each bin,
but unlike `Chi2DataVar`, the Gaussian approximation is not
used. This makes it more-suitable for use with low-count data.
The standard deviation for each bin is calculated using the
approximation from [1]_:
sigma(i,S) = 1 + sqrt(N(i,s) + 0.75)
where the higher-order terms have been dropped. This is accurate
to approximately one percent. For data where the background has
not been subtracted then the error term is:
sigma(i) = sigma(i,S)
whereas with background subtraction,
sigma(i)^2 = sigma(i,S)^2 + [A(S)/A(B)]^2 sigma(i,B)^2
A(B) is the off-source "area", which could be
the size of the region from which the background is extracted, or
the length of a background time segment, or a product of the two,
etc.; and A(S) is the on-source "area". These terms may be defined
for a particular type of data: for example, PHA data sets A(B) to
`BACKSCAL * EXPOSURE` from the background data set and A(S) to
`BACKSCAL * EXPOSURE` from the source data set.
See Also
--------
Chi2DataVar, Chi2ModVar, Chi2XspecVar
Notes
-----
The accuracy of the error term when the background has been
subtracted has not been determined. A preferable approach to
background subtraction is to model the background as well as the
source signal.
References
----------
.. [1] "Confidence limits for small numbers of events in
astrophysical data", Gehrels, N. 1986, ApJ, vol 303,
p. 336-346.
http://adsabs.harvard.edu/abs/1986ApJ...303..336G
"""
def __init__(self, name='chi2gehrels'):
Chi2.__init__(self, name)
[docs] @staticmethod
def calc_staterror(data):
return _statfcts.calc_chi2gehrels_errors(data)
[docs]class Chi2ConstVar(Chi2):
"""Chi Squared with constant variance.
The variance is the same in each bin, and set to be the mean
number of counts in the data:
sigma(i)^2 = (1/K) * (sum)_(j=1)^K N(j,S) + [A(S)/A(B)]^2 N(j,B)
where K is the number of on-source (and off-source) bins included
in the fit. The background term appears only if an estimate of the
background has been subtracted from the data.
"""
def __init__(self, name='chi2constvar'):
Chi2.__init__(self, name)
[docs] @staticmethod
def calc_staterror(data):
return _statfcts.calc_chi2constvar_errors(data)
[docs]class Chi2DataVar(Chi2):
"""Chi Squared with data variance.
The variance in each bin is estimated from the data value in that
bin.
If the number of counts in each bin is large, then the shape of
the Poisson distribution from which the counts are sampled tends
asymptotically towards that of a Gaussian distribution, with
variance
sigma(i)^2 = N(i,S) + [A(S)/A(B)]^2 N(i,B)
where N is the number of on-source (and off-source) bins included
in the fit. The background term appears only if an estimate of the
background has been subtracted from the data.
A(B) is the off-source "area", which could be
the size of the region from which the background is extracted, or
the length of a background time segment, or a product of the two,
etc.; and A(S) is the on-source "area". These terms may be defined
for a particular type of data: for example, PHA data sets A(B) to
`BACKSCAL * EXPOSURE` from the background data set and A(S) to
`BACKSCAL * EXPOSURE` from the source data set.
See Also
--------
Chi2Gehrels, Chi2ModVar, Chi2XspecVar
"""
def __init__(self, name='chi2datavar'):
Chi2.__init__(self, name)
[docs] @staticmethod
def calc_staterror(data):
return _statfcts.calc_chi2datavar_errors(data)
[docs]class Chi2ModVar(Chi2):
"""Chi Squared with model amplitude variance.
The variance in each bin is estimated from the *model* value in
that bin. This contrasts with other Chi-squared statics - such
as `Chi2DataVar`, Chi2XspecVar`, and `Chi2Gehrels` - which use
the data values. The variance is
sigma(i)^2 = S(i) + [A(S)/A(B)]^2 B(i,off)
where B(i,off) is the background model amplitude in bin i of the
off-source region.
A(B) is the off-source "area", which could be
the size of the region from which the background is extracted, or
the length of a background time segment, or a product of the two,
etc.; and A(S) is the on-source "area". These terms may be defined
for a particular type of data: for example, PHA data sets A(B) to
`BACKSCAL * EXPOSURE` from the background data set and A(S) to
`BACKSCAL * EXPOSURE` from the source data set.
See Also
--------
Chi2DataVar, Chi2Gehrels, Chi2XspecVar
Notes
-----
The background should not be subtracted from the data when this
statistic is used, as it underestimates the variance when fitting
background-subtracted data.
"""
_calc = _statfcts.calc_chi2modvar_stat
def __init__(self, name='chi2modvar'):
Chi2.__init__(self, name)
# Statistical errors are not used
[docs] @staticmethod
def calc_staterror(data):
return numpy.zeros_like(data)
[docs]class Chi2XspecVar(Chi2):
"""Chi Squared with data variance (XSPEC style).
The variance in each bin is estimated from the data value in that
bin.
The calculation of the variance is the same as `Chi2DataVar`
except that if the number of counts in a bin is less than 1 then
the variance for that bin is set to 1. Note that this does not
handle background subtraction, so to match XSPEC extra logic is
present in sherpa.astro.data.DataPHA.get_staterror to handle that
case.
See Also
--------
Chi2DataVar, Chi2Gehrels, Chi2ModVar
"""
def __init__(self, name='chi2xspecvar'):
Chi2.__init__(self, name)
[docs] @staticmethod
def calc_staterror(data):
return _statfcts.calc_chi2xspecvar_errors(data)
[docs]class UserStat(Stat):
"""Support simple user-supplied statistic calculations.
Notes
-----
This class is used by the `sherpa.ui.load_user_stat`
to provide a user-definable statistic calculation as
a function. For more complicated cases it is suggested that
users should write their own class instead of using
this one.
"""
def __init__(self, statfunc=None, errfunc=None, name='userstat'):
self._statfuncset = False
self.statfunc = (lambda x: None)
self._staterrfuncset = False
self.errfunc = (lambda x: None)
if statfunc is not None:
self.statfunc = statfunc
self._calc = statfunc
self._statfuncset = True
if errfunc is not None:
self.errfunc = errfunc
self._staterrfuncset = True
Stat.__init__(self, name)
def __getstate__(self):
state = self.__dict__.copy()
# Function pointers to methods of the class
# (of type 'instancemethod') are NOT picklable
# remove them and restore later with a coord init
del state['statfunc']
del state['errfunc']
return state
def __setstate__(self, state):
# Populate the function pointers we deleted at pickle time with
# no-ops.
self.__dict__['statfunc'] = (lambda x: None)
self.__dict__['errfunc'] = (lambda x: None)
self.__dict__.update(state)
[docs] def set_statfunc(self, func):
self.statfunc = func
self._statfuncset = True
[docs] def set_errfunc(self, func):
self.errfunc = func
self._staterrfuncset = True
[docs] def calc_staterror(self, data):
if not self._staterrfuncset:
raise StatErr('nostat', self.name, 'calc_staterror()')
return self.errfunc(data)
[docs] def calc_stat(self, data, model):
if not self._statfuncset:
raise StatErr('nostat', self.name, 'calc_stat()')
fitdata, modeldata = self._get_fit_model_data(data, model)
return self.statfunc(fitdata[0],
modeldata,
staterror=fitdata[1],
syserror=fitdata[2],
weight=None) # TODO weights
[docs]class WStat(Likelihood):
"""Poisson Log-likelihood function including background (XSPEC style).
This is equivalent to the XSPEC implementation of the
W statistic for CStat [1]_, and includes the background data in
the fit statistic. If a model is being fit to the background then
the CStat statistic should be used.
The following description is taken from [1]_.
Suppose that each bin in the background spectrum is given its own
parameter so that the background model is b_i = f_i. A standard fit
for all these parameters would be impractical; however there is an
analytical solution for the best-fit f_i in terms of the other
variables which can be derived by using the fact that the derivative
of the likelihood (L) will be zero at the best fit. Solving for the
f_i and substituting gives the profile likelihood::
W = 2 sum_(i=1)^N t_s m_i + (t_s + t_b) f_i -
S_i ln(t_s m_i + t_s f_i) - B_i ln(t_b f_i) -
S_i (1- ln(S_i)) - B_i (1 - ln(B_i))
where::
f_i = (S_i + B_i - (t_s + t_b) m_i + d_i) / (2 (t_s + t_b))
d_i = sqrt([(t_s + t_b) m_i - S_i - B_i]^2 +
4(t_s + t_b) B_i m_i)
If any bin has S_i and/or B_i zero then its contribution to W (W_i)
is calculated as a special case. So, if S_i is zero then::
W_i = t_s m_i - B_i ln(t_b / (t_s + t_b))
If B_i is zero then there are two special cases. If
m_i < S_i / (t_s + t_b) then::
W_i = - t_b m_i - S_i ln(t_s / (t_s + t_b))
otherwise::
W_i = t_s m_i + S_i (ln(S_i) - ln(t_s m_i) - 1)
In practice, it works well for many cases but for weak sources can
generate an obviously wrong best fit. It is not clear why this happens
although binning to ensure that every bin contains at least one count
often seems to fix the problem. In the limit of large numbers of counts
per spectrum bin a second-order Taylor expansion shows that W tends to::
sum_(i=1)^N ( [S_i - t_s m_i - t_s f_i]^2 / (t_s (m_i + f_i)) +
[B_i - t_b f_i]^2 / (t_b f_i) )
which is distributed as chi^2 with N - M degrees of freedom, where
the model m_i has M parameters (include the normalization).
References
----------
.. [1] The description of the W statistic (`wstat`) in
https://heasarc.gsfc.nasa.gov/xanadu/xspec/manual/XSappendixStatistics.html
"""
_calc = _statfcts.calc_wstat_stat
_can_calculate_rstat = True
def __init__(self, name='wstat'):
Likelihood.__init__(self, name)
[docs] def calc_stat(self, data, model):
data, model = self._validate_inputs(data, model)
# Need access to backscal values and background data filtered
# and grouped in the same manner as the data. There is no
# easy access to this via the Data API (in part because the
# Data class has no knowledge of grouping or backscale values).
#
# An alternative approach would be to just calculate the
# statistic for each dataset individually and then
# sum the statistics for the return value, but the
# original code used this approach.
#
data_src = []
data_model = data.eval_model_to_fit(model)
data_bkg = []
nelems = []
exp_src = []
exp_bkg = []
backscales = []
# Why are we looping over model.parts as it isn't used?
# Is it just a way to restrict to only use those
# datasets for which we have a model?
#
for dset, mexpr in zip(data.datasets, model.parts):
y = dset.to_fit(staterrfunc=None)[0]
data_src.append(y)
nelems.append(y.size)
try:
bids = dset.background_ids
except AttributeError:
raise StatErr('usecstat') from None
nbkg = len(bids)
if nbkg == 0:
raise StatErr('usecstat')
elif nbkg > 1:
# TODO: improve warning
warnings.warn("Only using first background component for data set {}".format(dset.name))
bid = bids[0]
bset = dset.get_background(bid)
# TODO: the following should be reviewed to see what
# happens if optional information is missing (e.g. if
# BACKSCAL is not set we should default to all 1's,
# but does this code handle it?)
#
data_bkg.append(dset.apply_filter(bset.get_dep(False),
groupfunc=numpy.sum))
# The assumption is that the source and background datasets
# have the same number of channels (before any grouping or
# filtering is applied).
#
# Since the backscal values can be a scalar or array, it is
# easiest just to convert everything to an array.
#
dummy = numpy.ones(dset.get_dep(False).size)
# Combine the BACKSCAL values (use the default _middle
# scheme as this is used elsewhere when combining
# BACKSCAL values; perhaps there should be an API call
# for this?).
#
src_backscal = dset.apply_filter(dset.backscal * dummy,
groupfunc=dset._middle)
bkg_backscal = dset.apply_filter(bset.backscal * dummy,
groupfunc=dset._middle)
backscales.append(bkg_backscal / src_backscal)
# The AREASCAL values are applied to the exposure
# times, since this is how XSPEC handles this (at
# least that's my undertanding of a conversation with
# Keith Arnaud, for XSPEC ~ version 12.9). This requires
# turning an exposure into an array if there's no
# AREASCAl value
#
# For now we follow the same approach as the BACKSCAL
# values if the data is grouped.
#
#
if dset.areascal is None:
ascal = dummy[:dset.get_dep(True).size]
else:
ascal = dset.apply_filter(dset.areascal * dummy,
groupfunc=dset._middle)
exp_src.append(dset.exposure * ascal)
if bset.areascal is None:
ascal = dummy[:dset.get_dep(True).size]
else:
ascal = dset.apply_filter(bset.areascal * dummy,
groupfunc=dset._middle)
exp_bkg.append(bset.exposure * ascal)
data_src = numpy.concatenate(data_src)
exp_src = numpy.concatenate(exp_src)
exp_bkg = numpy.concatenate(exp_bkg)
data_bkg = numpy.concatenate(data_bkg)
backscales = numpy.concatenate(backscales)
return self._calc(data_src, data_model, nelems,
exp_src, exp_bkg,
data_bkg, backscales,
truncation_value)