Approximations based on the normal distribution

Non-central chi-squared distribution: cdf and sf (Penev)

ctx.chi2_nc_penev_cdf(x, nu, nc)

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

Returns an approximation to the cdf of the (non-)central chi-squared distribution.

[Penev2000` give the following first and second order Wiener germ approximation:

\[F_{\chi^2}\left(n, x; \lambda\right) \thickapprox \Phi \left(\text{sgn}(s) \sqrt{n(s-1)^2(1/(2s) + m- h(1-s)/s) - \log(A(s)) + 2B(s)/n} \right)\]

\[\text{where } m = \lambda/n; \quad h(y) = \frac{(1-y) \log(1-y)+y- \tfrac{1}{2}y^2}{y^2} ; \quad s= \frac{\sqrt{1+4xm/n}-1}{2m}\]

\[A(s) = \frac{1}{s} - \frac{2}{s} \cdot \frac{h(1-s)}{1+2ms}; \quad B(s) = \frac{(1+3m)^2}{9(1+2m)^3}\]

where \(\Phi(\cdot)\) denotes the cdf of the normal distribution (see ).

An example (CDF):

>>> from mpfunlab import fpm, mpm
>>> mpm.dps = 40; x = 12; nu = 10; nc = 30
>>> dx = dec.cdisn_penev(x, nu, nc); mx = mpm.cdisn_penev(x, nu, nc)
>>> ix = ipm.cdisn_penev(x, nu, nc); fx = fpm.cdisn_penev(x, nu, nc)
>>> gx = gmp.cdisn_penev(x, nu, nc); ax = apm.cdisn_penev(x, nu, nc)
>>> mpm.show([dx, mx, ix, fx, gx, ax])
dec:  3.113877423416836055714616090074444149943E-1
mpm:  3.113877423416836055714616090074444149943e-1
ipm:  3.113877423416836055714616090074444149944e-1 (4.147e-38%)
fpm:  3.11387742341684E-01
gmp:  3.113877423416836055714616090074444149943E-01
ipm:  3.113877423416836055714616090074444149943e-1 (4.055e-38%)

(Non-central) chi-squared distribution: qtf and isf (Canal)

ctx.chi2_nc_canal_qtf(LeftTail, RightTail, n)

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

Returns an approximation to the quantile of the (non-)central chi-squared distribution.

Let \(z_\alpha\) and \(\chi^2_{n,\alpha}\) be the \(\alpha\)-quantiles of the standard normal distribution and central chi-square distribution with \(n\) the degrees of freedom. For \(n=1\) and \(n=2\), the following closed form expressions can be used:

\[\chi^2_{1,\alpha} = z^2_{\alpha}, \quad \chi^2_{2,\alpha} = 2 \log(1 - \alpha)\]

At the extreme left tail of the distribution, for small \(x\), the CDF of a \(\chi^2\) variable with \(n\) degrees of freedom can be approximated by the density of a \(\chi^2\) variable with \(n+2\) degrees of freedom:

\[F_{\chi^2}(n,x) \thickapprox 2 f_{\chi^2}(n+2,x).\]

The density of a \(\chi^2\) variable with \(n+2\) degrees of freedom can be inverted in closed form using the Lambert \(W\) function, which leads to the following approximation:

\[\chi^2_{n,\alpha} \thickapprox f^{-1}_{\chi^2}(n+2,\alpha) = -2 W(t)/a , \quad \text{where}\]

\[a=\frac{1}{(n+2)/2-1}, \quad k=\log(\Gamma((n+2)/2), \quad d=a-\log(1-\alpha)+k, \quad t=-a e^{p+d}\]

This approximation is used for \(|t|<0.1\), and the Lambert \(W\) function is approximated as

\[W(x) \thickapprox x - x^2 + \tfrac{3}{2} x^3 - \tfrac{8}{3} x^4 - \tfrac{125}{24} x^5.\]

Otherwise, the quantile is approximated by inverting a formula proposed by Canal [177]:

\[\chi^2_{n,\alpha} \thickapprox n\left( \frac{1}{2}+ \frac{t}{2}- \frac{3}{2t}\right)^6, \quad \text{where}\]

\[t = \left({-5+2L + 2 \sqrt{13-5L+L^2}} \right)^{1/3} , \quad L = 6 \left(m + s \left(z_{\alpha} + a (z^2_{\alpha} - 1) - a^2 (2 z^3_{\alpha} - 5 z_{\alpha}) \right) \right)\]

\[m = \frac{5}{6} - \frac{1}{9n} - \frac{7}{648n^2} - \frac{25}{2187n^3}, \quad s^2 = \frac{1}{18n} + \frac{1}{162n^2} - \frac{37}{11664n^3}, \quad a = \frac{1}{162 \sqrt{2n^3}} = \gamma_1/6.\]

The noncentral quantile is approximated as

\[\chi^2_{n,\lambda,\alpha} \thickapprox (1+b) \chi^2_{n_1,\alpha} , \quad \text{where } n_1= \frac{(n+\lambda)^2}{n+2\lambda} , \quad b = \frac{\lambda}{n+\lambda}\]

An example (CDF):

>>> from mpfunlab import fpm, mpm
>>> mpm.dps = 40; x = 12; nu = 10; nc = 30
>>> dx = dec.chi2_nc_canal_qtf(x, nu, nc); mx = mpm.chi2_nc_canal_qtf(x, nu, nc)
>>> ix = ipm.chi2_nc_canal_qtf(x, nu, nc); fx = fpm.chi2_nc_canal_qtf(x, nu, nc)
>>> gx = gmp.chi2_nc_canal_qtf(x, nu, nc); ax = apm.chi2_nc_canal_qtf(x, nu, nc)
>>> mpm.show([dx, mx, ix, fx, gx, ax])
dec:  3.113877423416836055714616090074444149943E-1
mpm:  3.113877423416836055714616090074444149943e-1
ipm:  3.113877423416836055714616090074444149944e-1 (4.147e-38%)
fpm:  3.11387742341684E-01
gmp:  3.113877423416836055714616090074444149943E-01
ipm:  3.113877423416836055714616090074444149943e-1 (4.055e-38%)

Gamma distribution: qtf and isf (Canal)

ctx.gamma_canal_qtf(LeftTail, RightTail, a)

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

Returns an approximation to the quantile of the gamma distribution.

An example (CDF):

>>> from mpfunlab import fpm, mpm
>>> mpm.dps = 40; x = 12; nu = 10; nc = 30
>>> dx = dec.gamma_canal_qtf(x, nu, nc); mx = mpm.gamma_canal_qtf(x, nu, nc)
>>> ix = ipm.gamma_canal_qtf(x, nu, nc); fx = fpm.gamma_canal_qtf(x, nu, nc)
>>> gx = gmp.gamma_canal_qtf(x, nu, nc); ax = apm.gamma_canal_qtf(x, nu, nc)
>>> mpm.show([dx, mx, ix, fx, gx, ax])
dec:  3.113877423416836055714616090074444149943E-1
mpm:  3.113877423416836055714616090074444149943e-1
ipm:  3.113877423416836055714616090074444149944e-1 (4.147e-38%)
fpm:  3.11387742341684E-01
gmp:  3.113877423416836055714616090074444149943E-01
ipm:  3.113877423416836055714616090074444149943e-1 (4.055e-38%)

F distribution: qtf and isf (Davis)

ctx.fisher_f_davis_qtf(LeftTail, RightTail, m, n)

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

Returns an approximation to the quantile of the F distribution.

An example (CDF):

>>> from mpfunlab import fpm, mpm
>>> mpm.dps = 40; x = 12; nu = 10; nc = 30
>>> dx = dec.fisher_f_davis_qtf(x, nu, nc); mx = mpm.fisher_f_davis_qtf(x, nu, nc)
>>> ix = ipm.fisher_f_davis_qtf(x, nu, nc); fx = fpm.fisher_f_davis_qtf(x, nu, nc)
>>> gx = gmp.fisher_f_davis_qtf(x, nu, nc); ax = apm.fisher_f_davis_qtf(x, nu, nc)
>>> mpm.show([dx, mx, ix, fx, gx, ax])
dec:  3.113877423416836055714616090074444149943E-1
mpm:  3.113877423416836055714616090074444149943e-1
ipm:  3.113877423416836055714616090074444149944e-1 (4.147e-38%)
fpm:  3.11387742341684E-01
gmp:  3.113877423416836055714616090074444149943E-01
ipm:  3.113877423416836055714616090074444149943e-1 (4.055e-38%)

Beta distribution: qtf and isf (Davis)

ctx.beta_davis_qtf(LeftTail, RightTail, a, b)

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

Returns an approximation to the quantile of the beta distribution.

An example (CDF):

>>> from mpfunlab import fpm, mpm
>>> mpm.dps = 40; x = 12; nu = 10; nc = 30
>>> dx = dec.beta_davis_qtf(x, nu, nc); mx = mpm.beta_davis_qtf(x, nu, nc)
>>> ix = ipm.beta_davis_qtf(x, nu, nc); fx = fpm.beta_davis_qtf(x, nu, nc)
>>> gx = gmp.beta_davis_qtf(x, nu, nc); ax = apm.beta_davis_qtf(x, nu, nc)
>>> mpm.show([dx, mx, ix, fx, gx, ax])
dec:  3.113877423416836055714616090074444149943E-1
mpm:  3.113877423416836055714616090074444149943e-1
ipm:  3.113877423416836055714616090074444149944e-1 (4.147e-38%)
fpm:  3.11387742341684E-01
gmp:  3.113877423416836055714616090074444149943E-01
ipm:  3.113877423416836055714616090074444149943e-1 (4.055e-38%)

Pearson’s rho distribution: pdf (Winterbottom)

ctx.pearson_rho_wb_pdf(x)

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

Returns an approximation to the pdf of Pearson’s rho distribution. See also Winterbottom [1605], Winterbottom [1606].

An example (CDF):

>>> from mpfunlab import fpm, mpm
>>> mpm.dps = 40; x = 12; nu = 10; nc = 30
>>> dx = dec.pearson_rho_wb_pdf(x, nu, nc); mx = mpm.pearson_rho_wb_pdf(x, nu, nc)
>>> ix = ipm.pearson_rho_wb_pdf(x, nu, nc); fx = fpm.pearson_rho_wb_pdf(x, nu, nc)
>>> gx = gmp.pearson_rho_wb_pdf(x, nu, nc); ax = apm.pearson_rho_wb_pdf(x, nu, nc)
>>> mpm.show([dx, mx, ix, fx, gx, ax])
dec:  3.113877423416836055714616090074444149943E-1
mpm:  3.113877423416836055714616090074444149943e-1
ipm:  3.113877423416836055714616090074444149944e-1 (4.147e-38%)
fpm:  3.11387742341684E-01
gmp:  3.113877423416836055714616090074444149943E-01
ipm:  3.113877423416836055714616090074444149943e-1 (4.055e-38%)

Pearson’s rho distribution: cdf and sf (Winterbottom)

ctx.pearson_rho_wb_cdf(x)

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

Returns an approximation to the cdf of Pearson’s rho distribution. See also Winterbottom [1605], Winterbottom [1606].

An improved approximation has been found by inverting Winterbottom’s approximation for the confidence limits of \(\rho\):

\[F_R(N, r,\rho) \approx \Phi \left(\frac{k}{6} - \frac{2c}{k} - \frac{b}{3}\right), \quad \text{where} \quad m = N - 1,\]

\[a = \frac{1}{12\sqrt{m^3}} + \frac{6 r^4 - 3r^2 + 2}{48\sqrt{m^5}}, \quad b = \frac{-r3}{6 a m^2}, \quad s=\frac{1}{a\sqrt{m}} + \frac{1+r^2}{4a\sqrt{m^3}} + \frac{11 r^4 - 2r^2 + 1}{32a\sqrt{m^5}}\]

\[t = \frac{Z(r)-Z(\rho)}{a} + \frac{r}{2 a m} + \frac{5 r^3 + 9r}{24 a m^2}, \quad d=t+ \frac{bs}{3} - \frac{2b^3}{27}\]

\[c = s-b^2/3, \quad p=\sqrt{|12c^3+81d^2|}, \quad k=(108d+12p)^{1/3}\]

An example (CDF):

>>> from mpfunlab import fpm, mpm
>>> mpm.dps = 40; x = 12; nu = 10; nc = 30
>>> dx = dec.pearson_rho_wb_cdf(x, nu, nc); mx = mpm.pearson_rho_wb_cdf(x, nu, nc)
>>> ix = ipm.pearson_rho_wb_cdf(x, nu, nc); fx = fpm.pearson_rho_wb_cdf(x, nu, nc)
>>> gx = gmp.pearson_rho_wb_cdf(x, nu, nc); ax = apm.pearson_rho_wb_cdf(x, nu, nc)
>>> mpm.show([dx, mx, ix, fx, gx, ax])
dec:  3.113877423416836055714616090074444149943E-1
mpm:  3.113877423416836055714616090074444149943e-1
ipm:  3.113877423416836055714616090074444149944e-1 (4.147e-38%)
fpm:  3.11387742341684E-01
gmp:  3.113877423416836055714616090074444149943E-01
ipm:  3.113877423416836055714616090074444149943e-1 (4.055e-38%)

Pearson’s rho distribution: qtf and isf (Winterbottom)

ctx.pearson_rho_wb_qtf(x)

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

Returns an approximation to the qtf of Pearson’s rho distribution. See also Winterbottom [1605], Winterbottom [1606].

Asymptotic expansions typically rely on the Fisher \(z\)-transform \(Z(a)= \text{atanh}(a)\) and its inverse \(Z^{-1}(a) = \tanh(a)\).

Let \(m = N - 1\) and let \(u_\alpha\) = \(\Phi^{-1}(\alpha)\) be the lower \(100\alpha\) percentage point of the standard normal distribution. An approximation to the \(100\alpha\) percentage point \(r_\alpha\) and the lower \(100(1-\alpha)\) confidence limit on \(\rho\), \(\rho_L\), is then obtained by \(r_\alpha \approx Z^{-1}(y_1)\) , where

\begin{eqnarray} y_1 & = & Z(\rho) + \frac{u_\alpha}{\sqrt{m}} + \frac{\rho}{2m} + \frac{u_\alpha^3+3(3-\rho^2)u_\alpha}{12\sqrt{m^3}} + \frac{4\rho^3 u_\alpha^2 + 15\rho-\rho^3}{24m^2} \\ && +\: \frac{u_\alpha^5+(-60\rho^4+30\rho^2+80)u_\alpha^3 + (45\rho^4-21\rho^2+375)u_\alpha}{480\sqrt{m^5}}, \quad \text{and} \nonumber \end{eqnarray}

An example (CDF):

>>> from mpfunlab import fpm, mpm
>>> mpm.dps = 40; x = 12; nu = 10; nc = 30
>>> dx = dec.pearson_rho_wb_qtf(x, nu, nc); mx = mpm.pearson_rho_wb_qtf(x, nu, nc)
>>> ix = ipm.pearson_rho_wb_qtf(x, nu, nc); fx = fpm.pearson_rho_wb_qtf(x, nu, nc)
>>> gx = gmp.pearson_rho_wb_qtf(x, nu, nc); ax = apm.pearson_rho_wb_qtf(x, nu, nc)
>>> mpm.show([dx, mx, ix, fx, gx, ax])
dec:  3.113877423416836055714616090074444149943E-1
mpm:  3.113877423416836055714616090074444149943e-1
ipm:  3.113877423416836055714616090074444149944e-1 (4.147e-38%)
fpm:  3.11387742341684E-01
gmp:  3.113877423416836055714616090074444149943E-01
ipm:  3.113877423416836055714616090074444149943e-1 (4.055e-38%)

Pearson’s rho distribution: confidence limit for \(\rho\) (Winterbottom)

ctx.pearson_rho_wb_cl(x)

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

Returns an approximation to the confidence interval (CI) for the noncentrality parameter \(\rho\) of Pearson’s rho distribution, as decribed in Winterbottom [1606].

We use the Fisher \(z\)-transform \(Z(a)= \text{atanh}(a)\) and its inverse \(Z^{-1}(a) = \tanh(a)\).

Let \(m = N - 1\) and let \(u_\alpha\) = \(\Phi^{-1}(\alpha)\) be the lower \(100\alpha\) percentage point of the standard normal distribution. An approximation to the \(100\alpha\) percentage point \(r_\alpha\) and the lower \(100(1-\alpha)\) confidence limit on \(\rho\), \(\rho_L\), is then obtained by \(\rho_L \approx Z^{-1}(y_2)\), where

\begin{eqnarray} y_2 & = & Z(r) + \frac{u_\alpha}{\sqrt{m}} - \frac{r}{2m} + \frac{u_\alpha^3+3(1+r^2)u_\alpha}{12\sqrt{m^3}} - \frac{4r^3 u_\alpha^2 + 5r^3+9r}{24m^2} \\ && +\: \frac{u_\alpha^5+(60r^4-30r^2+20)u_\alpha^3 + (165r^4+30r^2+15)u_\alpha}{480\sqrt{m^5}}. \nonumber \end{eqnarray}

An example (CDF):

>>> from mpfunlab import fpm, mpm
>>> mpm.dps = 40; x = 12; nu = 10; nc = 30
>>> dx = dec.pearson_rho_wb_cl(x, nu, nc); mx = mpm.pearson_rho_wb_cl(x, nu, nc)
>>> ix = ipm.pearson_rho_wb_cl(x, nu, nc); fx = fpm.pearson_rho_wb_cl(x, nu, nc)
>>> gx = gmp.pearson_rho_wb_cl(x, nu, nc); ax = apm.pearson_rho_wb_cl(x, nu, nc)
>>> mpm.show([dx, mx, ix, fx, gx, ax])
dec:  3.113877423416836055714616090074444149943E-1
mpm:  3.113877423416836055714616090074444149943e-1
ipm:  3.113877423416836055714616090074444149944e-1 (4.147e-38%)
fpm:  3.11387742341684E-01
gmp:  3.113877423416836055714616090074444149943E-01
ipm:  3.113877423416836055714616090074444149943e-1 (4.055e-38%)

Singly noncentral t: pdf (Broda)

ctx.student_t_nc_broda_pdf(x)

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

Returns an approximation to the pdf of the singly noncentral t distribution, as decribed in Witkovský [1612], using the approximation to the cdf of the singly noncentral t distribution as decribed Broda and Paolella [167].

\[f_{\text{StudentT}}\left(x, n, \delta\right) = \frac{n}{x} \left( F_{\text{StudentT}}\left(x \sqrt{(n+2)/n}, n+2, \delta\right) - F_{\text{StudentT}}\left(x, n, \delta\right) \right), \quad x \ne 0,\]

\[f_{\text{StudentT}}\left(x, n, \delta\right) = \frac{\Gamma((n+1)/2)}{\sqrt{n \pi} \Gamma(n/2)} e^{- \frac{1}{2} \delta^2}\]

An example (CDF):

>>> from mpfunlab import fpm, mpm
>>> mpm.dps = 40; x = 12; nu = 10; nc = 30
>>> dx = dec.student_t_nc_broda_pdf(x, nu, nc); mx = mpm.student_t_nc_broda_pdf(x, nu, nc)
>>> ix = ipm.student_t_nc_broda_pdf(x, nu, nc); fx = fpm.student_t_nc_broda_pdf(x, nu, nc)
>>> gx = gmp.student_t_nc_broda_pdf(x, nu, nc); ax = apm.student_t_nc_broda_pdf(x, nu, nc)
>>> mpm.show([dx, mx, ix, fx, gx, ax])
dec:  3.113877423416836055714616090074444149943E-1
mpm:  3.113877423416836055714616090074444149943e-1
ipm:  3.113877423416836055714616090074444149944e-1 (4.147e-38%)
fpm:  3.11387742341684E-01
gmp:  3.113877423416836055714616090074444149943E-01
ipm:  3.113877423416836055714616090074444149943e-1 (4.055e-38%)

Singly noncentral t: cdf, sf (Broda)

ctx.student_t_nc_broda_cdf(x)

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

Returns an approximation to the cdf of the singly noncentral t distribution, as decribed in Broda and Paolella [167].

Singly noncentral t: A saddlepoint approximation based on the joint cumulant generating function of \((x_1, x_2)\), with \(x_1 \sim N(\mu, 1)\) and \(x_2 \sim \chi^2(n)\), is possible (see Broda (2007)), and takes the following form:

\[F_{t'}(x;n,\mu) = \Phi(w)+\phi(w)\left( \frac{1}{w} - \frac{d}{u} \right) + O(n^{-3/2}) , \quad \text{where}\]

\[s = \left(\mu x + \sqrt{4n(x^2+n)+\mu^2 x^2}\right) / (2x^2+2n),\]

\[t_1 = -\mu + x s, \quad t_2 = -x t_1 / (2 n s), \quad d = 1 / (t_1 s),\]

\[u = \sqrt{(\mu x s + 2n ) / (2n)} / s, \quad w = \text{sgn} (x - \mu) \sqrt{-\mu t_1 - 2n \log(s)},\]

and \(\Phi(\cdot)\) and \(\phi(\cdot)\) denote the cdf and pdf of the normal distribution, respectively.

For \(\mu = 0\), these equations simplify to

\[s = \sqrt{n / (x^2+n)}, \quad d = 1 / (x s), \quad u = 1 / s, \quad w = \text{sgn} (x) \sqrt{- 2n \log(s)},\]

An example (CDF):

>>> from mpfunlab import fpm, mpm
>>> mpm.dps = 40; x = 12; nu = 10; nc = 30
>>> dx = dec.student_t_nc_broda_cdf(x, nu, nc); mx = mpm.student_t_nc_broda_cdf(x, nu, nc)
>>> ix = ipm.student_t_nc_broda_cdf(x, nu, nc); fx = fpm.student_t_nc_broda_cdf(x, nu, nc)
>>> gx = gmp.student_t_nc_broda_cdf(x, nu, nc); ax = apm.student_t_nc_broda_cdf(x, nu, nc)
>>> mpm.show([dx, mx, ix, fx, gx, ax])
dec:  3.113877423416836055714616090074444149943E-1
mpm:  3.113877423416836055714616090074444149943e-1
ipm:  3.113877423416836055714616090074444149944e-1 (4.147e-38%)
fpm:  3.11387742341684E-01
gmp:  3.113877423416836055714616090074444149943E-01
ipm:  3.113877423416836055714616090074444149943e-1 (4.055e-38%)

Singly noncentral t: qtf, isf (Harley)

ctx.student_t_nc_harley_qtf(x)

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

Returns an approximation to the qtf of the singly noncentral t distribution, as decribed in Harley [380].

The qtf of the noncentral t-distribution can be approximated by the qtf of Pearson’s rho as follows: Let \(N=f+2\) and

\[\rho = \Re \left( \delta \sqrt{\frac{2}{2N-3 + \delta^2}} \right).\]

Determine \(r_{\alpha;\rho;N} = r_{\alpha}\), using Winterbottom’s approximation. Calculate \(t_{\alpha;\delta;f} = t_{\alpha}\) as

\[t_{\alpha} = \Re \left(\frac{r_{\alpha}}{1-r_{\alpha}^2} \sqrt{\frac{2f(1-\rho^2)}{1-\rho^2}} \right)\]

Note that in these equations both \(\rho\) and \(r_{\alpha}\) can assume values \(>1\). Luckily, Winterbottom’s approximation can handle this; we always take only the real part of the functions results to avoid issues with the branch cuts of the complex square root.

An example (CDF):

>>> from mpfunlab import fpm, mpm
>>> mpm.dps = 40; x = 12; nu = 10; nc = 30
>>> dx = dec.student_t_nc_broda_qtf(x, nu, nc); mx = mpm.student_t_nc_broda_qtf(x, nu, nc)
>>> ix = ipm.student_t_nc_broda_qtf(x, nu, nc); fx = fpm.student_t_nc_broda_qtf(x, nu, nc)
>>> gx = gmp.student_t_nc_broda_qtf(x, nu, nc); ax = apm.student_t_nc_broda_qtf(x, nu, nc)
>>> mpm.show([dx, mx, ix, fx, gx, ax])
dec:  3.113877423416836055714616090074444149943E-1
mpm:  3.113877423416836055714616090074444149943e-1
ipm:  3.113877423416836055714616090074444149944e-1 (4.147e-38%)
fpm:  3.11387742341684E-01
gmp:  3.113877423416836055714616090074444149943E-01
ipm:  3.113877423416836055714616090074444149943e-1 (4.055e-38%)

Singly noncentral t: confidence limit for \(\delta\) (Akahira)

ctx.student_t_nc_akahira_cl(x)

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

Returns an approximation to the confidence interval (CI) for the noncentrality parameter of the singly noncentral t distribution, as decribed in Akahira et al. [6].

Let \(T\) be a statistic according to the non-central \(t\)-distribution with \(n\) degrees of freedom and a non-centrality parameter \(\delta\). Then the lower confidence limit \(\widehat{\delta}\) of level \(1-\alpha\) and the two-sided confidence interval \(( \underline{\delta},\overline{\delta})\) of the non-centrality parameter \(\delta\) of level \(1-\alpha\) are given by Akahira et al. [6]:

\[\widehat{\delta} = bT - z_\alpha \sqrt{k} + h T^3 (z_\alpha^2 - 1)/k,\]

\[\underline{\delta} = bT - z_{\alpha/2} \sqrt{k} + h T^3 (z_{\alpha/2}^2 - 1)/k,\]

\[\overline{\delta} = bT + z_{\alpha/2} \sqrt{k} - h T^3 (z_{\alpha/2}^2 - 1)/k,\]

\[h=\frac{1}{24}\left(\frac{1}{n^2}+ \frac{1}{4n^3}\right) \quad b=\sqrt{\frac{2}{n}}\frac{\Gamma(\tfrac{1}{2} n+\tfrac{1}{2})}{\Gamma(\tfrac{1}{2} n)}\]

where \(k=1+(1-b^2)T^2\), and \(z_\alpha\) denotes the \(\alpha\)-quantile of the normal distribution.

An example (CDF):

>>> from mpfunlab import fpm, mpm
>>> mpm.dps = 40; x = 12; nu = 10; nc = 30
>>> dx = dec.student_t_nc_akahira_cl(x, nu, nc); mx = mpm.student_t_nc_akahira_cl(x, nu, nc)
>>> ix = ipm.student_t_nc_akahira_cl(x, nu, nc); fx = fpm.student_t_nc_akahira_cl(x, nu, nc)
>>> gx = gmp.student_t_nc_akahira_cl(x, nu, nc); ax = apm.student_t_nc_akahira_cl(x, nu, nc)
>>> mpm.show([dx, mx, ix, fx, gx, ax])
dec:  3.113877423416836055714616090074444149943E-1
mpm:  3.113877423416836055714616090074444149943e-1
ipm:  3.113877423416836055714616090074444149944e-1 (4.147e-38%)
fpm:  3.11387742341684E-01
gmp:  3.113877423416836055714616090074444149943E-01
ipm:  3.113877423416836055714616090074444149943e-1 (4.055e-38%)

Doubly noncentral t: cdf, sf (Broda)

ctx.student_t_nc2_broda_cdf(x)

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

Returns an approximation to the cdf of the doubly noncentral t distribution, as decribed in Broda and Paolella [167].

A saddlepoint approximation based on the joint cumulant generating function of \((x_1, x_2)\), with \(x_1 \sim N(\mu, 1)\) and \(x_2 \sim \chi^2(n, \theta)\), is possible (see Broda (2007)), and takes the following form:

\[F_{t''}(x;n;\mu,\theta) = \Phi(w)+\phi(w)\left( \frac{1}{w} - \frac{d}{u} \right) + O(n^{-3/2}) , \quad \text{where}\]

\[a = x^4 + 2 n x^2 + n^2, \quad c_2 = (-2 x^3 \mu - 2 x n \mu)/a, \quad c_1 = (x^2 \mu^2 - n x^2 - n^2 - \theta n) / a, \quad c_0 = (x n \mu)/a,\]

\[q= \tfrac{1}{3}c_1 - \tfrac{1}{9}c_2^2, \quad r= \tfrac{1}{6}(c_1c_2-3c_0)- \tfrac{1}{27}c_2^3,\]

\[s = \sqrt{-4q} \cos \left(\tfrac{1}{3} \arccos\left(r/ \sqrt{-q^3}\right) \right) - \tfrac{1}{3} c_2,\]

\[t_1 = -\mu + x s, t_2 = -x t_1 / (2 n s), \quad d = 1 / (t_1 s), \quad \nu = 1 / (1 - 2 t_2), \quad \alpha = \mu / \sqrt{1 + \theta / n)},\]

\[u = \sqrt{(x^2 + 2 n t_2) (2 n \nu^2 + 4 \theta \nu^3) + 4 n^2 s^2} / (2 n s^2), \quad w = \text{sgn} (x - \alpha) \sqrt{-\mu t_1 - n \log(\nu) - 2 \theta \nu t_2},\]

and \(\Phi(\cdot)\) and \(\phi(\cdot)\) denote the cdf and pdf of the normal distribution, respectively.

An example (CDF):

>>> from mpfunlab import fpm, mpm
>>> mpm.dps = 40; x = 12; nu = 10; nc = 30
>>> dx = dec.student_t_nc2_broda_cdf(x, nu, nc); mx = mpm.student_t_nc2_broda_cdf(x, nu, nc)
>>> ix = ipm.student_t_nc2_broda_cdf(x, nu, nc); fx = fpm.student_t_nc2_broda_cdf(x, nu, nc)
>>> gx = gmp.student_t_nc2_broda_cdf(x, nu, nc); ax = apm.student_t_nc2_broda_cdf(x, nu, nc)
>>> mpm.show([dx, mx, ix, fx, gx, ax])
dec:  3.113877423416836055714616090074444149943E-1
mpm:  3.113877423416836055714616090074444149943e-1
ipm:  3.113877423416836055714616090074444149944e-1 (4.147e-38%)
fpm:  3.11387742341684E-01
gmp:  3.113877423416836055714616090074444149943E-01
ipm:  3.113877423416836055714616090074444149943e-1 (4.055e-38%)

Doubly noncentral t: qtf, isf (Broda)

ctx.student_t_nc2_broda_qtf(x)

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

Returns an approximation to the qtf of the doubly noncentral t distribution, as decribed in Broda and Paolella [167].

Broda_2007_WorkingPaper proposes the following approximation:

\[F_{t''}(t;n;\mu,\theta) \approx F_{t'}(t/c;f;\mu), \quad \text{where}\]

\[f=\frac{7}{2} \left(-1+\sqrt{15-7g^2/h}\right)^{-1}, \quad c=\sqrt{h(1-2/f)}, \quad g=m_1/\sqrt{\mu^2/2}, h=m_2/(1+\mu^2),\]

\[m_1 = \mu \binom{n}{2}^{1/2} \frac{\Gamma((n-1)/2)}{\Gamma(n/2)} {}_1F_1 \left(\frac{1}{2}, \frac{n}{2}, -\frac{\theta}{2} \right)\]

\[m_2 = (1+\mu^2) \frac{n}{n-2} {}_1F_1 \left(\frac{2}{2}, \frac{n}{2}, -\frac{\theta}{2} \right)\]

An example (CDF):

>>> from mpfunlab import fpm, mpm
>>> mpm.dps = 40; x = 12; nu = 10; nc = 30
>>> dx = dec.student_t_nc2_broda_qtf(x, nu, nc); mx = mpm.student_t_nc2_broda_qtf(x, nu, nc)
>>> ix = ipm.student_t_nc2_broda_qtf(x, nu, nc); fx = fpm.student_t_nc2_broda_qtf(x, nu, nc)
>>> gx = gmp.student_t_nc2_broda_qtf(x, nu, nc); ax = apm.student_t_nc2_broda_qtf(x, nu, nc)
>>> mpm.show([dx, mx, ix, fx, gx, ax])
dec:  3.113877423416836055714616090074444149943E-1
mpm:  3.113877423416836055714616090074444149943e-1
ipm:  3.113877423416836055714616090074444149944e-1 (4.147e-38%)
fpm:  3.11387742341684E-01
gmp:  3.113877423416836055714616090074444149943E-01
ipm:  3.113877423416836055714616090074444149943e-1 (4.055e-38%)

Spearman’s rho, first 8 cumulants (David)

ctx.spearman_mu8(N)

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

Returns an approximation to the central distribution of Spearman’s \(\rho\), based on the first 8 moments.

Spearman’s \(\rho\) is calculated like Pearsons’s correlation coefficient, using the rank-transform on the \(x_i\) and \(y_i\). The exact test for \(H_0\) is the permutation test. Techniques for obtaining the exact distribution and approximations are given below

The first 8 cumulants are given by David et al. [206]

\[\kappa_2 = \frac{1}{n-1}\]

\[\kappa_4 = \frac{-6(19n^2+5n-36)}{25n(n+1)(n-1)^3}\]

\[\kappa_6 = \frac{48(583n^6+723n^5-2603n^4-2637n^3+4054n^2+2760n-1800)}{(245n^3(n-1)^5n1^3)}\]

\begin{eqnarray} \kappa_8 & = &\frac{144(41939n^{10}-83709n^9+304254n^8+578442n^7-1012323n^6- 1690125n^5)}{875n^5 (n+1)^5 (n-1)^7} \quad \\ & + &\frac{144(1800776n^4+2358048n^3-1616688n^2- 1080567n+846720)}{875n^5 (n+1)^5 (n-1)^7} \nonumber \end{eqnarray}

Note: See PermCumulants.SpearmanCum for further details.

An example (CDF):

>>> from mpfunlab import fpm, mpm
>>> mpm.dps = 40; x = 12; nu = 10; nc = 30
>>> dx = dec.spearman_mu8(x, nu, nc); mx = mpm.spearman_mu8(x, nu, nc)
>>> ix = ipm.spearman_mu8(x, nu, nc); fx = fpm.spearman_mu8(x, nu, nc)
>>> gx = gmp.spearman_mu8(x, nu, nc); ax = apm.spearman_mu8(x, nu, nc)
>>> mpm.show([dx, mx, ix, fx, gx, ax])
dec:  3.113877423416836055714616090074444149943E-1
mpm:  3.113877423416836055714616090074444149943e-1
ipm:  3.113877423416836055714616090074444149944e-1 (4.147e-38%)
fpm:  3.11387742341684E-01
gmp:  3.113877423416836055714616090074444149943E-01
ipm:  3.113877423416836055714616090074444149943e-1 (4.055e-38%)

Mann-Whitney U distribution: general alternatives specified by rank order probabilities (Sundrum)

ctx.mannwhitney_nc_mu4(N)

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

Returns an approximation to the noncentral distribution of Mann-Whitney’s \(U\), based on the first 4 moments.

The first 4 moments under the alternative are given by Sundrum [536] (note that \(\kappa_4 = \mu_4 - 3\mu2^3\)):

\[\mu_1 = p\]

\[\mu_2 = (p^2 + p - q - r)/mn + (q - p^2 )/n + (r - p^2 )/m + p^2\]

\begin{eqnarray} \mu_3 & = & 6(p^3 + u - pq - pr)/mn + (2p^3 + s - 3pq)/n^2 + (2p^3 + t - 3pr)/m^2 \\ \nonumber & + & 3(3pq + 2(pr - u - p^3) + q - s - p^2)/mn^2 + 3(3pr + 2(pq - u - p^3) \\ \nonumber & + & r - t - p^2)/m^2n +(4p^3 + 3p^2 + p + 6u + 2(s + t) - 3(q + r)(1 + 2p))/m^2n^2 \nonumber \end{eqnarray}

\begin{eqnarray} \mu_4 & = & 3(q - p^2)2/n^2 + 6(q - p^2)(r - p^2)/mn + 3(r - p^2)2/m^2 \\ \nonumber & + & (12qp^2 + a - 4sp - 3q^2 - 6p^4)/n^3 + (12rp^2 + b - 4tp - 3r2 - 6p^4)/m^3 \\ \nonumber & + & 6(7(rp^2 - p^4) + 12qp^2 + qp + 2(w-sp) - 3(q^2 + qr) - 8up - p^3)/mn^2 \\ \nonumber & + & (42qp^2 + 72rp^2 + 6rp+ 12x - 42p^4 - 18r2 - 18qr - 12tp - 48up - 6p^3)/m^2n \\ \nonumber & + & 6(6(p^4 - rp^2) + 3(q^2 - qp) - 12qp^2 + 4sp + 8up + 2(p^3 + qr - w) + s - a)/mn^3 \\ \nonumber & + & (36p^4 + 18r2 + 12qr - 72rp^2 - 36qp^2 + 24tp - 6b + 48up \\ \nonumber & - & 12x + 12p^3 - 18rp+ 6t)/m^3 n \\ \nonumber & + & (105p^4 + 42p^3 + 3p^2 + 33q^2 + 33r2 + 54qr - 174qp^2 - 174rp^2 - 42pq \\ \nonumber & - & 42pr + 36sp + 36tp + 192up - 36w - 36x + 6v + 36u)/ m^2n^2 \\ \nonumber & + & (132qp^2 + 108rp^2 - 66p^4 - 33q^2 - 36qr - 18r2 - 44sp - 24tp + 11a \\ \nonumber & - & 144up + 36w + 24x - 6v - 36p^3 - 36u - 7p^2 + 54pq + 36pr - 18s + 7q)/m^2n^3 \\ \nonumber & + & (132rp^2 + 108qp^2 - 66p^4 - 33r2 - 36qr - 18q^2 - 44tp - 24sp + 11b \\ \nonumber & - & 144up + 24w + 36x - 6v - 36p^3 - 36u - 7p^2 + 54pr + 36pq - 18t + 7r)/m^3n^2 \\ \nonumber & + & (6(3(q^2 + r2)- 12((q + r)p^2) + 4(p^3 + qr + sp + tp - w - x) - (a + b - v) + 16up \\ \nonumber & - & 6(p^4 + pq + pr + u) + 2(s + t)) - 7(q + r - p^2) + p)/ m^3n^3 \nonumber \end{eqnarray}

The parameters \(p, q, r, s, t, v, u, a, b, w\) and \(x\) can be calculated from the following rank order probabilities:

\begin{eqnarray} p & = & P_{1,1}(0,0) \\ \nonumber q & = & P_{2,1}(0,0,1) \\ \nonumber r & = & P_{1,2}(0,1,1) \\ \nonumber s & = & P_{3,1}(0,0,0,1) \\ \nonumber t & = & P_{1,3}(0,1,1,1) \\ \nonumber v & = & P_{2,2}(0,0,1,1) \\ \nonumber u & = & v + (1/4) P_{2,2}(0,1,0,1) \\ \nonumber a & = & P_{4,1}(0,0,0,0,1) \\ \nonumber b & = & P_{1,4}(0,1,1,1,1) \\ \nonumber w & = & 2a + (2/3) P_{3,2}(0,0,1,0,1) + (1/6) P_{3,2}(0,1,0,0,1) \\ \nonumber x & = & 2b + (2/3) P_{2,3}(0,1,0,1,1) + (1/6) P_{2,3}(0,1,1,0,1) \\ \nonumber \end{eqnarray}

The rank order probabilities which are required for the calculation of the first 4 moments can be estimated from the sample as follows: Let \(U_i\) be the number of \(Y`s in the sample greater than `X_{(i)}\), where \(X_{(i)}\) is the ith ordered values of the \(X\) amongst themselves. Then the probabilities are calculated according to the following scheme :

\[P_{1,1}(0,1) = \frac{1}{mn} \sum_{i=1}^m U_i\]

\[P_{2,1}(0,0,1) = \frac{2}{mn(m-1)} \sum_{i=1}^m (i-1)U_i\]

\[P_{3,1}(0,0,0,1) = \frac{3}{mn(m-1)(m-2)} \sum_{i=1}^m (i-1)(i-2)U_i\]

\[P_{4,1}(0,0,0,0,1) = \frac{4}{mn(m-1)(m-2)(m-3)} \sum_{i=1}^m (i-1)(i-2)(i-3)U_i\]

\[P_{1,2}(0,1,1) = \frac{1}{mn(n-1)} \sum_{i=1}^m U_i(U_i-1)\]

\[P_{1,3}(0,1,1,1) = \frac{1}{mn(n-1)(n-2)} \sum_{i=1}^m U_i(U_i-1)(U_i-2)\]

\[P_{1,4}(0,1,1,1,1) = \frac{1}{mn(n-1)(n-2)(n-3)} \sum_{i=1}^m U_i(U_i-1)(U_i-2)(U_i-3)\]

\[P_{2,2}(0,0,1,1) = \frac{2}{mn(m-1)(n-1)} \sum_{i=1}^m (i-1)U_i(U_i-1)\]

\[P_{3,2}(0,0,0,1,1) = \frac{3}{mn(m-1)(m-2)(n-1)} \sum_{i=1}^m (i-1)(i-2)U_i(U_i-1)\]

\[P_{2,3}(0,0,1,1,1) = \frac{2}{mn(m-1)(n-1)(n-2)} \sum_{i=1}^m (i-1)U_i(U_i-1)(U_i-2)\]

\[P_{2,2}(0,1,0,1) = \frac{4}{mn(m-1)(n-1)} \sum_{i=1}^m \sum_{j=i+1}^m (U_i - U_j) U_j\]

\[P_{3,2}(0,0,1,0,1) = \frac{12}{mn(m-1)(m-2)(n-1)} \sum_{i=1}^m \sum_{j=i+1}^m (i-1)(U_i - U_j) U_j\]

\[P_{3,2}(0,1,0,0,1) = \frac{12}{mn(m-1)(m-2)(n-1)} \sum_{i=1}^m \sum_{j=i+1}^m (i-j-1)(U_i - U_j) U_j\]

\[P_{2,3}(0,1,0,1,1) = \frac{6}{mn(m-1)(n-1)(n-2)} \sum_{i=1}^m \sum_{j=i+1}^m (U_i - U_j) U_j (U_j-1)\]

\[P_{2,3}(0,1,0,1,1) = \frac{6}{mn(m-1)(n-1)(n-2)} \sum_{i=1}^m \sum_{j=i+1}^m (U_i - U_j) U_j (U_j - 1)\]

\[P_{2,3}(0,1,1,0,1) = \frac{6}{mn(m-1)(n-1)(n-2)} \sum_{i=1}^m \sum_{j=i+1}^m (U_i - U_j) (U_i - U_j - 1) U_j\]

An example (CDF):

>>> from mpfunlab import fpm, mpm
>>> mpm.dps = 40; x = 12; nu = 10; nc = 30
>>> dx = dec.mannwhitney_nc_mu4(x, nu, nc); mx = mpm.mannwhitney_nc_mu4(x, nu, nc)
>>> ix = ipm.mannwhitney_nc_mu4(x, nu, nc); fx = fpm.mannwhitney_nc_mu4(x, nu, nc)
>>> gx = gmp.mannwhitney_nc_mu4(x, nu, nc); ax = apm.mannwhitney_nc_mu4(x, nu, nc)
>>> mpm.show([dx, mx, ix, fx, gx, ax])
dec:  3.113877423416836055714616090074444149943E-1
mpm:  3.113877423416836055714616090074444149943e-1
ipm:  3.113877423416836055714616090074444149944e-1 (4.147e-38%)
fpm:  3.11387742341684E-01
gmp:  3.113877423416836055714616090074444149943E-01
ipm:  3.113877423416836055714616090074444149943e-1 (4.055e-38%)

First 4 moments of Kendalls \(\tau\) in the general case (Sundrum)

ctx.kendall_tau_nc_mu4(N)

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

Returns an approximation to the noncentral distribution of Kendall’s \(\tau\), based on the first 4 moments.

Consider \(k\) pairs of random variables \((x_1,y_1),\ldots,(x_k,y_k)\). Let \(r_1,\ldots,r_k\) be a permutaton of \(1,\ldots,k\) and let \(s_1,\ldots,s_k\) be its reciprocal. We define as in Snow [530].

\(P(r_1,\ldots,r_k)\) = Pr[if \(x\)’s are ranked in order, ranks of corresponding \(y\)’s are \(r_1,\ldots,r_k\)].

\(P(r_1,\ldots,r_k)\) = Pr[if \(y\)’s are ranked in order, ranks of corresponding \(x\)’s are \(s_1,\ldots,s_k\)].

Let \(H(t) = 1\) if \(t>0\) and \(H(t)=0\) otherwise. Then

\[P(r_1,\ldots,r_k) = k! P \left[ \prod_{i=2}^k H(x_i - x_{i-1}) H(y_{s_i} - y_{s_{i-1}}) \right] = k! P \left[ \prod_{i=2}^k H(y_i - y_{i-1}) H(x_{r_i} - x_{r_{i-1}}) \right]\]

The rank order probabilities can then be estimated from the sample as

\[\hat{P}_{(12)} = \frac{1}{n(n-1)} \sum_{i \ne j} I\left[ (x_i-x_j)(y_i-y_j) > 0 \right]\]

\[\hat{P}_{(123)} = \frac{1}{n(n-1)(n-2)} \sum_{i \ne j_1 \ne j_2} I\left[ (x_i-x_{j_1})(y_i-y_{j_1}) (x_i-x_{j_2})(y_i-y_{j_2}) > 0 \right]\]

The first 4 moments of \(T_N\) are given as Sundrum [535]:

\[\mu_1 = p\]

\[\mu_2 = 2 [p(1 - p) + 2(k - p^2)(n - 2)] / [n (n - 1)]\]

\begin{eqnarray} \mu_3 & = & 4 [(t - 18kp + 10p^3) n^2 + (6k + 2u - 6p^2 - 5t + 72kp - 34p^3) n \\ \nonumber & + & (p + 9p^2 + 30p^3 - 12k - 4u + 6t - 72kp)] / [n (n - 1)]^2 \nonumber \end{eqnarray}

\begin{eqnarray} \mu_4 & = & 8 [6(k - p^2)2 n^4 + (6kp + 2y - 6p^3 - 16tp - 84k2 + 270kp^2 - 108p4) n^3 \\ \nonumber & + & (1.5p^2 + 6t + 2b - 126kp - 24up - 18y + 75p^3 + 120tp + 426k2 - 1446kp^2 + 505.5p4) n^2 \\ \nonumber & + & (14k + 12u - 15.5p^2 - 30t - 10b + 444kp + 96up + 52y - 213p^3 - 296tp - 924k2 \\ \nonumber & + & 2988kp^2 - 943.5p4)n \\ \nonumber & + & (p - 28k - 24u + 21p^2 + 36t + 12b - 432kp - 96up - 48y + 180p^3 + 240tp + 720k2 \\ \nonumber & - & 2160kp^2 + 630p4)] / [n (n - 1)]^3 \nonumber \end{eqnarray}

The parameters \(p, k, u, t, b, y\) can be calculated from the following rank order probabilities, writing \(P_{(1423)}\) for the probility that the ranks of the \(Y\) occur in this order when the \(X\) are sorted in ascending order:

\begin{eqnarray} p & = & P_{(12)} \\ \nonumber u & = & P_{(123)} \\ \nonumber k & = & P_{(123)} + \tfrac{1}{4}P_{(132)} \\ \nonumber t & = & 8(P_{(1234)} + P_{(1243))} + 6P_{(1342)} + 4P_{(1324)} + 2P_{(2143)} + P_{(2413)} + P_{(1432)} \\ \nonumber b & = & 15P_{(1234)} + 10P_{(1243)} + 5P_{(1324)} + 4P_{(1342)} + P_{(2143)} \\ \nonumber y & = & (160P_{(12453)} + 150(P_{(12435)} + P_{(12354))} + 125P_{(12345)} + 96P_{(21453)} + 90P_{(13254)} + 84P_{(13524)} \\ \nonumber &+& 80P_{(13425)} + 64P_{(13452)} + 45P_{(21354)} + 40P_{(12543)} + 32(P_{(23514)} + P_{(14352)} + P_{(13542)) } \\ \nonumber &+& 24P_{(21543)} + 22P_{(24153)} + 20P_{(14325)} + 18P_{(14523)} + 16P_{(25143)} + 12(P_{(24513)} + P_{(14532))} \\ \nonumber &+& 6(P_{(25314)} + P_{(15342))} + 4P_{(25413)} + 2(P_{(15432)} + P_{(35142)}))/5 \nonumber \end{eqnarray}

When \((x_i,y_i)\), \(i=1,2,\ldots,k\), are normal variables with correlation \(\rho\), then closed form expressions exist for \(p, k\) and \(u\) (writing \(a = \arcsin\left(\rho\right)/\pi\) and \(b = 2\arcsin\left(\tfrac{1}{2}\rho\right)/\pi\):

\begin{eqnarray} p & = & \tfrac{1}{2} + a \\ \nonumber k & = & \tfrac{1}{4} (\tfrac{10}{9} + 4a + 4a^2 - b^2) \\ \nonumber u & = & \tfrac{3}{8} (\tfrac{4}{9} + 4a -2b + 4a^2 - b^2) \nonumber \end{eqnarray}

The parameters \(t, b\) and \(y\) have to be evaluated by numerical intergration of the rank order probabilities

An example (CDF):

>>> from mpfunlab import fpm, mpm
>>> mpm.dps = 40; x = 12; nu = 10; nc = 30
>>> dx = dec.kendall_tau_nc_mu4(x, nu, nc); mx = mpm.kendall_tau_nc_mu4(x, nu, nc)
>>> ix = ipm.kendall_tau_nc_mu4(x, nu, nc); fx = fpm.kendall_tau_nc_mu4(x, nu, nc)
>>> gx = gmp.kendall_tau_nc_mu4(x, nu, nc); ax = apm.kendall_tau_nc_mu4(x, nu, nc)
>>> mpm.show([dx, mx, ix, fx, gx, ax])
dec:  3.113877423416836055714616090074444149943E-1
mpm:  3.113877423416836055714616090074444149943e-1
ipm:  3.113877423416836055714616090074444149944e-1 (4.147e-38%)
fpm:  3.11387742341684E-01
gmp:  3.113877423416836055714616090074444149943E-01
ipm:  3.113877423416836055714616090074444149943e-1 (4.055e-38%)