Distribution of the product of independent beta variables#

class ctx.dist_beta_product(p, bi, ci)#

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

The betaproduct distribution is a continuous probability distribution, defined as the product of \(p \ge 1\) independent beta variables, which have a beta distribution with \(b_j\) and \(c_j-b_j\) degrees of freedom, and the support interval \((0,1)\). See also Anderson [9], Muirhead [440], Butler [172], Ginzberg [364], pages 92-105, Marques et al. [433], Tang and Gupta [539].

Text explaining the concepts.

See also Olkin and Press [455], Wilks [1604], Geisser [360], Marques et al. [433], Mauchly [434], Coelho et al. [187].

dist_beta_product.pdf(x)#

Returns \(\text{pdf}_X(x)\), the probability density function (pdf) of a random variable \(X\), following the distribution of the product of independent beta variables. The pdf can expressed as a Meijer G function as follows (PhamGia_2008, p. 1702, Ginzberg_2013, p. 97 + 98):

\[\begin{split}\text{pdf}_X(x) = \frac{K e^{-x/2}}{2 \alpha} G^{m,0}_{m,m} \left( \left. \begin{matrix} a_1, \dots, a_n ; a_{n+1} \dots a_p \\ b_1, \dots, b_m ; b_{m+1} \dots b_q \end{matrix}\; \right| \; e^{\frac{-x}{2 \alpha}} \right)\end{split}\]
>>> from mpfunlab import *
>>> mp.dps = 30
>>> mu = 0; sigma = 1; x = 3;
>>> print ("pdf: ", beta_product(mu, sigma).pdf(x))
6.3563523462564525615615615614561356E-20

dist_beta_product.cdf(x)#

Returns \(\text{cdf}_X(x)\), the cumulative distribution function (cdf) of a random variable \(X\), following the distribution of the product of independent beta variables. The cdf can expressed as a Meijer G function as follows (PhamGia_2008, p. 1702, Ginzberg_2013, p. 97 + 98):

\[\begin{split}\text{cdf}_X(x) = = K e^{-x/2} G^{m+1,0}_{m+1,m+1} \left( \left. \begin{matrix} a_1, \dots, a_n ; a_{n+1} \dots a_p \\ b_1, \dots, b_m ; b_{m+1} \dots b_q \end{matrix}\; \right| \; e^{\frac{-x}{2 \alpha}} \right)\end{split}\]
>>> from mpfunlab import *
>>> mp.dps = 30
>>> mu = 0; sigma = 1; x = 3;
>>> print ("cdf: ", beta_product(mu, sigma).pdf(x))
6.3563523462564525615615615614561356E-20

dist_beta_product.sf(x)#

Returns \(\text{sf}_X(x)\), the survival function (sf) of a random variable \(X\), following the distribution of the product of independent beta variables. The cdf can expressed as a Meijer G function as follows (PhamGia_2008, p. 1702, Ginzberg_2013, p. 97 + 98):

\[\begin{split}\text{sf}_X(x) = = 1 - K e^{-x/2} G^{m+1,0}_{m+1,m+1} \left( \left. \begin{matrix} a_1, \dots, a_n ; a_{n+1} \dots a_p \\ b_1, \dots, b_m ; b_{m+1} \dots b_q \end{matrix}\; \right| \; e^{\frac{-x}{2 \alpha}} \right)\end{split}\]
>>> from mpfunlab import *
>>> mp.dps = 30
>>> mu = 0; sigma = 1; x = 3;
>>> print (" sf: ", beta_product(mu, sigma).pdf(x))
sf: 6.3563523462564525615615615614561356E-20

dist_beta_product.qtf(q)#

Returns \(\text{qtf}_X(x)\), the quantile function (qtf) of a random variable \(X\), following the distribution of the product of independent beta variables:

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 Nagarsenker’s approximation.

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

dist_beta_product.isf(q)#

Returns \(\text{isf}_X(q)\), the inverse survival function (isf) of a random variable \(X\), following the distribution of the product of independent beta variables:

There is no known closed form for the quantile function \(\text{isf}^{-1}_X(q)\): It is computed with Newton iterations where the starting values are from Nagarsenker’s approximation.

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

dist_beta_product.c_x(t)#

Returns \(C_X(t)\), the characteristic function of a random variable \(X\), following the distribution of the product of independent beta variables:

\[C_X(t) = \int_{0}^{\infty} e^{i tx} \text{pdf}_X(x) \mathrm{d} x\]
>>> from mpfunlab import *
>>> mp.dps = 30
>>> mu = 0; sigma = 1; t = 0.3;
>>> print ("c_x: ", beta_product(mu, sigma).c_x(t))
6.3563523462564525615615615614561356E+00

dist_beta_product.m_x(t)#

Returns \(M_X(t)\), the moment generating function of a random variable \(X\), following the distribution of the product of independent beta variables:

\[M_X(t) = \int_{0}^{\infty} e^{tx} \text{pdf}_X(x) \mathrm{d} x\]
>>> from mpfunlab import *
>>> mp.dps = 30
>>> mu = 0; sigma = 1; t = 0.3;
>>> print ("c_x: ", beta_product(mu, sigma).c_x(t))
6.3563523462564525615615615614561356E+00

dist_beta_product.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 the product of independent beta variables:

\[K_X(s) = \log\left(M_X(t)\right)\]

where \(\Gamma_p(\cdot)\) is the multivariate gamma function and \(\psi(\cdot)\) is the digamma function.

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

dist_beta_product.moments(k)#

Returns the first \(j\) raw moments, \(\mu_j, j = 1 \ldots k\), of a random variable \(X\), following the distribution of the product of independent beta variables:

\[\mu'_X(r) = \prod_{j=1}^p \frac{\Gamma(c_j)\Gamma(b_j +h)}{\Gamma(b_j)\Gamma(c_j+h)},.\]
>>> from mpfunlab import *
>>> mp.dps = 30
>>> mu = 0; sigma = 1; k = 6;
>>> print ("saddlepoint: ", beta_product(mu, sigma).moments(k))
6.3563523462564525615615615614561356E+00

dist_beta_product.cumulants(k)#

Returns the first \(j\) cumulants, \(\kappa_j, j = 1 \ldots k\), of a random variable \(X\), following the distribution of the product of independent beta variables. The cumulants are calculated from the moments.

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

Additional information

ctx.log_beta_prod_gp(k)#

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

Calculates the pdf, cdf and sf from the characteristic function using the procedure of Gil-Pelaez (see gil_pelaez_pdf() and gil_pelaez_cdf()).

This uses \(U = \log 2 W\).

For a large subset of the random variables \(W\) defined in in the previous section there exists a power transformation \(U=W^{2/N}\), such that

\[\operatorname{E}[U^h] = \prod_{j=1}^p \frac{\Gamma(c_j)\Gamma(b_j +h)}{\Gamma(b_j)\Gamma(c_j+h)},\]

i.e. \(U\) is distributed as the product of \(p\) independent random variables which have a beta distribution with \(b_j\) and \(c_j-b_j\) degrees of freedom. This allows to specify these distributions in a simpler form, just giving the \(b_j\) and \(c_j\). Its density defined on \((0,1)\)

These are functions for which the l.r.t. \(M/2\) is distributed as the product of \(p\) independent random variables which have a beta distribution with \(b_j\) and \(c_j-b_j\) degrees of freedom.

The characteristic function has the general form:

\[\phi_{M/2}(t) = \prod_{j=1}^p \frac{\Gamma(c_j)\Gamma(b_j - it)}{\Gamma(b_j)\Gamma(c_j -it)}\]

The distribution of the product of independent beta random variables can also be expressed as a Box-Davis distribution of the Box type with the following parameter choices for \(a, b, x_k, \xi_k, y_j\) and \(\eta_j\):

\[a = b = p, \quad x_k = y_j = N/2, \quad \xi_k = b_k - N/2, \quad \eta_j = c_j - N/2, \quad j,k= 1 \ldots p,\]
ctx.beta_prod_spa(x, n, results='c')#

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

Calculates the Luggannini-Rice saddlepoint approximation of the pdf, cdf and sf.

This uses \(2\log W\).

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`2log W` of central Wilks \(W\):

\[K_X(s) = \log \left[ \frac{\Gamma_p(n/2 + s)\Gamma_p((n + m)/2)}{\Gamma_p(n/2)\Gamma_p((n + m)/2 + s)} \right].\]
\[K'(s) = \sum_{i=1}^p \left[\psi \left( \tfrac{1}{2}n+s - \tfrac{1}{2}(i-1)\right) - \psi \left( \tfrac{1}{2}(n+m)+s - \tfrac{1}{2}(i-1)\right) \right]\]

where \(\Gamma_p(\cdot)\) is the multivariate gamma function and \(\psi(\cdot)\) is the digamma function.

ctx.beta_prod_spa_inv(x, n, results='qtf')#

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

Calculates the inverse Jensen saddlepoint approximation of the qtf and isf.

This uses \(2\log W\).

ctx.beta_product_bd(x, f, rho, omega)#

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

Calculates the Box-Davis approximation to the pdf, cdf and sf.

\[f = -2 \sum_{j=1}^p (b_j -c_j); \quad \theta = \frac{N}{2};\]
\[\rho = 1 - \frac{1}{f \cdot \theta} \sum_{j=1}^p \bigl( B_2(b_j - \theta) - B_2(c_j - \theta) \bigr)\]
\[\omega_r = \frac{(-1)^{r+1}}{r(r+1) (\rho \theta)^r} \sum_{j=1}^p \bigl( B_{r+1}(b_j - \rho \theta) - B_{r+1}(c_j - \rho \theta) \bigr)\]
ctx.beta_product_bd_inv(q, f, rho, omega)#

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

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

ctx.beta_product_bd_cdf_old(x, f, rho, omega)#

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

Calculates the Box-Davis approximation to the pdf, cdf and sf.

For the product-of-independent-beta-variates distributions, the parameters of the Box-Davis expansion are given in equations ….

The distribution of the product of independent beta random variables can also be expressed as a Box-Davis distribution of the Box type with the following parameter choices for \(a, b, x_k, \xi_k, y_j\) and \(\eta_j\):

\[a = b = p, \quad x_k = y_j = N/2, \quad \xi_k = b_k - N/2, \quad \eta_j = c_j - N/2, \quad j,k= 1 \ldots p,\]
\[\theta = \min_{i,j}(x_i, y_j),\]
\[f = -2 \left[\sum_{k=1}^a \xi_k - \sum_{j=1}^b \eta_j -\tfrac{1}{2}(a-b) \right],\]
\[\rho = 1-\frac{1}{f} \left[\sum_{k=1}^a x_k^{-1}\left(\xi_k^2-\xi_k+\tfrac{1}{6}\right) - \sum_{j=1}^b y_j^{-1}\left(\eta_j^2-\eta_j+\tfrac{1}{6}\right) \right],\]
\[\omega_r = \frac{(-1)^{r+1}}{r(r+1)} \left[\sum_{k=1}^a \frac{B_{r+1}((1-\rho)x_k+\xi_k)}{(\rho x_k)^r} - \sum_{j=1}^b \frac{B_{r+1}((1-\rho)y_j+\eta_j)}{(\rho y_j)^r} \right],\]
\begin{eqnarray} \log(K_B) & = & \log(K) + \frac{1}{2}(m-k) \log(2\pi) - \frac{f}{2}\log(\rho) \\ && +\: \sum_{i=1}^{m} \left( x_i + \xi_i - \frac{1}{2} \right)\log(x_i) - \sum_{i=1}^{k} \left( y_i + \eta_i - \frac{1}{2} \right)\log(x_i) \nonumber \end{eqnarray}