Lmo reference¶

L-kurtosis coefficient; the 4th sample L-moment ratio.

Alias for lmo.l_ratio(a, 4, 2, *, **).

Examples:

>>> import lmo, numpy as np
>>> x = np.random.default_rng(12345).standard_t(2, 99)
>>> lmo.l_kurtosis(x)
0.28912787
>>> lmo.l_kurtosis(x, trim=(1, 1))
0.19928182

L-location (or L-loc): unbiased estimator of the first L-moment, \(\lambda^{(s, t)}_1\).

Alias for lmo.l_moment(a, 1, *, **).

Examples:

The first moment (i.e. the mean) of the Cauchy distribution does not exist. This means that estimating the location of a Cauchy distribution from its samples, cannot be done using the traditional average (i.e. the arithmetic mean). Instead, a robust location measure should be used, e.g. the median, or the TL-location.

To illustrate, let’s start by drawing some samples from the standard Cauchy distribution, which is centered around the origin.

>>> import lmo
>>> import numpy as np
>>> rng = np.random.default_rng(12345)
>>> x = rng.standard_cauchy(200)

The mean and the untrimmed L-location (which are equivalent) give wrong results, so don’t do this:

>>> np.mean(x)
-3.6805
>>> lmo.l_loc(x)
-3.6805

Usually, the answer to this problem is to use the median. However, the median only considers one or two samples (depending on whether the amount of samples is odd or even, respectively). So the median ignores most of the available information.

>>> np.median(x)
0.096825
>>> lmo.l_loc(x, trim=(len(x) - 1) // 2)
0.096825

Luckily for us, Lmo knows how to deal with longs tails, as well – trimming them (specifically, by skipping the first \(s\) and last \(t\) expected order statistics).

Let’s try the TL-location (which is equivalent to the median)

>>> lmo.l_loc(x, trim=1)  # equivalent to `trim=(1, 1)`
0.06522

Estimates the generalized trimmed L-moment \(\lambda^{(s, t)}_r\) from the samples along the specified axis. By default, this will be the regular L-moment, \(\lambda_r = \lambda^{(0, 0)}_r\).

Examples:

Calculate the L-location and L-scale from student-T(2) samples, for different (symmetric) trim-lengths.

>>> import lmo, numpy as np
>>> x = np.random.default_rng(12345).standard_t(2, 99)
>>> lmo.l_moment(x, [1, 2], trim=(0, 0))
array([-0.01412282,  0.94063132])
>>> lmo.l_moment(x, [1, 2], trim=(1/2, 1/2))
array([-0.02158858,  0.5796519 ])
>>> lmo.l_moment(x, [1, 2], trim=(1, 1))
array([-0.0124483 ,  0.40120115])

The theoretical L-locations are all 0, and the the L-scale are 1.1107, 0.6002 and 0.4165, respectively.

Non-parmateric auto-covariance matrix of the generalized trimmed L-moment point estimates with orders r = 1, ..., r_max.

Examples:

Fitting of the cauchy distribution with TL-moments. The location is equal to the TL-location, and scale should be \(0.698\) times the TL(1)-scale, see Elamir & Seheult (2003).

>>> import lmo, numpy as np
>>> rng = np.random.default_rng(12345)
>>> x = rng.standard_cauchy(1337)
>>> lmo.l_moment(x, [1, 2], trim=(1, 1))
array([0.08142405, 0.68884917])

The L-moment estimates seem to make sense. Let’s check their standard errors, by taking the square root of the variances (the diagonal of the covariance matrix):

>>> lmo.l_moment_cov(x, 2, trim=(1, 1))
array([[ 4.89407076e-03, -4.26419310e-05],
       [-4.26419310e-05,  1.30898414e-03]])
>>> np.sqrt(_.diagonal())
array([0.06995764, 0.03617989])

Empirical Influence Function (EIF) of a sample L-moment.

Estimates the generalized L-moment ratio.

Equivalent to lmo.l_moment(a, r, *, **) / lmo.l_moment(a, s, *, **).

The L-moment with r=0 is 1, so the l_ratio(a, r, 0, *, **) is equivalent to l_moment(a, r, *, **).

Examples:

Estimate the L-location, L-scale, L-skewness and L-kurtosis simultaneously:

>>> import lmo
>>> import numpy as np
>>> rng = np.random.default_rng(12345)
>>> x = rng.lognormal(size=99)
>>> lmo.l_ratio(x, [1, 2, 3, 4], [0, 0, 2, 2])
array([1.53196368, 0.77549561, 0.4463163 , 0.29752178])
>>> lmo.l_ratio(x, [1, 2, 3, 4], [0, 0, 2, 2], trim=(0, 1))
array([0.75646807, 0.32203446, 0.23887609, 0.07917904])

Empirical Influence Function (EIF) of a sample L-moment ratio.

Non-parametric estimates of the Standard Error (SE) in the L-ratio estimates from lmo.l_ratio.

Examples:

Estimate the values and errors of the TL-loc, scale, skew and kurtosis for Cauchy-distributed samples. The theoretical values are [0.0, 0.698, 0.0, 0.343] (Elamir & Seheult, 2003), respectively.

>>> import lmo, numpy as np
>>> rng = np.random.default_rng(12345)
>>> x = rng.standard_cauchy(42)
>>> lmo.l_ratio(x, [1, 2, 3, 4], [0, 0, 2, 2], trim=(1, 1))
array([-0.25830513,  0.61738638, -0.03069701,  0.25550176])
>>> lmo.l_ratio_se(x, [1, 2, 3, 4], [0, 0, 2, 2], trim=(1, 1))
array([0.32857302, 0.12896501, 0.13835403, 0.07188138])

L-scale unbiased estimator for the second L-moment, \(\lambda^{(s, t)}_2\).

Alias for lmo.l_moment(a, 2, *, **).

Examples:

>>> import lmo, numpy as np
>>> x = np.random.default_rng(12345).standard_cauchy(99)
>>> x.std()
72.87715244
>>> lmo.l_scale(x)
9.501123995
>>> lmo.l_scale(x, trim=(1, 1))
0.658993279

Unbiased sample estimator for the L-skewness coefficient.

Alias for lmo.l_ratio(a, 3, 2, *, **).

Examples:

>>> import lmo, numpy as np
>>> x = np.random.default_rng(12345).standard_exponential(99)
>>> lmo.l_skew(x)
0.38524343
>>> lmo.l_skew(x, trim=(0, 1))
0.27116139

Calculates the L-loc(ation), L-scale, L-skew(ness) and L-kurtosis.

Equivalent to lmo.l_ratio(a, [1, 2, 3, 4], [0, 0, 2, 2], *, **) by default.

Examples:

>>> import lmo, scipy.stats
>>> x = scipy.stats.gumbel_r.rvs(size=99, random_state=12345)
>>> lmo.l_stats(x)
array([0.79014773, 0.68346357, 0.12207413, 0.12829047])

The theoretical L-stats of the standard Gumbel distribution are [0.577, 0.693, 0.170, 0.150].

Calculates the standard errors (SE’s) of the L-stats.

Equivalent to lmo.l_ratio_se(a, [1, 2, 3, 4], [0, 0, 2, 2], *, **) by default.

Examples:

>>> import lmo, scipy.stats
>>> x = scipy.stats.gumbel_r.rvs(size=99, random_state=12345)
>>> lmo.l_stats(x)
array([0.79014773, 0.68346357, 0.12207413, 0.12829047])
>>> lmo.l_stats_se(x)
array([0.12305147, 0.05348839, 0.04472984, 0.03408495])

The theoretical L-stats of the standard Gumbel distribution are [0.577, 0.693, 0.170, 0.150]. The corresponding relative z-scores are [-1.730, 0.181, 1.070, 0.648].

The coefficient of L-variation (or L-CV) unbiased sample estimator:

Alias for lmo.l_ratio(a, 2, 1, *, **).

Examples:

>>> import lmo, numpy as np
>>> x = np.random.default_rng(12345).pareto(4.2, 99)
>>> x.std() / x.mean()
1.32161112
>>> lmo.l_variation(x)
0.59073639
>>> lmo.l_variation(x, trim=(0, 1))
0.55395044

Sample L-cokurtosis coefficient matrix \(\tilde\Lambda^{(t_1, t_2)}_4\).

Alias for lmo.l_coratio(a, 4, 2, *, **).

L-colocation matrix of 1st L-comoment estimates, \(\Lambda^{(t_1, t_2)}_1\).

Alias for lmo.l_comoment(a, 1, *, **).

Examples:

Without trimming, the L-colocation only provides marginal information:

>>> import lmo, numpy as np
>>> rng = np.random.default_rng(12345)
>>> x = rng.multivariate_normal([0, 0], [[6, -3], [-3, 3.5]], 99).T
>>> lmo.l_loc(x, axis=-1)
array([-0.02678225,  0.03008309])
>>> lmo.l_coloc(x)
array([[-0.02678225, -0.02678225],
       [ 0.03008309,  0.03008309]])

But the trimmed L-locations are a different story…

>>> lmo.l_loc(x, trim=(1, 1), axis=-1)
array([-0.10488868, -0.00625729])
>>> lmo.l_coloc(x, trim=(1, 1))
array([[-0.10488868, -0.03797989],
       [ 0.03325074, -0.00625729]])

What this tells us, is somewhat of a mystery: trimmed L-comoments have been only been briefly mentioned once or twice in the literature.

Multivariate extension of lmo.l_moment.

Estimates the L-comoment matrix:

Whereas the L-moments are calculated using the order statistics of the observations, i.e. by sorting, the L-comoment sorts \(x_i\) using the order of \(x_j\). This means that in general, \(\lambda_{r [ij]}^{(t_1, t_2)} \neq \lambda_{r [ji]}^{(t_1, t_2)}\), i.e. \(\Lambda_{r}^{(t_1, t_2)}\) is not symmetric.

The \(r\)-th L-comoment \(\lambda_{r [ij]}^{(t_1, t_2)}\) reduces to the L-moment if \(i=j\), and can therefore be seen as a generalization of the (univariate) L-moments. Similar to how the diagonal of a covariance matrix contains the variances, the diagonal of the L-comoment matrix contains the L-moments.

Based on the proposed definition by Serfling & Xiao (2007) for L-comoments. Extended to allow for generalized trimming.

Examples:

Estimation of the second L-comoment (the L-coscale) from biviariate normal samples:

>>> import lmo, numpy as np
>>> rng = np.random.default_rng(12345)
>>> x = rng.multivariate_normal([0, 0], [[6, -3], [-3, 3.5]], 99).T
>>> lmo.l_comoment(x, 2)
array([[ 1.2766793 , -0.83299947],
       [-0.71547941,  1.05990727]])

The diagonal contains the univariate L-moments:

>>> lmo.l_moment(x, 2, axis=-1)
array([1.2766793 , 1.05990727])

Estimate the generalized matrix of L-comoment ratio’s.

Sample L-correlation coefficient matrix \(\tilde\Lambda^{(t_1, t_2)}_2\); the ratio of the L-coscale matrix over the L-scale column-vectors.

Alias for lmo.l_coratio(a, 2, 2, *, **).

The diagonal consists of all 1’s.

Where the pearson correlation coefficient measures linearity, the (T)L-correlation coefficient measures monotonicity.

Examples:

>>> import lmo, numpy as np
>>> rng = np.random.default_rng(12345)
>>> cov = np.array([[6, -3], [-3, 3.5]])
>>> x = rng.multivariate_normal([0, 0], [[6, -3], [-3, 3.5]], 99).T
>>> lmo.l_corr(x)
array([[ 1.        , -0.65247355],
       [-0.67503962,  1.        ]])

Let’s compare this with the theoretical correlation

>>> cov[0, 1] / np.sqrt(cov[0, 0] * cov[1, 1])
-0.6546536707079772

and the (Pearson) correlation coefficient matrix:

>>> np.corrcoef(x)
array([[ 1.        , -0.66383285],
       [-0.66383285,  1.        ]])

L-coscale matrix of 2nd L-comoment estimates, \(\Lambda^{(t_1, t_2)}_2\).

Alias for lmo.l_comoment(a, 2, *, **).

Analogous to the (auto-) variance-covariance matrix, the L-coscale matrix is positive semi-definite, and its main diagonal contains the L-scale’s. conversely, the L-coscale matrix is inherently asymmetric, thus yielding more information.

Examples:

>>> import lmo, numpy as np
>>> rng = np.random.default_rng(12345)
>>> x = rng.multivariate_normal([0, 0], [[6, -3], [-3, 3.5]], 99).T
>>> lmo.l_scale(x, trim=(1, 1), axis=-1)
array([0.66698774, 0.54440895])
>>> lmo.l_coscale(x, trim=(1, 1))
array([[ 0.66698774, -0.41025416],
       [-0.37918065,  0.54440895]])

Sample L-coskewness coefficient matrix \(\tilde\Lambda^{(t_1, t_2)}_3\).

Alias for lmo.l_coratio(a, 3, 2, *, **).

Calculates the L-coscale, L-corr(elation), L-coskew(ness) and L-cokurtosis.

Equivalent to lmo.l_coratio(a, [2, 2, 3, 4], [0, 2, 2, 2], *, **).

Statistical hypothesis test for non-normality, using the L-skewness and L-kurtosis coefficients on the sample data..

Adapted from Harri & Coble (2011), and includes Hosking’s correction.

Examples:

Compare the testing power with scipy.stats.normaltest given 10.000 samples from a contaminated normal distribution.

>>> import numpy as np
>>> from lmo.diagnostic import normaltest
>>> from scipy.stats import normaltest as normaltest_scipy
>>> rng = np.random.default_rng(12345)
>>> n = 10_000
>>> x = 0.9 * rng.normal(0, 1, n) + 0.1 * rng.normal(0, 9, n)
>>> normaltest(x)[1]
0.04806618
>>> normaltest_scipy(x)[1]
0.08435627

At a 5% significance level, Lmo’s test is significant (i.e. indicating non-normality), whereas scipy’s test isn’t (i.e. inconclusive).

Goodness-of-fit (GOF) hypothesis test for the null hypothesis that the observed L-moments come from a distribution with the given scipy.stats distribution or cumulative distribution function (CDF).

• H0: The theoretical probability distribution, with the given CDF, is a good fit for the observed L-moments.
• H1: The distribution is not a good fit for the observed L-moments.

The test statistic is the squared Mahalanobis distance between the \(n\) observed L-moments, and the theoretical L-moments. It asymptically (in sample size) follows the \(\chi^2\) distribution, with \(n\) degrees of freedom.

The sample L-moments are expected to be of consecutive orders \(r = 1, 2, \dots, n\). Generally, the amount of L-moments \(n\) should not be less than the amount of parameters of the distribution, including the location and scale parameters. Therefore, it is required to have \(n \ge 2\).

Examples:

Test if the samples are drawn from a normal distribution.

>>> import lmo
>>> import numpy as np
>>> from lmo.diagnostic import l_moment_gof
>>> from scipy.stats import norm
>>> rng = np.random.default_rng(12345)
>>> X = norm(13.12, 1.66)
>>> n = 1_000
>>> x = X.rvs(n, random_state=rng)
>>> x_lm = lmo.l_moment(x, [1, 2, 3, 4])
>>> l_moment_gof(X, x_lm, n).pvalue
0.82597

Contaminated samples:

>>> y = 0.9 * x + 0.1 * rng.normal(X.mean(), X.std() * 10, n)
>>> y_lm = lmo.l_moment(y, [1, 2, 3, 4])
>>> y_lm.round(3)
array([13.193, 1.286, 0.006, 0.168])
>>> l_moment_gof(X, y_lm, n).pvalue
0.0

Analogous to lmo.diagnostic.l_moment_gof, but using the L-stats (see lmo.l_stats).

Returns the absolute upper bounds \(L^{(s,t)}_r\) on L-moments \(\lambda^{(s,t)}_r\), proportional to the scale \(\sigma_X\) (standard deviation) of the probability distribution of random variable \(X\). So \(\left| \lambda^{(s,t)}_r(X) \right| \le \sigma_X \, L^{(s,t)}_r\), given that standard deviation \(\sigma_X\) of \(X\) exists and is finite.

The bounds are derived by applying the Cauchy-Schwarz inequality to the covariance-based definition of generalized trimmed L-moment, for \(r > 1\):

where \(B\) is the Beta function, \(P^{(\alpha, \beta)}_r\) the Jacobi polynomial, and \(F\) the cumulative distribution function of random variable \(X\).

After a lot of work, one can (and one did) derive the closed-form inequality:

for \(r \in \mathbb{N}_{\ge 2}\) and \(s, t \in \mathbb{R}_{\ge 0}\), where \(\Gamma\) is the Gamma function, and \(B\) the multivariate Beta function

For the untrimmed L-moments, this simplifies to

Unlike the standardized product-moments, the L-moment ratio’s with \( r \ge 2 \) are bounded above and below.

Specifically, Hosking derived in 2007 that

But this derivation relies on unnecessarily loose Jacobi polynomial bounds. If the actual min and max of the Jacobi polynomials are used instead, the following (tighter) inequality is obtained:

where

Examples:

Without trim, the lower- and upper-bounds of the L-skewness and L-kurtosis are:

>>> l_ratio_bounds(3)
(-1.0, 1.0)
>>> l_ratio_bounds(4)
(-0.25, 1.0)

For the L-kurtosis, the “legacy” bounds by Hosking (2007) are clearly looser:

>>> l_ratio_bounds(4, legacy=True)
(-1.0, 1.0)

For the symmetrically trimmed TL-moment ratio’s:

>>> l_ratio_bounds(3, trim=3)
(-1.2, 1.2)
>>> l_ratio_bounds(4, trim=3)
(-0.15, 1.5)

Similarly, those of the LL-ratio’s are

>>> l_ratio_bounds(3, trim=(0, 3))
(-0.8, 2.0)
>>> l_ratio_bounds(4, trim=(0, 3))
(-0.233333, 3.5)

The LH-skewness bounds are “flipped” w.r.t to the LL-skewness, but they are the same for the L*-kurtosis:

>>> l_ratio_bounds(3, trim=(3, 0))
(-2.0, 0.8)
>>> l_ratio_bounds(4, trim=(3, 0))
(-0.233333, 3.5)

The bounds of multiple L-ratio’s can be calculated in one shot:

>>> np.stack(l_ratio_bounds([3, 4, 5, 6], trim=(1, 2)))
array([[-1.        , -0.19444444, -1.12      , -0.14945848],
       [ 1.33333333,  1.75      ,  2.24      ,  2.8       ]])

Evaluate the approximate rejection point of an influence function \(\psi_{T|F}(x)\) given a statistical functional \(T\) (e.g. an L-moment) and cumulative distribution function \(F(x)\).

with a \(\epsilon\) a small positive number, corresponding to the tol param of e.g. l_moment_influence , which defaults to 1e-8.

Examples:

The untrimmed L-location isn’t robust, e.g. with the standard normal distribution:

>>> import numpy as np
>>> from scipy.stats import distributions as dists
>>> from lmo.diagnostic import rejection_point
>>> if_l_loc_norm = dists.norm.l_moment_influence(1, trim=0)
>>> if_l_loc_norm(np.inf)
inf
>>> rejection_point(if_l_loc_norm)
nan

For the TL-location of the Gaussian distribution, and even for the Student’s t distribution with 4 degrees of freedom (3 also works, but is very slow), they exist.

>>> influence_norm = dists.norm.l_moment_influence(1, trim=1)
>>> influence_t4 = dists.t(4).l_moment_influence(1, trim=1)
>>> influence_norm(np.inf), influence_t4(np.inf)
(0.0, 0.0)
>>> rejection_point(influence_norm), rejection_point(influence_t4)
(6.0, 206.0)

Evaluate the gross-error sensitivity of an influence function \(\psi_{T|F}(x)\) given a statistical functional \(T\) (e.g. an L-moment) and cumulative distribution function \(F(x)\).

Examples:

Evaluate the gross-error sensitivity of the standard exponential distribution’s LL-skewness (\(\tau^{(0, 1)}_3\)) and LL-kurtosis (\(\tau^{(0, 1)}_4\)) coefficients:

>>> from lmo.diagnostic import error_sensitivity
>>> from scipy.stats import expon
>>> ll_skew_if = expon.l_ratio_influence(3, 2, trim=(0, 1))
>>> ll_kurt_if = expon.l_ratio_influence(4, 2, trim=(0, 1))
>>> error_sensitivity(ll_skew_if, domain=(0, float('inf')))
1.814657
>>> error_sensitivity(ll_kurt_if, domain=(0, float('inf')))
1.377743

Evaluate the local-shift sensitivity of an influence function \(\psi_{T|F}(x)\) given a statistical functional \(T\) (e.g. an L-moment) and cumulative distribution function \(F(x)\).

Represents the effect of shifting an observation slightly from \(x\), to a neighbouring point \(y\). For instance, adding an observation at \(y\) and removing one at \(x\).

Examples:

Evaluate the local-shift sensitivity of the standard exponential distribution’s LL-skewness (\(\tau^{(0, 1)}_3\)) and LL-kurtosis (\(\tau^{(0, 1)}_4\)) coefficients:

>>> from lmo.diagnostic import shift_sensitivity
>>> from scipy.stats import expon
>>> ll_skew_if = expon.l_ratio_influence(3, 2, trim=(0, 1))
>>> ll_kurt_if = expon.l_ratio_influence(4, 2, trim=(0, 1))
>>> domain = 0, float('inf')
>>> shift_sensitivity(ll_skew_if, domain)
0.837735
>>> shift_sensitivity(ll_kurt_if, domain)
1.442062

Let’s compare these with the untrimmed ones:

>>> shift_sensitivity(expon.l_ratio_influence(3, 2), domain)
1.920317
>>> shift_sensitivity(expon.l_ratio_influence(4, 2), domain)
1.047565