Some essential notes on C-statistics

Basic concepts

C-statistics is based on Poisson likelihood:

$$L(M)=\prod _i\frac{M_i^{D_i}}{D_i!}\exp (-M_i)$$

where $L(M)$ is the likelihood for the model $M$, $D_i$ is the data in $i$th bin, $M_i$ denotes the model prediction in $i$th bin. Taking its logarithm and multiplying by $-2$, we can get

$$-2\ln L(M)=2\sum _i(M_i-D_i\ln M_i+\ln (D_i!)).$$

Omitting the factorial term, we can get the Cash-statistics (Cash 1979):

$$\tilde{C}=2\sum _i(M_i-D_i\ln M_i).$$

Approximating the factorial term by Stirling's formula, that is

$$\ln (D_i!)\approx D_i\ln D_i-D_i,$$

a modification of the original Cash-statistic, C-statistics, can be obtained as follows:

$$C=2\sum _i(M_i-D_i\ln M_i+D_i\ln D_i-D_i),$$

which is implemented in some popular fitting packages like XSPEC (Arnaud 1996), SHERPA (Freeman et al. 2001), and SPEX (Kaastra et al. 1996). $\tilde{C}$ is the same as $C$, up to a constant $\sum _i(D_i\ln D_i-D_i)$. $C$ is non-negative, $C$ is equal to $0$ if and only if all the $M_i$ are equal to $D_i$. Since the count rate is usually low for X-ray observation, it is better to use C-statistics fitting than ${\chi }^2$ fitting, i.e., getting the best-fit parameters by minimizing $C$ instead of ${\chi }^2$.

$1\sigma$ confidence intervals for parameters

Assuming that there are $k$ parameters in the model, that is $p_1$, $p_2$, ..., $p_k$, a set of $\{p_i\}$ ($i$ ranges from $1$ to $k$) that results in a minimum value of $C$ is the best-fit parameter set, i.e., the best-fit model. However, merely getting the best-fit parameter set is not enough, we also need to evaluate the $1\sigma$ confidence intervals for the parameters.

A simple schematic

For simplicity, we may as well take a model with just two parameters $p_1$ and $p_2$ as an example. $\{p_{1,\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}},p_{2,\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}\}$ is the best-fit parameter set, which generates the minimum $C$, i.e., $C_{\mathrm{m}\mathrm{i}\mathrm{n}}$. Considering the 2D parameter space $(p_1,p_2)$, $C$ will be larger than $C_{\mathrm{m}\mathrm{i}\mathrm{n}}$ in the neighborhood of $(p_{1,\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}},p_{2,\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}})$. $C\le C_{\mathrm{m}\mathrm{i}\mathrm{n}}+1$ defines a close region in the 2D parameter space, which is called the $1\sigma$ confidence region. $C=C_{\mathrm{m}\mathrm{i}\mathrm{n}}+1$ defines the corresponding boundary of the $1\sigma$ confidence region. The figure below can serve as a schematic of the $1\sigma$ confidence region in the 2D parameter space $(p_1,p_2)$, which assumes that $p_{1,\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}=0$ and $p_{2,\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}=0$, or equivalently, the origin of the 2D parameter space is set to be the best-fit parameter set.

The line segment AD represents the $1\sigma$ confidence interval for $p_1$. The line segment BC is not the $1\sigma$ confidence interval for $p_1$. If we fix the parameter $p_2=p_{2,\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}=0$ and calculate the error of $p_1$ in SPEX, we will get the $1\sigma$ confidence interval for $p_1$, which is the line segment BC. However, it is not reasonable, because there will always be some correlation between parameters. You can also find that the $1\sigma$ confidence interval for $p_1$ obtained by setting the rest parameters to the best-fit value will always be smaller than the $1\sigma$ confidence interval for $p_1$ obtained by freeing all the parameters (BC<AD).

The $1\sigma$ confidence region in the 2D parameter space $(p_1,p_2)$ is represented by the oblique ellipse. Points A and D denote the minimum and maximum value of $p_1$ of the oblique ellipse, respectively. Points B and C are two points of intersection between the oblique ellipse and the horizontal line $p_2=p_{2,\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}=0$.

The algorithm of finding the $1\sigma$ confidence interval

As shown in the Figure, you may think that we need to get the boundary, namely, the oblique ellipse first and then get the $1\sigma$ confidence interval for each parameter by projecting the oblique ellipse to each parameter axis (by getting the minimum and maximum value of each parameter). However, calculating the boundary of the oblique ellipse is not efficient. To calculate the $1\sigma$ confidence interval for the parameter $p_1$, the algorithm implementd in SPEX goes as follows,

• 1. Get the best-fit parameter set
• 2. Fix the parameter $p_1$ to a value that is close to the $p_{1,\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}$ but free all the rest parameters.
• 3. Obtain the minimum $C$
• 4. If the minimum $C$ is larger than $C_{\mathrm{m}\mathrm{i}\mathrm{n}}+1$, choose a new $p_1$ that is closer to $p_{1,\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}$, if the minimum $C$ is smaller than $C_{\mathrm{m}\mathrm{i}\mathrm{n}}+1$, choose a new $p_1$ that is farther from $p_{1,\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}$. Go to step 2.
• 5. Return the two points $p_1$ that can result in the minimum $C$ being equal to $C_{\mathrm{m}\mathrm{i}\mathrm{n}}+1$ by setting the rest parameters to be free. Actually, it is not exactly $C_{\mathrm{m}\mathrm{i}\mathrm{n}}+1$, there is some tolerance, the algorithm in SPEX will stop at some value that is close to $C_{\mathrm{m}\mathrm{i}\mathrm{n}}+1$, like $C_{\mathrm{m}\mathrm{i}\mathrm{n}}+1.01$.

Why the $1\sigma$ region is defined as $C_{\mathrm{m}\mathrm{i}\mathrm{n}}+1$?

Why it is $C_{\mathrm{m}\mathrm{i}\mathrm{n}}+1$ not some value like $C_{\mathrm{m}\mathrm{i}\mathrm{n}}+1.2$ or $C_{\mathrm{m}\mathrm{i}\mathrm{n}}+1.3$? It should date back to the ${\chi }^2$ test.

$$\mathrm{\Delta }{\chi }^2={\chi }^2(N)-{\chi }_{\mathrm{m}\mathrm{i}\mathrm{n}}^2(N-1)$$

is distributed like a ${\chi }^2$ variable with one degree of freedom, where N is the number of data points and ${\chi }^2(N-1{)}_{\mathrm{m}\mathrm{i}\mathrm{n}}$ is the ${\chi }^2$ the usual minimum fit statistic Bonamente (2019). Namely, $\mathrm{\Delta }{\chi }^2$ is distributed as ${\chi }^2(1)$. According to the ${\chi }^2(1)$ distribution, $\mathrm{\Delta }{\chi }^2\le 1$ will lead to a probability of 68%, which is the so-called $1\sigma$ confidence level. For C-statistics,

$$\mathrm{\Delta }C=C-C_{\mathrm{m}\mathrm{i}\mathrm{n}}$$

is approximately distributed like a ${\chi }^2(1)$ distribution. Therefore, we also use $C_{\mathrm{m}\mathrm{i}\mathrm{n}}+1$ to define the $1\sigma$ confidence region.