Multivariate inverted beta distribution is an alternative expression of multivariate F distribution and is a special case of multivariate Lomax distribution (Balakrishnan and Lai, 2009). Its probability density is given as
$$f(x_1, \cdots, x_p) = \frac{\Gamma(\sum_{i=1}^{p} l_i + a) \prod_{i=1}^{p} x_i^{l_i-1}}{\Gamma(a) [\prod_{i=1}^{p} \Gamma(l_i)] (1+\sum_{i=1}^{p} x_i)^{\sum_{i=1}^{p} l_i + a}},$$
where \(x_i>0, a>0, l_i>0, i=1,\cdots, p\).
Cumulative distribution function \(F(x_1, \dots, x_k)\) is obtained by multiple integral
$$F(x_1, \dots, x_k) = \int_{0}^{x_1} \cdots \int_{0}^{x_k} f(y_1, \cdots, y_k) dy_k \cdots dy_1.$$
This multiple integral is calculated by either adaptive multivariate integration using hcubature in package cubature (Narasimhan et al., 2018) or via Monte Carlo method.
Equicoordinate quantile is obtained by solving the following equation for \(q\) through the built-in one dimension root finding function uniroot:
$$\int_{0}^{q} \cdots \int_{0}^{q} f(x_1, \cdots, x_k) dx_k \cdots dx_1 = p,$$
where \(p\) is the given cumulative probability.
The survival function \(\bar{F}(x_1, \cdots, x_k)\) is obtained either by the following formula related to cumulative distribution function \(F(x_1, \dots, x_k)\) (Joe, 1997)
$$\bar{F}(x_1, \cdots, x_k) = 1 + \sum_{S \in \mathcal{S}} (-1)^{|S|} F_S(x_j, j \in S),$$
or via Monte Carlo method.
Random numbers \(X_1, \cdots, X_k\) from multivariate inverted beta distribution can be generated through parameter substitutions in simulation of generalized multivariate Lomax distribution by letting \(\theta_i = 1, i = 1, \cdots, k\).