gend()
    >>> ax.set_xlabel("$k$")
    >>> ax.set_title("Negative binomial cumulative distribution function")
    >>> plt.show()

    The negative binomial distribution is also available as
    `scipy.stats.nbinom`. Using `nbdtr` directly can be much faster than
    calling the ``cdf`` method of `scipy.stats.nbinom`, especially for small
    arrays or individual values. To get the same results one must use the
    following parametrization: ``nbinom(n, p).cdf(k)=nbdtr(k, n, p)``.

    >>> from scipy.stats import nbinom
    >>> k, n, p = 5, 3, 0.5
    >>> nbdtr_res = nbdtr(k, n, p)  # this will often be faster than below
    >>> stats_res = nbinom(n, p).cdf(k)
    >>> stats_res, nbdtr_res  # test that results are equal
    (0.85546875, 0.85546875)

    `nbdtr` can evaluate different parameter sets by providing arrays with
    shapes compatible for broadcasting for `k`, `n` and `p`. Here we compute
    the function for three different `k` at four locations `p`, resulting in
    a 3x4 array.

    >>> k = np.array([[5], [10], [15]])
    >>> p = np.array([0.3, 0.5, 0.7, 0.9])
    >>> k.shape, p.shape
    ((3, 1), (4,))

    >>> nbdtr(k, 5, p)
    array([[0.15026833, 0.62304687, 0.95265101, 0.9998531 ],
           [0.48450894, 0.94076538, 0.99932777, 0.99999999],
           [0.76249222, 0.99409103, 0.99999445, 1.        ]])
    