// TLC pour des vecteurs uniformes sur un triangle
//
// Charger au préalable la fonction fq2d.sci
//
// 1)Dessine un triangle, génère N échantillons de taille n de la loi uniforme
//   sur ce triangle. 
//
// 2)Affiche le nuage de points des N sommes normalisées des échantillons. 
// 
//
// 3)Affiche les lignes de niveau de la densité gaussienne limite
//
// suggestion d'essai M=4.*rand(2,3), n=500, N=800
M=input('Coordonnées des sommets du triangle (matrice 2x3)?')
n=input('Taille d''échantillon ?'); 
N=input('nombre d''échantillons');
xsom=M(1,:);ysom=M(2,:); 
SX=zeros(1,N);SY=zeros(1,N);
for i=1:N
U=rand(2,n); // n colonnes de 2-échantillons de U[0,1]
so=gsort(U,'r','i');// tri ascendant par colonne pour statistiques d'ordre
esp=[so; ones(1:n)]-[zeros(1:n); so];// n colonnes d'espacements 
X=xsom*esp-sum(xsom)./3; //abscisses de l'échantillon centré
Y=ysom*esp-sum(ysom)./3; //ordonnées de l'échantillon centré
SX(i)=n.^{-0.5}.*sum(X); // somme normalisée de l'échantillon no i
SY(i)=n.^{-0.5}.*sum(Y);
end
//
// Affichages triangle et nuage
//
xbasc();
// On veut que tout tienne dans le cadre: triangle et nuage
a=min([SX xsom]);b=max([SX xsom]);
c=min([SY ysom]);d=max([SY ysom]);
xmin=a-0.2*(b-a); 
xmax=b+0.2*(b-a); 
ymin=c-0.2*(d-c); 
ymax=d+0.2*(d-c); 
plot2d([xmin xmax],[ymin ymax],[-1 -1],"021");// cadre vide
xpoly(xsom,ysom,"lines",1);// tracé du triangle
xpoly(SX,SY,'marks');// tracé du nuage de points

// Calcul de la covariance par P*cov(V)*P'
// où V est un vecteur de loi uniforme sur le triangle
// associe a la base orthonormee canonique
// et P la matrice de transformation de cette base en (AB,AC)
// 
P=[M(:,2)-M(:,1), M(:,3)-M(:,1)];
CV=P*[1/18, -1/36; -1/36, 1/18]*P';
// Lignes de niveau de la densité gaussienne de la loi limite
x=linspace(xmin,xmax,101);y=linspace(ymin,ymax,101);
z=1/(2*%pi*sqrt(det(CV))).*exp(-0.5.*fq2d(x,y,inv(CV)));
contour2d(x,y,z,6,2:7,"000");