Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add MCMC tutorial using emcee #2077

Merged
merged 13 commits into from Mar 25, 2019

Conversation

5 participants
@facero
Copy link
Contributor

facero commented Mar 11, 2019

A first PR for the MCMC tutorial with emcee mentioned in #1923. This is in fact a very nice showcase of the power of the recent dataset class !

The notebook simulates a dataset with gammapy.cube.simulate.simulate_dataset, runs a chain of walkers, and plots results.
The interface needed with emcee was to define a prior function and a Likelihood function. This was relatively simple thanks to Dataset (once you figure out whether you need a + or - for the likelihood ...).

I have many questions and there are many things to improve but the first thing that we should decide is:
What is the best spectral model to test here ?
I started with a bright source (Norm=5e-12) and Ecut=50 TeV for 20h. Should we try a fainter source, a higher Ecut (100 TeV ?). It seems to reconstruct the parameters well but I want to try PeVatrons and/or fainter sources but maybe this is not the point of this tutorial. Any thoughts ?

Then there are several things that could be improved :

  • How do I update the sky model of an existing dataset ?
    This dataset.model = sky_model updates the model but the dataset gives always the same likelihood, independently of the new model like if its not reading it.

  • Once I understand this, it is simple to compute the Npred for each sampler chain and overlay the sample models on top of the data.

  • Currently bkg parameters can be included in the Likelihood as they do not appear in dataset.parameters.free_parameter. This is normal are these parameters are in dataset.background_model.parameters. Not sure how to access them easily. Freezing Bkg parameter for now.

facero added some commits Feb 26, 2019

@adonath adonath self-assigned this Mar 12, 2019

@adonath adonath added the docs label Mar 12, 2019

@adonath adonath added this to To Do in Documentation via automation Mar 12, 2019

@adonath adonath added this to the 0.11 milestone Mar 12, 2019

@adonath

This comment has been minimized.

Copy link
Member

adonath commented Mar 12, 2019

Thank you very much @facero. This is indeed a very nice show-case for the MapDataset.

Just a few preliminary comments:

  • In 1ae6f22 I implemented a setter for MapDataset.model. Because of #2071, the model were cached and not updated, when the model was changed. This is solved now.

  • Can you run gammapy jupyter --src tutorials/MCMC_fitting.ipynb strip to remove the output cells. The notebook will be executed and rendered as part of our docs build. And maybe also run gammapy jupyter --src tutorials/MCMC_fitting.ipynb black to apply a uniform code style to the code cells.

  • A nitpicky moment: maybe rename the notebook to mcmc_sampling.ipynb. I think (strictly
    speaking) running an MCMC sampler is not a fitting in a classical sense.

The problem concerning the background parameters I didn't understand yet. If the background parameters are free they should be included in MapDataset.parameters.free_parameters and are available for sampling as well.

Concerning which spectral model to use I don't have a strong opinion. I think this tutorial does not have to show a scientifically super meaningful analysis. I think assuming any plausible spectral model is probably fine.

@facero

This comment has been minimized.

Copy link
Contributor Author

facero commented Mar 12, 2019

Thanks ! I will implement those changes.
I've just retried the bkg stuff and indeed it now works. I can now have norm and tilt in MCMC. Maybe somthing has changed recenlty. The norm and tilt param were not showing in free_param during the sprint.

One thing still is when doing print(dataset.model) the bkg param do not appear as they are not part of a Skymodel. From a user point of view I would like to see all the param from my model (SKy or not Sky).

@registerrier , @cdeil any other comments ?

@cdeil

This comment has been minimized.

Copy link
Member

cdeil commented Mar 12, 2019

@cdeil any other comments ?

I don't have time to check this or other Gammapy additions out this week, sorry. Please go ahead and merge.

@adonath

This comment has been minimized.

Copy link
Member

adonath commented Mar 12, 2019

@facero You can alway print the full parameter list using:

print(dataset.parameters)

Or print the background model separately:

print(dataset.background_model)

One thing that is currently missing in Gammapy is the back-reference of parameters to their models and datasets they belong to. Once we have this we can add a nicer API to work with the global model parameter list.

@facero

This comment has been minimized.

Copy link
Contributor Author

facero commented Mar 12, 2019

Indeed the print(dataset.parameters) works but it doesn't have the nice print that print(dataset.model) has but it might be related to what you mention.

@facero

This comment has been minimized.

Copy link
Contributor Author

facero commented Mar 13, 2019

@adonath I'e updated the notebook (strip, black) and added a cool plot showing the SED of 100s of samples:
MCMC

@adonath adonath self-requested a review Mar 13, 2019

@adonath
Copy link
Member

adonath left a comment

Thanks @facero! Right now, I don't have any further comments.

@facero facero closed this Mar 13, 2019

Documentation automation moved this from To Do to Done Mar 13, 2019

@registerrier registerrier reopened this Mar 14, 2019

Documentation automation moved this from Done to In Progress Mar 14, 2019

@registerrier

This comment has been minimized.

Copy link
Contributor

registerrier commented Mar 14, 2019

Hi @facero . We should merge the PR first.

@adonath

This comment has been minimized.

Copy link
Member

adonath commented Mar 14, 2019

@registerrier Do you have any further comments? From my side it's OK to merge.

@registerrier

This comment has been minimized.

Copy link
Contributor

registerrier commented Mar 14, 2019

OK to merge from my side too. There is an unrelated Travis fail. But as far as I can see the notebook introduced by this PR is not tested in the tutorial check. Is it expected?

@adonath

This comment has been minimized.

Copy link
Member

adonath commented Mar 14, 2019

Thanks @registerrier, you're right. I just added the notebook to the index file. Let's wait for Travis-CI to finish the testing and then merge...

@Bultako

This comment has been minimized.

Copy link
Member

Bultako commented Mar 14, 2019

Just a small reminder, if you want to make it visible in the docs, we also need to provide a link somewhere in the RST files. For example in docs/tutorials.rst.

Thanks @registerrier, you're right. I just added the notebook to the index file.

@adonath

This comment has been minimized.

Copy link
Member

adonath commented Mar 14, 2019

Thanks @Bultako, I just linked the notebook in the docs.

@facero

This comment has been minimized.

Copy link
Contributor Author

facero commented Mar 14, 2019

apparently Travis CI failed because of fermipy related tests. how do you handle this usually ?

@adonath

This comment has been minimized.

Copy link
Member

adonath commented Mar 14, 2019

@cdeil already opened PR: fermiPy/fermipy#273

@facero

This comment has been minimized.

Copy link
Contributor Author

facero commented Mar 19, 2019

is this PR blocked because of the fermipy issue ?

@adonath

This comment has been minimized.

Copy link
Member

adonath commented Mar 25, 2019

@facero No, it was blocked because of busy maintainers :-). I'll merge now...

@adonath adonath merged commit 72909e6 into gammapy:master Mar 25, 2019

2 of 3 checks passed

continuous-integration/travis-ci/pr The Travis CI build failed
Details
Codacy/PR Quality Review Up to standards. A positive pull request.
Details
Scrutinizer Analysis: No new issues – Tests: passed
Details

Documentation automation moved this from In Progress to Done Mar 25, 2019

@adonath adonath changed the title MCMC tutorial with emcee Add MCMC tutorial using emcee Mar 26, 2019

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.