1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71
|
% 1/ préparation des données à simuler
% 1.a) X et Y sont tirés dans [0 1[
taille=10000;
X=rand(taille,1) ;
Y=rand(taille,1) ;
% 1.b) on simule Z = X^a * Y^b après avoir ajouté du bruit à X et Y
bruit = 0.1 ;
a=1.75 ;
b=-0.9 ;
Z=(X+bruit*rand(taille,1)).^a .* (Y+bruit*rand(taille,1)).^b ;
% 2/ définition du modèle à résoudre
px=13 ; % degré max en X
py=13 ; % degré max en Y
ptot=13 ; % degré max total
[i,j]=meshgrid(0:px,0:py) ; % génération de toutes les possibilités
u=find(i+j<=ptot) ; % écarter les degrés trop élevés
i=i(u) ; j=j(u) ;
nk=length(u) ; % Nombre de paramètres à calculer ?
% 3/ Résolution du problème
% 3.a) on calcule 4 matrices de dimention [taille * nk]
nx = length(X) ; % combien de données ?
X=repmat(X,1,nk) ;
Y=repmat(Y,1,nk) ;
i=repmat(i',nx,1) ;
j=repmat(j',nx,1) ;
% 3.b) on les assemble dans la matrice M
M=X.^i.*Y.^j ;
% 3.c1) résolution normale par le \ de Matlab
k_matlab = M\Z ;
% 3.c2) pour compraison, résolution par inv(X'X)
k_manuel = inv(M' * M) * M' * Z ;
% 3.d) comparaison des coefficients du polynome par les deux méthodes
% C'est quasiment la meme chose, mais pas tout à fait.
% Dans certains cas dégénérés, l'inversion peut etre fausse et l'anti-divison contnuer à trouver le bon résultat.
disp('Réultats de chaque calcul') ;
disp ([k_matlab k_manuel]) ;
disp('différence relative entre les deux méthodes') ;
disp(k_matlab-k_manuel) ;
% 4. Calcul des Z simulés
Z_matlab = M*k_matlab ;
Z_manuel = M*k_manuel ;
differences = (Z_matlab-Z_manuel) ;
disp ('l''ecart type des differences en Z entre chaque version est') ;
disp (std(differences)) ;
% 5. calcul du coeficient de correlation
residus = Z-Z_matlab ;
r2 = var(Z_matlab)/var(Z) ;
disp(sprintf('coefficient de correlation r2 = %f', r2)) ;
% 6. calcul des z pour des valeurs de X et Y differentes de celles
% utilisées pour le calcul de l'approximation ?
taille2 = 10 ; % valeur forcement differente de "taille" utilisé en 1.a)
xx1 = rand(taille2,1) ;
yy1 = rand(taille2,1) ; |
Partager