Modèle
polynomial
Données synthétiques pour l'estimation
non paramétrique avec la validation croisée.
Description
On observe des paires (x,y) suivant le modèle
y=polynôme(x)+bruit
et l'on cherche à estimer le polynôme, ainsi
que son ordre.
Les données sont ici.
L'objectif est le suivant
1/ Estimer le polynôme pour les ordres allant (par
exemple) de 1 a 9
et tracer sur un même graphique
pour tous les ordres
la variance estimée non-biaisée,
la variance MV et CV.
2/ Localiser les trois minimas. Proposer un ordre raisonnable pour le polynôme.
Donner ses coefficients et l'erreur standard de chacun, ainsi que la variance du bruit.
3/ Superposer sur un même graphique les points
(x,y) et la représentation du polynôme.
4/ Faire le graphique de la question précédente
pour un ordre de 20.
On estimera les 10 modèles dans une boucle et on récupérera
à chaque étape des trois sigmas.
Un exemple préliminaire
Les données pour cette partie sont la population des États-Unis tous les 10 ans de 1900 à 2010.
La paire (x,y) est ici (année,population).
On se propose de faire simplement le tracé qui
superpose les points et le polynôme d'ordre ord prédit.
Noter l'utilisation de poly dans la formule du modèle
(R fabrique les variables polynomiales), et l'utilisation de data.frame dans
predict. Le nom de la variable de prédiction est an.
dat=read.table('popus.txt',h=T)
ord=2
mod=lm(pop~poly(an,ord,raw=T),data=dat)
ann=seq(1900,2090,1)
popn=predict(mod,data.frame(an=ann))
plot(c(ann,dat$an),c(popn,dat$pop),type='n')
lines(ann,popn,col='blue')
points(dat$an,dat$pop,col='red')
Faire augmenter ord j'usqu'à 5 et commenter le résultat.
Analyse
On reprend les données du début.
(A) Etude à l'ordre 3 seulement.
Il s'agit d'estimer le modèle d'ordre trois et calculer les trois
écart-types sig, sigcv et sigmv (associés à
la variance estimée non-biaisée, CV et MV).
On superposera ensuite sur un même graphique les points
(x,y) et la représentation du polynôme.
On utilisera les remarques suivantes concernant un modèle
linéaire estimé par la commande
mod=lm(y~......):
1/ Les résidus sont en attributs de
mod
2/ Le sigma non biaisé est en attribut de
summary(mod), mais on pourra trouver
plus simple de le calculer en même temps que le reste.
3/ Les hi sont en attributs de
lm.influence(mod):h=lm.influence(mod)$hat.
(B) Pour maintenant traiter tous les ordres, il faut faire une boucle.
S'inspirer
du schéma de programme suivant:
deg=9;
##### initialisations des tableaux de résultat
sig=(1:deg);
sigcv=sig;
sigmv=sig;
n=...; à compléter
##### boucle
for (i in 1:deg)
{p=i+1;
mod=lm(...,data=pol);
sig[i]=...; à compléter
sigcv[i]=...; à compléter
sigmv[i]=...; à compléter
}