# How can I plot the F distribution

I wish to plot the F distribution by declaring its density function. However with the following code I get an error saying that fdst has not been declared:

\documentclass[border=5mm]{standalone}
\usepackage{pgfplots}

\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={fdst(\x,\n,\m))= (gamma((\n+\m)/2)/(gamma(\n/2) *gamma(\m/2))) *(\n/\m)^(\n/2) *((\x *((\n-2)/2))/((1+(\n *\x)/\m)^((\n+\m)/2)));}
]

\begin{axis}[
axis lines=left,
enlargelimits=upper,
samples=50
]

\end{axis}
\end{tikzpicture}
\end{document}


Any help will be much appreciated.

• There are various parenthesis missing or additional in your input. I'm not sure what the equations should read, but if I at least make things match then the example will compile to something. – Joseph Wright Sep 25 '14 at 14:27
• As per Joseph Wright's comment if you fix the mismatching round brackets than a plot gets produced. One example is fdst(\x,\n,\m))= where you have an extra ). With this change things compile even though there are still other mismatches lingering. For example I am unable to locate the opening ( to match the closing one in )*exp. – Peter Grill Sep 25 '14 at 15:01
• Both are right with the mismatches. The p.d.f. is the standard definition with the Gamma function, but the output looks strange. The option below with the Beta function looks great. Thank you for your answers. – Toño Sep 25 '14 at 19:34

You can use the approximation of the gamma function used in Plot the probability density function of the gamma distribution, use that to define the beta function and then define the f distribution in terms of that:

\documentclass[border=5mm]{standalone}
\usepackage{pgfplots}

\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={
beta(\x,\y)=gamma(\x)*gamma(\y)/gamma(\x+\y);
},
declare function={
fdst(\x,\a,\b) = 1 / beta(\a/2, \b/2) * (\a/\b)^(\a/2) * \x^(\a/2-1) * (1 + \a/\b*\x)^(-(\a + \b)/2);
}
]

\begin{axis}[
axis lines=left,
enlargelimits=upper,
samples=100,
xmin=0, ymin=0,
domain=0.01:4,
legend cell align=left
]

\addplot [very thick,blue] {fdst(x,1,1)}; \addlegendentry{$d_1=1,\hphantom{00} d_2=1$}
\addplot [very thick,orange] {fdst(x,100,100)}; \addlegendentry{$d_1=100, d_2=100$}
\addplot [very thick,purple] {fdst(x,5,2)}; \addlegendentry{$d_1=5,\hphantom{00} d_2=2$}

\end{axis}
\end{tikzpicture}
\end{document}

• Excellent solution using the Beta function as well. It is for sure more elegant, and it looks like the F distribution. – Toño Sep 25 '14 at 19:37

Run with xelatex:

\documentclass{article}
\usepackage{pst-func}
\begin{document}

\psset{xunit=2cm,yunit=10cm,plotpoints=100}
\begin{pspicture*}(-0.5,-0.07)(5.5,0.8)
\psline[linestyle=dashed](0.5,0)(0.5,0.75)
\psline[linestyle=dashed](! 2 7 div 0)(! 2 7 div 0.75)
\psset{linewidth=1pt}
\psFDist{0.1}{5}
\psFDist[linecolor=red,nue=3,mue=12]{0.01}{5}
\psFDist[linecolor=blue,nue=12,mue=3]{0.01}{5}
\psaxes[Dy=0.1,ticksize=-4pt 0]{->}(0,0)(5,0.75)
\end{pspicture*}

\end{document}


• Oh well that was easy! :) – azetina Sep 25 '14 at 15:01
• Nice output. And how label the axes with pspicture option? Thanks. – Toño Sep 25 '14 at 19:29

Just an extension to Jake's answer by implementing the functions in lua. There is some fussing around to convert to and from pgfplots numerical representation (which probably isn't quite correct). Also a different approximation for the beta function is used:

\documentclass[tikz,border=5]{standalone}
\usepackage{pgfplots}
{\catcode\~=11 \gdef\tilde{~}}
{\catcode\%=11 \gdef\percent{%}}

\directlua{%

function lua2pgfplots(v)
if v \tilde= v then
return "3Y0e0]"
end
if v == math.huge then
return "4Y0e0]"
end
if v == -math.huge then
return "5Y0e0]"
end
if v == 0 then
return "0Y0e0]"
end
if v > 0 then
return string.format("1Y\percent e]", v)
else
return string.format("2Y\percent e]", v)
end
end

function pgfplots2lua(v)
local f, x
f, x = string.match(v, "(\percent d)Y(.*).")
f = tonumber(f)
x = tonumber(x)
if f == 1 then
return x
end
if f == 2 then
return -x
end
if f == 3 then
return 0/0
end
if f == 4 then
return math.huge
end
if f == 5 then
return -math.huge
end
return x
end
}

\directlua{
function beta(x,y)
return math.sqrt(2*math.pi)*(math.pow(x,x-.5)*math.pow(y,y-.5)) /
math.pow(x+y, x+y-.5)
end

function fdist(x, d1, d2)
return 1/beta(d1/2,d2/2)*math.pow(d1/d2, d1/2) *
math.pow(x, d1/2-1) *
math.pow(1+d1/d2*x, -(d1+d2)/2)
end
}

\pgfmathdeclarefunction{fdist}{3}{%
\edef\pgfmathresult{%
\directlua{%
local f
f = fdist(pgfplots2lua("#1"),pgfplots2lua("#2"),pgfplots2lua("#3"))
f = lua2pgfplots(f)
tex.print(f)
}%
}%
}
\begin{document}

\begin{tikzpicture}
\begin{axis}[
axis lines=left,
enlargelimits=upper,
samples=100,
xmin=0, ymin=0,
domain=0.01:4,
legend cell align=left
]

\addplot [very thick, red]   {fdist(x,1,1)};  \addlegendentry{$d_1=1,\hphantom{00} d_2=1$}
\addplot [very thick, black] {fdist(x,2,1)};  \addlegendentry{$d_1=2,\hphantom{00} d_2=1$}
\addplot [very thick, blue]  {fdist(x,5,2)};  \addlegendentry{$d_1=5,\hphantom{00} d_2=2$}
\addplot [very thick, green] {fdist(x,10,1)}; \addlegendentry{$d_1=10,\hphantom{0} d_2=1$}
\addplot [very thick, gray]  {fdist(x,100,100)}; \addlegendentry{$d_1=100, d_2=100$}

\end{axis}
\end{tikzpicture}
\end{document}