New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add flux points upper limit estimation #1084

Merged
merged 4 commits into from Jul 7, 2017
Jump to file or symbol
Failed to load files and symbols.
+138 −11
Diff settings

Always

Just for now

Copy path View file
@@ -345,7 +345,7 @@ def _info_flux_points(self):
flux_points = self.flux_points.table.copy()
energy_cols = ['e_ref', 'e_min', 'e_max']
flux_cols = ['dnde', 'dnde_errn', 'dnde_errp']
flux_cols = ['dnde', 'dnde_errn', 'dnde_errp', 'dnde_ul']
flux_points = flux_points[energy_cols + flux_cols]
for _ in energy_cols:
Copy path View file
@@ -557,7 +557,6 @@ def compute_flux_point(self, energy_group):
global_model=self.model,
energy_range=energy_group.energy_range,
)
energy_ref = self.compute_energy_ref(energy_group)
return self.fit_point(
@@ -617,9 +616,110 @@ def compute_approx_model(global_model, energy_range):
par.frozen = True
return approx_model
def fit_point(self, model, energy_group, energy_ref):
def compute_flux_point_ul(self, fit, best_fit, delta_ts=4, negative=False):
"""
Compute upper limits for flux point values.
Parameters
----------
fit : `SpectrumFit`
Instance of spectrum fit.
best_fit : `SpectrumFitResult`
Best fit result.
delta_ts : float (4)
Difference in log-likelihood for given confidence interval.
See Example below.
negative : bool
Compute limit in negative direction.
Examples
--------
To compute ~95% confidence upper limits (or 2 sigma) you can use:
from scipy.stats import chi2, norm
sigma = 2
cl = 1 - 2 * norm.sf(sigma) # using two sided p-value
delta_ts = chi2.isf(1 - cl, df=1)
Returns
-------
dnde_ul : `~astropy.units.Quantity`
Flux point upper limit.
"""
from scipy.optimize import brentq
# this is a prototype for fast flux point upper limit
# calculation using brentq
stat_best_fit = best_fit.statval
amplitude = best_fit.model.parameters['amplitude'].value / 1E-12
amplitude_err = best_fit.model.parameters.error('amplitude') / 1E-12
if negative:
amplitude_max = amplitude
amplitude_min = amplitude_max - 1E3 * amplitude_err
else:
amplitude_max = amplitude + 1E3 * amplitude_err
amplitude_min = amplitude
def ts_diff(x):
fit.model.parameters['amplitude'].value = x * 1E-12
fit.predict_counts()
fit.calc_statval()
return (stat_best_fit + delta_ts) - fit.total_stat
try:
result = brentq(ts_diff, amplitude_min, amplitude_max,
maxiter=100, rtol=1e-2)
return 1E-12 * result * fit.model.parameters['amplitude'].unit
except (RuntimeError, ValueError):
# Where the root finding fails NaN is set as amplitude
log.debug('Flux point upper limit computation failed.')
return np.nan * fit.model.parameters['amplitude'].unit
def compute_flux_point_sqrt_ts(self, fit, best_fit):
"""
Compute sqrt(TS) for flux point.
Parameters
----------
fit : `SpectrumFit`
Instance of spectrum fit.
best_fit : `SpectrumFitResult`
Best fit result.
Returns
-------
sqrt_ts : float
Sqrt(TS) for flux point.
"""
amplitude = best_fit.model.parameters['amplitude'].value
stat_best_fit = best_fit.statval
fit.model.parameters['amplitude'].value = 0
fit.predict_counts()
fit.calc_statval()
stat_null = fit.total_stat
fit.model.parameters['amplitude'].value = amplitude
ts = np.abs(stat_null - stat_best_fit)
return np.sign(amplitude) * np.sqrt(ts)
def fit_point(self, model, energy_group, energy_ref, sqrt_ts_threshold=1):
from .fit import SpectrumFit
# Set reference and remove min amplitude
model.parameters['reference'].value = energy_ref.to('TeV').value
model.parameters['amplitude'].parmin = None
fit = SpectrumFit(self.obs, model)
erange = energy_group.energy_range
@@ -640,13 +740,23 @@ def fit_point(self, model, energy_group, energy_ref):
e_max = energy_group.energy_range.max
e_min = energy_group.energy_range.min
dnde, dnde_err = res.model.evaluate_error(energy_ref)
sqrt_ts = self.compute_flux_point_sqrt_ts(fit, best_fit=res)
dnde_ul = self.compute_flux_point_ul(fit, best_fit=res)
dnde_errp = self.compute_flux_point_ul(fit, best_fit=res, delta_ts=1.) - dnde
dnde_errn = dnde - self.compute_flux_point_ul(fit, best_fit=res, delta_ts=1., negative=True)
return OrderedDict([
('e_ref', energy_ref),
('e_min', e_min),
('e_max', e_max),
('dnde', dnde.to(DEFAULT_UNIT['dnde'])),
('dnde_err', dnde_err.to(DEFAULT_UNIT['dnde'])),
('dnde_ul', dnde_ul.to(DEFAULT_UNIT['dnde'])),
('is_ul', sqrt_ts < sqrt_ts_threshold),
('sqrt_ts', sqrt_ts),
('dnde_errp', dnde_errp),
('dnde_errn', dnde_errn)
])
@@ -164,10 +164,13 @@ def test_flux_points(config):
),
obs=obs(),
seg=seg(obs()),
dnde=2.7465439050126e-11 * u.Unit('cm-2 s-1 TeV-1'),
dnde_err=4.755502901867284e-12 * u.Unit('cm-2 s-1 TeV-1'),
res=-0.11262182922477647,
res_err=0.1536450758523701,
dnde=2.7465e-11 * u.Unit('cm-2 s-1 TeV-1'),
dnde_err=4.7555e-12 * u.Unit('cm-2 s-1 TeV-1'),
dnde_errn=4.5333e-12 * u.Unit('cm-2 s-1 TeV-1'),
dnde_errp=5.0050e-12 * u.Unit('cm-2 s-1 TeV-1'),
dnde_ul=3.7998e-11 * u.Unit('cm-2 s-1 TeV-1'),
res=-0.1126,
res_err=0.1536,
)
tester = FluxPointTester(config)
@@ -181,6 +184,7 @@ def test_flux_points(config):
class FluxPointTester:
def __init__(self, config):
self.config = config
self.rtol = 0.5E-2 # accuracy of 0.5%
self.setup()
def setup(self):
@@ -217,11 +221,24 @@ def test_values(self):
actual = flux_points.table['dnde'].quantity[0]
desired = self.config['dnde']
assert_quantity_allclose(actual, desired)
assert_quantity_allclose(actual, desired, rtol=self.rtol)
actual = flux_points.table['dnde_err'].quantity[0]
desired = self.config['dnde_err']
assert_quantity_allclose(actual, desired)
assert_quantity_allclose(actual, desired, rtol=self.rtol)
actual = flux_points.table['dnde_ul'].quantity[0]
desired = self.config['dnde_ul']
assert_quantity_allclose(actual, desired, rtol=self.rtol)
actual = flux_points.table['dnde_errn'].quantity[0]
desired = self.config['dnde_errn']
assert_quantity_allclose(actual, desired, rtol=self.rtol)
actual = flux_points.table['dnde_errp'].quantity[0]
desired = self.config['dnde_errp']
assert_quantity_allclose(actual, desired, rtol=self.rtol)
def test_spectrum_result(self):
result = SpectrumResult(
@@ -231,11 +248,11 @@ def test_spectrum_result(self):
actual = result.flux_point_residuals[0][0]
desired = self.config['res']
assert_quantity_allclose(actual, desired)
assert_quantity_allclose(actual, desired, rtol=self.rtol)
actual = result.flux_point_residuals[1][0]
desired = self.config['res_err']
assert_quantity_allclose(actual, desired)
assert_quantity_allclose(actual, desired, rtol=self.rtol)
result.plot(energy_range=[1, 10] * u.TeV)
ProTip! Use n and p to navigate between commits in a pull request.