In [1]:
from threeML import *
from threeML.analysis_results import *
from threeML.io.progress_bar import progress_bar
%matplotlib notebook
import matplotlib.pyplot as plt
In [2]:
# These are the same simulated dataset we use in the test of the XY plugin
x = np.linspace(0, 10, 50)
poiss_sig = [44, 43, 38, 25, 51, 37, 46, 47, 55, 36, 40, 32, 46, 37, 44, 42, 50, 48, 52, 47, 39, 55, 80, 93, 123, 135,
96, 74, 43, 49, 43, 51, 27, 32, 35, 42, 43, 49, 38, 43, 59, 54, 50, 40, 50, 57, 55, 47, 38, 64]
In [4]:
# A simple test of the likelihood minimization
y = np.array(poiss_sig)
xy = XYLike("test", x, y, poisson_data=True)
fitfun = Line() + Gaussian()
fitfun.a_1.bounds = (-10, 10.0)
fitfun.b_1.bounds = (-100, 100.0)
fitfun.F_2 = 60.0
fitfun.F_2.bounds = (1e-3, 200.0)
fitfun.mu_2 = 5.0
fitfun.mu_2.bounds = (0.0, 100.0)
fitfun.sigma_2.bounds = (1e-3, 10.0)
model = Model(PointSource('fake',0.0, 0.0, fitfun))
data = DataList(xy)
jl = JointLikelihood(model, DataList(xy))
_ = jl.fit()
In [5]:
# This functionality is still there, although it cannot be propagated
# (error propagation assume Gaussian errors, i.e., symmetric errors)
# In this case though errors are pretty symmetric, so we are likely in the case
# where the MLE is actually normally distributed
_ = jl.get_errors()
In [6]:
# This is the novelty. This access the new class of the results
ar = jl.results
In [7]:
ar.write_to("test_mle.fits", overwrite=True)
In [8]:
ar_reloaded = load_analysis_results("test_mle.fits")
In [9]:
ar_reloaded.get_statistic_frame()
Out[9]:
In [10]:
# You can get a DataFrame with the results
ar_reloaded.get_data_frame()
Out[10]:
In [10]:
from threeML.analysis_results import AnalysisResultsSet
analysis_set = AnalysisResultsSet([ar, ar_reloaded])
analysis_set.set_bins("testing", [-1, 1], [3, 5], unit = 's')
analysis_set.write_to("analysis_set_test.fits", overwrite=True)
In [11]:
analysis_set = load_analysis_results("analysis_set_test.fits")
In [12]:
analysis_set[0].display()
In [13]:
# or display again the results. However, this time I choose to display
# hpd intervals
ar.display("hpd")
In [14]:
# You can use the results for propagating errors non-linearly for analytical functions
p1 = ar.get_variates("fake.spectrum.main.composite.a_1")
p2 = ar.get_variates("fake.spectrum.main.composite.b_1")
print("Propagating a+b, with a and b respectively:")
print(p1)
print(p2)
print("\nThis is the result (with errors):")
res = p1 + p2
print(res)
print res.equal_tail_confidence_interval()
# The propagation accounts for covariances, so for example this
# has error of zero (of course) since there is perfect covariance
print("\nThis is 50 * a/a:")
print(50 * p1/p1)
# You can use arbitrary (np) functions
print("\nThis is arcsinh(a + 5*b) / np.log10(b) (why not?)")
print(np.arcsinh(p1 + 5*p2) / np.log10(p2))
# Errors can become asymmetric. For example, the ratio of two gaussians is
# asymmetric notoriously:
print("\nRatio a/b:")
print(p2/p1)
# You can always use it with arbitrary functions like:
def my_function(x, a, b):
return b*x**a
print("\nPropagating using a custom function:")
print(my_function(2.3, p1, p2))
In [15]:
# This is an example of an error propagation to get the plot of the model with its errors
# (which are propagated without assuming linearity on parameters)
def go(fitfun, ar, model):
# Gather the parameter variates
arguments = {}
for par in fitfun.parameters.values():
if par.free:
this_name = par.name
this_variate = ar.get_variates(par.path)
# Do not use more than 1000 values (would make computation too slow for nothing)
if len(this_variate) > 1000:
this_variate = np.random.choice(this_variate, size=1000)
arguments[this_name] = this_variate
# Prepare the error propagator function
pp = ar.propagate(ar.optimized_model.fake.spectrum.main.shape.evaluate_at, **arguments)
# You can just use it as:
print(pp(5.0))
#Make the plot
energies = np.linspace(0, 10, 100)
low_curve = np.zeros_like(energies)
middle_curve = np.zeros_like(energies)
hi_curve = np.zeros_like(energies)
free_parameters = model.free_parameters
with progress_bar(len(energies), title="Propagating errors") as p:
with use_astromodels_memoization(False):
for i, e in enumerate(energies):
this_flux = pp(e)
low_bound, hi_bound = this_flux.equal_tail_confidence_interval()
low_curve[i], middle_curve[i], hi_curve[i] = (low_bound, this_flux.median, hi_bound)
p.increase()
plt.plot(energies, middle_curve, '--', color='black')
plt.fill_between(energies, low_curve, hi_curve, alpha=0.5, color='blue')
go(fitfun, ar, model)
In [16]:
# Exactly the same can be done for a Bayesian analysis
# Let's run it first
for parameter in ar.optimized_model:
model[parameter.path].value = parameter.value
model.fake.spectrum.main.composite.a_1.set_uninformative_prior(Uniform_prior)
model.fake.spectrum.main.composite.b_1.set_uninformative_prior(Uniform_prior)
model.fake.spectrum.main.composite.F_2.set_uninformative_prior(Log_uniform_prior)
model.fake.spectrum.main.composite.mu_2.set_uninformative_prior(Uniform_prior)
model.fake.spectrum.main.composite.sigma_2.set_uninformative_prior(Log_uniform_prior)
bs = BayesianAnalysis(model, data)
samples = bs.sample(20, 100, 1000)
In [17]:
ar2 = bs.results
In [18]:
ar2.write_to("test_bayes.fits", overwrite=True)
In [19]:
ar2_reloaded = load_analysis_results("test_bayes.fits")
In [20]:
np.allclose(ar2_reloaded.samples, ar2.samples)
Out[20]:
In [21]:
ar2.get_data_frame("equal tail")
Out[21]:
In [22]:
p1 = ar2.get_variates("fake.spectrum.main.composite.a_1")
p2 = ar2.get_variates("fake.spectrum.main.composite.b_1")
print(p1)
print(p2)
res = p1 + p2
print(res)
In [23]:
go(fitfun, ar2, model)
In [ ]: