How to plot 2D probability density function of a (2x2) gaussian mixture model?

Hello, I am new to using pymc and currently learning it. My goal is to build a (2X2) Gaussian mixture model that can effectively fit a 2D dataset comprising two distinct components. My data looks like this:

Following the example 3 in the page about pymc.miture, I define the model:

k_comp = 2 ## number of components 
k_dim = 2 ## number of dimensions
with pm.Model() as model:
    # w is a in shape of '(k_comp,)' 
    # the weights would be shared across the k_dim replication dimensions
    w = pm.Dirichlet('w', a=np.ones(k_comp), shape=(k_comp,))
    
    # mu and sigma is in shape of (k_dim, k_comp), each component at each dimension has an independent mu and sigma
    mu = pm.Normal('mu', mu=np.arange(k_comp), sigma=1, shape=(k_dim, k_comp))
    sigma = pm.HalfNormal('sigma', sigma=np.zeros((k_comp))+1, shape=(k_dim, k_comp))

 
    components = pm.Normal.dist(mu=mu, sigma=sigma, shape=(k_dim, k_comp))

    # The mixture is an array of k_dim elements
    # Each element can be thought of as an independent scalar mixture of k_comp
    # components with different means
    like = pm.Mixture('like', w=w, comp_dists=components, observed=data)
    
    idata = pm.sample()

I can run the model successfully and get the output. I use az.plot_trace(idata) to plot the posterior and chain for my results. It’s like

However, I am confused about how to plot the 2D histogram distribution of the model, similar to the one I have for my data. Additionally, I came across an example that shows how to plot the Probability Density Function (PDF) and group membership for a 1D Gaussian model. I am seeking guidance on how to create similar plots for the 2D Gaussian models in my context. Any assistance would be greatly appreciated.