2012-05-01 2 views
0

Salut je ne suis pas sûr si mon algorithme est correct, je suis en train de reproduire la fonction mvnrnd de matlab mais en gsl. J'ai trouvé un algorithme dans quelques articles de journal qui produit un vecteur de normale multivariée, mais j'ai besoin d'une matrice de nombres aléatoires normaux multivariésmatlab mvnrnd en gsl

disons que la distribution est Z ~ (mu, sigma);

en supposant que sigma est une matrice déjà définie positive.

un algorithme je l'ai trouvé sur le web dit

1. cholskey(sigma) = A 
2. generate uniform gaussian vector r 
3. matrix vector triangular product with gsl_blas_dtrmv A * r 
4. add mu to Ar and that will be a vector of multivariate normal random numbers 

ma méthode ci-dessous

sont les modifications suivantes belowcorrect au produit une matrice de variables aléatoires

1. cholskey(sigma) = A 
    2. generate uniform gaussian matrix R 
    3. matrix matrix scalar product AR 
    4. add mu to AR and that will be a matrix of multivariate normal random numbers 

Répondre

2

Oui, ça est correct. Voir par exemple this Wikipedia entry on multivariate normal RNGs qui présente cette section:

valeurs de dessin à partir de la distribution

Un procédé largement utilisé pour l'élaboration d'un vecteur x aléatoire à partir de la N dimensions distribution normale multivariée avec le vecteur moyen μ et matrice de covariance fonctionne Σ comme suit:

  1. Trouver une matrice réelle A tel que A AT = Σ. Lorsque Σ est positif-défini, la décomposition de Cholesky est typiquement utilisée. [...]

  2. Soit z = (z1, ..., zn) T un vecteur dont les composantes sont des variables aléatoires normales N standards indépendants (qui peuvent être générés, par exemple, en utilisant la transformée de Box-Muller).

  3. Soit x soit μ + Az. Ceci a la distribution désirée due à la propriété de transformation affine.

qui décrit le même algorithme.

Il existe également plusieurs implémentations pour R, par exemple mvrnorm dans le package MASS fourni avec chaque installation R.

+0

merci! Je ne pensais pas que j'aurais quelqu'un célèbre comme toi pour répondre à ma question – pyCthon