WStat

class sherpa.stats.WStat(name='wstat')

Bases: sherpa.stats.Likelihood

Maximum 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]	(1, 2) The description of the W statistic (wstat) in https://heasarc.gsfc.nasa.gov/xanadu/xspec/manual/XSappendixStatistics.html

Methods Summary

calc_stat(data, model)
calc_staterror(data)

Methods Documentation

calc_stat(data, model)

calc_staterror(data)