Es natural utilizar una cópula gaussiana para esta construcción. Esto equivale a transformar las distribuciones marginales de un $d$ -de una variable aleatoria gaussiana en marginales Beta especificados. Los detalles se dan a continuación.
En realidad, la pregunta sólo describe $2d + d(d-1)/2$ parámetros: dos parámetros $a_i, b_i$ para cada distribución marginal Beta y el $d(d-1)/2$ coeficientes de correlación $\rho_{i,j}=\rho_{j,i},$ $1\le i \lt j \le d.$ Estas últimas determinan la matriz de covarianza $\Sigma$ de la variable aleatoria gaussiana $Z$ (que también podría tener marginales estandarizados y, por tanto, tiene varianzas unitarias en la diagonal). Es convencional escribir
$$\mathbf Z \sim \mathcal{N}(\mathbf 0, \Sigma).$$
Por lo tanto, escribir $\Phi$ para la función de distribución Normal estándar (su cdf) y $F_{a,b}^{-1}$ para la Beta $(a,b)$ función cuantílica, definir
$$X_i = F_{a_i, b_i}^{-1}\left(\Phi(Z_i)\right).$$
Por construcción, el $X_i$ tienen los marginales Beta deseados y su matriz de correlación está determinada por la $d(d-1)/2$ entradas independientes en $\Sigma.$
Aquí, para ilustrar, hay un R aplicación de una función para generar $n$ iid valores Beta multivariados según esta receta. La persona que llama especifica los parámetros Beta $(a)$ y $(b)$ como vectores junto con la matriz de correlación $\Sigma,$ seguido de cualquier parámetro opcional que se pase al generador gaussiano multivariante MASS::mvrnorm . La salida es un $n\times d$ cuyas filas son las realizaciones.
rmbeta <- function(n, a, b, Sigma, ...) {
require(MASS)
Z <- mvrnorm(n, rep(0, nrow(Sigma)), Sigma, ...)
t(qbeta(t(pnorm(Z, 0, sqrt(diag(Sigma)))), a, b))
}
En este código, qbeta implementa $F^{-1}_{a,b}$ y pnorm implementa $\Phi.$ rmbeta debería ser igualmente sencillo de escribir en casi cualquier plataforma de cálculo estadístico.
Una prueba de este algoritmo establecerá que las columnas de la salida (1) están correlacionadas y (2) tienen las distribuciones Beta previstas. En el siguiente ejemplo con $d=4,$ en el que se generaron un millón de valores multivariantes (que tardaron aproximadamente un segundo), las correlaciones de entrada fueron $\rho_{i,j} = (-0.8)^{|i-j|}.$ Esto especifica fuertes correlaciones negativas entre columnas vecinas, correlaciones positivas bastante fuertes para $(Z_1,Z_3)$ y $(Z_2,Z_4),$ y una débil correlación negativa para $(Z_1,Z_4).$ Estos patrones de correlación persistirán, al menos cualitativamente, cuando el $Z_i$ se transforman en el $X_i.$ Estas correlaciones son evidentes en esta matriz de dispersión de los dos primeros miles de $X_i$ (trazar millones de puntos sería superfluo y llevaría demasiado tiempo):
![Figure 1]()
Las correlaciones son evidentes en las nubes de puntos. Las correlaciones observadas entre estos valores generados se obtienen a partir del cor y se obtiene la matriz
X.1 X.2 X.3 X.4
X.1 1.000 -0.794 0.610 -0.509
X.2 -0.794 1.000 -0.739 0.628
X.3 0.610 -0.739 1.000 -0.778
X.4 -0.509 0.628 -0.778 1.000
Se acercan notablemente a las correlaciones de $-0.800,$ $0.640,$ y $-0.512$ especificada para la distribución gaussiana matriz del $Z_i.$ Sin embargo, hay que tener en cuenta que (a diferencia de los gaussianos multivariantes) estas nubes tienden a seguir tendencias curvilíneas: esto se debe a las diferentes formas de las distribuciones marginales.
Aquí están las distribuciones marginales. (Los títulos del histograma muestran los parámetros Beta).
![Figure 2]()
Las curvas rojas son las funciones de densidad previstas: los datos se ajustan bien a ellas.
No he analizado si todas esas distribuciones Beta multivariantes son unimodales. No deberíamos esperar que lo sean a menos que todos los pares de $(a_i,b_i)$ Los parámetros incluyen al menos un valor de $1$ o mayor. (En estos casos, creo que se mantendrá la unimodalidad, pero no ofrezco ninguna prueba de ello. (Los ejemplos ofrecidos en https://stats.stackexchange.com/a/91944/919. )
Finalmente, poco de esta construcción es específica de las distribuciones Beta: puede crear distribuciones multivariadas correlacionadas con cualquier distribución marginal prevista mediante esta construcción de cópula gaussiana. Basta con sustituir las funciones cuantílicas $F_{a,b}^{-1}$ por cualquier función cuantílica (y cualquier parámetro) que desee.
Anexo
Aquí está el R código necesario para reproducir los datos de las figuras. Genera los parámetros Beta a y b al azar. Funciones pairs y hist se utilizaron posteriormente para las parcelas.
set.seed(17)
d <- 4
rho <- -0.8
a <- 1 + rexp(d, 1/5)
b <- 1 + rexp(d, 1/5)
Sigma <- outer(1:d, 1:d, function(i,j) rho^abs(i-j))
X <- rmbeta(1e6, a, b, Sigma)
