We would like to draw a bivariate normal distribution and show where the means from the two variables meet in the space.
There are a couple of things in the code that might be useful:
You can define mathematical functions using declare function={<name>(<argument macros>)=<function>;}, which will help to keep the code clean and avoid repetitions.
You can define a new colormap using \pgfplotsset{colormap={<name>}{<color model>(<distance>)=(<value1>); <color model>(<distance 2>)=(<value2>)} }. This is a very powerful feature, so definitely worth reading in the pgfplots manual.
The legend created by colorbar is a whole new plot, so you can configure it with all the usual axis options.
There are different ways for defining 3D functions: \addplot3 {<function>}; will evaluate <function> at every point on a grid and assume the result to be a z-value. \addplot3 ({<x>},{<y>},{<z>}); defines a parametric function in 3D space, which allows you to draw three-dimensional lines, among other things.
Edit and compile if you like:\documentclass[border=10pt]{standalone} \usepackage{pgfplots} \pgfplotsset{width=7cm,compat=1.8} \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$}, 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}