% TP5ex2 Régression linéaire gaussienne simple

clear
disp ('Régression linéaire gaussienne simple');
n=input('Valeur de n : ');
a=input('Valeur de a : ');
b=input('Valeur de b : ');
sigcar=input('Valeur de la variance : ');

x=randn(n,1); % Choix d'une loi ...
epsilon=sqrt(sigcar).*randn(n,1);
Y=a+b.*x+epsilon;
X=[ones(n,1) x];

tetac=inv(X'*X)*(X')*Y;
ac=tetac(1,1);
bc=tetac(2,1);

Ychap=X*tetac;
sig2c=sum((Y-Ychap).^2)/(n-2); % sigma carré chapeau

% intervalles de confiance
deg=n-2;
t=qt(0.975,deg);
mx=mean(x);
mx2=mean(x.*x);
vx=(n-1)/n*var(x);
va=mx2/(n*vx);
vb=1/(n*vx);
acm=ac-t*sqrt(sig2c*va);
acp=ac+t*sqrt(sig2c*va);
Ia=[acm acp];
bcm=bc-t*sqrt(sig2c*vb);
bcp=bc+t*sqrt(sig2c*vb);
Ib=[bcm bcp];

disp ('Intervalle de confiance de ac');
Ia
disp ('Intervalle de confiance de bc');
Ib

% intervalle de confiance de sig2c

if (n>200) % Utilisation du théorème limite central pour n grand
   k1=qnorm(0.025)*sqrt(2*n)+n;
   k2=qnorm(0.975)*sqrt(2*n)+n;
else
   k1=qchisq(0.023,deg);
   k2=qchisq(0.975,deg);
end

disp('estimateur de la variance');
sig2c
disp('Intervalle de confiance de sig2c');
Isig2c=[deg*sig2c/k2 deg*sig2c/k1]

% Tracé
clf;
x=sort(x);
hold on;

plot(x,Y,'*');


vxx=(mx2-2*mx*x+x.*x)/(n*vx);
plot(x,ac+bc*x+t*sqrt(sig2c*vxx),'g') %int. de conf de la droite de régréssion linéaire
plot(x,ac+bc*x-t*sqrt(sig2c*vxx),'g') %int. de conf de la droite de régréssion linéaire
plot(x,ac+bc*x,'r');
title('Régression linéaire gaussienne simple');
xlabel('Observations et droite des moindres carrés');
hold off;

%n=20 a=1 b=10 var=100 loi norm : joli graphique

%Ellipsoiide de confiance
ff=qf(.95,2,n-2);
fs=2*sigc*ff;
tea=[0:0.01:2*pi];
deno=n*(cos(tea).*cos(tea)+2*mx*cos(tea).*sin(tea)+mx2*sin(tea).*sin(tea));
r=sqrt(fs./deno);
plot(ac-r.cos(tea),bc-r.sin(tea))
hold on
plot(a,b,'r*')
xlabel('Ellipsoide de confiance')
grid
hold off