[Tex/LaTex] How to plot the F distribution

pgfplotsplottikz-pgf

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.

Best Answer

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}

Related Solutions

PGFPlots – How to Plot Two Time Series with Bounds Using PGFPlots

This happens because PGFPlots only uses one "stack" per axis: You're stacking the second confidence interval on top of the first. The easiest way to fix this is probably to use the approach described in "Is there an easy way of using line thickness as error indicator in a plot?": After plotting the first confidence interval, stack the upper bound on top again, using stack dir=minus. That way, the stack will be reset to zero, and you can draw the second confidence interval in the same fashion as the first:

\documentclass{standalone}
\usepackage{pgfplots, tikz}

\usepackage{pgfplotstable}

\pgfplotstableread{
temps   y_h y_h__inf    y_h__sup    y_f y_f__inf    y_f__sup    

1   0.237340    0.135170    0.339511    0.237653    0.135482    0.339823    
2   0.561320    0.422007    0.700633    0.165871    0.026558    0.305184    
3   0.694760    0.534205    0.855314    0.074856    -0.085698   0.235411    
4   0.728306    0.560179    0.896432    0.003361    -0.164765   0.171487    
5   0.711710    0.544944    0.878477    -0.044582   -0.211349   0.122184    
6   0.671241    0.511191    0.831291    -0.073347   -0.233397   0.086703    
7   0.621177    0.471219    0.771135    -0.088418   -0.238376   0.061540    
8   0.569354    0.431826    0.706882    -0.094382   -0.231910   0.043146    
9   0.519973    0.396571    0.643376    -0.094619   -0.218022   0.028783    
10  0.475121    0.366990    0.583251    -0.091467   -0.199598   0.016664    
}{\table}

\begin{document}
\begin{tikzpicture}

    \begin{axis}
    % y_h confidence interval
    \addplot [stack plots=y, fill=none, draw=none, forget plot]   table [x=temps, y=y_h__inf]   {\table} \closedcycle;
    \addplot [stack plots=y, fill=gray!50, opacity=0.4, draw opacity=0, area legend]   table [x=temps, y expr=\thisrow{y_h__sup}-\thisrow{y_h__inf}]   {\table} \closedcycle;
    % subtract the upper bound so our stack is back at zero
    \addplot [stack plots=y, stack dir=minus, forget plot, draw=none] table [x=temps, y=y_h__sup] {\table};

    % y_f confidence interval
    \addplot [stack plots=y, fill=none, draw=none, forget plot]   table [x=temps, y=y_f__inf]   {\table} \closedcycle;
    \addplot [stack plots=y, fill=gray!50, opacity=0.4, draw opacity=0, area legend]   table [x=temps, y expr=\thisrow{y_f__sup}-\thisrow{y_f__inf}]   {\table} \closedcycle;

    % the line plots (y_h and y_f)    
    \addplot [stack plots=false, very thick,smooth,blue]  table [x=temps, y=y_h]   {\table};
    \addplot [stack plots=false, very thick,smooth,blue]  table [x=temps, y=y_f]   {\table};
    \end{axis}

\end{tikzpicture}
\end{document}

PGFPlots – Plotting the Probability Density Function of the Gamma Distribution

Here's the gamma PDF using the gamma function approximation that Wolfram Alpha spits out:

\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={gammapdf(\x,\k,\theta) = \x^(\k-1)*exp(-\x/\theta) / (\theta^\k*gamma(\k));}
]

\begin{axis}[
    axis lines=left,
    enlargelimits=upper,
    samples=50,
    legend entries={$k=1\quad \theta=2$,$k=2\quad \theta=2$, $k=9\quad \theta=0.5$}
]
\addplot [smooth, domain=0:20] {gammapdf(x,1,2)};
\addplot [smooth, domain=0:20, red] {gammapdf(x,2,2)};
\addplot [smooth, domain=0:20, blue] {gammapdf(x,9,0.5)};
\end{axis}
\end{tikzpicture}
\end{document}

Best Answer

Related Solutions

PGFPlots – How to Plot Two Time Series with Bounds Using PGFPlots

PGFPlots – Plotting the Probability Density Function of the Gamma Distribution

Related Question