Distribution of the modified likelihood ratio test (LRT) for a given covariance matrix and mean vector

class ctx.dist_lrt_x0_s0(p, n)

where ctx is fpm, mpm, ipm, dec, gmp or apm.

The distribution of the modified (LRT) for a given covariance matrix and mean vector is a continuous probability distribution with \(p \ge 1\) predictor variables, error degress of freedom \(n \ge 1\), and the support interval \((0,1)\). See also Anderson [9], Davis [211]. For tables see Nagarsenker 1974

The likelihood criterion \(\lambda_2\) for testing the hypothesis that a sample of size \(N\) is drawn from a \(p\)-variate normal population with a given covariance matrix \(\boldsymbol{\Sigma}_0\) and mean vector \(\boldsymbol{\mu}_0\) is given by

\[\lambda_2 = \left(\frac{e}{N} \right)^{\tfrac{1}{2}pN} \left|\boldsymbol{S \Sigma}_0^{-1}\right|^{\tfrac{1}{2}N} \exp \left(-\tfrac{1}{2}\left[\text{tr} \left(\boldsymbol{S \Sigma}_0^{-1}\right) +N(\boldsymbol{\bar{x}}-\boldsymbol{\mu}_0)' \boldsymbol{\Sigma}_0 (\boldsymbol{\bar{x}}-\boldsymbol{\mu}_0)\right]\right)\]

where \(\boldsymbol{S}\) is the sample covariance matrix and \(\boldsymbol{\bar{x}}\) is the sample mean vector.

dist_lrt_x0_s0.pdf(x)

Returns \(\text{pdf}_X(x)\), the probability density function (pdf) of a random variable \(X\), following the distribution of Box’s LRTs.

The pdf can be calculated (in principle in arbitrary precision) by numerical inversion of the characteristic function, using the algorithm by Gil-Pelaez. The PDF of Y is the inverse Fourier transform of its characteristic function,

\[\text{pdf}_X(x) = \frac{1}{\pi} \int_{0}^{\infty} \Re \left ( e^{-itx} C_X(t) \right ) \mathrm{d} t.\]

where \(\Re (z)\) denotes the real part of \(z\).

>>> from mpfunlab import *
>>> mp.dps = 30
>>> mu = 0; sigma = 1; x = 3;
>>> print ("pdf: ", lrt_x0_s0(mu, sigma).pdf(x))
6.3563523462564525615615615614561356E-20

dist_lrt_x0_s0.cdf(x)

Returns \(\text{cdf}_X(x)\), the cumulative distribution function (cdf) of a random variable \(X\), following the distribution of Box’s LRTs.

The cdf can be calculated (in principle in arbitrary precision) by numerical inversion of the characteristic function, using the algorithm by Gil-Pelaez. Gil-Pelaez derived the following inversion formula which requires integration of a real-valued function, only. In particular,

\[\text{cdf}_X(x) = \frac{1}{2} - \frac{1}{\pi} \int_{0}^{\infty} \Im \left ( \frac{ e^{-itx} C_X(t)}{t} \right ) \mathrm{d} t.\]

where \(\Im (z)\) denotes the imaginary part of \(z\).

>>> from mpfunlab import *
>>> mp.dps = 30
>>> mu = 0; sigma = 1; x = 3;
>>> print ("cdf: ", lrt_x0_s0(mu, sigma).pdf(x))
6.3563523462564525615615615614561356E-20

dist_lrt_x0_s0.sf(x)

Returns \(\text{sf}_X(x)\), the survival function (sf) of a random variable \(X\), following the distribution of Box’s LRTs.

The sf can be calculated (in principle in arbitrary precision) by numerical inversion of the characteristic function, using the algorithm by Gil-Pelaez. Gil-Pelaez derived the following inversion formula which requires integration of a real-valued function, only. In particular,

\[\text{sf}_X(x) = \frac{1}{2} + \frac{1}{\pi} \int_{0}^{\infty} \Im \left ( \frac{ e^{-itx} C_X(t)}{t} \right ) \mathrm{d} t.\]

where \(\Im (z)\) denotes the imaginary part of \(z\).

>>> from mpfunlab import *
>>> mp.dps = 30
>>> mu = 0; sigma = 1; x = 3;
>>> print (" sf: ", lrt_x0_s0(mu, sigma).pdf(x))
sf: 6.3563523462564525615615615614561356E-20

dist_lrt_x0_s0.qtf(q)

Returns \(\text{qtf}_X(x)\), the quantile function (qtf) of a random variable \(X\), following the distribution of Box’s LRTs:

There is no known closed form for the quantile function \(\text{cdf}^{-1}_X(q)\): It is computed with Newton iterations where the starting values are from a central chi-square approximation.

>>> from mpfunlab import *
>>> mp.dps = 30
>>> mu = 0; sigma = 1; q = 0.3;
>>> print ("qtf: ", lrt_x0_s0(mu, sigma).qtf(q))
qtf: 6.3563523462564525615615615614561356E+00

dist_lrt_x0_s0.isf(q)

Returns \(\text{isf}_X(q)\), the inverse survival function (isf) of a random variable \(X\), following the distribution of Box’s LRTs:

\[\text{isf}_X(q) = \text{qtf}_X(1-q) = \text{cdf}^{-1}_X(1-q).\]

>>> from mpfunlab import *
>>> mp.dps = 30
>>> mu = 0; sigma = 1; q = 0.3;
>>> print ("isf: ", lrt_x0_s0(mu, sigma).isf(q))
6.3563523462564525615615615614561356E+00

dist_lrt_x0_s0.c_x(t)

Returns \(C_X(t)\), the characteristic function of a random variable \(X\), following the distribution of Box’s LRTs:

The characteristic function of \(M = -2 \log (\lambda_2)\) has the form (Nagarsenker 1974, equation 2.3)

\[C_X(t) = \frac{\left(2e/N\right)^{-it \cdot pN}}{(1-2it)^{-\tfrac{1}{2}pN(1-2it)}} \prod_{j=1}^p \frac{ \Gamma\left[\tfrac{1}{2}N(1-2it)- j/2\right]}{ \Gamma \left[\tfrac{1}{2}(N-j)\right]}\]

>>> from mpfunlab import *
>>> mp.dps = 30
>>> mu = 0; sigma = 1; t = 0.3;
>>> print ("c_x: ", lrt_x0_s0(mu, sigma).c_x(t))
6.3563523462564525615615615614561356E+00

dist_lrt_x0_s0.m_x(t)

Returns \(M_X(t)\), the moment generating function of a random variable \(X\), following the distribution of Box’s LRTs:

\[M_X(t) = \frac{\left(2e/N\right)^{-t \cdot pN}}{(1-2t)^{-\tfrac{1}{2}pN(1-2t)}} \prod_{j=1}^p \frac{ \Gamma\left[\tfrac{1}{2}N(1-2t)- j/2\right]}{ \Gamma \left[\tfrac{1}{2}(N-j)\right]}\]

>>> from mpfunlab import *
>>> mp.dps = 30
>>> mu = 0; sigma = 1; t = 0.3;
>>> print ("c_x: ", lrt_x0_s0(mu, sigma).c_x(t))
6.3563523462564525615615615614561356E+00

dist_lrt_x0_s0.k_x(s, k=0)

Returns \(K_X(s)\), the cumulant generating function, and its \(j^{\text{th}}\) derivatives, \(K_X^{(j)}(s), j = 1 \ldots k\), of a random variable \(X\), following the distribution of Box’s LRTs:

\[M_{X}(t) = \log(M_{X}(t))\]

>>> from mpfunlab import *
>>> mp.dps = 30
>>> mu = 0; sigma = 1; t = 0.3; k = 6;
>>> print ("c_x: ", lrt_x0_s0(mu, sigma).k_x(t, k))
6.3563523462564525615615615614561356E+00

dist_lrt_x0_s0.moments(k)

Returns the first \(j\) raw moments, \(\mu_j, j = 1 \ldots k\), of a random variable \(X\), following the distribution of Box’s LRTs. The moments are calculated from the cumulants.

>>> from mpfunlab import *
>>> mp.dps = 30
>>> mu = 0; sigma = 1; k = 6;
>>> print ("saddlepoint: ", lrt_x0_s0(mu, sigma).moments(k))
6.3563523462564525615615615614561356E+00

dist_lrt_x0_s0.cumulants(k)

Returns the first \(j\) cumulants, \(\kappa_j, j = 1 \ldots k\), of a random variable \(X\), following the distribution of Box’s LRTs

The cumulants \(\kappa_j\) of \(M = -2 \log(W)\) are given by

\[\kappa_1 = tbd\]

\[\kappa_j = tbd\]

where \(\psi(\cdot)\) is the digamma function, and its derivatives \(\psi^{(j)}(\cdot)\) are polygamma functions.

>>> from mpfunlab import *
>>> mp.dps = 30
>>> mu = 0; sigma = 1; k = 6;
>>> print ("saddlepoint: ", lrt_x0_s0(mu, sigma).cumulants(k))
6.3563523462564525615615615614561356E+00

Additional information

We have for the Box-Davis expansion (see Davis 1971, equations 23-24):

\[f=\tfrac{1}{2} p(p+3); \quad \rho=1-\frac{2p^2+9p+11}{6N(p+3)}\]

\[\omega_r = \frac{2(-1)^r}{r(r+1)(r+2) \rho^r} \sum_{s=1}^{r+1} \binom{r+2}{s+1} (1-\rho-N^{-1})^{r+1-s} \frac{\delta_s+p(s+1)/2^{s+1}}{(\tfrac{1}{2}N)^{s-1}},\]

where \(\delta\) is defined in Box (1949). For tables see Nagarsenker 1974

Approximations

ctx.lrt_x0_vc0_bd(q, f, rho, omega)

where ctx is ipm, dec, mpm, or gmp.

We have for the Box-Davis expansion (see Davis 1971, equations 23-24):

\[f=\tfrac{1}{2} p(p+3); \quad \rho=1-\frac{2p^2+9p+11}{6N(p+3)}; \quad \delta_{s,p} = -\frac{s+1}{2} \sum_{j=0}^{p-1} B_s(-j/2)\]

\[\omega_r = \frac{2(-1)^r}{r(r+1)(r+2) \rho^r} \sum_{s=1}^{r+1} \binom{r+2}{s+1} (1-\rho-N^{-1})^{r+1-s} \frac{\delta_{s,p} + p(s+1)/2^{s+1}}{(\tfrac{1}{2}N)^{s-1}},\]

where \(B_s(\cdot)\) is the Bernoulli polynomial of degree \(s\). . For tables see Nagarsenker 1974

ctx.lrt_x0_vc0_bd_inv(q, f, rho, omega)

where ctx is ipm, dec, mpm, or gmp.

Calculates the Box-Davis approximation to the qtf and isf.