Arviz: contour plot

[Hossein-1990]

Hi everybody
I’m going to plot a 95% probability contours using Arviz. I know that plot_pair is for constructing contours, but there isn't any option to adjust the HDI interval. Is there any other function to approach this problem?

Thanks

[ahartikainen]

Our plot_pair basically calculates density with _fast_kde_2d function and then pushes that density on a grid to matplotlib contour function.

https://github.com/arviz-devs/arviz/blob/d10e5b3f1d71152f06e6d138807118c36cbd122b/arviz/stats/density_utils.py#L852

I'm not sure if there is a way to set specific values for contours, but if there is, then it should be possible through kde_kwargs -> contour_kwargs

az.plot_pair(..., kde_kwargs={"contour_kwargs":  {"levels": [hdi_025, hdi_975]}})

But then, what are the correct hdi_025 and hdi_975 values.

@OriolAbril @aloctavodia 2D density -> user manually calls the _fast_kde_2d and then finds the correct density levels from that grid? What values should user look for?

[OriolAbril]

I guess it is possible provided you trust the kde and don't mind computing things only within the grid it's defined in. As Ari suggested, you can choose which levels to plot when calling matplotlib contour function, but I don't think this is of any use in practice. There are too many things to check and that will probably depend on each specific pair of variables (kwargs are common for all of them).

Here is a pseudo list of steps I would follow:

1. compute the 2d kde with _fast_kde_2d _and ensure it is normalized_ (I don't remember if we do normalize it, I don't think this is visible in the countour plots)
2. define a function to find the levels you want, that will probably have to be a numerical root finding algorithm. You choose a kde level, find the area in the 2d space with kde values larger than that, and integrate the kde in that area only, compare to your probability and find a level close enough. This will most probably give different levels for each pair of variables
   
   
   
   • You can also define a fine enough grid between 0 and max(kde_2d) and compute the integral for each point in the grid, then choose the level closer to your desired value. If done carefully, I believe you'll be able to evaluate some simple quadrature in a vectorized way along all the grid at once.
   
   
3. use those levels as levels in the matplotlib contourplot.

[ahartikainen]

One more thing to consider, you can probably calculate the area with the calculated contour (see bokeh code for calculating the contour locations with mpl functions manually).

(but then need to integrate over that area still)

[Hossein-1990]

I guess it is possible provided you trust the kde and don't mind computing things only within the grid it's defined in. As Ari suggested, you can choose which levels to plot when calling matplotlib contour function, but I don't think this is of any use in practice. There are too many things to check and that will probably depend on each specific pair of variables (kwargs are common for all of them).

Here is a pseudo list of steps I would follow:

1. compute the 2d kde with _fast_kde_2d _and ensure it is normalized_ (I don't remember if we do normalize it, I don't think this is visible in the countour plots)

2. define a function to find the levels you want, that will probably have to be a numerical root finding algorithm. You choose a kde level, find the area in the 2d space with kde values larger than that, and integrate the kde in that area only, compare to your probability and find a level close enough. This will most probably give different levels for each pair of variables

   • You can also define a fine enough grid between 0 and max(kde_2d) and compute the integral for each point in the grid, then choose the level closer to your desired value. If done carefully, I believe you'll be able to evaluate some simple quadrature in a vectorized way along all the grid at once.
3. use those levels as levels in the matplotlib contourplot.

In my opinion, it isn't logical too (to specify certain contour), but apparently this is common in some areas! this is the one of those:

Note that, the original programming language for this graph is R! (maybe its code in R is easier than python!)

[Hossein-1990]

One more thing to consider, you can probably calculate the area with the calculated contour (see bokeh code for calculating the contour locations with mpl functions manually).

(but then need to integrate over that area still)

Thank you. For being more clarified, I added an image to show the case which I'm going to plot. But as I mentioned above, seemingly in R it is more common to plot those kinds of plots (and easy to write too)

[OriolAbril]

You should take a look at the actual source code for the plot if you can, the idea should be the same, and the main problem is conceptual.

I also remembered that corner.py supports this already and also supports InferenceData inputs. It uses a similar apporach to the idea I outlined in the sub-bullet point on how to compute it, in this case, using a midpoint quadrature in a clever and elegant way, you can see the code here

[aloctavodia]

I will try to implement a general solution, unless someone else is interested.

[Hossein-1990]

Thanks to all responses, I could solve the problem with corner library. These are codes for an artificial example:

you can extend the method to any function.

[wm1995]

I actually also had this issue, and had been in the process of writing a pull request to fix it. I've never written a pull request before, but I do think I have working code to solve this, if it would be helpful.

The algorithm I used was to take the grid output by _fast_kde_2d, flatten and sort it, and then take the cumulative sum. I could then calculate the HPD thresholds from that and pass them as arguments to plt.contourf.

The code I have seems to work - should I submit it as a pull request? Are there any guidelines for this?

Incidentally, I found that passing arguments for the levels of contours doesn't work with the bokeh backend as it stands, but I wasn't sure whether I should look at that as well (perhaps I should file a separate bug report for that?).

[OriolAbril]

I'll let @aloctavodia guide you on this as he may have already done some work on that.

We do have contributing guidelines at CONTRIBUTING.md.

[wm1995]

Thank you for the guidance! I'll wait to hear from @aloctavodia before doing anything further.

[aloctavodia]

Thanks @wm1995 for your offer to solve this. Let me know if you need help. If possible I would recommend you make a PR and then we can discuss over that, but let me know if you need help prior to that.