# Ce fichier "R" permet d'engendrer des v.a. i.i.d. bivariees de loi normale, de les afficher, ainsi que de representer sur le meme graphe une ou plusieurs zone(s) de tolerance...
#---------------------------

n=200; # nombre d'observations
a=5; # parametre qui determine l'etendue du graphe
beta=seq(.70,.90,.10); # determine les zones de tolerance considerees (beta=1-alpha)

#---------------------------
# choix de la moyenne de la distribution normale

#mu=c(0,0);
mu=c(4,3);

#---------------------------
# choix de la matrice de variance-covariance de la distribution normale

#Sigma=diag(k);A=chol(Sigma);
#Sigma=.6*array(c(4,0,0,1),c(2,2));A=t(chol(Sigma));
#theta=pi/6;ratio=3;scale=.4;A=scale*array(c(ratio*cos(theta),ratio*sin(theta),-sin(theta),cos(theta)),c(2,2));Sigma=A%*%t(A);
Sigma=.4*array(c(5,3,3,2.25),c(2,2));A=t(chol(Sigma));
#A=array(c(1.5,2,0,1),c(2,2));Sigma=A%*%t(A);

#---------------------------

k=2;
Z=array(rnorm(k*n),c(n,k));
Mu=t(array(mu,c(k,n)));
X=t(A%*%t(Z))+Mu;

chibeta=sqrt(qchisq(beta,k));
d=beta;

plot(X[,1],X[,2],xlim=c(mu[1]-a,mu[1]+a),ylim=c(mu[2]-a,mu[2]+a),xlab="x_1",ylab="x_2",main="Normale bivariee", col="blue",asp=1);
points(mu[1],mu[2],col='green');

for (i in 1:length(beta)) {
t=seq(0,2*pi,.01);
Sun=array(c(cos(t),sin(t)),c(length(t),2));
MuSun=t(array(mu,c(k,length(t))));
ellipse=t(A%*%t(chibeta[i]*Sun))+MuSun;
d[i]=sum(sqrt(Z[,1]^2+Z[,2]^2)<=chibeta[i])/n;
abline(v=0,lwd=.5);abline(h=0,lwd=.5);

# La commande suivante affiche les zones de tolerance considerees
polygon(ellipse);
}

# Les commandes suivantes affichent le nombre d'observations et les parametres de la loi consideree
print("Nombres d'observations:");
print(n);
print("Moyenne:");
print(mu);
print("Variance-covariance:");
print(Sigma);

# Les commandes suivantes affichent les proportions d'observations effectivement obtenues dans chaque zone de tolerance consideree
print("Proportion d'observations dans la zone de tolerance a x%:");
print(array(c(beta,d),c(length(d),2)));

# Les commandes suivantes affichent les points d'intersection du bord de la plus grande zone de tolerance consideree avec ses axes principaux
#w=eigen(Sigma)$vectors;
#lam=eigen(Sigma)$values;
#vp=t(mu+chibeta[length(beta)]*w%*%diag(sqrt(lam)));
#vm=t(mu-chibeta[length(beta)]*w%*%diag(sqrt(lam)));
#points(vp,col='red');
#points(vm,col='red');