在 TikZ 中使用 \pgfmathdeclarefunction 和阶乘运算符进行绘图

Question

我正在尝试绘制一个包含负二项分布的图表。我正尝试在 Tikz 中执行此操作。以下是该代码的最小工作版本：

\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}

运行该程序会产生以下结果：

这是错误的。为了进行比较，下面是使用相同公式在 R 中绘制的图：

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')

得出以下结论：

...哪个是对的。

我非常确定这个问题是由于我在使用 \pgfmathdeclarefunction 调用的“nbinom”函数中使用了阶乘运算符而引起的。例如，如果我设置了具有正态分布的相同最小工作版本...

\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}

...它工作正常：

是的，我可以在 R 中拼凑这个图表。但我知道我原则上可以在 Tikz 中制作一个更好的版本。这真是让我头疼。有人遇到过这个问题吗？他们有什么想法吗？

Answer 1

你分析的是对的，factorial或者!只支持整数，实数的小数部分被截断了。

下面的例子定义了一个函数facreal，它使用 Gamma 函数以实数作为参数来计算阶乘，请参阅阶乘的扩展维基百科以及 Gamma 函数的近似值（维基百科）。

\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}

具有更多采样点的更快版本：

\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}

Answer 2

一个可能感兴趣的人会感兴趣的 MetaPost 解决方案，插入到 LuaLaTeX 程序中。

实数阶乘函数通常用 Gamma 函数定义。我从 Anthony Phan 的 MetaPost 中找到了该函数的定义mps包裹，这是我刚刚发现的。（CTAN 上尚未提供。）另请参阅更详细的介绍的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}