Low-level API¶

Calculates the “sandwich” matrix product (A @ X @ A.T) along the specified X axis.

Parameters:

Returns:

Construct the lower-diagonal Pascal matrix \(L_{k \times k\)}$, or its matrix inverse \(L^{-1}\).

Implemented using recursion, unlike the slow naive implementation from the equivalent scipy.linalg.pascal and scipy.linalg.invpascal functions using kind='lower'. By using the binomial recurrence relation, assuming \(0 < j < i\), \(\binom{i}{j} = \frac{i}{j} \binom{i-1}{j-1}\), the following recursive definition is obtained:

Examples:

>>> import numpy as np
>>> pascal(4, dtype=np.int_)
array([[1, 0, 0, 0],
       [1, 1, 0, 0],
       [1, 2, 1, 0],
       [1, 3, 3, 1]])
>>> pascal(4, dtype=np.int_, inv=True)
array([[ 1,  0,  0,  0],
       [-1,  1,  0,  0],
       [ 1, -2,  1,  0],
       [-1,  3, -3,  1]])
>>> np.rint(np.linalg.inv(pascal(4))).astype(np.int_)
array([[ 1,  0,  0,  0],
       [-1,  1,  0,  0],
       [ 1, -2,  1,  0],
       [-1,  3, -3,  1]])

Now, let’s compare with scipy:

>>> import scipy.linalg
>>> scipy.linalg.invpascal(4, kind='lower').astype(np.int_)
array([[ 1,  0,  0,  0],
       [-1,  1,  0,  0],
       [ 1, -2,  1,  0],
       [-1,  3, -3,  1]])

Inverse regulatized lower-diagonal Pascal matrix, \(\bar{L}_{ij} = L^{-1}_ij / i\).

Used to linearly combine order statistics order statistics into L-moments.

Shifted Legendre polynomial coefficient matrix \(\widetilde{P}\) of shape (k, k).

The \(j\)-th coefficient of the shifted Legendre polynomial of degree \(k\) is at \((k, j)\):

Useful for transforming probability-weighted moments into L-moments.

Parameters:

Returns:

Examples:

Calculate \(\widetilde{P}_{4 \times 4}\):

>>> from lmo.linalg import sh_legendre
>>> sh_legendre(4, dtype=int)
array([[  1,   0,   0,   0],
       [ -1,   2,   0,   0],
       [  1,  -6,   6,   0],
       [ -1,  12, -30,  20]])

Shifted Jacobi polynomial coefficient matrix \(\widetilde{P}^{(a,b)}\) of shape (k, k).

The \(j\)-th coefficient of the shifted Jacobi polynomial of degree \(k\) is at \((k, j)\):

The “shift” refers to the change of variables \(x \mapsto 2x - 1\) in the (unshifted) Jacobi (or hypergeometric) polynomials.

The (shifted) Jacobi polynomials \(\widetilde{P}^{(a,b)}\) generalize the (shifted) Legendre polynomials \(\widetilde{P}\): \(\widetilde{P}^{(0, 0)} = \widetilde{P}\)

Parameters:

Returns:

Examples:

Calculate \(\widetilde{P}^{(1, 1)}_{4 \times 4}\):

>>> from lmo.linalg import sh_jacobi
>>> sh_jacobi(4, 1, 1, dtype=int)
array([[  1,   0,   0,   0],
       [ -2,   4,   0,   0],
       [  3, -15,  15,   0],
       [ -4,  36, -84,  56]])

Let’s compare \(\widetilde{P}^(1, \pi)_3\) with the scipy Jacobi poly1d. This requires manual shifting \(x \mapsto f(x)\), with \(f(x) = 2x - 1\):

>>> import numpy as np
>>> import scipy.special as sc
>>> f_x = np.poly1d([2, -1])  # f(x) = 2*x + 1
>>> sc.jacobi(3, 1, np.pi)(f_x)
poly1d([ 125.80159497, -228.55053774,  128.54584648,  -21.79690371])
>>> sh_jacobi(4, 1, np.pi)[3]
array([ -21.79690371,  128.54584648, -228.55053774,  125.80159497])

Apart from the reversed coefficients of numpy.poly1d (an awkward design choice, but it’s fixed in the new numpy.polynomial module.)

A toeplitz-like transformation matrix construction, that prepends \(i\) zeroes to \(i\)-th row, so that the input shape is mapped from (n, k) to (n, k + n).

So all values \(i > j \vee i + j \ge k\) are zero in the succession matrix.

Parameters:

Returns:

Examples:

>>> from lmo.linalg import succession_matrix
>>> c = np.arange(1, 9).reshape(4, 2)
>>> c
array([[1, 2],
       [3, 4],
       [5, 6],
       [7, 8]])
>>> succession_matrix(c)
array([[1, 2, 0, 0, 0],
       [0, 3, 4, 0, 0],
       [0, 0, 5, 6, 0],
       [0, 0, 0, 7, 8]])

Linearization of the trimmed L-moment recurrence relations, following the (corrected) derivation by Hosking (2007) from the (shifted) Jacobi Polynomials.

This constructs a \(r \times r + t_1 + t_2\) matrix \(T^{(t_1, t_2)}\) that “trims” conventional L-moments. E.g. the first 3 \((1, 1)\) trimmed L-moments can be obtained from the first \(3+1+1=5\) (untrimmed) L-moments (assuming they exist) with trim_matrix(3, (1, 1)) @ l_moment(x, np.ogrid[:5] + 1).

The big “L” in “L-moment”, referring to it being a Linear combination of order statistics, has been prominently put in the name by Hosking (1990) for a good reason. It means that transforming order statistics to a bunch of L-moments, can be done using a single matrix multiplication (see lmo.linalg.sh_legendre). By exploiting liniarity, it can easily be chained with this trim matrix, to obtain a reusable order-statistics -> trimmed L-moments transformation (matrix).

Note that these linear transformations can be used in exactly the same way to e.g. calculate several population TL-moments of some random varianble, using nothing but its theoretical probablity-weighted moments (PWMs).

Parameters:

Returns:

Examples:

>>> from lmo.linalg import trim_matrix
>>> trim_matrix(3, (0, 1))
array([[ 1.        , -1.        ,  0.        ,  0.        ],
       [ 0.        ,  0.75      , -0.75      ,  0.        ],
       [ 0.        ,  0.        ,  0.66666667, -0.66666667]])
>>> trim_matrix(3, (1, 0))
array([[1.        , 1.        , 0.        , 0.        ],
       [0.        , 0.75      , 0.75      , 0.        ],
       [0.        , 0.        , 0.66666667, 0.66666667]])

Factorial power, or falling factorial.

It is defined as

Parameters:

Returns:

Incomplete (upper) gamma function.

It is defined as

for \( a \ge 0 \) and \( x \ge 0 \).

Parameters:

Returns:

Harmonic number \( H_n = \sum_{k=1}^{n} 1 / k \), extended for real and complex argument via analytic contunuation.

Examples:

>>> harmonic(0)
0.0
>>> harmonic(1)
1.0
>>> harmonic(2)
1.5
>>> harmonic(42)
4.32674
>>> harmonic(np.pi)
1.87274
>>> harmonic(-1 / 12)
-0.146106

Parameters:

Returns:

Evaluate the (weighted) \( L^2 \)-norm of a shifted Jacobi polynomial.

Specifically,

with

the normalized Jacobi polynomial on \( [0, 1] \).

Evaluate the Fourier-Jacobi series, using the Clenshaw summation algorithm.

If \( c \) is of length \( n + 1 \), this function returns the value:

Here, \( \jacobi{n}{a}{b}{x} \) is a Jacobi polynomial of degree \( n = |\vec{c}| \), which is orthogonal iff \( (a, b) \in (-1,\ \infty)^2 \) and \( x \in [0,\ 1] \).

Parameters:

Returns:

Numerical inversion of the PPF.

Evaluate the (differential / continuous) entropy \( H(X) \) of a univariate random variable \( X \), from its quantile density function (QDF), \( q(u) = \frac{\mathrm{d} F^{-1}(u)}{\mathrm{d} u} \), with \( F^{-1} \) the inverse of the CDF, i.e. the PPF / quantile function.

The derivation follows from the identity \( f(x) = 1 / q(F(x)) \) of PDF \( f \), specifically:

Parameters:

Returns:

Evaluate the theoretical L-comoment matrix of a multivariate probability distribution, using the joint PDF \(f(\vec x) \equiv f(x_1, x_2, \ldots, x_n)\) of random vector \(\vec{X}\), and the marginal CDFs \(F_k\) of its \(k\)-th random variable.

The L-comoment matrix is defined as

with elements

where \(U_j = F_j(X_j)\) and \(u_j = F_j(x_j)\) denote the (marginal) probability integral transform of \(X_j\) and \(x_j \sim X_j\). Furthermore, \(\widetilde{P}^{(\alpha, \beta)}_k\) is a shifted Jacobi polynomial, and

a positive constant.

For \(r \ge 2\), it can also be expressed as

and without trim (\(s = t = 0\)), this simplifies to

with \(\tilde{P}_n = \tilde{P}^{(0, 0)}_n\) the shifted Legendre polynomial. This last form is precisely the definition introduced by Serfling & Xiao (2007).

Note that the L-comoments along the diagonal, are equivalent to the (univariate) L-moments, i.e.

Examples:

Find the L-coscale and TL-coscale matrices of the multivariate Student’s t distribution with 4 degrees of freedom:

>>> from scipy.stats import multivariate_t
>>> df = 4
>>> loc = np.array([0.5, -0.2])
>>> cov = np.array([[2.0, 0.3], [0.3, 0.5]])
>>> X = multivariate_t(loc=loc, shape=cov, df=df)

>>> from scipy.special import stdtr
>>> std = np.sqrt(np.diag(cov))
>>> cdf0 = lambda x: stdtr(df, (x - loc[0]) / std[0])
>>> cdf1 = lambda x: stdtr(df, (x - loc[1]) / std[1])

>>> l_cov = l_comoment_from_pdf(X.pdf, (cdf0, cdf1), 2)
>>> l_cov.round(4)
array([[1.0413, 0.3124],
       [0.1562, 0.5207]])
>>> tl_cov = l_comoment_from_pdf(X.pdf, (cdf0, cdf1), 2, trim=1)
>>> tl_cov.round(4)
array([[0.4893, 0.1468],
       [0.0734, 0.2447]])

The (Pearson) correlation coefficient can be recovered in several ways:

>>> cov[0, 1] / np.sqrt(cov[0, 0] * cov[1, 1])  # "true" correlation
0.3
>>> l_cov[0, 1] / l_cov[0, 0]
0.3
>>> l_cov[1, 0] / l_cov[1, 1]
0.3
>>> tl_cov[0, 1] / tl_cov[0, 0]
0.3
>>> tl_cov[1, 0] / tl_cov[1, 1]
0.3

Parameters:

Other Parameters:

Returns:

Evaluate the theoretical L-comoment ratio matrix of a multivariate probability distribution, using the joint PDF \(f_{\vec{X}}(\vec{x})\) and \(n\) marginal CDFs \(F_X(x)\) of random vector \(\vec{X}\).

Evaluate the population L-moment of a continuous probability distribution, using its Cumulative Distribution Function (CDF) \(F_X(x) = P(X \le x)\).

where,

\(\widetilde{P}^{(\alpha, \beta)}_k(x)\) the shifted (\(x \mapsto 2x-1\)) Jacobi polynomial, \(H(x)\) the Heaviside step function, and \(I_x(\alpha, \beta)\) the regularized incomplete gamma function, and \(u = F_X(x)\) the probability integral transform of \(x \sim X\).

Examples:

Evaluate the first 4 L- and TL-moments of the standard normal distribution:

>>> from scipy.special import ndtr  # standard normal CDF
>>> l_moment_from_cdf(ndtr, [1, 2, 3, 4])
array([0.        , 0.56418958, 0.        , 0.06917061])
>>> l_moment_from_cdf(ndtr, [1, 2, 3, 4], trim=1)
array([0.        , 0.29701138, 0.        , 0.01855727])

Evaluate the first 4 TL-moments of the standard Cauchy distribution:

>>> def cdf_cauchy(x: float) -> float:
...     return np.arctan(x) / np.pi + 1 / 2
>>> l_moment_from_cdf(cdf_cauchy, [1, 2, 3, 4], trim=1)
array([0.        , 0.69782723, 0.        , 0.23922105])

Parameters:

Other Parameters:

Raises:

Returns:

Evaluate the population L-moment of a univariate probability distribution, using its Percentile Function (PPF), \(x(F)\), also commonly known as the quantile function, which is the inverse of the Cumulative Distribution Function (CDF).

where

and \(\widetilde{P}^{(\alpha, \beta)}_k(x)\) the shifted (\(x \mapsto 2x-1\)) Jacobi polynomial.

Examples:

Evaluate the L- and TL-location and -scale of the standard normal distribution:

>>> from scipy.special import ndtri  # standard normal inverse CDF
>>> l_moment_from_ppf(ndtri, [1, 2])
array([0.        , 0.56418958])
>>> l_moment_from_ppf(ndtri, [1, 2], trim=1)
array([0.        , 0.29701138])

Parameters:

Other Parameters:

Raises:

Returns:

Evaluate the population L-moments \( \tlmoment{s, t}{r} \) for \( r > 1 \) from the quantile distribution function (QDF), which is the derivative of the PPF (quantile function).

Population L-ratio’s from a CDF.

Population L-ratio’s from a PPF.

Calculates the theoretical- / population- L-moments (for \(r \le 2\)) and L-ratio’s (for \(r > 2\)) of a distribution, from its CDF.

By default, the first num = 4 population L-stats are calculated:

• \(\lambda^{(s,t)}_1\) - L-location
• \(\lambda^{(s,t)}_2\) - L-scale
• \(\tau^{(s,t)}_3\) - L-skewness coefficient
• \(\tau^{(s,t)}_4\) - L-kurtosis coefficient

This function is equivalent to l_ratio_from_cdf(cdf, [1, 2, 3, 4], [0, 0, 2, 2], *, **).

Calculates the theoretical- / population- L-moments (for \(r \le 2\)) and L-ratio’s (for \(r > 2\)) of a distribution, from its quantile function.

By default, the first num = 4 population L-stats are calculated:

• \(\lambda^{(s,t)}_1\) - L-location
• \(\lambda^{(s,t)}_2\) - L-scale
• \(\tau^{(s,t)}_3\) - L-skewness coefficient
• \(\tau^{(s,t)}_4\) - L-kurtosis coefficient

This function is equivalent to l_ratio_from_cdf(cdf, [1, 2, 3, 4], [0, 0, 2, 2], *, **).

L-moments that are estimated from \(n\) samples of a distribution with CDF \(F\), converge to the multivariate normal distribution as the sample size \(n \rightarrow \infty\).

Here, \(\vec{l}^{(s, t)} = \left[l^{(s, t)}_r, \dots, l^{(s, t)}_{r_{max}} \right]^T\) is a vector of estimated sample L-moments, and \(\vec{\lambda}^{(s, t)}\) its theoretical (“true”) counterpart.

This function calculates the covariance matrix

where \(u = F_X(x)\) and \(v = F_Y(y)\) (marginal) probability integral transforms, and

the shifted Jacobi polynomial \(p^{(s, t)}_n(u) = P^{(t, s)}_{n-1}(2u - 1)\), \(P^{(t, s)}_m\), and \(w^{(s, t)}(u) = u^s (1 - u)^t\) its weight function.

Parameters:

Other Parameters:

Returns:

Raises:

Similar to l_moment_from_cdf, but for the lmo.l_stats.

As the sample size \(n \rightarrow \infty\), the L-moment ratio’s are also distributed (multivariate) normally. The L-stats are defined to be L-moments for \(r\le 2\), and L-ratio coefficients otherwise.

The corresponding covariance matrix has been found to be

where \(\bf{\Lambda}^{(s, t)}\) is the covariance matrix of the L-moments from l_moment_cov_from_cdf, and \(\tau^{(s,t)}_r = \lambda^{(s,t)}_r / \lambda^{(s,t)}_2\) the population L-ratio.

Parameters:

Other Parameters:

Influence Function (IF) of a theoretical L-moment.

with \(F\) the CDF, \(\tilde{P}^{(s,t)}_{r-1}\) the shifted Jacobi polynomial, and

where \(B\) is the (complete) Beta function.

The proof is trivial, because population L-moments are linear functionals.

Parameters:

Other Parameters:

Returns:

Construct the influence function of a theoretical L-moment ratio.

where the generalized L-moment ratio is defined as

Because IF’s are a special case of the general Gâteuax derivative, the L-ratio IF is derived by applying the chain rule to the L-moment IF.

Parameters:

Other Parameters:

Returns:

Return a PPF (quantile function, or inverse CDF), with the specified. L-moments \( \tlmoment{s, t}{1}, \tlmoment{s, t}{2}, \ldots, \tlmoment{s, t}{R} \). Other L-moments are considered zero.

For \( R \) L-moments, this function returns

where \( \shjacobi{n}{a}{b}{x} \) is an \( n \)-th degree shifted Jacobi polynomial, which is orthogonal for \( (a, b) \in (-1, \infty)^2 \) on \( u \in [0, 1] \).

This nonparametric quantile function estimation method was first described by J.R.M. Hosking in 2007. However, his derivation contains a small, but obvious error, resulting in zero-division for \( r = 1 \). So Lmo derived this correct version himself, by using the fact that L-moments are the disguised coefficients of the PPF’s generalized Fourier-Jacobi series expansion.

With Parseval’s theorem it can be shown that, if the probability-weighted moment \( M_{2,s,t} \) (which is the variance if \( s = t = 0 \)) is finite, then \( \hat{Q}_R(u) = Q(u) \) as \( R \to \infty \).

Parameters:

Returns:

Return the QDF (quantile density function, the derivative of the PPF), with the specified L-moments \( \tlmoment{s, t}{1}, \tlmoment{s, t}{2}, \ldots, \tlmoment{s, t}{R} \). Other L-moments are considered zero.

This function returns

where \( \shjacobi{n}{a}{b}{x} \) is an \( n \)-th degree shifted Jacobi polynomial, which is orthogonal for \( (a, b) \in (-1, \infty)^2 \) on \( u \in [0, 1] \).

See ppf_from_l_moments for options.

Fit the distribution parameters using the (Generalized) Method of L-Moments (L-(G)MM).

The goal is to find the “true” parameters \(\bm{\theta^*}\) of the distribution. In practise, this is done using a reasonably close estimate, \(\bm{\hat\theta}\).

In the (non-Generalized) Method of L-moments (L-MM), this is done by solving the system of equations \(\ell^{(s, t)}_r = \lambda^{(s, t)}_r\), for \(r = 1, \dots, n\), with \(n\) the number of free parameters. Because the amount of parameters matches the amount of L-moment conditions, the solution is point-defined, and can be found using simple least squares.

L-GMM extends L-MM by allowing more L-moment conditions than there are free parameters, \(m > n\). This requires solving an over-identified system of \(m\) equations:

where \(W_m\) is a \(m \times m\) weight matrix.

The weight matrix is initially chosen as the matrix inverse of the non-parametric L-moment covariance matrix, see lmo.l_moment_cov. These weights are then plugged into the the equation above, and fed into scipy.optimize.minimize, to obtain the initial parameter estimates.

In the next step(s), Monte-Carlo sampling is used to draw samples from the distribution (using the current parameter estimates), with sample sizes matching that of the data. The L-moments of these samples are consequently used to to calculate the new weight matrix.

Parameters:

Other Parameters:

Raises:

Returns: