Fit multiple background components from simulations

[israelmcmc]

Currently, we only fit a single background component. The basic code is relatively stragithforward, we just need to add more nuisance parameters. The hard part though is to make it converge.

[israelmcmc]

For now, it would be fine to add the capabilities, in principle, to fit an arbitrary number of components. But work on making it converge for only 2. See also #207 for many many components.

[GallegoSav]

@israelmcmc @ckarwin @eneights If I understood well the fit is not even converging when using the total background model right ? I have indeed a error message when I am trying with total_bg_3months_binned_for_continuum.hdf5 :

---------------------------------------------------------------------------
FitFailed                                 Traceback (most recent call last)
Cell In[30], line 5
      1 plugins = DataList(cosi) # If we had multiple instruments, we would do e.g. DataList(cosi, lat, hawc, ...)
      3 like = JointLikelihood(model, plugins, verbose = False)
----> 5 like.fit()

File ~/software/COSItools/python-env/lib/python3.9/site-packages/threeML/classicMLE/joint_likelihood.py:324, in JointLikelihood.fit(self, quiet, compute_covariance, n_samples)
    320 # Perform the fit, but first flush stdout (so if we have verbose=True the messages there will follow
    321 # what is already in the buffer)
    322 sys.stdout.flush()
--> 324 xs, log_likelihood_minimum = self._minimizer.minimize(
    325     compute_covar=compute_covariance
    326 )
    328 if log_likelihood_minimum == minimization.FIT_FAILED:
    329     log.error("The fit failed to converge.")

File ~/software/COSItools/python-env/lib/python3.9/site-packages/threeML/minimizer/minimization.py:632, in Minimizer.minimize(self, compute_covar)
    628 # Gather the best fit values from the minimizer and the covariance matrix (if provided)
    630 try:
--> 632     internal_best_fit_values, function_minimum = self._minimize()
    634 except FitFailed:
    636     raise

File ~/software/COSItools/python-env/lib/python3.9/site-packages/threeML/minimizer/minuit_minimizer.py:210, in MinuitMinimizer._minimize(self)
    206 if not self.minuit.valid:
    208     self._print_current_status()
--> 210     raise FitFailed(
    211         "MIGRAD call failed. This is usually due to unconstrained parameters."
    212     )
    214 else:
    215 
    216     # Gather the optimized values for all parameters from the internal
    217     # iminuit dictionary
    219     best_fit_values = []

FitFailed: MIGRAD call failed. This is usually due to unconstrained parameters.

[israelmcmc]

@GallegoSav Yes, that's correct. It converged for one single component, but not when summing up all of them. It's not clear to me if adding more background was what caused this, or if it's simply that the fit is not very robust and we just got lucky before. In any case, this would be the starting point indeed, but any trick to make it convert should also help the multiple components fit.

[israelmcmc]

~@GallegoSav I was asking Yong to add the tests for this PR, but you'll need to do it as well at some point hehe. Here are some instructions: https://github.com/cositools/cosipy/blob/c8eb53631f686621d93a8b56395ce66ae3d71825/docs/dev/unit_tests.rst~ Wrong issue! I meant to post this in #171

[GallegoSav]

@GallegoSav Yes, that's correct. It converged for one single component, but not when summing up all of them. It's not clear to me if adding more background was what caused this, or if it's simply that the fit is not very robust and we just got lucky before. In any case, this would be the starting point indeed, but any trick to make it convert should also help the multiple components fit.

@israelmcmc It seems to be due to high background level indeed I have similar error when I am using just the albedo photon scale by 10 or the total background .

[israelmcmc]

@GallegoSav I see. I'd check the significance of the Crab when using the total background, to see if this is below the sensitivity. Using the 3ML, you can get the log-likelihood with and without the source. Keep all parameters fixed to the known values and get the likelihood (source hypothesis L1). Then set the source flux normalization to 0 and get the likelihood again (null hypothesis L0). The TS = 2*log(L1/L0). I remember @ckarwin did this, but I don't remember exactly for which case.

[GallegoSav]

@GallegoSav I see. I'd check the significance of the Crab when using the total background, to see if this is below the sensitivity. Using the 3ML, you can get the log-likelihood with and without the source. Keep all parameters fixed to the known values and get the likelihood (source hypothesis L1). Then set the source flux normalization to 0 and get the likelihood again (null hypothesis L0). The TS = 2*log(L1/L0). I remember @ckarwin did this, but I don't remember exactly for which case.

I obtained TS=1e-4 , assuming Wilks' theorem and 4 degrees of freedom that would give pvalue close to one so no possibility to exclude the null hypothesis.

I think it is not really surprising because there are no cuts at all on the CDS or pointing cuts so this is very difficult to fit directly the crab

[israelmcmc]

Thanks @GallegoSav. I'll open the discussion for this on Slack, since this doesn't seem like an issue with the code itself.