Message d'erreur sur Matlab

Merci Makitoch,

J'ai un fichier tex qui contient les données suivants:

Y X1 X2
9.95000000000000 2 50
24.45000000000000 8 110
31.75000000000000 11 120
35.00000000000000 10 550
25.02000000000000 8 295
16.86000000000000 4 200
14.38000000000000 2 375
9.60000000000000 2 52
24.35000000000000 9 100
27.50000000000000 8 300
17.08000000000000 4 412
37.00000000000000 11 400
41.95000000000000 12 500
11.66000000000000 2 360
21.65000000000000 4 205
17.89000000000000 4 400
69.00000000000000 20 600
10.30000000000000 1 585
34.93000000000000 10 540
46.59000000000000 15 250
44.88000000000000 15 290
54.12000000000000 16 510
56.63000000000000 17 590
22.13000000000000 6 100
21.15000000000000 5 400

Y= réponse et X=facteurs

QUESTION
-----------------

1. Définir les matrices X et Y dans Matlab
2. calculer les estimés des paramètres par régression multilinéaire

le script suivant nous a été donné pour y arriver:

% Exemple 12_1 Montgomery


load ex12_1Montgomery.txt % ouvrir le fichier de données

X=[ones(size(ex12_1Montgomery,1),1) ex12_1Montgomery(:,2:3)] % matrice X
Y=ex12_1Montgomery(:,1); % matrice Y

Beta=inv(X'*X)*X'*Y; % calcul des paramètres de régression (estimés)

% Paramètres de régression normalisés
B0star = mean(Y);
B1star = Beta(2,1)*std(X(:,2));
B2star = Beta(3,1)*std(X(:,3));

% ANOVA - MODEL SIGNIFICATIF? ************
[n p]=size(X); % n=nombre de points expérimentaux, p=nombre de paramètres de régression
Yhat=X*Beta; % calcul des prédictions du modèle
SSE=(Y-Yhat)'*(Y-Yhat); % Somme des carrés des résidus
MSE=SSE/(n-p); % estimé de la variance des résidus du modèle

Ybar=mean(Y)*ones(n,1); % calcul de la moyenne de Y
SSR=(Yhat-Ybar)'*(Yhat-Ybar); % Somme des carrés expliquée par le modèle
MSR=SSR/(p-1); % estimé de la variance de Y expliquée par le modèle

SST=(Y-Ybar)'*(Y-Ybar); % Somme des carrés totale
MST=SST/(n-1); % estimé de la variance totale de Y

f=MSR/MSE; % statistique f
pvalue=1-fcdf(f,p-1,n-p); % calcul du p-value à partir de la loi F (conclusion, modèle significatif)

R2=SSR/SST; % calcul de R2
R2adj=1-((SSE/(n-p))/(SST/(n-1))); %calcul de R2 ajusté
plot(Yhat,Y,'.',[0 70],[0 70],'-k'),xlabel('Prédictions'),ylabel('Observations') % graphique des observations vs prédictions

% INTERVALLES DE CONFIANCE INDIVIDUELS SUR LES ESTIMÉS ****
C=inv(X'*X); % Calcul de la matrice C
[Beta-tinv(0.975,n-p)*sqrt(diag(C)*MSE) Beta Beta+tinv(0.975,n-p)*sqrt(diag(C)*MSE)] % intervalles de confiance à 95%

% FONCTION MATLAB DE RÉGRESSION MULTILINÉAIRE
[B,BINT,R,RINT,STATS] = regress(Y,X,0.05); % régression multilinéaire
STATS % [R2 MSR/MSE P-VALUE MSE]
[BINT(:,1) B BINT(:,2)] % intervalles de confiance à 95%

% NORMALISATION ********************************************
X2=zscore(X(:,2:3)); % normalisation des colonnes 2 à 3 de X
X3=[ones(n,1) X2]; % ajout d'une colonne de 1 à X2
[B3,BINT3,R3,RINT3,STATS3] = regress(Y,X3,0.05); % régression multilinéaire Y vs X3
e3=BINT3(:,2)-BINT3(:,1); % calcul de la longueur des intervalles
e3=e3./2; % demi-intervalles
errorbar([1:3]',B3,e3) % barres d'erreurs sur les estimés des paramètres

%tinv(P,v)=inverse de la distribution T cumulative - calcul de t_alpha
%finv(P,v1,v2)=inverse de la distribution F cumulative - calcul de f_alpha
%fcdf(X,v1,v2)=distribution F cumulative - calcul du p-value

% ANOVA (MSR/MSE=STATS(1,2))
falpha=finv(0.95,3-1,25-3)
p_value=1-fcdf(STATS(1,2),3-1,25-3) % STATS(1,3)

% Analyse des résidus
subplot(3,2,1),rcoplot(R,RINT)
subplot(3,2,2),plot(X*B,R,'.b',[0 70],[0 0],'-k'),xlabel('yhat'),ylabel('Résidus')
subplot(3,2,3),plot(X(:,2),R,'.b',[0 20],[0 0],'-k'),xlabel('x_1'),ylabel('Résidus')
subplot(3,2,4),plot(X(:,3),R,'.b',[0 600],[0 0],'-k'),xlabel('x_2'),ylabel('Résidus')
subplot(3,2,5),normplot(R)

% ajout de x1^2 au modèle car LOF (voir analyse des résidus)
Xsq=[X X(:,2).^2];
[Bsq,BINTsq,Rsq,RINTsq,STATSsq] = regress(Y,Xsq,0.05);
[BINTsq(:,1) Bsq BINTsq(:,2)]
subplot(3,2,1),rcoplot(R,RINT)
subplot(3,2,2),plot(Xsq*Bsq,Rsq,'.b',[0 70],[0 0],'-k'),xlabel('yhat'),ylabel('Résidus')
subplot(3,2,3),plot(Xsq(:,2),Rsq,'.b',[0 20],[0 0],'-k'),xlabel('x_1'),ylabel('Résidus')
subplot(3,2,4),plot(Xsq(:,3),Rsq,'.b',[0 600],[0 0],'-k'),xlabel('x_2'),ylabel('Résidus')
subplot(3,2,5),plot(Xsq(:,4),Rsq,'.b',[0 400],[0 0],'-k'),xlabel('x_2^2'),ylabel('Résidus')
subplot(3,2,6),normplot(Rsq)

% Prédiction de valeurs futures
xk=[1 19 230 19^2]; % veut prédire y pour x1=19 et x2=230
yhatk=xk*Bsq % prédiction de la réponse
INTyhatk=tinv(0.975,n-p-1)*sqrt(xk*inv(Xsq'*Xsq)*xk'*MSE); %intervalle de confiance à 95% sur yk
[yhatk-INTyhatk yhatk yhatk+INTyhatk]