How to extract the overlap of two samples from Cobaya and draw the contour plot?

[Biagio De Simone]

Good afternoon,

Thank you for welcoming me into the CosmoCoffee forum! I have a question about Cobaya.

I have to check the distribution that results from the overlap of two or more samples drawn from Cobaya package (for simplicity, let's consider only two for the moment). For the first distribution, the code is the following:

Code: Select all

gdsamples1of3 = MCSamplesFromCobaya(updated_info, products.products()["sample"],ignore_rows=0.3)
gdplot = gdplt.getSubplotPlotter(width_inch=5)
p1 = gdsamples1of3.getParams()
O_m1=p1.O_m
H01=p1.H0
gdsamples1of3.addDerived(O_m1, name='O_m1', label=r"\Omega_{0m}")
gdsamples1of3.addDerived(H01, name='H01', label=r"H_0(Km \: s^{-1} \: Mpc^{-1})")
gdplot.triangle_plot(gdsamples1of3, ["O_m1","H01"], filled=True)

while for the second the code is similar

Code: Select all

gdsamples2of3 = MCSamplesFromCobaya(updated_info, products.products()["sample"],ignore_rows=0.3)
gdplot = gdplt.getSubplotPlotter(width_inch=5)
p1 = gdsamples2of3.getParams()
O_m1=p1.O_m
H01=p1.H0
gdsamples2of3.addDerived(O_m1, name='O_m1', label=r"\Omega_{0m}")
gdsamples2of3.addDerived(H01, name='H01', label=r"H_0(Km \: s^{-1} \: Mpc^{-1})")
gdplot.triangle_plot(gdsamples2of3, ["O_m1","H01"], filled=True)

When I run the codes and plot the distributions I obtain the overlap of the whole contours, e.g. in the following case with multiple distributions:

Code: Select all

plot=g.plot_2d([gdsamples1of3,gdsamples2of3,gdsamples3of3,gdsamplesfull],'O_m1','H01',filled=False,colors=['cyan','orange','lime','red'],lims=[0, 1, 64, 76])

What I would like to obtain is only the sub-distribution given by the overlap of all the contours so that I can plot it alone. I draw it in black in the following image.



Can you please clarify me how to do it? Thank you in advance!

[Aruna Harikant]

Hi Simone,

I have tried to solve your problem by following https://getdist.readthedocs.io/en/latest/plot_gallery.html.
This is one example plot from Plot Gallery.


I have followed some steps here to plot only the overlap region between the two distributions...
I have created grid of points and checked for those lie within both contours. Then. created a mask that shows only the overlapping points.

Code: Select all

%matplotlib inline
g = plots.get_single_plotter(width_inch=3, ratio=1)
g.plot_2d([samples, samples2], 'x1', 'x2', filled=True)
# Get the paths from the contour plots
paths1 = g.subplots[0][0].collections[0].get_paths()
paths2 = g.subplots[0][0].collections[1].get_paths()

# Create a grid of points in the range of your plots
x = np.linspace(-5, 5, 500)
y = np.linspace(-5, 5, 500)
X, Y = np.meshgrid(x, y)
points = np.c_[X.ravel(), Y.ravel()]

for path1 in paths1:
    for path2 in paths2:
        # Check which points are inside each contour
        mask1 = path1.contains_points(points).reshape(X.shape)
        mask2 = path2.contains_points(points).reshape(X.shape)
        
        # Find overlapping region and plot it
        overlap_mask = np.logical_and(mask1, mask2)
        g.subplots[0][0].scatter(X[overlap_mask], Y[overlap_mask], color='purple', s=1, alpha=0.6)

# Remove the original unfilled contours
for collection in g.subplots[0][0].collections:
    collection.remove()

g.add_legend(['Overlap Region'], colored_text=True);
plt.show()

Please let me know if this approach solves your problem!

[Antony Lewis]

I assume the contour is supposed to be the joint constraint. In general you cannot do this without generating new samples by sampling from the joint likelihood (using all three bins). (if you have more than 2 parameters, the overlap of the 2D projected contours is in general not the same as two the 2D projection of the ND-intersection of the ND distributions)

[Biagio De Simone]

Dear All,

Thank you for your support. The graphical solution of Aruna Harikant is useful for the graphical showing of the intersection and surely I will leverage it but my idea was to consider the joint distribution of the posteriors.

As Antony Lewis was pointing out, I should rewrite the likelihood considering the joint contribution of the three or more bins. But in this case it is already given since it is the one in red called "Full Pantheon".

Thank you again for your kind support.

Biagio

[Aruna Harikant]

Hi Simone,

I'm grateful for Prof. Lewis's insights and based on his suggestions, I am addressing this issue. Accordingly, the focus of this problem is on separating the joint constraint. For eg. I think we can separate it out, one method is to apply kernel density estimation.

example.png (40.42 KiB) Viewed 11104 times

[/img].

Please let me know if you still face any problems separating it out.

[Biagio De Simone]

Dear Aruna Harikant,

This is a very good idea! So, what should I do with the code? Extract the posteriors for the different bins as list of values, join together the values and run a KDE method in Python?

Thanks to you and Prof. Lewis for your kind support!

Biagio

[Jesus Torrado]

Hi both,

As stated by Antony a few posts above, this is not correct: the joint constraint is imposed by the product of the likelihoods, which intuitively is closer to the intersection than the union (what you are doing). As Antony said, you have to generate samples from the joint posterior, containing the two likelihoods.

Alternatively, you can use post-processing (https://cobaya.readthedocs.io/en/latest/post.html) to reweight ("add") one of the samples with the other likelihood. Since you have samples from each individual likelihood, you can reweight each with the other likelihood to check that you get consistent results (if the result differs significantly between the two approaches, you cannot reweight, and need to sample from the joint posterior).

[Biagio De Simone]

Dear Prof. Torrado,

Thank you for the further clarification. I will proceed with the joint likelihoods or the alternative approach you suggested.

With best regards,

Biagio