How can I plot the F distribution

Question

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
]
 \addplot [very thick,yellow!80!black] {fdst(x,3,3)};

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

Any help will be much appreciated.

Answer 1

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}

Answer 2

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}

Answer 3

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}