Special functions

Clean-room implementations of the numerics the families need. Every one is validated against mpmath at 60 digits as well as against R.

rcopula.special.stable

Samplers for the frailty (mixing) distributions of Archimedean copulas.

Every Archimedean copula whose generator is completely monotone can be written as a frailty mixture: draw :math:V \sim F, draw iid :math:E_j \sim \mathrm{Exp}(1), and set :math:U_j = \psi(E_j / V). The resulting :math:U has copula :math:C with generator :math:\psi = \mathcal{L}[F]. This is the Marshall-Olkin algorithm, and it is how rcopula samples every Archimedean family in :math:d > 2.

Marshall, A. W. and Olkin, I. (1988). Families of multivariate
distributions. *JASA* 83(403), 834-841.

Each family needs its own :math:V_0:

======== ==================================== ============================== Family :math:V_0 Implemented by ======== ==================================== ============================== Clayton Gamma(1/theta, 1) numpy directly Gumbel positive stable S(1/theta, 1, ...) :func:rstable_positive Frank logarithmic series, p = 1 - e^-theta :func:rlog_series Joe Sibuya(1/theta) :func:rsibuya AMH Geometric(1 - theta) numpy directly ======== ==================================== ==============================

References

Chambers, J. M., Mallows, C. L. and Stuck, B. W. (1976). A method for simulating stable random variables. JASA 71(354), 340-344. The CMS algorithm used by :func:rstable_positive. Kemp, A. W. (1981). Efficient generation of logarithmically distributed pseudo-random variables. Applied Statistics 30(3), 249-253. The "LK" algorithm used by :func:rlog_series. Hofert, M. (2008). Sampling Archimedean copulas. Computational Statistics & Data Analysis 52(12), 5163-5174. Sibuya sampling and the family-to-frailty correspondence.

sinc

sinc(x: NDArray[float64] | float) -> NDArray[np.float64]

Compute sin(x) / x safely, returning 1 at x = 0 instead of dividing by zero.

This is the unnormalised sinc function (note: numpy.sinc is the normalised sin(pi x) / (pi x)). It appears in Zolotarev's representation of the stable density. Near zero a Taylor approximation 1 - x**2 / 6 is used.

Parameters:

Name	Type	Description	Default
x	float or numpy.ndarray of float, any shape	Arguments, in radians.	required

Returns:

Type	Description
numpy.ndarray of float64, same shape as ``x``	sin(x) / x, elementwise (0-d for a scalar input).

Examples:

>>> import numpy as np
>>> float(sinc(0.0)), bool(np.isclose(sinc(np.pi / 2), 2 / np.pi))
(1.0, True)

rstable_positive

rstable_positive(size: int, alpha: float, rng: Generator) -> NDArray[np.float64]

Draw random positive numbers from a stable law; used to simulate the Gumbel copula.

Precisely, draws from the positive stable law :math:S(\alpha, 1, \gamma, 0; 1) with :math:\gamma = \cos(\pi\alpha/2)^{1/\alpha}. A "frailty" is a shared random factor: every coordinate of a Gumbel sample is generated from the same draw of it, which is what makes them dependent. You only need this if you are writing your own Archimedean sampler.

This is the frailty distribution of the Gumbel copula with :math:\alpha = 1/\theta. Its Laplace transform is :math:\mathbb{E}[e^{-tV}] = e^{-t^{\alpha}}, which is exactly the Gumbel generator.

Parameters:

Name	Type	Description	Default
size	int	Number of draws.	required
alpha	float	Stability index in (0, 1]. alpha = 1 degenerates to the constant 1 (the independence case, theta = 1).	required
rng	Generator	Source of randomness, e.g. numpy.random.default_rng(0).	required

Returns:

Type	Description
numpy.ndarray of float64, shape (size,)	Strictly positive draws.

Raises:

Type	Description
ValueError	If alpha is not in (0, 1].

Notes

Implements the Chambers-Mallows-Stuck transformation directly rather than going through scipy.stats.levy_stable: the generic scipy path is far slower and its S0/S1 parameterisation mapping is easy to get subtly wrong for the totally-skewed case needed here.

With :math:W \sim \mathrm{Exp}(1) and :math:U \sim \mathrm{Unif}(0, \pi),

.. math::

V = \frac{\sin(\alpha U)}{\sin(U)^{1/\alpha}}
    \left(\frac{\sin((1-\alpha)U)}{W}\right)^{(1-\alpha)/\alpha}

Examples:

The defining Laplace transform, checked by Monte Carlo:

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> v = rstable_positive(200_000, 0.5, rng)
>>> t = 1.3
>>> bool(abs(np.mean(np.exp(-t * v)) - np.exp(-t ** 0.5)) < 5e-3)
True

rlog_series

rlog_series(size: int, p: float, rng: Generator, log1mp: float | None = None) -> NDArray[np.float64]

Draw random positive integers from the logarithmic series law; used to simulate Frank.

The logarithmic series distribution with parameter p puts probability :math:P(V = k) = -p^k / (k \log(1 - p)) for :math:k = 1, 2, \dots. This is the frailty of the Frank copula with :math:p = 1 - e^{-\theta}.

Implements Kemp's (1981) "LK" algorithm: a chop-down search from the mode for small p, and an inversion-based shortcut otherwise.

Parameters:

Name	Type	Description	Default
size	int	Number of draws.	required
p	float	Series parameter in (0, 1].	required
rng	Generator	Source of randomness.	required
log1mp	float or None	A separately computed :math:\log(1 - p). Supply it whenever p was formed as 1 - e^{-\theta}: past theta around 37 that expression rounds to exactly 1, so log1p(-p) becomes -inf and the sampler fails outright -- even though the quantity it actually needs is just -theta, known exactly. The Frank copula reaches that at tau = 0.92, well inside the range people fit.	None

Returns:

Type	Description
numpy.ndarray of float64, shape (size,)	Integer-valued draws, each >= 1, stored as floats.

Raises:

Type	Description
ValueError	If p is not in (0, 1], or if log(1 - p) is not finite and negative (for example p has rounded to 1 and log1mp was not given).

Examples:

Mean of the logarithmic distribution is -p / ((1-p) log(1-p)):

>>> import numpy as np
>>> rng = np.random.default_rng(1)
>>> p = 0.6
>>> v = rlog_series(200_000, p, rng)
>>> expected = -p / ((1 - p) * np.log1p(-p))
>>> bool(abs(v.mean() - expected) / expected < 0.02)
True

Passing log1mp keeps it working where p itself has saturated:

>>> theta = 60.0
>>> p = -np.expm1(-theta)
>>> bool(p == 1.0)
True
>>> v = rlog_series(1000, p, rng, log1mp=-theta)
>>> bool(np.all(v >= 1))
True

rsibuya

rsibuya(size: int, alpha: float, rng: Generator) -> NDArray[np.float64]

Draw random positive integers from the Sibuya law; used to simulate the Joe copula.

The Sibuya distribution with parameter alpha in (0, 1] has :math:P(V = k) = \binom{k - 1 - \alpha}{k - 1}\frac{\alpha}{k}, with generating function :math:1 - (1 - t)^{\alpha}. This is the frailty of the Joe copula with :math:\alpha = 1/\theta.

Sampled by inversion (Hofert 2008): with :math:U \sim \mathrm{Unif}(0,1), return 1 if :math:U \le \alpha, else find the largest k whose exact survival probability :math:P(V \ge k) is at least :math:1 - U, starting from the tail asymptote and correcting with a short monotone search.

Parameters:

Name	Type	Description	Default
size	int	Number of draws.	required
alpha	float	Parameter in (0, 1]. alpha = 1 gives the constant 1.	required
rng	Generator	Source of randomness.	required

Returns:

Type	Description
numpy.ndarray of float64, shape (size,)	Integer-valued draws, each >= 1, stored as floats.

Raises:

Type	Description
ValueError	If alpha is not in (0, 1].

Examples:

Sibuya is heavy-tailed with infinite mean for alpha < 1; the sampler must therefore produce occasional very large values.

>>> import numpy as np
>>> rng = np.random.default_rng(2)
>>> v = rsibuya(100_000, 0.5, rng)
>>> bool((v >= 1).all() and v.max() > 100)
True

alpha = 1 degenerates to the point mass at 1:

>>> bool((rsibuya(1000, 1.0, rng) == 1).all())
True

The survival function is matched, not just the shape:

>>> from scipy.special import gammaln
>>> v = rsibuya(400_000, 0.7, np.random.default_rng(3))
>>> k = 5
>>> exact = np.exp(gammaln(k - 0.7) - gammaln(1 - 0.7) - gammaln(k))
>>> bool(abs((v >= k).mean() - exact) < 5e-3)
True

retstable

retstable(size: int, alpha: float, v0: NDArray[float64] | float, h: float, rng: Generator) -> NDArray[np.float64]

Draw random positive numbers from a tilted stable law; used to simulate nested Clayton.

This is the exponentially tilted positive stable distribution. It is defined by its Laplace transform,

.. math:: \mathbb{E}\bigl[e^{-tV}\bigr] = \exp\bigl(-V_0\,\bigl[(t + h)^{\alpha} - h^{\alpha}\bigr]\bigr),

which is the inner frailty of a nested Clayton copula and the reason nested Archimedean copulas have never been available in Python: this distribution has no closed-form density, no quantile function, and no implementation in NumPy, SciPy or any copula package.

Parameters:

Name	Type	Description	Default
size	int	Number of draws. Ignored when v0 is an array, whose length wins.	required
alpha	float	Stability index in (0, 1]. At alpha = 1 the law is degenerate at v0.	required
v0	float or numpy.ndarray of float, shape (m,)	The multiplier in the exponent, >= 0. May vary per draw, which is what the nested sampler needs -- each observation carries its own outer frailty.	required
h	float	Tilt, >= 0. h = 0 gives the untilted stable.	required
rng	Generator	Source of randomness.	required

Returns:

Type	Description
numpy.ndarray of float64, shape (size,) or (m,)	Non-negative draws; (m,) when v0 is an array (it is flattened first).

Raises:

Type	Description
ValueError	If alpha is not in (0, 1], h < 0, or any v0 < 0.

Notes

Method. Rejection from the untilted stable accepts with probability exactly :math:e^{-V_0 h^\alpha}, which is fine for a small exponent and hopeless otherwise -- worse, for V0 ~ Gamma(k, 1) with k >= 1 the expected number of attempts is infinite, so a naive implementation does not merely run slowly, it occasionally does not finish.

Infinite divisibility fixes it. The law is the sum of k independent draws with V0/k in place of V0 -- the Laplace transforms multiply to give back the original -- and each of those accepts with probability :math:e^{-V_0 h^\alpha / k}. Taking :math:k = \lceil V_0 h^\alpha\rceil puts every piece's acceptance at about :math:e^{-1}, so the work becomes linear in :math:V_0 h^\alpha rather than exponential. The result is exact: no approximation enters at any point.

Examples:

The defining Laplace transform, checked by Monte Carlo:

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> alpha, v0, h, t = 0.5, 2.0, 1.0, 0.8
>>> v = retstable(200_000, alpha, v0, h, rng)
>>> want = np.exp(-v0 * ((t + h) ** alpha - h ** alpha))
>>> bool(abs(np.mean(np.exp(-t * v)) - want) < 3e-3)
True

The mean is alpha * v0 * h**(alpha - 1), by differentiating it:

>>> bool(abs(v.mean() - alpha * v0 * h ** (alpha - 1)) < 0.02)
True

A zero tilt gives the ordinary positive stable back:

>>> untilted = retstable(5, 0.5, 1.0, 0.0, np.random.default_rng(1))
>>> bool(np.all(untilted > 0))
True

rcopula.special.mvtnorm

Multivariate normal and Student-t distribution functions.

scipy.stats.multivariate_normal.cdf exists but is unsuitable here for two reasons. It is not deterministic — two identical calls on a 10-dimensional problem differ by ~5e-6 — and it accepts no random_state, so the non-determinism cannot even be pinned down. A copula library needs a CDF that returns the same number every time, both for reproducible research and so that golden fixtures mean anything.

This module therefore provides its own, with the accuracy split by dimension:

• d = 2: exact, via Owen's T function. No quadrature at all.
• d = 3: a one-dimensional Gauss-Legendre integral of the exact bivariate CDF, accurate to about 1e-14 (:func:tvn_cdf).
• d >= 4: the Genz separation-of-variables transformation evaluated on a scrambled Sobol sequence with a fixed seed. Randomised QMC converges far faster than plain Monte Carlo on this integrand and, with the seed pinned, is reproducible to the last bit.

References

Genz, A. (1992). Numerical computation of multivariate normal probabilities. Journal of Computational and Graphical Statistics 1(2), 141-149. Genz, A. and Bretz, F. (2009). Computation of Multivariate Normal and t Probabilities. Springer Lecture Notes in Statistics 195. Chapter 4.2 gives the Student-t radial mixture used in :func:mvt_cdf. Owen, D. B. (1956). Tables for computing bivariate normal probabilities. Annals of Mathematical Statistics 27(4), 1075-1090.

bvn_cdf

bvn_cdf(h: ArrayLike, k: ArrayLike, rho: float) -> NDArray[np.float64]

Return the probability that two correlated standard normals both fall below given limits.

This is the standard bivariate normal CDF :math:\Phi_2(h, k; \rho) = :math:P(X \le h,\ Y \le k) for standard normal :math:X, Y with correlation :math:\rho. The Gaussian copula's CDF in two dimensions is exactly this function evaluated at normal quantiles, so it has to be exact and repeatable; it is computed in closed form, with no random numbers.

Uses the Owen (1956) decomposition

.. math::

\Phi_2(h,k;\rho) = \tfrac{1}{2}\bigl(\Phi(h) + \Phi(k)\bigr)
                   - T(h, a_h) - T(k, a_k) - \beta

with :math:a_h = (k/h - \rho)/\sqrt{1-\rho^2}, and :math:\beta = 1/2 when :math:hk < 0 (or :math:hk = 0 with :math:h + k < 0), else 0.

Parameters:

Name	Type	Description	Default
h	array_like of float	Upper limits for :math:X and :math:Y, broadcast against each other. +-inf is allowed.	required
k	array_like of float	Upper limits for :math:X and :math:Y, broadcast against each other. +-inf is allowed.	required
rho	float	Correlation in [-1, 1]. The endpoints -1 and 1 use their closed forms.	required

Returns:

Type	Description
numpy.ndarray of float64, broadcast shape of ``h`` and ``k``	:math:P(X \le h,\ Y \le k).

Raises:

Type	Description
ValueError	If rho lies outside [-1, 1].

Examples:

Independence factorises:

>>> import numpy as np
>>> from rcopula.special.mvtnorm import bvn_cdf
>>> float(np.round(bvn_cdf(0.5, 0.3, 0.0), 12))
0.427262552836
>>> float(np.round(bvn_cdf(0.5, 0.3, 0.0) - ndtr_(0.5) * ndtr_(0.3), 15))
0.0

At the origin there is a closed form, 1/4 + arcsin(rho)/(2*pi):

>>> rho = 0.6
>>> exact = 0.25 + np.arcsin(rho) / (2 * np.pi)
>>> bool(np.isclose(bvn_cdf(0.0, 0.0, rho), exact, rtol=1e-14))
True

ndtr_

ndtr_(x: float) -> float

Return the standard normal CDF at x as a Python float.

A doctest helper, exposed so the :func:bvn_cdf example reads cleanly; it is not part of __all__.

Parameters:

Name	Type	Description	Default
x	float	Point at which to evaluate.	required

Returns:

Type	Description
float	:math:\Phi(x) = P(Z \le x) for a standard normal :math:Z.

tvn_cdf

tvn_cdf(upper: ArrayLike, corr: ArrayLike) -> NDArray[np.float64]

Return the probability that three correlated standard normals all fall below given limits.

This is the trivariate normal CDF :math:\Phi_3(x_1, x_2, x_3; R), the three-dimensional Gaussian-copula CDF on the normal scale. It is computed to near machine precision and returns the same number on every call. :func:mvn_cdf dispatches here when d = 3.

Conditioning on the third coordinate reduces the problem to a one dimensional integral of the exact bivariate CDF:

.. math::

\Phi_3(x_1,x_2,x_3; R) = \int_{-\infty}^{x_3} \phi(t)\,
    \Phi_2\!\left(\frac{x_1 - \rho_{13}t}{\sigma_1},
                    \frac{x_2 - \rho_{23}t}{\sigma_2};
                    \rho_{12\cdot 3}\right) \mathrm{d}t,

with :math:\sigma_i = \sqrt{1-\rho_{i3}^2} and the partial correlation :math:\rho_{12\cdot3} = (\rho_{12} - \rho_{13}\rho_{23})/(\sigma_1\sigma_2).

The integrand is analytic, so a fixed Gauss-Legendre rule reaches ~1e-14 -- seven orders of magnitude better than the quasi-Monte-Carlo path this replaces at d = 3, and comparable to R's TVPACK.

Parameters:

Name	Type	Description	Default
upper	array_like of float, shape (n, 3) or (3,)	Upper integration limits, one row per point. Must be finite; use :func:mvn_cdf for infinite limits.	required
corr	array_like of float, shape (3, 3)	Correlation matrix (unit diagonal). Only the three off-diagonal entries corr[0, 1], corr[0, 2] and corr[1, 2] are read.	required

Returns:

Type	Description
numpy.ndarray of float64, shape (n,)	:math:P(X_1 \le x_1, X_2 \le x_2, X_3 \le x_3) for each row.

Notes

If a coordinate is perfectly correlated with the third (:math:|\rho_{13}| = 1 or :math:|\rho_{23}| = 1), the conditioning integral is undefined and the Genz quasi-Monte-Carlo integrator is used instead.

Examples:

>>> import numpy as np
>>> from rcopula.special.mvtnorm import tvn_cdf
>>> corr = np.array([[1.0, 0.5, 0.3], [0.5, 1.0, 0.4], [0.3, 0.4, 1.0]])
>>> float(np.round(tvn_cdf([[0.0, 0.0, 0.0]], corr)[0], 12))
0.22366080778

Independence factorises exactly:

>>> vals = tvn_cdf([[0.4, -0.2, 1.1]], np.eye(3))[0]
>>> from scipy.special import ndtr
>>> bool(np.isclose(vals, ndtr(0.4) * ndtr(-0.2) * ndtr(1.1), rtol=1e-13))
True

mvn_cdf

mvn_cdf(upper: ArrayLike, corr: ArrayLike, *, n_points: int | None = None) -> NDArray[np.float64]

Return the probability that d correlated standard normals all fall below given limits.

This is the multivariate normal CDF with unit variances, :math:P(X_1 \le x_1, \dots, X_d \le x_d) for :math:X \sim N(0, R). The Gaussian copula's CDF is this function evaluated at normal quantiles. Unlike scipy.stats.multivariate_normal.cdf it returns the same number on every call, which reproducible research and golden test fixtures need.

How it is computed depends on the dimension: d = 1 is the normal CDF, d = 2 is exact (:func:bvn_cdf), d = 3 is a near-exact quadrature (:func:tvn_cdf), and d >= 4 uses Genz's method on a fixed-seed scrambled Sobol sequence.

Parameters:

Name	Type	Description	Default
upper	array_like of float, shape (n, d) or (d,)	Upper integration limits, one row per point. +inf removes that coordinate's constraint; -inf (or nan) gives probability 0.	required
corr	array_like of float, shape (d, d)	Correlation matrix (unit diagonal, positive definite for d >= 4).	required
n_points	int or None	Number of quasi-Monte-Carlo points, used only when d >= 4. None means 2**16.	None

Returns:

Type	Description
numpy.ndarray of float64, shape (n,)	P(X_1 <= upper_1, ..., X_d <= upper_d) for each row.

Raises:

Type	Description
ValueError	If upper does not have d columns.
LinAlgError	If d >= 4 and corr is not positive definite.

Notes

Deterministic: repeated calls with the same arguments return bit-identical results. d = 2 is exact and d = 3 is accurate to about 1e-14; d >= 4 carries a QMC error of roughly 1e-8 at the default point count.

Examples:

>>> import numpy as np
>>> from rcopula.special.mvtnorm import mvn_cdf
>>> corr = np.array([[1.0, 0.5], [0.5, 1.0]])
>>> float(np.round(mvn_cdf([[0.0, 0.0]], corr)[0], 12))
0.333333333333

Determinism, which scipy's version does not offer:

>>> c = np.eye(5) * 0.5 + 0.5
>>> a = mvn_cdf([[0.3] * 5], c)
>>> b = mvn_cdf([[0.3] * 5], c)
>>> bool(a[0] == b[0])
True

mvt_cdf

mvt_cdf(upper: ArrayLike, corr: ArrayLike, df: float, *, n_points: int | None = None) -> NDArray[np.float64]

Return the probability that a d-dimensional Student-t vector falls below given limits.

This is the multivariate Student-t CDF with unit scale, :math:P(X_1 \le x_1, \dots, X_d \le x_d) for :math:X \sim t_\nu(0, R). The Student-t copula's CDF is this function evaluated at t quantiles. Like :func:mvn_cdf it is deterministic.

Uses the radial mixture representation: if :math:X \sim t_\nu(0, R) then :math:X = Z / \sqrt{W/\nu} with :math:Z \sim N(0, R) and :math:W \sim \chi^2_\nu, so the t probability is the normal probability at rescaled limits, averaged over :math:W (Genz & Bretz 2009, §4.2).

Unlike R's mvtnorm::pmvt, non-integer degrees of freedom are supported — the radial variable is continuous, so there is no reason to restrict it, and t copulas are routinely fitted with fractional df.

Parameters:

Name	Type	Description	Default
upper	array_like of float, shape (n, d) or (d,)	Upper integration limits, one row per point. +inf is allowed when d <= 3.	required
corr	array_like of float, shape (d, d)	Correlation matrix (unit diagonal).	required
df	float	Degrees of freedom, > 0. May be non-integer.	required
n_points	int or None	Number of quasi-Monte-Carlo points, used only when d >= 4 (d <= 3 uses a fixed Gauss-Legendre rule). None means 2**16.	None

Returns:

Type	Description
numpy.ndarray of float64, shape (n,)	P(X_1 <= upper_1, ..., X_d <= upper_d) for each row.

Raises:

Type	Description
ValueError	If df <= 0, or if upper does not have d columns.

Notes

Accuracy is about 1e-14 for d <= 3 and df >= 2, about 1e-7 for df < 2 with extreme limits, and limited by the QMC error (roughly 1e-7) for d >= 4. Very small df degrades gracefully rather than failing: about 4e-4 at df = 0.02 and 4e-3 at df = 0.01.

Examples:

Large df converges to the normal case:

>>> import numpy as np
>>> from rcopula.special.mvtnorm import mvn_cdf, mvt_cdf
>>> corr = np.array([[1.0, 0.4], [0.4, 1.0]])
>>> t_big = mvt_cdf([[0.5, -0.2]], corr, df=1e6)[0]
>>> normal = mvn_cdf([[0.5, -0.2]], corr)[0]
>>> bool(abs(t_big - normal) < 1e-4)
True

Non-integer degrees of freedom, which mvtnorm::pmvt rejects:

>>> float(mvt_cdf([[0.0, 0.0]], corr, df=3.5)[0]) > 0.3
True

rcopula.special.debye

Debye functions :math:D_n(x).

.. math::

D_n(x) = \frac{n}{x^n} \int_0^x \frac{t^n}{e^t - 1}\,\mathrm{d}t

scipy does not provide these, yet they are exactly what the Frank copula's dependence measures are made of:

.. math::

\tau(\theta)  &= 1 - \frac{4}{\theta}\bigl(1 - D_1(\theta)\bigr) \\
\rho(\theta)  &= 1 - \frac{12}{\theta}\bigl(D_1(\theta) - D_2(\theta)\bigr)

Implemented from the standard series representations:

Abramowitz, M. and Stegun, I. A. (1964). *Handbook of Mathematical Functions*,
§27.1. NIST DLMF §25.12 gives the same expansions.

Two regimes are used, switching at :math:|x| = 2:

• Small :math:|x| \le 2 — the Bernoulli series, convergent for :math:|x| < 2\pi:

.. math::

  D_n(x) = 1 - \frac{n x}{2(n+1)}
         + n \sum_{k \ge 1} \frac{B_{2k}\, x^{2k}}{(2k+n)\,(2k)!}

Direct evaluation of the integral form here would cancel catastrophically.

• Large :math:|x| > 2 — expand :math:(e^t-1)^{-1} = \sum_{k\ge1} e^{-kt} and integrate termwise:

.. math::

  D_n(x) = \frac{n\,n!}{x^n} \sum_{k \ge 1} \frac{1}{k^{n+1}}
           \left[1 - e^{-kx} \sum_{j=0}^{n} \frac{(kx)^j}{j!}\right]

The bracket tends to :math:1, so the terms decay only like :math:k^{-(n+1)} and truncating the sum directly would leave a zeta-tail error of order :math:K^{-n} — about 1% at :math:K = 60, :math:n = 1. Instead the constant part is summed in closed form, :math:\sum_{k\ge1} k^{-(n+1)} = \zeta(n+1), leaving

.. math::

  D_n(x) = \frac{n\,n!}{x^n} \left[ \zeta(n+1)
         - \sum_{k \ge 1} \frac{e^{-kx}}{k^{n+1}}
           \sum_{j=0}^{n} \frac{(kx)^j}{j!} \right]

whose remaining sum decays geometrically like :math:e^{-kx}. This also makes the correct asymptote :math:D_n(x) \to n\,n!\,\zeta(n+1)/x^n exact.

Negative arguments use the exact reflection :math:D_n(-x) = D_n(x) + \frac{n x}{n + 1}, which follows from the series (only even powers of :math:x appear beyond the linear term).

even_bernoulli_over_factorial

even_bernoulli_over_factorial(count: int) -> NDArray[np.float64]

B_{2k} / (2k)! for k = 1..count, from exact rational arithmetic.

scipy.special.bernoulli is not accurate enough for a series meant to reach machine precision: its B_4 is off by 1.7e-12 relative, which capped the Debye functions (and so Frank's tau and rho) at ~4e-14. The numbers are computed here as exact fractions by the standard recurrence :math:\sum_{j=0}^{m} \binom{m+1}{j} B_j = 0 and divided by the factorial before rounding, so each value is correctly rounded.

debye_n

debye_n(x: ArrayLike, n: int) -> NDArray[np.float64]

Evaluate the Debye function of order n, a smooth integral the Frank copula needs.

The Debye function :math:D_n(x) is a normalised integral (defined in the module docstring) that starts at 1 when x = 0 and decays towards 0 as x grows. The package needs it because the Frank copula's Kendall's tau and Spearman's rho are written in terms of :math:D_1 and :math:D_2, and scipy does not provide it. Use :func:debye1 or :func:debye2 for the two orders the copulas actually use.

Parameters:

Name	Type	Description	Default
x	array_like of float, any shape	Real arguments; a Python float works too. x = 0 returns 1 (the limiting value); negative arguments are handled exactly by reflection; nan gives nan.	required
n	int	Order, n >= 1.	required

Returns:

Type	Description
numpy.ndarray of float64, same shape as ``x``	:math:D_n(x), elementwise. A scalar input gives a 0-d result (a numpy.float64).

Raises:

Type	Description
ValueError	If n < 1.

Notes

:math:D_n is decreasing on :math:(0, \infty) with :math:D_n(0) = 1 and :math:D_n(x) \to n\,n!\,\zeta(n+1) / x^n as :math:x \to \infty.

Examples:

>>> import numpy as np
>>> from rcopula.special import debye_n
>>> float(debye_n(0.0, 1))
1.0
>>> float(np.round(debye_n(1.0, 1), 10))
0.7775046341

The reflection identity holds exactly:

>>> x = 3.7
>>> lhs = float(debye_n(-x, 2))
>>> rhs = float(debye_n(x, 2)) + 2 * x / 3
>>> bool(np.isclose(lhs, rhs, rtol=1e-14))
True

debye1

debye1(x: ArrayLike) -> NDArray[np.float64]

Evaluate the first-order Debye function :math:D_1(x), used for Frank's Kendall tau.

Shorthand for debye_n(x, 1). It appears in the Frank copula's Kendall tau, :math:\tau(\theta) = 1 - 4\bigl(1 - D_1(\theta)\bigr)/\theta, which is how the package converts between Frank's theta and tau.

Parameters:

Name	Type	Description	Default
x	array_like of float, any shape	Real arguments. x = 0 returns 1.	required

Returns:

Type	Description
numpy.ndarray of float64, same shape as ``x``	:math:D_1(x), elementwise (0-d for a scalar input).

Examples:

Frank's tau at theta = 5:

>>> import numpy as np
>>> from rcopula.special import debye1
>>> theta = 5.0
>>> tau = 1 - 4 * (1 - float(debye1(theta))) / theta
>>> float(np.round(tau, 10))
0.4567009582

debye2

debye2(x: ArrayLike) -> NDArray[np.float64]

Evaluate the second-order Debye function :math:D_2(x), used for Frank's Spearman rho.

Shorthand for debye_n(x, 2). It appears in the Frank copula's Spearman rho, :math:\rho(\theta) = 1 - 12\bigl(D_1(\theta) - D_2(\theta)\bigr)/\theta.

Parameters:

Name	Type	Description	Default
x	array_like of float, any shape	Real arguments. x = 0 returns 1.	required

Returns:

Type	Description
numpy.ndarray of float64, same shape as ``x``	:math:D_2(x), elementwise (0-d for a scalar input).

Examples:

Frank's rho at theta = 5:

>>> import numpy as np
>>> from rcopula.special import debye1, debye2
>>> theta = 5.0
>>> rho = 1 - 12 * (float(debye1(theta)) - float(debye2(theta))) / theta
>>> float(np.round(rho, 10))
0.6434871081

rcopula.special.logexp

Accurate evaluation of log(1 - exp(-a)) and log(1 + exp(x)).

Both expressions are numerically treacherous: each has a regime where the naive formula cancels catastrophically and a regime where it overflows. The piecewise-optimal cutoffs used here are derived in

Mächler, M. (2012). *Accurately Computing log(1 - exp(-|a|))*.
Assessed by the Rmpfr package vignette.
https://cran.r-project.org/package=Rmpfr/vignettes/log1mexp-note.pdf

For log1mexp the optimal switch is at a = log 2: below it log(-expm1(-a)) is accurate, above it log1p(-exp(-a)) is. For log1pexp the note gives the four-regime cutoffs -37 / 18 / 33.3, chosen so that the discarded term falls below the double-precision epsilon.

Also provided is a sign-aware log-sum-exp, needed because the d-dimensional Archimedean densities are alternating sums whose terms span many orders of magnitude.

log1mexp

log1mexp(a: ArrayLike) -> NDArray[np.float64]

Compute log(1 - exp(-a)) without the rounding errors of the naive formula.

Typing np.log(1 - np.exp(-a)) loses every digit when a is tiny (1 - exp(-a) rounds to 0) and is inaccurate in other ranges too. Copula log-densities and log-CDFs are full of this expression, so the package routes it through here. Use it whenever you need the log of a probability written as 1 - exp(-a).

Parameters:

Name	Type	Description	Default
a	array_like of float, any shape	Non-negative values; a Python float works too. a = 0 yields -inf (since log 0 = -inf); negative values yield nan, because 1 - exp(-a) < 0 there. No error is raised for either.	required

Returns:

Type	Description
numpy.ndarray of float64, same shape as ``a``	log(1 - exp(-a)), elementwise. A scalar input gives a numpy.float64.

Notes

Uses log(-expm1(-a)) for a <= log 2 and log1p(-exp(-a)) above it, the switch point that minimises the worst-case relative error (Mächler 2012).

Examples:

>>> import numpy as np
>>> from rcopula.special import log1mexp
>>> float(log1mexp(1e-10))          # naive log(1 - exp(-a)) loses all precision here
-23.025850929990458
>>> bool(np.isinf(log1mexp(0.0)))
True

log1pexp

log1pexp(x: ArrayLike) -> NDArray[np.float64]

Compute log(1 + exp(x)) (the softplus function) without overflow or rounding loss.

The naive formula overflows for large x (exp(1000) is inf) and loses accuracy for very negative x. This version is accurate for every real input, which matters inside log-likelihoods that are evaluated far into the tails.

Parameters:

Name	Type	Description	Default
x	array_like of float, any shape	Any real values; a Python float works too. nan gives nan.	required

Returns:

Type	Description
numpy.ndarray of float64, same shape as ``x``	log(1 + exp(x)), elementwise. A scalar input gives a numpy.float64.

Notes

Four regimes, with cutoffs from Mächler (2012):

================== ============================== Range Formula ================== ============================== x <= -37 exp(x) -37 < x <= 18 log1p(exp(x)) 18 < x <= 33.3 x + exp(-x) x > 33.3 x ================== ==============================

Examples:

>>> import numpy as np
>>> from rcopula.special import log1pexp
>>> float(log1pexp(0.0))
0.6931471805599453
>>> float(log1pexp(1000.0))         # naive exp(1000) would overflow
1000.0

signed_logsumexp

signed_logsumexp(log_abs: ArrayLike, signs: ArrayLike, axis: int | None = None) -> tuple[NDArray[np.float64], NDArray[np.float64]]

Add up positive and negative numbers that are stored as logarithms, without overflow.

Computes log|S| and sign(S) where S = sum(signs * exp(log_abs)). It is scipy.special.logsumexp extended to terms that may be negative, and it is what you need when the individual terms are far too large or too small to hold as ordinary floats but their (possibly cancelling) sum is wanted.

This is the workhorse for d-dimensional Archimedean densities, whose generator derivatives are alternating sums with terms spanning hundreds of orders of magnitude: forming them directly overflows, and dropping the sign information gives the wrong answer.

Parameters:

Name	Type	Description	Default
log_abs	array_like of float, any shape	Logarithms of the absolute values of the terms. -inf stands for a zero term.	required
signs	array_like of float, broadcastable to ``log_abs.shape``	Signs (+1 or -1, or 0 to ignore a term) of the terms.	required
axis	int or None	Axis to sum over. None sums the flattened input.	None

Returns:

Name	Type	Description
log_abs_sum	numpy.ndarray of float64	log|S|, with axis removed from the shape (a numpy.float64 when axis=None). Equals -inf where the terms cancel exactly.
sign_sum	numpy.ndarray of float64	sign(S), one of -1.0, 0.0, 1.0; same shape as log_abs_sum.

Examples:

Cancellation that would be lost in ordinary arithmetic:

>>> import numpy as np
>>> from rcopula.special import signed_logsumexp
>>> lv, sg = signed_logsumexp([0.0, 0.0], [1.0, -1.0])
>>> bool(np.isneginf(lv)), float(sg)
(True, 0.0)

A plain sum, exp(1) - exp(0) = 1.718...:

>>> lv, sg = signed_logsumexp([1.0, 0.0], [1.0, -1.0])
>>> float(sg) * float(np.exp(lv))
1.7182818284590453

rcopula.special.combinatorics

Stirling and Eulerian numbers.

These are the coefficient machinery behind the d-dimensional Archimedean densities. The Gumbel generator's d-th derivative, in particular, is a polynomial whose coefficients are built from Stirling numbers of both kinds:

.. math::

a_{d,k}(\theta) = (-1)^{d-k} \sum_{j=k}^{d}
    \theta^{-j}\, s(d, j)\, S(j, k)

as given in

Hofert, M., Mächler, M. and McNeil, A. J. (2012). Likelihood inference for
Archimedean copulas in high dimensions under known margins.
*Journal of Multivariate Analysis* 110, 133-150.

Exact integer recurrences are used and results are cached, because these are evaluated inside likelihood loops and the values grow large enough that recomputing in floating point would both cost time and lose digits.

stirling1

stirling1(n: int, k: int) -> float

Return the signed Stirling number of the first kind, s(n, k).

A whole number from combinatorics. In plain terms, |s(n, k)| counts the ways to arrange n labelled objects into k cycles; the sign is (-1)**(n - k). Together with :func:stirling2 it gives the polynomial coefficients of the Gumbel generator's high-order derivatives, which the d-dimensional Gumbel density needs. You will rarely call it directly.

Parameters:

Name	Type	Description	Default
n	int	Number of objects, n >= 0.	required
k	int	Number of cycles, k >= 0.	required

Returns:

Type	Description
float	s(n, k), computed exactly in integer arithmetic and returned as a float. 0.0 when k > n.

Raises:

Type	Description
ValueError	If n or k is negative.

Examples:

>>> from rcopula.special.combinatorics import stirling1
>>> [stirling1(4, k) for k in range(5)]
[0.0, -6.0, 11.0, -6.0, 1.0]

stirling1_all

stirling1_all(n: int) -> NDArray[np.float64]

Return a whole row of signed Stirling numbers of the first kind at once.

Equivalent to [stirling1(n, k) for k in range(1, n + 1)] but taken from the same cached row in one go.

Parameters:

Name	Type	Description	Default
n	int	Number of objects, n >= 0.	required

Returns:

Type	Description
numpy.ndarray of float64, shape (n,)	s(n, k) for k = 1, ..., n. Empty when n = 0.

Examples:

>>> from rcopula.special.combinatorics import stirling1_all
>>> stirling1_all(4).tolist()
[-6.0, 11.0, -6.0, 1.0]

stirling2

stirling2(n: int, k: int) -> float

Count the ways to split n labelled items into k non-empty groups.

This is the Stirling number of the second kind S(n, k): the number of ways to partition n labelled objects into k non-empty unlabelled subsets. Together with :func:stirling1 it gives the polynomial coefficients of the Gumbel generator's high-order derivatives, which the d-dimensional Gumbel density needs.

Parameters:

Name	Type	Description	Default
n	int	Number of objects, n >= 0.	required
k	int	Number of groups, k >= 0.	required

Returns:

Type	Description
float	S(n, k), computed exactly in integer arithmetic and returned as a float. 0.0 when k > n.

Raises:

Type	Description
ValueError	If n or k is negative.

Examples:

>>> from rcopula.special.combinatorics import stirling2
>>> [stirling2(4, k) for k in range(5)]
[0.0, 1.0, 7.0, 6.0, 1.0]

stirling2_all

stirling2_all(n: int) -> NDArray[np.float64]

Return a whole row of Stirling numbers of the second kind at once.

Equivalent to [stirling2(n, k) for k in range(1, n + 1)] but taken from the same cached row in one go.

Parameters:

Name	Type	Description	Default
n	int	Number of objects, n >= 0.	required

Returns:

Type	Description
numpy.ndarray of float64, shape (n,)	S(n, k) for k = 1, ..., n. Empty when n = 0.

Examples:

>>> from rcopula.special.combinatorics import stirling2_all
>>> stirling2_all(4).tolist()
[1.0, 7.0, 6.0, 1.0]

eulerian

eulerian(n: int, k: int) -> float

Count the orderings of n items that go up exactly k times.

This is the Eulerian number A(n, k): the number of permutations of 1..n with exactly k ascents (positions where the next entry is larger). The package uses them as the coefficients of the polylogarithm of negative integer order, which the d-dimensional densities of several Archimedean families are built from.

Parameters:

Name	Type	Description	Default
n	int	Permutation length, n >= 0.	required
k	int	Number of ascents, k >= 0.	required

Returns:

Type	Description
float	A(n, k), computed exactly in integer arithmetic and returned as a float. 0.0 when k > n; note that A(n, n) = 0 for n >= 1.

Raises:

Type	Description
ValueError	If n or k is negative.

Examples:

>>> from rcopula.special.combinatorics import eulerian
>>> [eulerian(4, k) for k in range(4)]
[1.0, 11.0, 11.0, 1.0]

eulerian_all

eulerian_all(n: int) -> NDArray[np.float64]

Return a whole row of Eulerian numbers at once.

Equivalent to [eulerian(n, k) for k in range(n)] (the entries that can be non-zero) but taken from the same cached row in one go.

Parameters:

Name	Type	Description	Default
n	int	Permutation length, n >= 0.	required

Returns:

Type	Description
numpy.ndarray of float64, shape (max(n, 1),)	A(n, k) for k = 0, ..., n - 1. For n = 0 this is [1.0].

Examples:

>>> from rcopula.special.combinatorics import eulerian_all
>>> eulerian_all(4).tolist()
[1.0, 11.0, 11.0, 1.0]