6

I am trying to draw a diagram that has a negative-binomial distribution in it. I am trying to do this in Tikz. Here is a minimal working version of that code:

\documentclass[landscape]{article}
\usepackage{pgfplots}
\pagestyle{empty}
\begin{document}

\pgfmathdeclarefunction{nbinom}{2}{%
  \pgfmathparse{%
((x+#1-1)!/(x!*(#1-1)!))*(1-#2)^#1*(#2^x)}%
}

\begin{tikzpicture}

\begin{axis}[
  no markers, domain=0:15, samples=100,
  axis lines*=left, xlabel=Foo,
  height=5cm, width=12cm,
  ytick=\empty,
  enlargelimits=false, clip=false, axis on top
  ]
  \addplot [very thick,black] {nbinom(2,.6)};
\end{axis}

\end{tikzpicture}
\end{document}

Running this produces the following: enter image description here

This is wrong. For comparison, here is a plot made in R, using the same formula:

xn <- 0:15
rn <- 2
pn <- .6
hxn <- (factorial(xn+rn-1)/(factorial(xn)*factorial(rn-1)))*(1-pn)^rn*pn^xn
plot(xn, hxn, type="l", ylab="", xlab="Foo", yaxt='n', xaxt='n')

This yields the following:

enter image description here

...which is correct.

I am pretty sure that this issue is stemming from the use of the factorial operator in the "nbinom" function I call with \pgfmathdeclarefunction. For example, if I set up the same minimal working version with a normal distribution...

\documentclass[landscape]{article}
\usepackage{pgfplots}
\pagestyle{empty}
\begin{document}

\pgfmathdeclarefunction{gauss}{2}{%
  \pgfmathparse{1/(#2*sqrt(2*pi))*exp(-((x-#1)^2)/(2*#2^2))}%
}

\begin{tikzpicture}

\begin{axis}[
  no markers, domain=0:10, samples=100,
  axis lines*=left, xlabel=Foo,
  height=5cm, width=12cm,
  ytick=\empty,
  enlargelimits=false, clip=false, axis on top
  ]
  \addplot [very thick,black] {gauss(5,1)};
\end{axis}

\end{tikzpicture}
\end{document}

...it works fine:

enter image description here

Yes, I could just hack this diagram together in R. But I know I can do a better version, in principle, in Tikz. This is just killing me. Has anyone encountered this problem, and do they have any ideas?

0

2 Answers 2

10

Your analysis is right, factorial or ! only support integer numbers, the decimal part of real numbers is truncated.

The following example defines a function facreal, which uses the Gamma function to calculate the factorials with real numbers as argument, see extensions of the factorial in Wikipedia and approximations for the Gamma function (Wikipedia).

\documentclass[landscape]{article}
\usepackage{pgfplots}
%\pgfplotsset{compat=1.12}
\pagestyle{empty}

\pgfmathsetmacro\twopi{2*pi}

\pgfmathdeclarefunction{lngamma}{1}{%
  \pgfmathsetmacro\lngammatmp{#1*#1*#1}%
  \pgfmathparse{%
    #1*ln(#1) - #1 - .5*ln(#1/\twopi)
    + 1/12/#1 - 1/360/\lngammatmp + 1/1260/\lngammatmp/#1/#1
  }%
}   
\pgfmathdeclarefunction{facreal}{1}{%
  \pgfmathparse{exp(lngamma(#1+1))}% 
}

\pgfmathdeclarefunction{nbinom}{2}{%
  \pgfmathparse{%
    facreal((x+#1-1)/(facreal(x)*facreal(#1-1)))*(1-#2)^#1*(#2^x)%
  }%
}

\begin{document}
\begin{tikzpicture}

\begin{axis}[
  no markers, domain=0:15, samples=100,
  axis lines*=left, xlabel=Foo,
  height=5cm, width=12cm,
  ytick=\empty,
  enlargelimits=false, clip=false, axis on top
  ]
  \addplot [very thick,black] {nbinom(2,.6)};
\end{axis}

\end{tikzpicture}
\end{document}

Result

A faster version with more sample points:

\documentclass[landscape]{article}
\usepackage{pgfplots}
%\pgfplotsset{compat=1.12}
\pagestyle{empty}

% Calculation, see
% https://en.wikipedia.org/wiki/Gamma_function#Approximations
\makeatletter
\pgfmath@def{lngamma}{A}{0.159154943092}% 1/(2*pi)
\pgfmath@def{lngamma}{B}{0.0833333333333}% 1/12
\pgfmath@def{lngamma}{C}{0.00277777777778}% 1/360
\pgfmath@def{lngamma}{D}{0.000793650793651}% 1/1260
\pgfmathdeclarefunction{lngamma}{1}{%
  \pgfmathmultiply{#1}{#1}%
  \let\pgfmath@lngamma@tmp\pgfmathresult
  % tmp = x^2
  \pgfmathdivide\pgfmath@lngamma@D\pgfmath@lngamma@tmp
  % result = 1/(1260 x^2)
  \pgfmathsubtract\pgfmathresult\pgfmath@lngamma@C
  % result = -1/360 + 1/(1260 x^2)
  \pgfmathdivide\pgfmathresult\pgfmath@lngamma@tmp
  % result = -1/(360 x^2) + 1/(1260 x^4)
  \pgfmathadd\pgfmathresult\pgfmath@lngamma@B
  % result = 1/12 - 1/(360 x^2) + 1/(1260 x^4)
  \pgfmathdivide\pgfmathresult{#1}%
  \let\pgfmath@lngamma@sum\pgfmathresult
  % sum = 1/(12 x) - 1/(360 x^3) + 1/(1260 x^5)
  \pgfmathmultiply{#1}\pgfmath@lngamma@A
  % result = x/(2 pi)
  \pgfmathln\pgfmathresult
  % result = ln(x/(2 pi))
  \pgfmathmultiply\pgfmathresult{.5}%
  % result = (1/2) ln(x/(2 pi))
  \pgfmathadd\pgfmathresult{#1}%
  \let\pgfmath@lngamma@tmp\pgfmathresult
  % tmp = x + (1/2) ln(x/(2 pi))
  \pgfmathln{#1}%
  % result = ln(x)
  \pgfmathmultiply\pgfmathresult{#1}%
  % result = x ln(x)
  \pgfmathsubtract\pgfmathresult\pgfmath@lngamma@tmp
  % result = x ln(x) - x - (1/2) ln(x/(2 pi))
  \pgfmathadd\pgfmathresult\pgfmath@lngamma@sum
  % result = x ln(x) - x - (1/2) ln(x/(2 pi))
  %          + 1/(12 x) - 1/(360 x^3) + 1/(1260 x^5)
}
\makeatother

\pgfmathdeclarefunction{facreal}{1}{%
  \pgfmathadd{#1}{1}%
  \pgfmathlngamma\pgfmathresult
  \pgfmathexp\pgfmathresult
}

\pgfmathdeclarefunction{nbinom}{2}{%
  \pgfmathparse{%
    facreal((x+#1-1)/(facreal(x)*facreal(#1-1)))*(1-#2)^#1*(#2^x)%
  }%
}
\begin{document}    
\begin{tikzpicture}

\begin{axis}[
  no markers, domain=0:15, samples=500,
  axis lines*=left, xlabel=Foo, 
  height=5cm, width=12cm,  
  ytick=\empty,
  enlargelimits=false, clip=false, axis on top
  ]
  \addplot [very thick,black] {nbinom(2,.6)};
\end{axis}

\end{tikzpicture}
\end{document}

Result with 500 samples

3
  • Impressive math knowledge :) Commented May 23, 2015 at 11:52
  • 3
    I have made a feature request for the functions gamma and factorial for real numbers. Commented May 23, 2015 at 12:30
  • Oh, bravo! My suspicion was that the problem lay in using the factorial operator on non-integers. I say that, obviously, because the spikes were at the integers. It seemed that the function was using a constant value for the factorial, even as it increased x in the exponent elsewhere in the function. And I vaguely remember learning that the gamma could be used for a real-number generalization of factorials, but hadn't put the pieces together. Much obliged! Commented May 23, 2015 at 17:24
1

A MetaPost solution, for whom it may interest, inserted in a LuaLaTeX program.

The factorial function for real values is defined, as usual, with the Gamma function. I've picked up the MetaPost definition of this function from Anthony Phan's mps package, which I've just discovered. (Not available on CTAN yet.) See also this more detailed presentation of mps.

\documentclass[border=2mm]{standalone}
\usepackage{luamplib}
  \mplibsetformat{metafun}
  \mplibtextextlabel{enable}
  \mplibnumbersystem{double}
\begin{document}
  \begin{mplibcode}
    % Definitions of lnGamma and Gamma
    % from Anthony Phan's mps-math package
    vardef lnGamma(expr x) =
      mlog(x+4.5)/256*(x-.5)-x-4.5
      +mlog(2.50662827465*(1+76.18009173/x
        -86.50532033/(x+1)
        +24.01409822/(x+2)
        -1.231739516/(x+3)
        +0.00120858003/(x+4)
        -0.00000536382/(x+5)))/256
    enddef;

    vardef Gamma(expr x) = mexp(256*lnGamma(x)) enddef;

    vardef factorial(expr x) = Gamma(x+1) enddef;  

    vardef nbinom(expr a, b)(expr x) =
      factorial(x+a-1) / (factorial(x)*factorial(b-1)) * ((1-b)**a)*(b**x)
    enddef;

    u = .75cm; v = 40cm;
    beginfig(1);
      draw function(2, "x", "nbinom(2,.6)(x)", 0, 15, .1) xyscaled (u, v)
        withpen pencircle scaled 1.5bp;
      draw origin -- (15u, 0);
      draw origin -- (0, .15v);
      label.bot("$0$", origin);
      labeloffset := 6bp;
      for i = 1 upto 15:
        draw (i*u, -3bp) -- (i*u, 3bp);
        label.bot("$" & decimal i & "$", (i*u, 0));
      endfor
    endfig;
  \end{mplibcode}
\end{document}

enter image description here

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .