Add different options to compute stat_array on FluxPointsDatasets

[QRemy]

Add different options to compute stat_array on FluxPointsDatasets :

• chi2 : etimate from chi2 statistics.
• profile : estimate from interpolation of the likelihood profile.
• distrib : estimate from probability distributions,
  assumes that flux points correspond to asymmetric gaussians and upper limits complemantary error functions.

Default is chi2, in that case upper limits are ignored and the mean of asymetrics error is used.
The distrib case provides an approximation if the profile is not available
which allows to take into accounts upper limits and asymetrics errors.

Updated /tutorials/analysis-1d/spectral_analysis.py to present the different cases:

Without ul the results are the same:

[codecov]

Codecov Report

Attention: Patch coverage is 97.36842% with 3 lines in your changes are missing coverage. Please review.

Project coverage is 94.46%. Comparing base (a048535) to head (a55f23c).
Report is 79 commits behind head on main.

Files	Patch %	Lines
gammapy/datasets/flux_points.py	97.46%	2 Missing ⚠️
gammapy/estimators/points/tests/test_lightcurve.py	95.45%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #5135      +/-   ##
==========================================
+ Coverage   94.25%   94.46%   +0.20%     
==========================================
  Files         234      235       +1     
  Lines       35226    35493     +267     
==========================================
+ Hits        33204    33528     +324     
+ Misses       2022     1965      -57     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[AtreyeeS]

Thanks @QRemy ! The profile is an elegant solution for considering upper limits during the fit!

I don't understand the implementation of the distrib option. Maybe you can add some information in gammapy.stat docs page?

Otherwise, only minor comments inline

[QRemy]

Thanks @AtreyeeS I implemented most of the comments.

I don't understand the implementation of the distrib option. Maybe you can add some information in gammapy.stat docs page?

Assuming gaussian errors the likelihood is given by the probability density function of the normal distribution.
For the upper limit case it is necessary to marginalize over the unknown measurement so we integrate the normal distribution up to the upper limit value which gives the complementary error function.

see eq. C7 from https://iopscience.iop.org/article/10.1088/0004-637X/773/2/168/pdf
(in eq C6 they are doing something more complicated as they also marginalised over the unknown true value)

[adonath]

Thanks @QRemy, the statistical approach to handle the ULs looks reasonable to me. I was wonderign whether we can maybe find a more descriptive name. I think even "chi2-marginalize-ul" would be better. We should definitely give the paper reference in addition.

[AtreyeeS]

Thank you @QRemy ! Elegant indeed.
Maybe you can add a little info the documentation (what you wrote in the comment, and the paper reference).
And add stat_type option in the example here

gammapy/examples/tutorials/analysis-time/light_curve_simulation.py

Line 295 in fdd05df

dataset_fp.models = model

[AtreyeeS]

The distrib option is often returning a nan for some bins, leading to the total stat becoming nan and a failure of the minimisation. Any idea why?

[registerrier]

Thanks @QRemy ! I have left some comments inline.

My main comment is that you could compute most of the necessary arrays in advance. I suppose that for large SEDs this might be much more efficient.

[bkhelifi]

Thanks @QRemy.. I have some questions for you.

[bkhelifi]

Thanks @QRemy . All good for me;)

[registerrier]

Thanks @QRemy . This looks good to me.

I have left a simple comment regarding not providing the dictionary of stat functions but simply their keys.

Otherwise, this would be much cleaner with a stat functions being defined as external objects. This can be a further PR though.