How to create a complex Gaussian Joint PDF plot? (Example figure)

Question

I am essentially trying to replicate this figure as close as possible:

Especially the one on the right with the red curve.

These are figures showing a curve generated by taking products of powers of the 2 Gaussian PDFs $\Theta_i$, $\Theta_j$ in the first figure. This needs 2 main things, the first is plotting the red curve the points and the labels, which I can do do easily in metapost. But the second portion requires plotting the combined density PDF of the 2 Gaussian PDFs as a color gradient and I don't think metapost can do this. I then tried using tkiz to plot this instead, but A) I am not very profficient with it as I tend to use metapost more often, and B) It seems really hard to define custom color shading for the background on tkiz.

Following a prior answer I managed to get a gaussian ellipsoid with a black dot at its mean, but any further manipulation started to break the figure. For example trying to add a second gaussian ellipsoid would render it on top of the first and occlude it. Moreover defining multidimensional functions in tkiz seems tricky (then again the same is true for metapost).

Any suggestions? I am ok with using tkiz, metapost or any other tool.

I stole this from the linked question:

\documentclass{standalone}

\usepackage{pgfplots}
\usepgfplotslibrary{patchplots}

\begin{document}

\pgfplotsset{
colormap={whitered}{color(0cm)=(white); color(1cm)=(orange!75!red)}
}

\begin{tikzpicture}[
    rotate=30,
    declare function={mu1=1;},
    declare function={mu2=2;},
    declare function={sigma1=0.5;},
    declare function={sigma2=1;},
    declare function={normal(\m,\s)=1/(2*\s*sqrt(pi))*exp(-(x-\m)^2/(2*\s^2));},
    declare function={bivar(\ma,\sa,\mb,\sb)=
        1/(2*pi*\sa*\sb) * exp(-((x-\ma)^2/\sa^2 + (y-\mb)^2/\sb^2))/2;}]
\begin{axis}[
    colormap name=whitered,
    width=15cm,
    view={0}{90},
    enlargelimits=false,
    domain=0:2,
    y domain=0:4,
    samples=17,
    axis equal image,
    hide axis,
]
\addplot3 [surf, draw=none, shader=interp, patch type=bilinear] {bivar(mu1,sigma1,mu2,sigma2)};
\end{axis}
\end{tikzpicture}
\end{document}

Could you please post some example such that others have at least the data to play with? — , May 14 '20 at 02:21
This is just some random Gaussian blob, isn't it? So you just want to have two of those in one plot? — , May 14 '20 at 02:41
Almost (and that would be good enough). However if you look closely at the picture, the 2 gaussian blobs interact with each other where they touch. The density of a given point in that figure is equal to the common mixture of both PDFs. But that;s icing on the cake if I can at minimum get the 2 gaussian ellipsoids with the curve I am happy. — Makogan, May 14 '20 at 02:44

score 4 · Accepted Answer · answered May 14 '20 at 07:18

This does not really produce the figures you want. I am not even sure if you should do a plot of that type when you neither have real data nor a function at hand. Rather, you can just draw the ellipses as you need them.

I am posting this because in your link the nontrivial principle axes of the Gaussian were obtained by rotating the picture. Here, on the other hand, a rotation angle is introduced that allow you to rotate the principle axes. So you can plot two relatively rotated Gaussians, for instance. In the code the rotation angle are alpha1 and alpha2. The peaks are given by the mu parameters, so you have the full 2D Euclidean group at your disposal to move the Gaussians around.

\documentclass[tikz,border=3mm]{standalone}
\usepackage{pgfplots}
\pgfplotsset{compat=1.17} 

\begin{document}
% 
% \pgfplotsset{
% colormap={whitered}{color(0cm)=(white); color(1cm)=(orange!75!red)}
% }

\begin{tikzpicture}[
    declare function={mu1=1;mu2=2;sigma1=0.5;sigma2=1.5;alpha1=-60;
        mu3=2.25;mu4=1;sigma3=0.3;sigma4=1;alpha2=0;},
    declare function={normal(\m,\s)=1/(2*\s*sqrt(pi))*exp(-(x-\m)^2/(2*\s^2));},
    declare function={bivar(\ma,\sa,\mb,\sb,\a)=
        1/(2*pi*\sa*\sb) * 
        exp(-(pow((x-\ma)*cos(\a)+(y-\mb)*sin(\a),2)/\sa^2 
            + pow(-(x-\ma)*sin(\a)+(y-\mb)*cos(\a),2)/\sb^2))/2;}]
\begin{axis}[%x=1cm,y=1cm,z=1cm,
    %colormap name=whitered,
    width=15cm,
    view={0}{90},
    enlargelimits=false,
    domain=0:4,
    y domain=0:4,
    samples=51,
    axis equal image,
    hide axis,
    colormap/hot
]
\addplot3[surf,% draw=none, 
    shader=interp,
    %smooth,contour gnuplot={number=14,labels=false}
    ]   
    {bivar(mu1,sigma1,mu2,sigma2,alpha1)+bivar(mu3,sigma3,mu4,sigma4,alpha2)
    +10*bivar(mu1,sigma1,mu2,sigma2,alpha1)*bivar(mu3,sigma3,mu4,sigma4,alpha2)};
\end{axis}
\end{tikzpicture}
\end{document}

This does not look like realistic data? Because this isn't realistic data.

That image looks perfect actually. I am not dealing with data, this is actually for geometric modelling. But this image is awesome! Thank you — Makogan, May 14 '20 at 07:21

How to create a complex Gaussian Joint PDF plot? (Example figure)

1 Answers1