<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Analysis Results
3ML stores the results of a fit in a container we call an "Analysis Result" (AR). The structure of this object is designed to be useable in a live sense within an active analysis (python script, ipython interactive shell, jupyter notebook) as well as storable as a FITS file for saving results for later.
The structure is nearly the same between MLE and Bayesian analyses in order to make a seamless functionality between all analyses.
Step1:  Let's take a look at what we can do with an AR. First, we will simulate some data.
Step2:  MLE Results
First we will demonstrate how AR's work for an MLE analysis on our synthetic data. As we will see, most of the functionality exists in the Bayesian AR's as well. 
Let's do a simple likelihood maximization of our data and model.
Step3:  We can get our errors as always, but the results cannot be propagated (error propagation assumes 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.
Step4:  We need to get the AnalysisResults object that is created after a fit is performed. The AR object is a member of the JointLikelihood object
Step5:  We can display the results of the analysis. Note, when a fit is performed, the post display is actaully from the internal AR.
Step6:  By default, the equal tail intervals are displayed. We can instead display highest posterior densities (equal in the MLE case)
Step7:  The AR stores several properties from the analysis
Step8:  Saving results to disk
The beauty of the analysis result is that all of this information can be written to disk and restored at a later time. The statistical parameters, best-fit model, etc. can all be recovered.
AR's are stored as a structured FITS file. We write the AR like this
Step9:  The FITS file can be examines with any normal FITS reader.
Step10:  However, to easily pull the results back into the 3ML framework, we use the ${\tt load_analysis_results}$ function
Step11:  You can get a DataFrame with the saved results
Step12:  Analysis Result Sets
When doing time-resolved analysis or analysing a several objects, we can save several AR's is a set. This is achieved with the analysis result set. We can pass an array of AR's to the set and even set up descriptions for the different entries.
Step13:  Error propagation
In 3ML, we propagate errors for MLE reults via sampling of the covariance matrix instead of Taylor exanding around the maximum of the likelihood and computing a jacobain. Thus, we can achieve non-linear error propagation.
You can use the results for propagating errors non-linearly for analytical functions
Step14:  The propagation accounts for covariances. For example this
has error of zero (of course) since there is perfect covariance.
Step15:  You can use arbitrary (np) functions
Step16:  Errors can become asymmetric. For example, the ratio of two gaussians is
asymmetric notoriously
Step17:  You can always use it with arbitrary functions
Step18:  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)
Step19:  Bayesian Analysis Results
Analysis Results work exactly the same under Bayesian analysis. 
Let's run the analysis first.
Step20:  Again, we grab the results from the BayesianAnalysis object
Step21:  We can write and read the results to/from a file
Step22:  The AR holds the posterior samples from the analysis. We can see the saved and live reults are the same
Step23:  NOTE
Step24:  We can return pandas DataFrames with equal tail or HPD results.
Step25:  Error propagation operates the same way. Internally, the process is the same as the MLE results, however, the samples are those of the posterior rather than the (assumed) covariance matrix.
Step26:  To demonstrate how the two objects (MLE and Bayes) are the same, we see that our plotting function written for the MLE result works on our Bayesian results seamlessly.
<ASSISTANT_TASK:>
Python Code:
from threeML import *
from threeML.analysis_results import *
from threeML.io.progress_bar import progress_bar
from jupyterthemes import jtplot
%matplotlib inline
jtplot.style(context="talk", fscale=1, ticks=True, grid=False)
import matplotlib.pyplot as plt
plt.style.use("mike")
import astropy.units as u
Explanation: Analysis Results
3ML stores the results of a fit in a container we call an "Analysis Result" (AR). The structure of this object is designed to be useable in a live sense within an active analysis (python script, ipython interactive shell, jupyter notebook) as well as storable as a FITS file for saving results for later.
The structure is nearly the same between MLE and Bayesian analyses in order to make a seamless functionality between all analyses.
End of explanation
gen_function = Line(a=0.0, b=2) + Gaussian(F=30.0, mu=25.0, sigma=1)
# Generate a dataset using the line and a gaussian.
# constant 20% error
x = np.linspace(0, 50, 50)
xy = XYLike.from_function(
    "sim_data", function=gen_function, x=x, yerr=0.2 * gen_function(x)
)
xy.plot()
Explanation: Let's take a look at what we can do with an AR. First, we will simulate some data.
End of explanation
fitfun = Line() + Gaussian()
fitfun.a_1.bounds = (-10, 10.0)
fitfun.b_1.bounds = (-100, 100.0)
fitfun.F_2 = 25.0
fitfun.F_2.bounds = (1e-3, 200.0)
fitfun.mu_2 = 25.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()
Explanation: MLE Results
First we will demonstrate how AR's work for an MLE analysis on our synthetic data. As we will see, most of the functionality exists in the Bayesian AR's as well. 
Let's do a simple likelihood maximization of our data and model.
End of explanation
jl.get_errors();
Explanation: We can get our errors as always, but the results cannot be propagated (error propagation assumes 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.
End of explanation
ar = jl.results
Explanation: We need to get the AnalysisResults object that is created after a fit is performed. The AR object is a member of the JointLikelihood object
End of explanation
ar.display()
Explanation: We can display the results of the analysis. Note, when a fit is performed, the post display is actaully from the internal AR.
End of explanation
ar.display("hpd")
Explanation: By default, the equal tail intervals are displayed. We can instead display highest posterior densities (equal in the MLE case)
End of explanation
ar.analysis_type
ar.covariance_matrix
ar.get_point_source_flux(1*u.keV, .1*u.MeV)
ar.optimized_model
Explanation: The AR stores several properties from the analysis:
End of explanation
ar.write_to("test_mle.fits", overwrite=True)
Explanation: Saving results to disk
The beauty of the analysis result is that all of this information can be written to disk and restored at a later time. The statistical parameters, best-fit model, etc. can all be recovered.
AR's are stored as a structured FITS file. We write the AR like this:
End of explanation
import astropy.io.fits as fits
ar_fits = fits.open('test_mle.fits')
ar_fits.info()
Explanation: The FITS file can be examines with any normal FITS reader.
End of explanation
ar_reloaded = load_analysis_results("test_mle.fits")
ar_reloaded.get_statistic_frame()
Explanation: However, to easily pull the results back into the 3ML framework, we use the ${\tt load_analysis_results}$ function:
End of explanation
ar_reloaded.get_data_frame()
Explanation: You can get a DataFrame with the saved results:
End of explanation
from threeML.analysis_results import AnalysisResultsSet
analysis_set = AnalysisResultsSet([ar, ar_reloaded])
# index as time bins
analysis_set.set_bins("testing", [-1, 1], [3, 5], unit="s")
# write to disk
analysis_set.write_to("analysis_set_test.fits", overwrite=True)
analysis_set = load_analysis_results("analysis_set_test.fits")
analysis_set[0].display()
Explanation: Analysis Result Sets
When doing time-resolved analysis or analysing a several objects, we can save several AR's is a set. This is achieved with the analysis result set. We can pass an array of AR's to the set and even set up descriptions for the different entries.
End of explanation
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_interval())
Explanation: Error propagation
In 3ML, we propagate errors for MLE reults via sampling of the covariance matrix instead of Taylor exanding around the maximum of the likelihood and computing a jacobain. Thus, we can achieve non-linear error propagation.
You can use the results for propagating errors non-linearly for analytical functions:
End of explanation
print("\nThis is 50 * a/a:")
print(50 * p1/p1)
Explanation: The propagation accounts for covariances. For example this
has error of zero (of course) since there is perfect covariance.
End of explanation
print("\nThis is arcsinh(a + 5*b) / np.log10(b) (why not?)")
print(np.arcsinh(p1 + 5 * p2) / np.log10(p2))
Explanation: You can use arbitrary (np) functions
End of explanation
print("\nRatio a/b:")
print(p2 / p1)
Explanation: Errors can become asymmetric. For example, the ratio of two gaussians is
asymmetric notoriously:
End of explanation
def my_function(x, a, b):
    return b * x ** a
print("\nPropagating using a custom function:")
print(my_function(2.3, p1, p2))
Explanation: You can always use it with arbitrary functions:
End of explanation
def go(fitfun, ar, model):
    fig, ax = plt.subplots()
    # 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, 50, 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_interval()
                low_curve[i], middle_curve[i], hi_curve[i] = (
                    low_bound,
                    this_flux.median,
                    hi_bound,
                )
                p.increase()
    ax.plot(energies, middle_curve, "--", color="black")
    ax.fill_between(energies, low_curve, hi_curve, alpha=0.5, color="blue")
go(fitfun, ar, model)
Explanation: 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)
End of explanation
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)
Explanation: Bayesian Analysis Results
Analysis Results work exactly the same under Bayesian analysis. 
Let's run the analysis first.
End of explanation
ar2 = bs.results
Explanation: Again, we grab the results from the BayesianAnalysis object:
End of explanation
ar2.write_to("test_bayes.fits", overwrite=True)
ar2_reloaded = load_analysis_results("test_bayes.fits")
Explanation: We can write and read the results to/from a file:
End of explanation
np.allclose(ar2_reloaded.samples, ar2.samples)
Explanation: The AR holds the posterior samples from the analysis. We can see the saved and live reults are the same:
End of explanation
#ar2.corner_plot();
# with chain consumer (pretty!)
ar2.corner_plot_cc();
Explanation: NOTE: MLE AR's store samples as well. These are the samples from the covariance matrix
We can examine the marginal distributions of the parameters:
End of explanation
ar2.get_data_frame("equal tail")
ar2.get_data_frame("hpd")
Explanation: We can return pandas DataFrames with equal tail or HPD results.
End of explanation
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)
Explanation: Error propagation operates the same way. Internally, the process is the same as the MLE results, however, the samples are those of the posterior rather than the (assumed) covariance matrix.
End of explanation
go(fitfun, ar2, model)
Explanation: To demonstrate how the two objects (MLE and Bayes) are the same, we see that our plotting function written for the MLE result works on our Bayesian results seamlessly.
End of explanation
<END_TASK>