Skew multivarié Normal dans Stan

Question

J'essaie d'estimer un modèle hiérarchique bayésien dans Rstan et je voudrais inclure dans mon modèle une distribution normale asymétrique multivariée. Ce n'est pas une distribution qui est déjà définie dans Stan, mais la documentation semble suggérer que l'on pourrait l'implémenter en utilisant des facteurs de Cholesky. Par exemple, la documentation de Stan 2.15.0 indique aux pp. 333-334:

"Cette reparameterization d'une distribution normale multivariée en termes de variables normales standard peut être étendue à d'autres distributions multivariées qui peuvent être conceptualisées en tant que contaminations du normal multivarié, tel que le Student t multivarié et la distribution normale multidimensionnelle asymétrique. "

Est-ce que quelqu'un a une idée de comment faire cela? J'ai envisagé d'implémenter le skew multivariate normal moi-même dans Stan, mais il ne semble pas qu'il y ait une forme fermée agréable pour la distribution qui se prêterait à une implémentation simple ...

Answer 1

Donc le paquet sn a la fonction rmsn montre clairement la section suivante du code:

function (n = 1, xi = rep(0, length(alpha)), Omega, alpha, tau = 0, dp = NULL) { 
//.... 
    lot <- dp2cpMv(dp = dp0, family = "SN", aux = TRUE) 
    d <- length(dp0$alpha) 
    y <- matrix(rnorm(n * d), n, d) %*% chol(lot$aux$Psi) 
//.... 
} 

Ceci est très similaire à la façon dont les mvtnorm::rmvn fonctionne exactement la matrice chol() méthode vient de la fonction dp2cpMv dans la bibliothèque de base here. Vous pouvez le porter dans votre bloc functions{}. Et dans votre programme stan vous allez (je suppose la matrice de CoV dans chol() est un paramètre)

functions { 
    // Contains code for your ported cholesky factor 
} 
transformed parameters { 
    matrix[K, J] z; 
    cholesky_factor_corr[K] L_tri = dp2cvClone(...); // Cholesky factor from your function{} block.. 
    beta = foo + (L_tri * z)'; // Assuming foo is baseline parameter representing the mean of dimensionality J*K. 
.... 
} 

parameters { 
    matrix[K, J] beta; //# J levels/groups and K dimensional parameters 
    to_vector(z) ~ normal_pdf(0, 1); 
    .... 
} 

Aussi, vous pouvez passer des choses à partir des paramètres transformés bloc aux données ou un bloc de données transformées si les entrées à dp2cvClone() sont Les données. Mais vous avez compris. Cette dernière section de code est extraite du manuel stan section 8.15 et paraphrasée par souci de brièveté, et espère capturer les parties importantes dont vous avez besoin pour le faire fonctionner.