Approximation polynomiale

Version imprimable

Bonjour à tous,

Je cherche à "fitter" les points d'un vecteur y avec un polynôme dont la puissance la plus faible est d'ordre 1, mathématiquement on écrit :
y = sigma(de k=1 à k=N) pk * x^k, pk étant les coefficients du polynôme.

J'ai essayé de coder ça en utilisant la fonction polyfit :

Code:

1
2
3
4
5
6
7
8
9
 
x=[0.00447213600000000,0.0180949800000000,0.0317178200000000,0.0453406500000000,0.0589634900000000,0.0725863300000000,0.0862091700000000,0.0998320100000000,0.113454900000000,0.127077700000000,0.140700500000000,0.154323400000000,0.167946200000000,0.181569100000000,0.195191900000000,0.208814700000000,0.222437600000000,0.236060400000000,0.249683200000000,0.263306100000000,0.276928900000000,0.290551800000000,0.304174600000000,0.317797400000000,0.331420300000000,0.345043100000000,0.358666000000000,0.372288800000000,0.385911600000000,0.399534500000000,0.413157300000000,0.426780200000000,0.440403000000000,0.454025800000000,0.467648700000000,0.481271500000000,0.494894400000000,0.508517200000000,0.522140000000000,0.535762900000000,0.549385700000000,0.563008600000000];
 
y = [0.279589606042921 + 0.372724610293237i,0.150469709434042 + 0.0158359760483071i,0.259755509082023 - 0.0729330482967179i,0.327093967340840 - 0.148265993630941i,0.382124430436839 - 0.227073704701423i,0.492016975825251 - 0.312928300760763i,0.641994080914315 - 0.417584163020943i,0.666709855436498 - 0.504390295652359i,0.664065756966530 - 0.587147827075405i,0.657417499164075 - 0.663826644416700i,0.596039702508291 - 0.712067994510609i,0.528534501906463 - 0.748494576693688i,0.450680529607647 - 0.765959068203925i,0.373077606639998 - 0.770424645400167i,0.299200965989157 - 0.767932948676417i,0.230394973592732 - 0.758918218892205i,0.169072274821718 - 0.733312669572486i,0.118153499891171 - 0.695587042595892i,0.0799067917615191 - 0.653113203185179i,0.0540346383423679 - 0.609980335377763i,0.0338318463496849 - 0.567688187282820i,0.0145381142932576 - 0.526284674220806i,-0.00185810788924279 - 0.486391889668856i,-0.0137909136559662 - 0.449709133366984i,-0.0222113897878540 - 0.416854042552800i,-0.0303962043346147 - 0.387534938934198i,-0.0400026770482179 - 0.361489685635558i,-0.0499921255229547 - 0.338669564981779i,-0.0585566138010606 - 0.318655521163127i,-0.0652913414660232 - 0.301295337059639i,-0.0701167327554450 - 0.285978627287176i,-0.0733278339141154 - 0.272512230318287i,-0.0755008063941609 - 0.260718351576263i,-0.0769712615202275 - 0.250007683569512i,-0.0776370578248482 - 0.240528727619081i,-0.0781741029017742 - 0.232268094346295i,-0.0786774950938140 - 0.225094139128329i,-0.0785018846107360 - 0.218719490505836i,-0.0771874116299962 - 0.212511197914517i,-0.0747255611087394 - 0.206113121181531i,-0.0719820602354802 - 0.199918468936144i,-0.0696807510115635 - 0.194252761102507i ];
 
p= polyfit(x,y,16) ;
p(1,end)=0 ;   % j'annule le coefficient correspondant à l'ordre 0
f=polyval(p,x) ;
plot(x,abs(y),'o',x,abs(f),'-')

Malheureusement, cette solution n'aboutit pas à un bon fit, parce que j'ai du annuler le dernier coef d'ordre 0 puisque la fonction polyfit calcule les coefficients d'un polynôme de degré allant de 0 à N, or moi je veux qu'elle interpole mes points avec un polynôme de degré allant de 1 à N.

Auriez vous une autre idée pour résoudre ce problème ?

En vous remerciant d'avance.

10/12/2013, 19h13
shaiHulud

Bonjour,

Dans ton cas, tu peux fitter y./x = P(x). (ce n'est pas exactement la mm chose)
Sinon, il faut ajouter des contraintes dans ton optimisation.
11/12/2013, 09h26
le fab

plus efficace que polyfit, utilise l’opérateur \
exemple

là tu fais ce que tu veux (genre interdire les dégrées 0 et 12)

Fab
11/12/2013, 11h15
dptmt

Citation:

Envoyé par VV33D

Bonjour,

Dans ton cas, tu peux fitter y./x = P(x). (ce n'est pas exactement la mm chose)
Sinon, il faut ajouter des contraintes dans ton optimisation.

mais pour la suite j'ai besoin des coefficients du polynôme qui fit bien y et non pas y./x .

Merci en tout cas.
11/12/2013, 11h21
dptmt

Citation:

Envoyé par le fab

plus efficace que polyfit, utilise l’opérateur \
exemple

là tu fais ce que tu veux (genre interdire les dégrées 0 et 12)

Fab

Est ce que tu peux me dire à quoi correspondent y et x dans l'exemple 8 du lien ?

J'ai fait un fit polynomial en utilisant l'opérateur mldivide \ .
voici ce que ça donne au niveau du code

Code:

1
2
3
4
5
6
7
8
9
10
 
x=[0.00447213600000000,0.0180949800000000,0.0317178200000000,0.0453406500000000,0.0589634900000000,0.0725863300000000,0.0862091700000000,0.0998320100000000,0.113454900000000,0.127077700000000,0.140700500000000,0.154323400000000,0.167946200000000,0.181569100000000,0.195191900000000,0.208814700000000,0.222437600000000,0.236060400000000,0.249683200000000,0.263306100000000,0.276928900000000,0.290551800000000,0.304174600000000,0.317797400000000,0.331420300000000,0.345043100000000,0.358666000000000,0.372288800000000,0.385911600000000,0.399534500000000,0.413157300000000,0.426780200000000,0.440403000000000,0.454025800000000,0.467648700000000,0.481271500000000,0.494894400000000,0.508517200000000,0.522140000000000,0.535762900000000,0.549385700000000,0.563008600000000];
 
y = [0.279589606042921 + 0.372724610293237i,0.150469709434042 + 0.0158359760483071i,0.259755509082023 - 0.0729330482967179i,0.327093967340840 - 0.148265993630941i,0.382124430436839 - 0.227073704701423i,0.492016975825251 - 0.312928300760763i,0.641994080914315 - 0.417584163020943i,0.666709855436498 - 0.504390295652359i,0.664065756966530 - 0.587147827075405i,0.657417499164075 - 0.663826644416700i,0.596039702508291 - 0.712067994510609i,0.528534501906463 - 0.748494576693688i,0.450680529607647 - 0.765959068203925i,0.373077606639998 - 0.770424645400167i,0.299200965989157 - 0.767932948676417i,0.230394973592732 - 0.758918218892205i,0.169072274821718 - 0.733312669572486i,0.118153499891171 - 0.695587042595892i,0.0799067917615191 - 0.653113203185179i,0.0540346383423679 - 0.609980335377763i,0.0338318463496849 - 0.567688187282820i,0.0145381142932576 - 0.526284674220806i,-0.00185810788924279 - 0.486391889668856i,-0.0137909136559662 - 0.449709133366984i,-0.0222113897878540 - 0.416854042552800i,-0.0303962043346147 - 0.387534938934198i,-0.0400026770482179 - 0.361489685635558i,-0.0499921255229547 - 0.338669564981779i,-0.0585566138010606 - 0.318655521163127i,-0.0652913414660232 - 0.301295337059639i,-0.0701167327554450 - 0.285978627287176i,-0.0733278339141154 - 0.272512230318287i,-0.0755008063941609 - 0.260718351576263i,-0.0769712615202275 - 0.250007683569512i,-0.0776370578248482 - 0.240528727619081i,-0.0781741029017742 - 0.232268094346295i,-0.0786774950938140 - 0.225094139128329i,-0.0785018846107360 - 0.218719490505836i,-0.0771874116299962 - 0.212511197914517i,-0.0747255611087394 - 0.206113121181531i,-0.0719820602354802 - 0.199918468936144i,-0.0696807510115635 - 0.194252761102507i ];
 
X=[x, x.^2, x.^3 , x.^4 , x.^5 , x.^6];
p=zeros(size(X,2),1);
p = X \ y ;
f = X*p;
plot(x,abs(y),'o',x,abs(f),'-')

ça marche bien avec ce code mais je veux l'optimiser afin qu'il me calcule automatiquement le degré du polynôme (ici j'ai mis 6) assurant le meilleur fit possible de y ! Comment peut on s'y prendre ?

11/12/2013, 13h35
shaiHulud

Bonjour,

Aucune recette miracle pour optimiser m(n) sur un paramètre discret n. En général deux approches:
- Prolonge n par un paramètre régulier x + optimisation classique sur x, puis on regarde les deux entiers proches de l'optimum x*
- Exploration systématique.

Dans ton cas, n correspond à la complexité du modèle, donc plus n est grand, meilleur est le fit sur un échantillon donnée. Il faut donc soit
- Pénaliser la qualité du fit (l'erreur d'interpolation) par une pénalisation dépendant de n, ton nombre de paramètre.
- Faire une validation croisée quelconque (Nfold, jacknife, bootstrap, fenêtres glissantes ...) : tu estimes P(x) sur un sous échantillon, tu calcules la qualité du fit sur un autre sous échantillon (on appelle ca out-of-sample). Tu sélectionnes le N correspondant au meilleur fit out of sample
11/12/2013, 14h03
le fab

la réponse est 42 !

en effet tes données sont composées de 42 valeurs, donc un polynome d'ordre 42 devrait passer exactement par tous les points
désolé j'ai pas pu m’empêcher :mrgreen: