8

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.

  • 1
    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
8

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
7

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}

enter image description here

  • 1
    Oh well that was easy! :) – azetina Sep 25 '14 at 15:01
  • 1
    Nice output. And how label the axes with pspicture option? Thanks. – Toño Sep 25 '14 at 19:29
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}

enter image description here

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.