# Gaussian kernel density estimation with data from file

I am trying to draw a histogram next to a density function, both with data from a file. The histogram is already working:

\documentclass[tikz,border=3.14mm]{standalone}
\usepackage{pgfplots}

\begin{filecontents}{example.dat}
71
54
55
54
98
76
93
95
86
88
68
68
50
61
79
79
73
57
56
57
97
80
91
94
85
88
45
58
78
81
74
60
57
58
95
81
\end{filecontents}

\begin{document}
\begin{tikzpicture}
\begin{axis}[ybar, ymin=0]

hist={
density, % <-- EDIT
bins=11
}] table [y index=0] {example.dat};

\end{axis}
\end{tikzpicture}

\end{document}


My question is with the density function. I would like to draw it using this formula (kernel density estimation):

EDIT start

where:

• n: Number of datapoints

• sigma: Standard deviation. The value of it is chosen (by me), so that the resulting curve has a specific number of local maxima. This is not part of the question and you can just use some fixed number for it, so that it looks smooth.

• x_i: Datapoint i

• x: Function input variable (f(x))

EDIT end

I can't just plot it using \addplot... because f(x) depends on all datapoints x_i.

I was thinking about using something like that somewhere:

\pgfplotstableread{example.dat}\table
\pgfplotstablegetrowsof{\table}
\pgfmathsetmacro{\R}{\pgfplotsretval-1}

\pgfplotsinvokeforeach{0,...,\R}{
\pgfplotstablegetelem{#1}{0}\of{\table}
\pgfmathsetmacro \value {\pgfplotsretval}
% sum up all e^0.5(\value-x)/sigma somhow
}


But I couldn't find a way to define a variable where I can add values in every iteration.

Here is an image from Wikipedia on Kernel density estimation:

The blue curve on the right is kind of what I would like to draw. What's the best way of achieving that?

You can sum these things up as follow. I use \pgfplotsforeachungrouped in order to avoid making the variables global. The following uses your sigma and your normalized Gaussian, and there is a factor of 5 to account for the bar width.

\documentclass[tikz,border=3.14mm]{standalone}
\usepackage{pgfplots}
\pgfplotsset{compat=1.16}
\begin{filecontents*}{example.dat}
71
54
55
54
98
76
93
95
86
88
68
68
50
61
79
79
73
57
56
57
97
80
91
94
85
88
45
58
78
81
74
60
57
58
95
81
\end{filecontents*}

\begin{document}
\begin{tikzpicture}
\pgfplotstablegetrowsof{\datatable}
\pgfmathsetmacro{\R}{\pgfplotsretval-1}
\pgfmathsetmacro\mysum{0}
\pgfmathsetmacro\mysigma{8}
\pgfplotsforeachungrouped \X in {0,...,\R}{
\pgfplotstablegetelem{\X}{0}\of{\datatable}
\edef\mysum{\mysum+(5/(sqrt(2*pi)*\mysigma))*exp(-(x-\pgfplotsretval)^2/(2*\mysigma*\mysigma))}
}

\begin{axis}[ ymin=0]

hist={
bins=11
}] table [y index=0] {example.dat};
\end{axis}
\end{tikzpicture}

\end{document}


OLDER:

\documentclass[tikz,border=3.14mm]{standalone}
\usepackage{pgfplots}
\pgfplotsset{compat=1.16}
\begin{filecontents*}{example.dat}
71
54
55
54
98
76
93
95
86
88
68
68
50
61
79
79
73
57
56
57
97
80
91
94
85
88
45
58
78
81
74
60
57
58
95
81
\end{filecontents*}

\begin{document}
\begin{tikzpicture}
\pgfplotstablegetrowsof{\datatable}
\pgfmathsetmacro{\R}{\pgfplotsretval-1}
\pgfmathsetmacro\mysum{0}
\pgfplotsforeachungrouped \X in {0,...,\R}{
\pgfplotstablegetelem{\X}{0}\of{\datatable}
\edef\mysum{\mysum+2*exp(-(x-\pgfplotsretval-0.5)^2)}
% sum up all e^0.5(\value-x)/sigma somhow
}

\begin{axis}[ ymin=0]

hist={
bins=11
}] table [y index=0] {example.dat};
\end{axis}
\end{tikzpicture}

\end{document}


If you use

\edef\mysum{\mysum+sqrt(2)*exp(-0.25*(x-\pgfplotsretval-0.5)^2)}


OLD ANSWER: I am not sure I got the normalization of the Gaussian right.

\documentclass[tikz,border=3.14mm]{standalone}
\usepackage{pgfplots}
\pgfplotsset{compat=1.16}
\begin{filecontents*}{example.dat}
71
54
55
54
98
76
93
95
86
88
68
68
50
61
79
79
73
57
56
57
97
80
91
94
85
88
45
58
78
81
74
60
57
58
95
81
\end{filecontents*}

\begin{document}
\begin{tikzpicture}
\pgfplotstablegetrowsof{\datatable}
\pgfmathsetmacro{\R}{\pgfplotsretval-1}
\pgfmathsetmacro\mysum{0}
\pgfplotsforeachungrouped \X in {0,...,\R}{
\pgfplotstablegetelem{\X}{0}\of{\datatable}
\pgfmathsetmacro\mysum{\mysum+\pgfplotsretval}
% sum up all e^0.5(\value-x)/sigma somhow
}
\pgfmathsetmacro{\myaverage}{\mysum/\R}
\pgfmathsetmacro\mysigma{0}
\pgfplotsforeachungrouped \X in {0,...,\R}{
\pgfplotstablegetelem{\X}{0}\of{\datatable}
\pgfmathsetmacro\mysigma{\mysigma+pow(\pgfplotsretval-\myaverage,2)}

}
%\typeout{\mysum,\myaverage,\mysigma}

\begin{axis}[ ymin=0]

hist={
bins=11
}] table [y index=0] {example.dat};
\end{axis}
\end{tikzpicture}

\end{document}


• Thanks for your answer, but I think you misunderstood my question a bit. Maybe it wasn't clear enough. I want to draw a sum of gaussians (red in my last image) to get to the blue one. Sigma is a parameter I provide. But I think I can come up with a solution with the help of your answer! – Nathanael May 16 at 16:11
• @Nathanael I changed my answer. Do you have a concrete prescriptionor explanation for all the elements of your formula? Or is sigma just the standard variance? – marmot May 16 at 16:32
• Thanks. Looks good now. Will update my question when I get home in a couple of hours. Sigma is just the standard deviation of those gaussians. By changing it one can control the "smoothness" of the resulting curve. – Nathanael May 16 at 16:58
• Edited my q. and had a closer look at your edit. I added \pgfmathsetmacro{\sigma}{8} and changed your code a bit: \edef\mysum{\mysum+2*exp(-(x-\pgfplotsretval-0.5)^2)} -> \edef\mysum{\mysum+1/(sqrt(2*pi)*\sigma*\R)*exp(-1/(2*\sigma*\sigma)*(x-\pgfplotsretval)^2)}. By adding density to hist={...} both the histogram and the KDE-function fit on one scale and the KDE-function is an actual density function (Integral == 1). My original histogram (w.o. density) is probably the reason you omitted the 1/(sqrt(2*pi)*\sigma*\R) part, right? I hope my scaling is correct now... – Nathanael May 17 at 0:19
• @Nathanael Yes, I feel that there is not yet a clear prescription of what should be done. One approach could be to demand that the areas under the Gaussian and the bars should coincide. However, this does not yet tell us uniquely what the widths of the Gaussians should be. You could set it to some value that looks good, which is a nice option. There might be some other preferred choice, and it could be that different users will have their own prescription. If you tell me what a preferred prescription is I will be happy to update my answer. – marmot May 17 at 0:29
\documentclass{article}
\begin{filecontents}{example.dat}
71
54
.
.
.
95
81
\end{filecontents}
\begin{document}
<<echo=F,fig.cap="Histogram and density plot.">>=
hist(data$V1, freq=F, col="gray", main="", xlab="Example data") lines(density(data$V1),col="blue",lwd=3)
@
\end{document}


Of course, you can have some control over the density function, for instance with:

lines(density(data\$V1,adjust=.5, bw=8),col="blue",lwd=3)


The result will be ...

• Thanks for your answer. This doesn't look like normal TikZ/PGF though, are you calling Gnuplot or something? – Nathanael May 16 at 22:41
• @Nathanael Hot. It is R code in a noweb file (.Rnw). The R package "knitr" convert this in true LaTex (.tex) file than can be compiled by pdflatex (edited with RStudio, all the process is just click on "compile PDF"). With the tikzDevice R package installed and dev='tikz'in the chunk options (not showed), the graph will look just as a original TikZ plot, using the current LaTeX font (but imho look better with a different font, as showed here, using the R default font). – Fran May 17 at 9:26
• Here there are an making a R plots with Tikz device. – Fran May 17 at 9:38
• I am trying this approach atm, since my actual file has 16k datapoints and while marmots answer works, it takes ages to finish (Running for over 1h now). I was able to reproduce your output and changed the font as you suggested. I have a followup questions now: Is it possible to make the output look exactly like pgfplots output (Axes, rotated ticklabels etc...)? I would like all my figures to have the same stlye/feel. Thanks. – Nathanael May 17 at 11:44
• @Nathanael Not sure if you can reproduce every possible pgf plot 100% exactly in R but of course you can rotate axis labels and change almost any aspect of a R plots, beside that there are tons of plot types (a sample here). – Fran May 17 at 12:25