Multivariate uniform distribution of Cook and Johnson (1981) is a joint distribution of uniform variables over \((0,1]\) and its probability density is given as
$$f(x_1, \cdots, x_k) = \frac{\Gamma(a+k)}{\Gamma(a)a^k}\prod_{i=1}^{k} x_i^{(-1/a)-1} \left[\sum_{i=1}^{k} x_i^{-1/a} - k +1 \right]^{-(a+k)},$$
where \(0 < x_i <=1, a>0, i=1,\cdots, k\). In fact, Cook-Johnson's uniform distribution is also called Clayton copula (Nelsen, 2006).
Cumulative distribution function \(F(x_1, \dots, x_k)\) is given as
$$F(x_1, \cdots, x_k) = \left[ \sum_{i=1}^{k} x_i^{-1/a} - k + 1 \right]^{-a}.$$
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 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).$$
Random numbers \(X_1, \cdots, X_k\) from Cook-Johnson<U+2019>s multivariate uniform distribution can be generated through transformation of multivariate Lomax random variables \(Y_1, \cdots, Y_k\) by letting \(X_i = (1+\theta_i Y_i)^{-a}, i = 1, \cdots, k\); see Nayak (1987).