I want to combine a histogram with a probability density function. Following the advise of Jake

My plan now is to combine the two in one plot, like the following from an R example by Cai or the one by Dirk:

I don't need to fit the distribution, I already extracted the estimates using EasyFit. I only want to draw a histogram together with the fitted function gamma(1.7435, 21.263).

I read the gnuplot example, again from Jake, but couldn't figure out if this is the only solution. Especially if I have the function already.

Can you help me out?

MWE - What I have so far:

\documentclass{article}
\usepackage{pgfplots}
\usepackage{subcaption}
\begin{document}
\begin{figure}[htbp]
\centering
\begin{subfigure}{0.35\textwidth}
\centering
\begin{tikzpicture}
\usepgfplotslibrary{statistics}
\begin{axis}[
axis x line=middle,
axis y line=middle,
ybar,
ymin=0
]
\addplot +[hist] table [y index=0, col sep=comma,row sep=newline] {
data
48
51
48
51
66
0
0
25
31
60
9
43
49
15
22
1
};
\end{axis}
\end{tikzpicture}
\subcaption[subfigcapskip = 50pt]{Histogram}
\label{fig:histogram}
\end{subfigure}%
\hfill
\begin{subfigure}{0.35\textwidth}
\centering

\begin{tikzpicture}[
declare function={gamma(\z)=
2.506628274631*sqrt(1/\z)+ 0.20888568*(1/\z)^(1.5)+ 0.00870357*(1/\z)^(2.5)- (174.2106599*(1/\z)^(3.5))/25920- (715.6423511*(1/\z)^(4.5))/1244160)*exp((-ln(1/\z)-1)*\z;},
declare function={gammapdf(\x,\k,\theta) = 1/(\theta^\k)*1/(gamma(\k))*\x^(\k-1)*exp(-\x/\theta);}
]

\begin{axis}[
axis lines=left,
]
\addplot [smooth, domain=0:66, red] {gammapdf(x,1.7435, 21.263)};
\end{axis}
\end{tikzpicture}
\subcaption[subfigcapskip = 50pt]{Plot}
\label{fig:plot}
\end{subfigure}
\end{figure}
\end{document}


Edit - Final Solution, thanks to Enthusiastic Student!

\documentclass[margin=1mm]{standalone}
\usepackage{filecontents}
\usepackage{pgfplots}
\usepgfplotslibrary{statistics}

\begin{document}

\begin{tikzpicture}[
declare function={gamma(\z)= 2.506628274631*sqrt(1/\z)+ 0.20888568*(1/\z)^(1.5)+ 0.00870357*(1/\z)^(2.5)- (174.2106599*(1/\z)^(3.5))/25920- (715.6423511*(1/\z)^(4.5))/1244160)*exp((-ln(1/\z)-1)*\z;},
declare function={gammapdf(\x,\k,\theta) = 1/(\theta^\k)*1/(gamma(\k))*\x^(\k-1)*exp(-\x/\theta);},
declare function={cauchypdf(\x,\mu,\gamma) = 1/(pi*\gamma*(1+((\x-\mu)/\gamma)^2));}
]
\pgfplotsset{set layers}
\begin{axis}[
axis y line*=left,
ylabel=Histogram,
every axis y label/.style={rotate=90, blue, at={(-0.15,0.5)},},
ybar,
ymin=0,
yticklabel style={blue}
]

\addplot +[hist] table [y index=0, col sep=comma,row sep=newline] {
data
48
51
48
51
66
0
0
25
31
60
9
43
49
15
22
1
};
\end{axis}

\begin{axis}[
ylabel style={black},
axis y line*=right,
axis x line=none,
ylabel=PDF's,
every axis y label/.style={rotate=-90, black, at={(1.15,0.5)},},
/pgfplots/every y tick scale label/.style={
at={(1,1)},
above right,
inner sep=0pt,
yshift=0.3em,
},
yticklabel style={black},
legend entries={
$\textbf{Gamma}~\alpha=1.7435~ \beta=21.263$,
$~\textbf{Cauchy}~\sigma=16.1340~ \mu=39.554$
},
legend style={at={(1,-0.15)}},
ymin=0,
]

\addplot [thick, smooth, domain=0:66, black, mark=square] {gammapdf(x,1.7435, 21.263)};
\addplot [thick, smooth, domain=0:66, red, mark=triangle] {cauchypdf(x,39.554, 16.134)};
\end{axis}
\end{tikzpicture}
\end{document}

• How should the y-axis labels be? Because the value of the gamma function is about a hundred times smaller than the histogram and when you plot them both together, the plotted function will disappear. Commented Nov 10, 2014 at 14:12
• Hm .. good question. I would say it has to be normed. Sadly I don't know how to do this.
– lony
Commented Nov 10, 2014 at 14:18
• I am just working on it and will post some answer. Commented Nov 10, 2014 at 14:18
• You can justify center your legend by having this code: legend style={at={(0.5,-0.15)},anchor=north}, Commented Nov 10, 2014 at 17:14
• I asked this question based on this question of yours. tex.stackexchange.com/questions/211519/… Commented Nov 10, 2014 at 17:40

\documentclass[margin=1mm]{standalone}
\usepackage{filecontents}
\usepackage{pgfplots}
\usepgfplotslibrary{statistics}

\begin{filecontents}{data.csv}
data
48
51
48
51
66
0
0
25
31
60
9
43
49
15
22
1
\end{filecontents}

\begin{document}

\begin{tikzpicture}[
declare function={gamma(\z)=
2.506628274631*sqrt(1/\z)+ 0.20888568*(1/\z)^(1.5)+ 0.00870357*(1/\z)^(2.5)- (174.2106599*(1/\z)^(3.5))/25920- (715.6423511*(1/\z)^(4.5))/1244160)*exp((-ln(1/\z)-1)*\z;},
declare function={gammapdf(\x,\k,\theta) = 1/(\theta^\k)*1/(gamma(\k))*\x^(\k-1)*exp(-\x/\theta);},
]
\pgfplotsset{set layers}
\begin{axis}[
axis y line*=left,
ylabel=Histogram,
every axis y label/.style={rotate=90, blue, at={(-0.15,0.5)},},
ybar,
ymin=0,
yticklabel style={blue}
]
\addplot +[hist] table [y index=0, col sep=comma,row sep=newline] {data.csv};
\end{axis}
\begin{axis}[
ylabel style={red},
axis y line*=right,
axis x line=none,
ylabel=Function Plot,
every axis y label/.style={rotate=-90, red, at={(1.15,0.5)},},
/pgfplots/every y tick scale label/.style={
at={(1,1)},above right,inner sep=0pt,yshift=0.3em},
yticklabel style={red}
]
\addplot [smooth, domain=0:66, red] {gammapdf(x,1.7435, 21.263)};
\end{axis}
\end{tikzpicture}
\end{document}

• Really, really, really nice! I can even add more pdf's to compare them. Thank you!
– lony
Commented Nov 10, 2014 at 15:20
• @lony I edited to show the y labels. Commented Nov 10, 2014 at 15:20
• @lony If you are printing your plots in black and white, I recommend you to add some legends to the plot and use the line styles too. Commented Nov 10, 2014 at 15:24