3

I am trying to plot below function as 3d surface plot, where $\theta_1,\theta_2$ are the x and y axes respectively. The Z axis is supposed to be $L(\theta_1,\theta_2)$.

$$
L(\theta_1,\theta_2) = \prod_{i=1}^{m} \dfrac{1}{\sqrt{2\pi\theta_2}}{\text{exp}}{\Big[ -\dfrac{  (x_i-\theta_1)^2 }{2\theta_2}  \Big]} \\
=  \Big(   \dfrac{1}{\sqrt{2\pi\theta_2}} \Big)^{m}{\text{exp}}{\Big[ 
\dfrac{-\sum_{i=1}^{m}(x_i-\theta_1)^2}{2\theta_2} \Big]}
\tag{1}
$$

Just in case of above latex is not rendered this is my equation:
enter image description here

I managed to get to an extent, but unable to find how to compute the summation by passing a series of values for \m., and eventually plot proper 3d surface plot for the function above. Kindly help.

MWE:

\documentclass{article}
\usepackage{pgfplots}
\pgfplotsset{compat=1.15}
\pgfmathdeclarefunction{joint_normal}{3}{%
    \pgfmathparse{ (1/(2*pi*#3))^(#1)*exp( -(#1-#2)^2/(2*#3^2)  )}%
}
\begin{document}
\begin{tikzpicture}
\begin{axis}[
    grid=both,
    restrict z to domain*=0:1,
    zmin=0,
    colormap/hot,
    %point meta min=-0.2, 
    %point meta max=1,    
    view={20}{20}  %tune here to change viewing angle
    ]

    \def\m{5}
    \addplot3[surf,domain=-30:30,domain y=0:30, samples=25] { joint_normal(\m, x, y) };
\end{axis}
\end{tikzpicture}
\end{document}

Current output:
enter image description here

online editor to try out: here

Note: One could assume $X_i$ values varying from any range say 0 to 10, from a normal distribution N(5, 4). Any range would do, but for not spending too much time on thinking of range, I provide this one. I could later tinker it as needed.

  • You can do summations with the tikzmath library by using recursions. Other than that, I am not aware of any other way of doing the sum in an elegant way in this framework. – marmot Oct 18 '18 at 18:10
  • I am new to tikz, can you please show with an example for above problem? – Parthiban Rajendran Oct 18 '18 at 18:11
  • Well, to be able to do that, one would need to know what the x_i are. You sum over the x_i but as long as you do not specify what they are it is impossible to plot the function. – marmot Oct 18 '18 at 18:20
  • One could assume $X_i$ values varying from any range say 0 to 10, from a normal distribution N(5, 4). Any range would do, but for not spending too much time on thinking of range, I provide this one. I could later tinker it as needed. – Parthiban Rajendran Oct 18 '18 at 18:25
  • If you can make a concrete proposal for the x_i I'll be happy to give it a shot. (Sorry, was offline for a few hours and will be really online in another few hours) – marmot Oct 18 '18 at 23:50
4

Here is an MWE, which at the same time contains the explanation.

\documentclass{article}
\usepackage{amsmath}
\usepackage{pgfplots}
\pgfplotsset{compat=1.16}
\usetikzlibrary{math} 
\begin{document}
% based on https://tex.stackexchange.com/a/307032/121799
% and https://tex.stackexchange.com/a/451326/121799
\def\xvalues{{0,1,2,4,5,7}} % notice that the `0` is the 0th entry, which is not used here
\tikzset{evaluate={
        function myN(\x,\z,\k) { % \x = \theta_1 and \z=\theta_2 
            if \k == 1 then {
                return myn(\x,\xvalues[1],\z);
            } else {
                return myN(\x,\z,\k-1)
                +myn(\x,\xvalues[\k],\z);
            };
        };
    },
declare function={myn(\x,\y,\z)=(-(\x-\y)*(\x-\y))/(2*\z*\z) ;
L(\x,\z,\k)=pow(2*pi*\z,-\k/2)*exp(myN(\x,\z,\k));}}

\section*{How to plot sums in Ti\emph{k}Z/pgfplots}

We define the argument of the exponential as
\begin{equation}
 n_k(\theta_1,\theta_2)~=~-\frac{(\theta_1-x_k)^2}{2\theta_2^2}
\end{equation}
and their sum as
\begin{equation}
 N_k(\theta_1,\theta_2)~=~\sum\limits_{\ell=1}^k n_k(\theta_1,\theta_2)\;.
\end{equation}
This means that $N_k$ can be defined recursively as
\begin{equation}
 N_k(\theta_1,\theta_2)~=~N_{k-1}(\theta_1,\theta_2)+n_k(\theta_1,\theta_2)\;,
\end{equation}
and this is the point where the Ti\emph{k}Z library \texttt{math} comes into
play. It allows us to do the recursive deinition. Examples are shown in
Figure~\ref{fig:N_k}.

\begin{figure}[htb]
\centering
\begin{tikzpicture}
\begin{axis}[samples=101,
    use fpu=false,mark=none,
    xlabel=$x$,ylabel=$y$,
    xmin=0, xmax=10,
    domain=0:10,legend pos=south west
    ] 
  \addplot [mark=none] {myN(x,1,1)};
  \addlegendentry{$N_1$}
  \addplot+ [mark=none] {myN(x,1,2)};
  \addlegendentry{$N_2$}
  \addplot+ [mark=none] {myN(x,1,3)};
  \addlegendentry{$N_3$}
\end{axis}
\end{tikzpicture}
\caption{$N_1$, $N_2$ and $N_3$ for $\theta_2=1$ and $\{x_k\}=\{1,2,4\}$.}
\label{fig:N_k}
\end{figure}

\clearpage

Of course, one can then define functions of these sums,
\begin{equation}
L_k(\theta_1,\theta_2)~=~  
\Big(   \dfrac{1}{\sqrt{2\pi\theta_2}} \Big)^{m}\,\exp\Bigl[ 
\dfrac{-\sum_{i=1}^{k}(x_i-\theta_1)^2}{2\theta_2} \Bigr]\;.
\end{equation}
Examples are shown in Figure~\ref{fig:L_k}.

\begin{figure}[htb]
\centering
\begin{tikzpicture}
\begin{axis}[samples=101,
    use fpu=false,mark=none,
    xlabel=$x$,ylabel=$y$,
    xmin=0, xmax=10,
    domain=0:10,legend pos=north east
    ] 
  \addplot [mark=none] {L(x,1,1)};
  \addlegendentry{$L_1$}
  \addplot+ [mark=none] {L(x,1,2)};
  \addlegendentry{$L_2$}
  \addplot+ [mark=none] {L(x,1,3)};
  \addlegendentry{$L_3$}
\end{axis}
\end{tikzpicture}
\caption{$L_1$, $L_2$ and $L_3$ for $\theta_2=1$ and $\{x_k\}=\{1,2,4\}$.}
\label{fig:L_k}
\end{figure}

\end{document}

enter image description here

The second page contains (hopefully) what you are seeking for.

enter image description here

I'd also like to urge you not to confuse TikZ/pgfplots with a computer algegra system. You can do these things, but should not be too surprised if the performance is below the one of, say, Mathematica.

And here is a 3D example, similar to what you do in your MWE.

\documentclass[tikz,border=3.14mm]{standalone}
\usepackage{pgfplots}
\pgfplotsset{compat=1.16}
\usetikzlibrary{math} 
\begin{document}
% based on https://tex.stackexchange.com/a/307032/121799
% and https://tex.stackexchange.com/a/451326/121799
\def\xvalues{{0,1,2,4,5,7}} % notice that the `0` is the 0th entry, which is not used here
\tikzset{evaluate={
        function myN(\x,\z,\k) { % \x = \theta_1 and \z=\theta_2 
            if \k == 1 then {
                return myn(\x,\xvalues[1],\z);
            } else {
                return myN(\x,\z,\k-1)
                +myn(\x,\xvalues[\k],\z);
            };
        };
    },
declare function={myn(\x,\y,\z)=(-(\x-\y)*(\x-\y))/(2*\z*\z) ;
L(\x,\z,\k)=pow(2*pi*\z,-\k/2)*exp(myN(\x,\z,\k));}}

\begin{tikzpicture}
\begin{axis}[use fpu=false,
    grid=both,
    restrict z to domain*=0:1,
    zmin=0,
    colormap/hot,
    %point meta min=-0.2, 
    %point meta max=1,    
    view={20}{20}  %tune here to change viewing angle
    ]

    \addplot3[surf,domain=-1:9,domain y=1:4, samples=25] { L(x, y,4) };
\end{axis}
\end{tikzpicture}
\end{document}

enter image description here

  • 1
    @manooooh Good catch, I have just copied that from the OP. – marmot Oct 19 '18 at 3:54
  • wow, wow.. I am so grateful for you for this MWE with explanation. yes, as you said beyond a point, performance bothers me. So I am also trying via python in parallel. – Parthiban Rajendran Oct 19 '18 at 18:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.