How to draw a sample from a probability function and the show it in a tikz surf plot?

Question

I saw this post bivariate which plots a bivariate normal probability density function. What I want is to draw (statistically) a number of points based on this probability function and then show it in the same plot as the beautiful plot in the aforementioned question.

The following is a the source code from the above source:

\documentclass{standalone}

\usepackage{pgfplots}

\begin{document}

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

\begin{tikzpicture}[
    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={45}{65},
    enlargelimits=false,
    grid=major,
    domain=-1:4,
    y domain=-1:4,
    samples=26,
    xlabel=$x_1$,
    ylabel=$x_2$,
    zlabel={$P$},
    colorbar,
    colorbar style={
        at={(1,0)},
        anchor=south west,
        height=0.25*\pgfkeysvalueof{/pgfplots/parent axis height},
        title={$P(x_1,x_2)$}
    }
]
\addplot3 [surf] {bivar(mu1,sigma1,mu2,sigma2)};
\addplot3 [domain=-1:4,samples=31, samples y=0, thick, smooth] (x,4,{normal(mu1,sigma1)});
\addplot3 [domain=-1:4,samples=31, samples y=0, thick, smooth] (-1,x,{normal(mu2,sigma2)});

\draw [black!50] (axis cs:-1,0,0) -- (axis cs:4,0,0);
\draw [black!50] (axis cs:0,-1,0) -- (axis cs:0,4,0);

\node at (axis cs:-1,1,0.18) [pin=165:$P(x_1)$] {};
\node at (axis cs:1.5,4,0.32) [pin=-15:$P(x_2)$] {};
\end{axis}
\end{tikzpicture}
\end{document}

What I want is to draw a number of points (say 100) from the probability distribution, and then draw these points on the same or another plot.

I changed the link to the original source which draws a bivariate normal. — user85361, Jan 06 '19 at 11:14
Please help us help you and add a minimal working example (MWE) that illustrates your problem. Reproducing the problem and finding out what the issue is will be much easier when we see compilable code, starting with \documentclass{...} and ending with \end{document}. — DG', Jan 06 '19 at 12:11
@DG', thank you very much. I hope this makes the problem more clear. — user85361, Jan 06 '19 at 13:11

caverac · Accepted Answer · 2019-01-06T14:12:50.040

This is one option, I just used this other post to generate the numbers

\documentclass[border=10pt]{standalone}

\usepackage{pgfplots}

% https://tex.stackexchange.com/questions/148091/gaussian-random-numbers
\newcount\gaussF
\edef\gaussR{0}
\edef\gaussA{0}

\makeatletter
\pgfmathdeclarefunction{invgauss}{2}{%
  \pgfmathln{#1}% <- might need parsing
  \pgfmathmultiply@{\pgfmathresult}{-2}%
  \pgfmathsqrt@{\pgfmathresult}%
  \let\@radius=\pgfmathresult%
  \pgfmathmultiply{6.28318531}{#2}% <- might need parsing
  \pgfmathdeg@{\pgfmathresult}%
  \pgfmathcos@{\pgfmathresult}%
  \pgfmathmultiply@{\pgfmathresult}{\@radius}%
}

\pgfmathdeclarefunction{randnormal}{0}{%
  \pgfmathrnd@
  \ifdim\pgfmathresult pt=0.0pt\relax%
    \def\pgfmathresult{0.00001}%
  \fi%
  \let\@tmp=\pgfmathresult%
  \pgfmathrnd@%
  \ifdim\pgfmathresult pt=0.0pt\relax%
    \def\pgfmathresult{0.00001}%
  \fi
  \pgfmathinvgauss@{\pgfmathresult}{\@tmp}%
}

\makeatother

\pgfplotsset{width=7cm,compat=newest}
\pgfplotsset{%
  colormap={whitered}{color(0cm)=(white);
  color(1cm)=(orange!75!red)}
}
\begin{document}
\begin{tikzpicture}[
  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           = {45}{65},
    enlargelimits  = false,
    grid           = major,
    domain         = -1:4,
    y domain       = -1:4,
    samples        = 26,
    xlabel         = $x_1$,
    ylabel         = $x_2$,
    zlabel         = {$P$},
    xmin           = -1,
    xmax           = 4,
    ymin           = -1,
    ymax           = 4,
    colorbar,
    colorbar style = {
      at     = {(1,0)},
      anchor = south west,
      height = 0.25*\pgfkeysvalueof{/pgfplots/parent axis height},
      title  = {$P(x_1,x_2)$}
    }
  ]
    \addplot3 [surf] {bivar(mu1,sigma1,mu2,sigma2)};
    \addplot3 [domain=-1:4,samples=31, samples y=0, thick, smooth]
      (x,4,{normal(mu1,sigma1)});
    \addplot3 [domain=-1:4,samples=31, samples y=0, thick, smooth]
      (-1,x,{normal(mu2,sigma2)});

    \pgfplotsinvokeforeach{1,...,500} {
      \pgfmathsetmacro{\xx}{mu1 + sigma1 * randnormal}
      \pgfmathsetmacro{\yy}{mu2 + sigma2 * randnormal}
      \addplot3 [mark=*, mark size=1pt] coordinates {(\xx, \yy, 0)};
    }

    \draw [black!50] (axis cs:-1,0,0) -- (axis cs:4,0,0);
    \draw [black!50] (axis cs:0,-1,0) -- (axis cs:0,4,0);

    \node at (axis cs:-1,1,0.18) [pin=165:$P(x_1)$] {};
    \node at (axis cs:1.5,4,0.32) [pin=-15:$P(x_2)$] {};
  \end{axis}

\end{tikzpicture}
\end{document}

score 2 · Answer 2 · answered Jan 06 '19 at 15:55

Another proposal based on the inverse Gaussian trick.

\documentclass{standalone}

\usepackage{pgfplots}

\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;
    normal(\m,\s)=1/(2*\s*sqrt(pi))*exp(-(x-\m)^2/(2*\s^2));
    bivar(\ma,\sa,\mb,\sb)=
        1/(2*pi*\sa*\sb) * exp(-((x-\ma)^2/\sa^2 + (y-\mb)^2/\sb^2))/2;
    mybivar(\x,\y,\ma,\sa,\mb,\sb)=
        1/(2*pi*\sa*\sb) * exp(-((\x-\ma)^2/\sa^2 + (\y-\mb)^2/\sb^2))/2;   
    invgauss(\a,\b) = sqrt(-2*ln(\a))*cos(deg(2*pi*\b));}]
% invgauss based on https://tex.stackexchange.com/questions/254484/how-to-make-a-graph-of-heteroskedasticity-with-tikz-pgf/254494#254494    
\foreach\X in {1,...,100}   
{\pgfmathsetmacro{\myx}{mu1+0.5*sigma1*invgauss(rnd,rnd)}
\pgfmathsetmacro{\myy}{mu2+0.5*sigma2*invgauss(rnd,rnd)}
\pgfmathsetmacro{\myz}{mybivar(\myx,\myy,mu1,sigma1,mu2,sigma2)}
\ifnum\X=1
\xdef\lstOne{(\myx,\myy,\myz)}
\else
\xdef\lstOne{\lstOne (\myx,\myy,\myz)}
\fi}
\begin{axis}[
    colormap name=whitered,
    width=15cm,
    view={45}{65},
    enlargelimits=false,
    grid=major,
    domain=-1:4,
    y domain=-1:4,
    samples=26,
    xlabel=$x_1$,
    ylabel=$x_2$,
    zlabel={$P$},
    colorbar,
    colorbar style={
        at={(1,0)},
        anchor=south west,
        height=0.25*\pgfkeysvalueof{/pgfplots/parent axis height},
        title={$P(x_1,x_2)$}
    }
]
\addplot3 [surf] {bivar(mu1,sigma1,mu2,sigma2)};
\addplot3 [domain=-1:4,samples=31, samples y=0, thick, smooth] (x,4,{normal(mu1,sigma1)});
\addplot3 [domain=-1:4,samples=31, samples y=0, thick, smooth] (-1,x,{normal(mu2,sigma2)});

\addplot3 [only marks, mark=*, mark layer=like plot,mark size=1pt] coordinates {\lstOne};

\draw [black!50] (axis cs:-1,0,0) -- (axis cs:4,0,0);
\draw [black!50] (axis cs:0,-1,0) -- (axis cs:0,4,0);

\node at (axis cs:-1,1,0.18) [pin=165:$P(x_1)$] {};
\node at (axis cs:1.5,4,0.32) [pin=-15:$P(x_2)$] {};
\end{axis}
\end{tikzpicture}
\end{document}

How to draw a sample from a probability function and the show it in a tikz surf plot?

2 Answers2