Discussion :

[thecrazydonut]

Bien le bonjour,

J'ai un problème... En fait plutôt deux, un purement matlab et l'autre davantage d'algorithmique.

Situation : J'ai une fonction de densité de probabilité gaussienne avec une moyenne et un écart type qui sont tous les deux connus. Je voudrais pouvoir échantillonner cette densité avec un nombre restreint de tirages (nTirages) tout en m'assurant de préserver les statistiques de ma distribution (donc préserver la moyenne et l'écart type). Mon idée initiale est de faire un nombre de tirages suffisant et de calculer les quantiles de cet échantillons (calcul de N quantiles correspondant à mes nTirages tirages restreints). Entre chaque quantile devrait se trouver un nombre équivalent de points aléatoires. On pourrait donc remplacer ces groupes par leur quantiles respectifs.

Code :

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

1
2
3
4
5
6
7
8
nTirage=1e5; % Nombre tirage
nQuantile=1e4; % Nombre de quantiles
x=normrnd(0,1,1,nTirage); % Tirage aléatoire N(0,1)
q=quantile(x,nQuantile); % Calcul de nQuantile (1/(N + 1), 2/(N + 1), ..., N/(N + 1))
[muh,sigh]=normfit(q); % Vérification
fprintf('muh: %0.3f, sigh: %0.3f\n',muh,sigh)
 
>> muh: 0.001, sigh: 0.999 % yeah !

Mes questions:
1) La question purement Matlab. C'est totalement ridicule de devoir tirer un grand nombre de points pour approximer la vraie distribution dans le but de restreindre le nombre de tirages. Est-ce qu'il est possible de calculer ces quantiles directement sur la fonction normale exacte sans passer par le tirage aléatoire (l'échantillonnage) ?
2) La question algorithmique. Le code proposé fonctionne correctement tant que le nombre de quantile est important. Par contre, dès que le nombre de quantile décroît, l'écart type dévie de la valeur attendue (systématique trop faible). Lorsque je définie nQuantile=10 ou 20, alors sigh=0.83... Aïe, on est loin de la vérité.

Merci à vous !

[thecrazydonut]

OK... J'ai pris un peu le problème à l'envers. C'est la fonction norminv qui correspond à mes besoins.

Code :

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

1
2
3
4
5
6
nQuantile=5e1; % Nombre de quantiles échantillonés
q=norminv(linspace(0.01,0.99,nQuantile),0,1);
[muh,sigh]=normfit(q); % Vérification
fprintf('muh: %0.3f, sigh: %0.3f\n',muh,sigh)
 
>>muh: 0.000, sigh: 0.997