# 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 ellipse(s) ou rectangle(s) de confiance pour la moyenne...
#---------------------------

n=10; # nombre d'observations
a=5; # parametre qui determine l'etendue du graphe
beta=seq(.95,.95,.05); # 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=t(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;

muhat=c(mean(X[,1]),mean(X[,2]));
Sigmahat=var(X);
Ahat=t(chol(Sigmahat));

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", col="blue",asp=1);
points(mu[1],mu[2],col='green');
abline(v=mu[1],lwd=.5,col='green');abline(h=mu[2],lwd=.5,col='green');

# Ellipses de confiance exactes (en noir)
rayonbeta=sqrt((k*(n-1))/(n*(n-k))*qf(beta,k,n-k));
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(muhat,c(k,length(t))));
ellipse=t(Ahat%*%t(rayonbeta[i]*Sun))+MuSun;
abline(v=0,lwd=.5);abline(h=0,lwd=.5);
polygon(ellipse);
}

# Ellipses de confiance asymptotiques (en rouge)
rayonbeta=sqrt((1/n)*qchisq(beta,k));
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(muhat,c(k,length(t))));
ellipse=t(Ahat%*%t(rayonbeta[i]*Sun))+MuSun;
abline(v=0,lwd=.5);abline(h=0,lwd=.5);
polygon(ellipse,border='red');
}

# Rectangles de confiance exacts (en noir)
rayonbeta1=sqrt((k*(n-1))/(n*(n-k))*Sigmahat[1,1]*qf(beta,k,n-k));
rayonbeta2=sqrt((k*(n-1))/(n*(n-k))*Sigmahat[2,2]*qf(beta,k,n-k));
for (i in 1:length(beta)) {
rect(muhat[1]-rayonbeta1[i],muhat[2]-rayonbeta2[i],muhat[1]+rayonbeta1[i],muhat[2]+rayonbeta2[i])
}

# Produits des intervalles de confiance exacts individuels (en rouge)
rayonbeta1=sqrt(1/n*Sigmahat[1,1]*qf(beta,1,n-1));
rayonbeta2=sqrt(1/n*Sigmahat[2,2]*qf(beta,1,n-1));
for (i in 1:length(beta)) {
rect(muhat[1]-rayonbeta1[i],muhat[2]-rayonbeta2[i],muhat[1]+rayonbeta1[i],muhat[2]+rayonbeta2[i],border='red')
}