Fitting

Fast Norm Fit

class bctools.fitting.FastNormFit(max_iter=1000, conv_frac_tol=0.001, zero_ts_tol=1e-05, allow_negative=False)[source]

Bases: object

Perform a fast poisson maximum likelihood ratio fit of the normalization of a source over background.

The likelihood ratio as a function of the norm is computed as follow

\[TS(N) = 2 \sum_i \left( \frac{\log P(d_i; b_i+N e_i)}{\log P(d_i; b_i)} \right)\]

where \(P(d; \lambda)\) is the Poisson probability of obtaining \(d\) count where \(\lambda\) is expected on average; \(b\) is the estimated number of background counts; \(N\) is the normalization; and \(e\) is the expected excess -i.e. signal- per normalization unit -i.e. the number of excess counts equals \(N\).

It can be shown that \(TS(N)\) has analytic derivative of arbitrary order and that the Newton’s method is guaranteed to converge if initialized at \(N=0\).

Note

The background is not allowed to float. It is assumed the error on the estimation of the background is small compared to the fluctuation of the background itself (i.e. \(N_{B}/N_{off} \lll 1\)).

Note

Because of the Poisson probability, \(TS(N)\) is only well-defined for \(N \geq 1\). By default, underfluctuations are set to \(TS(N=0) = 0\). For cases when there is benefit in letting the normalization float to negative values, you can use allow_negative, but in that case the results are only valid in the Gaussian regime.

Parameters:

max_iter (int) – Maximum number of iteration
conv_frac_tol (float) – Convergence stopping condition, expressed as the ratio between the size of the last step and the current norm value.
zero_ts_tol (float) – If zero_ts_tol < TS < 0, then TS is set to 0 without failed flag status (analytically, TS < 0 should never happen)
allow_negative (bool) – Allow the normalization to float toward negative values

static dts(data, bkg, unit_excess, norm, order=1)[source]

Get the derivative of TS with respecto to the normalization, at given normalization.

Parameters:

data (array) – Observed counts
bkg (array) – Background estimation. Same size as data. Every element should be >0
unit_excess (array) – Expected excess per normalization unit. Same size as data.
norm (float or array) – Normalization value
order (int) – Derivative order

Returns:

d^n TS / dN^n, same shape as norm

Return type:

float or array

solve(data, bkg, unit_excess)[source]

Get the maximum TS, fitted normalization and normalization error (\(\Delta TS = 1\))

Note

The normalization error is obtained from approximating the TS function as as parabola (i.e. valid in the Gaussian regime). TS and norm are indeed valid in the Poisson regime.

Parameters:

data (array) – Observed counts
bkg (array) – Background estimation. Same size as data. Every element should be >0
unit_excess (array) – Expected excess per normalization unit. Same size as data.

Returns:

ts, norm, norm error and status (0 = good)

Return type:

(float, float, float, bool)

static ts(data, bkg, unit_excess, norm)[source]

Get TS for a given normalization.

Parameters:

data (array) – Observed counts
bkg (array) – Background estimation. Same size as data. Every element should be >0
unit_excess (array) – Expected excess per normalization unit. Same size as data.
norm (float or array) – Normalization value

Returns:

TS, same shape as norm

Return type:

float or array