D = [1958    294    634
1959 314 728
1961 383 819
1963 402 938
1965 475 1136
1967 786 1317];


%%

x = D(:,2);
y = D(:,3);

plot(x,y,'or');

X = [x ones(size(x))];


n = length(x);
sx = sum(x);
sy = sum(y);
sx2 = x'*x;
sy2 = y'*y;
sxy = x'*y;

a = (n*sxy - sx*sy)/(n*sx2 - sx*sx)
a = sum((x-mean(x)).*(y-mean(y)))/sum((x-mean(x)).^2)
a = ((x-mean(x))'*(y-mean(y)))/sum((x-mean(x)).^2)
cv = cov(x,y);
a = cv(1,2)/var(x)
b = mean(y) - a*mean(x)

A = X\y
A = (X'*X)\(X'*y)
A = inv(X'*X)*X'*y
[Q,R] = qr(X);
A = R\Q'*y

polyfit(x,y,1)

hold on
plot(x,a*x+b,'b');

%% question 2) Qualité

R = cv(1,2)/sqrt(var(x)*var(y))
R = corrcoef(x,y)
R2 = R^2
T = [x y a*x+b y-a*x-b (y-a*x-b).^2 ]
R2 = 1 - sum(T(:,5))/(n*var(y,1))

%% question 3) prediction

yp = a*800 + b
t  = icdf('t',0.975,4);
sig = sqrt(sum(T(:,5))/(n-2))
interval = t*sig*sqrt(1 + 1/n + (800 - mean(x))^2/(n*var(x,1)))

%%
Xr = X;
Xr(6,:) = [];
A = Xr\y(1:5);
plot(x,X*A,'m')

an = A(1)
bn = A(2)
T = [x(1:5) y(1:5) an*x(1:5)+bn y(1:5)-an*x(1:5)-bn (y(1:5)-an*x(1:5)-bn).^2 ]
R2 = 1 - sum(T(:,5))/(5*var(y(1:5),1))
Last modified: lundi, 15 juillet 2013, 12:52