Discussion :

[Invité]

Bonjour,

Pour m'exercer, j'essaye de générer un échantillon gaussien 2D de cette forme :



Ce que j'ai à disposition, c'est une fonction prédéfinie randn() qui me donne un nombre aléatoire selon la loi gaussienne d'espérance zéro et de variance un.

Pour l'instant, j'arrive à créer un échantillon gaussien "vertical" ou "horizontal" de cette façon :

Code pseudocode :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
coordonnée_x = randn() * écart_type_x + moyenne_x
coordonnée_y = randn() * écart_type_y + moyenne_y

Si "écart_type_x" est plus grand que "écart_type_y" alors l'échantillon a une forme horizontale :



Mais je ne sais pas comment introduire une covariance différente de zéro entre l'axe des x et l'axe des y.

Merci d'avance pour vos indications.

[plegat]

Salut

Mais je ne sais pas comment introduire une covariance différente de zéro entre l'axe des x et l'axe des y.

sans rentrer dans les détails de probabilités (étant donné que je suis une bille dans ce domaine...), ton besoin m'inspire une petite rotation...

Multiplies tes coordonnées par une matrice de rotation, ça devrait te donner le nuage que tu cherches.

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
 
coordonnée_x_final = coordonnée_x*cos(theta)-coordonnée_y*sin(theta)
coordonnée_y_final =coordonnée_x*cos(theta)+coordonnée_y*sin(theta)

avec theta entre -pi/2 et pi/2.

Je te laisse traduire ça en terme de covariance...

A adapter si les moyennes ne sont pas à zéro (en gros sortir le nuage avec 0 de moyenne, faire la rotation, et translater vers la moyenne).

[Invité]

Merci pour la solution géométrique, ça pourra me servir en dépannage, mais il faudrait que j'arrive à utiliser la solution "naturelle" avec la matrice de covariance. Parce qu'elle généralise bien quand on augmente le nombre de dimensions.

Et j'ai l'impression que ce genre de choses (génération d'échantillon) est la principale utilité d'une matrice de covariance. Donc j'aimerai maitriser le truc.

Apparemment, la formule de la densité de probabilité d'une gaussienne utilise la matrice de covariance (le sigma majuscule) :



Cependant, je ne vois pas comment cette densité de probabilité peut me servir pour créer mon échantillon gaussien oblique, car elle ne calcule que des probabilités entre zéro et un ...

[pseudocode]

Cependant, je ne vois pas comment cette densité de probabilité peut me servir pour créer mon échantillon gaussien oblique, car elle ne calcule que des probabilités entre zéro et un ...

Pour avoir une distribution "oblique", il faut que la distribution horizontale dépende de la distribution verticale. C'est à dire que les variables x et y soient corrélées.

La densité de proba s'écrit donc comme une fonction à deux variables.



avec,



Google : Bivariate Normal Distribution

[Invité]

Je suis pas sûr de comprendre : comment je fais pour instancier un point avec cette formule ?

D'abord je choisis le centre (mu) et la matrice de covariance (sigma) que je veux, je fixe la corrélation (rho) à 1 par exemple, mais après comment je fait pour trouver la valeur de x et y ? (pour instancier un point)

[pseudocode]

Je suis pas sûr de comprendre : comment je fais pour instancier un point avec cette formule ?

D'abord je choisis le centre (mu) et la matrice de covariance (sigma) que je veux, je fixe la corrélation (rho) à 1 par exemple, mais après comment je fait pour trouver la valeur de x et y ? (pour instancier un point)

rho fixé à "1", ce n'est pas une très bonne idée. Déjà parce que ca signifie que les 2 variables sont en fait totalement corrélées (=c'est 1 seule variable), et puis la division par zéro dans la formule ce n'est pas conseillé.

Pour instancier un point dans une loi normale, la méthode la plus usuelle c'est de passer par une transformation linéaire.

Sinon, il y a la méthode brutale du "sample rejection". Tu tires un candidat (x,y) suivant une loi uniforme. Tu calcules la pdf de x,y suivant la formule. Tu tires une valeur aléatoire uniforme u en 0 et 1. Si (u > pdf(x,y)) alors tu rejètes le candidats et tu recommences.

[Invité]

J'y avais pensé à la méthode de rejeter un point avec une certaine probabilité. Mais je trouve pas ça très bien :
- déjà il faut tirer un nombre entre plus l'infini et moins l'infini, ce que je ne peux pas faire avec des flottants (il faudrait faire une approximation en prennant un certain intervalle autour de la moyenne)
- ça ne me fait pas progresser

Mais cette méthode à néanmoins un avantage : comme elle est simple, on comprend bien ce qu'on fait.

Je vais utiliser la méthode avec la factorisation de Cholesky. Si ça intéresse quelqu'un, voici le code octave de cette méthode (ce n'est pas de moi : l'original est ici) :

Code Octave :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
function echantillon = echantillonGaussien(nb_observations, centre, matrice_covariance)
  nb_dimensions = numel(centre);
  echantillon = randn(nb_observations, nb_dimensions);
  if nb_observations > 1
    echantillon = echantillon - repmat(mean(echantillon), nb_observations, 1);
  end
 
  [factorisation_cholesky, definie_positive] = chol(matrice_covariance);
  if definie_positive == 0
    echantillon = echantillon * factorisation_cholesky;
  else
    echantillon = echantillon * sqrtm(matrice_covariance);
  end
 
  echantillon = echantillon + repmat(centre, nb_observations, 1);
end

Merci bien de m'avoir aidé