tomasstolker / species

Toolkit for atmospheric characterization of directly imaged exoplanets
https://species.readthedocs.io
MIT License
22 stars 10 forks source link

Adding priors to corner plot output in plot_posterior #78

Closed gabrielastro closed 7 months ago

gabrielastro commented 9 months ago

In plot_posterior (from species.plot.plot_mcmc), would it be possible to add a simple option to overplot or underplot in the diagonal panels the prior (say in red or as a dotted grey line)? Thanks!

gabrielastro commented 8 months ago
  1. Noticed in passing: max_prob=True and inc_luminosity=True|inc_mass=True seem incompatible:
    
    File […]/species/plot/plot_mcmc.py:860, in plot_posterior(tag, burnin, title, offset, title_fmt, limits, max_prob, vmr, inc_luminosity, inc_mass, inc_log_mass, inc_pt_param, inc_loglike, inc_abund, output, object_type, param_inc)
    858 if i > j:
    859     if max_prob:
    --> 860         ax.axhline(par_val[i], color="tomato")
    861         ax.plot(par_val[j], par_val[i], "s", color="tomato")
    863     if limits is not None:

IndexError: tuple index out of range



3. When the marginalised posterior does not have a peak (increasing or decreasing to an edge), is the median meaningful? Maybe it would make more sense to plot vertical lines for upper limits and to write "*Q* < *X* (1 sigma)" or so as the column title for quantity `Q`, and to place `max_prob` at the ege.
tomasstolker commented 7 months ago

Thanks @gabrielastro! I have implemented the first two suggestions 👍.

Robustly determining an upper limit from a cutoff posterior is not possible I think, so it is up to the user to decide if/how to interpret the results.

Note that there is the get_samples method of Database for any custom analysis of the samples.

tomasstolker commented 7 months ago

The plot_posterior has the show_priors parameter for overplotting the normal priors. Showing the uniform priors seems not useful since it is just a horizontal line...

gabrielastro commented 7 months ago

Thanks a lot! That works well. Thanks for the pointer to get_samples. a.) Would it be possible by default to under-plot the prior? b.) Would it still be possible to offer the option of showing the (top-hat) prior specifically for Teff? Alternatively, maybe plot by default the whole range considered for the fitting (through bounds by the user or automatically based on the models), not the region around the maximum posterior. A common outcome of high-SNR fitting is a very small errorbar on Teff and it is not clear from the plot what input range of models was considered. I get very different peaks (not up to the edges) when taking Teff = (1000,1500) or (1400,1800), for example… [I cannot afford (RAM-wise) to load all exorem-high-res spectra so have to try it out in blocks.]

gabrielastro commented 7 months ago

For a cut-off posterior, would it be more meaningful to write Q \geqslant X1 for a rising posterior cut off at X1, and not to plot the vertical dashed lines? However, if the cut-off is physical (mass or vsini towards zero, eccentricity towards the edges, etc.), there are well-defined upper limits that could be plotted automatically, no :thinking:?

tomasstolker commented 7 months ago

I get very different peaks (not up to the edges) when taking Teff = (1000,1500) or (1400,1800), for example…

Hm that's odd and should not happen I think. I am guessing your fits have not fully converged to the global minimum.

For any detailed adjustments, I would suggest manually adjusting the Figure object that is returned by plot_posterior to your own needs 😊.

tomasstolker commented 7 months ago

Ah but the overlap in these two Teff cases is only small. I assume that the solutions don't end up differently in the 1400-1500 K regime?

gabrielastro commented 7 months ago

Hmm, not sure anymore actually :no_mouth:… I did not write the Teff range in the file names and there is no easy way to check with database.list_content(), it seems; it just says in the - bounds section: - teff: <HDF5 dataset "teff": shape (2,), type "<i8"> but not what the two values are….

But the exo-rem-highres grid I was using has some apparently not monotonic behaviour over Teff=1100–1200 K; whether physical or numerical, I want to check at some point but it is probably linked to that. Stay tuned :radio:…

For any detailed adjustments, I would suggest manually adjusting the Figure object Ok!

So then maybe only 1.) underplotting and 3.) dealing with the upper or lower limits are left, if possible. Thanks in any case for the changes you already made!

tomasstolker commented 7 months ago

Maybe I have misunderstood, but you can simply check the posterior plot to compare the Teff constraints between the two fits. Or, you can get the posterior samples with get_samples and plot Teff.

But it is not so important 😉, my guess is that the solutions for these two fits will not occur at different Teff in the overlapping regime (1400-1500 K), but at different Teff outside that regime.

1.) Hm I think it would make the plot less clear so I will leave it like this. 3.) Mass and vsin(i) can not go to zero I think. If the mass does, then there is probably an issue with the parallax/radius. If it is a matter of plotting then you could use the inc_log_mass parameter in plot_posterior. The vsin(i) can only go to zero if the instrument resolution is infinite.

gabrielastro commented 7 months ago

Ah! thank you very much for the explanations. Right, that makes sense.

Just the prior behind the other curves would actually be more logical (it is a prior… :wink:) and would make those panels clearer. But it is a small thing indeed :microscope:

tomasstolker commented 7 months ago

I tried but it looks worse. You can easily adjust it with the returned Figure object!

gabrielastro commented 7 months ago

Ok, thank you very much! Good to know that these things can be changed ‘after the fact’ (I am not a python plotter :construction_worker_man: , just a species user… :chart_with_upwards_trend:)!