Spectral1D Simulation and Analysis with gammapy

[Gaia Verna]

Hi everyone,

I'm working with 1D Spectral Simulation and Analysis for CTA with gammapy.

I'm using SpectrumSimulation and SpectrumFit classes and I'm considering a simple case of an exponential cutoff power law (ECPL) simulated source and I'm performing the fit with the same spectral model:

I have found some strange behavior in the fit results and I'm looking for some help.

I decided to simulate 1000 observations of the same sources and to inspect the best spectral parameter values from the fit each time.

I noticed that the spectral parameters are found out only in some of these simulated observations: at a certain moment in the loop, the best fitted values starts to be -nan- and this behavior seems to "condition" the results for the parameters in all the next iterations (i.e. independent observations)

I'm trying to simulate very faint and/or soft sources...this probably conditions the instability of the fit results.

I attach a simple script.

thank you,

cheers,

Gaia

[Axel Donath]

Hi Gaia,

I have taken a look at your script and here are a few things to consider:

1. Move the creation of `ecpl` model inside the simulation loop. When fitting a model in Gammapy the model parameters

are changed in place. In case the fit fails and the parameters are set to NaN and the model cannot be reused. This should

fix the problem you see.

2. Use `sim.simulate_obs(seed=)` to use a different seed for each of the simulations, otherwise you'll always get the same result.

3. I would split the simulation and fitting part into two steps: first do the simulation, write the results to file, then read the data in

again and fit it to get your parameter distributions. 

4. For the fitting part it's probably useful to use an astropy.Table object to collect your fitting results. So define a table with the values and

errors of the model parameters as table columns and the add one fitted observation per row. Again you can then write the table to file

and read again e.g. to plot the parameter distributions.

5. The `FitSpectrum` class will be removed in the next Gammapy version `v.012` and you'll have to adapt your script, if you want to update 

Gammapy (it's just a small API change, but let us know, when you want to update and have problems with the transition). 

Hope this helps, Axel

[Gaia Verna]

Hi Axel,

your first suggestion (creation of a NEW spectral model to use in the fit for each observation) seems to fix the problem.

Now only for few observations the fit results are -nan- but this is now related (I think...) to the fact that for faint(/soft) sources observed for few time the statistic is not enough to reconstruct the spectral characteristics i.e. the fit doesn't converge.

For what concerns the seed I used the default initialization that is

def simulate_obs(self, obs_id, seed="random-seed"):

...so the seed is different for each observation and the results too.

Anyway I agree that if I would fix the seed to a constant value, nothing would have been change.

Thank you very much for your help,

Gaia

--
You received this message because you are subscribed to a topic in the Google Groups "gammapy" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/gammapy/KxSBYgB47-A/unsubscribe.
To unsubscribe from this group and all its topics, send an email to gammapy+u...@googlegroups.com.
Visit this group at https://groups.google.com/group/gammapy.
To view this discussion on the web visit https://groups.google.com/d/msgid/gammapy/f56e5e9b-c211-40bc-a1db-8205b089b266%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Axel Donath]

Am Dienstag, 21. Mai 2019 12:58:04 UTC+2 schrieb Gaia Verna:

Hi Axel,

your first suggestion (creation of a NEW spectral model to use in the fit for each observation) seems to fix the problem.

Now only for few observations the fit results are -nan- but this is now related (I think...) to the fact that for faint(/soft) sources observed for few time the statistic is not enough to reconstruct the spectral characteristics i.e. the fit doesn't converge.

For what concerns the seed I used the default initialization that is

def simulate_obs(self, obs_id, seed="random-seed"):

...so the seed is different for each observation and the results too.

But you also want reproducible results, that's why I suggested to seed the random number generation in a deterministic way.

 

Anyway I agree that if I would fix the seed to a constant value, nothing would have been change.

Thank you very much for your help,

I'm glad, I could help...

Axel 

 

Gaia

To unsubscribe from this group and all its topics, send an email to gammapy+unsubscribe@googlegroups.com.