ajustement quadratique (estimation des coefficients)

Postez ici vos questions, réponses, commentaires ou suggestions - Les sujets seront ultérieurement répartis dans les archives par les modérateurs

Modérateur : Groupe des modérateurs

Terry Gilles
Messages : 8
Enregistré le : 23 Avr 2008, 12:20

ajustement quadratique (estimation des coefficients)

Messagepar Terry Gilles » 07 Mai 2008, 08:59

Bonjour,

je suis tout nouveau sur ce forum, je ne suis qu'un débutant et j'aimerais beaucoup améliorer mes connaissances dans ce logiciel.
Voilà mon problème:
je dispose d'un tableau de données:

Jan. Feb. Mar. Apr. May Jun. Jul. Aug. Sep. Oct. Nov. Dec.
1984 8,844 8,882 8,920 8,958 8,996 9,033 9,071 9,109 9,147 9,185 9,223 9,260
1985 9,298 9,336 9,374 9,412 9,450 9,488 9,525 9,563 9,601 9,639 9,677 9,715
1986 9,752 9,790 9,828 9,866 9,904 9,942 9,980 10,017 10,055 10,093 10,131 10,169
1987 10,207 10,244 10,282 10,320 10,358 10,396 10,434 10,472 10,509 10,547 10,585 10,623
1988 10,661 10,699 10,736 10,774 10,812 10,850 10,888 10,926 10,964 11,001 11,039 11,077
1989 11,115 11,153 11,191 11,228 11,266 11,304 11,342 11,380 11,418 11,456 11,493 11,531
1990 11,569 11,607 11,645 11,683 11,720 11,758 11,796 11,834 11,872 11,910 11,947 11,985
1991 12,023 12,061 12,099 12,137 12,175 12,212 12,250 12,288 12,326 12,364 12,402 12,439
1992 12,477 12,515 12,553 12,591 12,629 12,667 12,704 12,742 12,780 12,818 12,856 12,894
1993 12,931 12,969 13,007 13,045 13,083 13,121 13,159 13,196 13,234 13,272 13,310 13,348
1994 13,386 13,423 13,461 13,499 13,537 13,575 13,613 13,651 13,688 13,726 13,764 13,802
1995 13,840 13,878 13,915 13,953 13,991 14,029 14,067 14,105 14,143 14,180 14,218 14,256
1996 14,294 14,332 14,370 14,407 14,445 14,483 14,521 14,559 14,597 14,635 14,672 14,710
1997 14,748 14,786 14,824 14,862 14,899 14,937 14,975 15,013 15,051 15,089 15,126 15,164
1998 15,202 15,240 15,278 15,316 15,354 15,391 15,429 15,467 15,505 15,543 15,581 15,618
1999 15,656 15,694 15,732 15,770 15,808 15,846 15,883 15,921 15,959 15,997 16,035 16,073
2000 16,110 16,148 16,186 16,224 16,262 16,300 16,338 16,375 16,413 16,451 16,489 16,527
2001 16,565 16,602 16,640 16,678 16,716 16,754 16,792 16,830 16,867 16,905 16,943 16,981
2002 17,019 17,057 17,094 17,132 17,170 17,208 17,246 17,284 17,322 17,359 17,397 17,435
2003 17,473 17,511 17,549 17,586 17,624 17,662 17,700 17,738 17,776 17,814 17,851 17,889
2004 17,927 17,965 18,003 18,041 18,078 18,116 18,154 18,192 18,230 18,268 18,305 18,343
2005 18,381 18,419 18,457 18,495 18,533 18,570 18,608 18,646 18,684 18,722 18,760 18,797
2006 18,835 18,873 18,911 18,949 18,987 19,025


Je dois faire un ajustement quadratique de la forme m(t)=at^2 + bt +c.
Ce modèle est linéaire.
Grâce à R, je suis censé trouver les estimations des coefficients:

a = -0,00002875
b = 0,04564
c = 8,453

je voudrais savoir comment retrouver ces valeurs sous R.

Merci!

Nicolas Péru
Messages : 1408
Enregistré le : 07 Aoû 2006, 08:13

Messagepar Nicolas Péru » 07 Mai 2008, 09:28

Pour les modèles linéaires on peut utiliser la fonction lm()
vous donnera la manière de rédiger le modèle en langage S.

Ce qui n'est pas spécifié dans l'aide c'est l'ajout d'un paramètre d'ordre 2 qui s'écrit pour une varaible y en fonction d'une variable x avec un élément quadratique :

Code : Tout sélectionner

monmodele <- lm(y~x+I(x)^2)


C'est la forme de base. Il vaut mieux être familier de la théorie du modèle linéaire.
Ce qui est utile ensuite ce sont les fonctions summary() et anova()...

Bon courage.

Nicolas

Renaud Lancelot
Messages : 2484
Enregistré le : 16 Déc 2004, 08:01
Contact :

Messagepar Renaud Lancelot » 07 Mai 2008, 09:30

En utilisant la fonction lm. Il va falloir transformer le tableau (probablement en utilisant la fonction reshape) pour lui donner une forme d'un tableau (data.frame) avec deux colonnes: temps (t) à partir de l'origine (à voir comment vous codez le temps: en jours, ou en mois,...) et valeur observée m. NB: éviter de donner le nom "t" au temps car c'est un nom de fct R (transposition).

Renaud

Terry Gilles
Messages : 8
Enregistré le : 23 Avr 2008, 12:20

Messagepar Terry Gilles » 07 Mai 2008, 09:39

Je vais tester...
J'ai pas l'habitude d'utiliser R, donc il se peut que je repose des questions.

Merci du coup de pouce!

Terry Gilles
Messages : 8
Enregistré le : 23 Avr 2008, 12:20

Messagepar Terry Gilles » 13 Mai 2008, 08:25

Re,

j'ai transformé le tableau de cette manière:

w z x y
1 Janvier 1982 7.29
2 Février 1982 6.60
3 Mars 1982 7.83
4 Avril 1982 9.14
5 Mai 1982 8.89
6 Juin 1982 9.37
...

je n'ai pas utilisé la fonction reshape vu que je ne sais pas comment écrire.
J'ai fait ça manuellement sous excel.

Puis, j'ai écrit le code suivant sous R:

chroniquedata1<-read.table("X:/chroniquedata1.txt", sep = " ", header = TRUE)
chroniquedata1

chroniquedata1[1 :12,]
serie=chroniquedata1[,4]

serie.ts=ts(serie,start=c(1982,1),frequency=12)

plot(serie.ts,main="titre",xlab="Axe des temps",ylab="axe des ordonnées")

chroniquedata1 <- lm(y~x+I(x)^2)



Alors là, avec la dernière ligne, je voudrais bien obtenir une équation du 2e degré à 3 coefficients mais forcemment ça ne marche pas.

ça affiche: "Erreur dans eval(expr, envir, enclos) : objet "y" non trouvé"

Pourtant j'ai bien mis "header = TRUE"

si quelqu'un pourrait donner suite à ce fil, ce serait sympa!

Eric Pagot
Messages : 195
Enregistré le : 15 Fév 2007, 17:10

Messagepar Eric Pagot » 13 Mai 2008, 08:27

peut-être qu'il faut préciser d'où viennent les données (data=...)
Vétérinaire CTPA

Terry Gilles
Messages : 8
Enregistré le : 23 Avr 2008, 12:20

Messagepar Terry Gilles » 13 Mai 2008, 08:34

comment ça? les données sont au début du post.

le nom de mon fichier et où il est localisé: chroniquedata1<-read.table("X:/...")

peut-être que j'ai mal compris...

Aurélien Madouasse
Messages : 352
Enregistré le : 26 Fév 2007, 11:23

Messagepar Aurélien Madouasse » 13 Mai 2008, 08:42

Bonjour,

Il faut dire à lm ou prendre les données:

Code : Tout sélectionner

chroniquedata2 <- lm(y~x+I(x)^2, data=chroniquedata1)


Aurélien

Terry Gilles
Messages : 8
Enregistré le : 23 Avr 2008, 12:20

Messagepar Terry Gilles » 13 Mai 2008, 08:51

Merci Aurélien,

je pouvais pas savoir, j'ai copié ton code et ça me sort:

Call:
lm(formula = y ~ x + I(x)^2, data = chroniquedata1)

Coefficients:
(Intercept) x I(x)
-868.4767 0.4424 NA



c'est pas vraiment ce que je veux obtenir...

Tillard
Messages : 87
Enregistré le : 17 Déc 2004, 10:32

Messagepar Tillard » 13 Mai 2008, 09:00

Bonjour

Code : Tout sélectionner

chroniquedata2 <- lm(y~x+I(x^2), data=chroniquedata1)

I(x^2) et pas I(x)^2
cordialement
Emmanuel Tillard
UMR ERRC (Elevage des Ruminants en Regions Chaudes)
CIRAD - St PIERRE (La Réunion)
tel: 02 62 49 92 54

Terry Gilles
Messages : 8
Enregistré le : 23 Avr 2008, 12:20

Messagepar Terry Gilles » 13 Mai 2008, 09:42

Merci Tillard,

en effet ça m'affiche 3 valeurs:

Call:
lm(formula = y ~ x + I(x^2), data = chroniquedata1)

Coefficients:
(Intercept) x I(x^2)
-6.589e+03 6.180e+00 -1.439e-03


mais je voudrais bien obtenir à la place:
a = −0,00002875
b = 0,04564
c = 8,453

bon, je donne mon fichier au complet si vous voulez bien essayer de votre côté:

w z x y
1 Janvier 1984 08,844
2 Février 1984 08,882
3 Mars 1984 08,92
4 Avril 1984 08,958
5 Mai 1984 08,996
6 Juin 1984 09,033
7 Juillet 1984 09,071
8 Août 1984 09,109
9 Septembre 1984 09,147
10 Octobre 1984 09,185
11 Novembre 1984 09,223
12 Décembre 1984 09,26
1 Janvier 1985 09,298
2 Février 1985 09,336
3 Mars 1985 09,374
4 Avril 1985 09,412
5 Mai 1985 09,45
6 Juin 1985 09,488
7 Juillet 1985 09,525
8 Août 1985 09,563
9 Septembre 1985 09,601
10 Octobre 1985 09,639
11 Novembre 1985 09,677
12 Décembre 1985 09,715
1 Janvier 1986 09,752
2 Février 1986 09,79
3 Mars 1986 09,828
4 Avril 1986 09,866
5 Mai 1986 09,904
6 Juin 1986 09,942
7 Juillet 1986 09,98
8 Août 1986 10,017
9 Septembre 1986 10,055
10 Octobre 1986 10,093
11 Novembre 1986 10,131
12 Décembre 1986 10,169
1 Janvier 1987 10,207
2 Février 1987 10,244
3 Mars 1987 10,282
4 Avril 1987 10,32
5 Mai 1987 10,358
6 Juin 1987 10,396
7 Juillet 1987 10,434
8 Août 1987 10,472
9 Septembre 1987 10,509
10 Octobre 1987 10,547
11 Novembre 1987 10,585
12 Décembre 1987 10,623
1 Janvier 1988 10,661
2 Février 1988 10,699
3 Mars 1988 10,736
4 Avril 1988 10,774
5 Mai 1988 10,812
6 Juin 1988 10,85
7 Juillet 1988 10,888
8 Août 1988 10,926
9 Septembre 1988 10,964
10 Octobre 1988 11,001
11 Novembre 1988 11,039
12 Décembre 1988 11,077
1 Janvier 1989 11,115
2 Février 1989 11,153
3 Mars 1989 11,191
4 Avril 1989 11,228
5 Mai 1989 11,266
6 Juin 1989 11,304
7 Juillet 1989 11,342
8 Août 1989 11,38
9 Septembre 1989 11,418
10 Octobre 1989 11,456
11 Novembre 1989 11,493
12 Décembre 1989 11,531
1 Janvier 1990 11,569
2 Février 1990 11,607
3 Mars 1990 11,645
4 Avril 1990 11,683
5 Mai 1990 11,72
6 Juin 1990 11,758
7 Juillet 1990 11,796
8 Août 1990 11,834
9 Septembre 1990 11,872
10 Octobre 1990 11,91
11 Novembre 1990 11,947
12 Décembre 1990 11,985
1 Janvier 1991 12,023
2 Février 1991 12,061
3 Mars 1991 12,099
4 Avril 1991 12,137
5 Mai 1991 12,175
6 Juin 1991 12,212
7 Juillet 1991 12,25
8 Août 1991 12,288
9 Septembre 1991 12,326
10 Octobre 1991 12,364
11 Novembre 1991 12,402
12 Décembre 1991 12,439
1 Janvier 1992 12,477
2 Février 1992 12,515
3 Mars 1992 12,553
4 Avril 1992 12,591
5 Mai 1992 12,629
6 Juin 1992 12,667
7 Juillet 1992 12,704
8 Août 1992 12,742
9 Septembre 1992 12,78
10 Octobre 1992 12,818
11 Novembre 1992 12,856
12 Décembre 1992 12,894
1 Janvier 1993 12,931
2 Février 1993 12,969
3 Mars 1993 13,007
4 Avril 1993 13,045
5 Mai 1993 13,083
6 Juin 1993 13,121
7 Juillet 1993 13,159
8 Août 1993 13,196
9 Septembre 1993 13,234
10 Octobre 1993 13,272
11 Novembre 1993 13,31
12 Décembre 1993 13,348
1 Janvier 1994 13,386
2 Février 1994 13,423
3 Mars 1994 13,461
4 Avril 1994 13,499
5 Mai 1994 13,537
6 Juin 1994 13,575
7 Juillet 1994 13,613
8 Août 1994 13,651
9 Septembre 1994 13,688
10 Octobre 1994 13,726
11 Novembre 1994 13,764
12 Décembre 1994 13,802
1 Janvier 1995 13,84
2 Février 1995 13,878
3 Mars 1995 13,915
4 Avril 1995 13,953
5 Mai 1995 13,991
6 Juin 1995 14,029
7 Juillet 1995 14,067
8 Août 1995 14,105
9 Septembre 1995 14,143
10 Octobre 1995 14,18
11 Novembre 1995 14,218
12 Décembre 1995 14,256
1 Janvier 1996 14,294
2 Février 1996 14,332
3 Mars 1996 14,37
4 Avril 1996 14,407
5 Mai 1996 14,445
6 Juin 1996 14,483
7 Juillet 1996 14,521
8 Août 1996 14,559
9 Septembre 1996 14,597
10 Octobre 1996 14,635
11 Novembre 1996 14,672
12 Décembre 1996 14,71
1 Janvier 1997 14,748
2 Février 1997 14,786
3 Mars 1997 14,824
4 Avril 1997 14,862
5 Mai 1997 14,899
6 Juin 1997 14,937
7 Juillet 1997 14,975
8 Août 1997 15,013
9 Septembre 1997 15,051
10 Octobre 1997 15,089
11 Novembre 1997 15,126
12 Décembre 1997 15,164
1 Janvier 1998 15,202
2 Février 1998 15,24
3 Mars 1998 15,278
4 Avril 1998 15,316
5 Mai 1998 15,354
6 Juin 1998 15,391
7 Juillet 1998 15,429
8 Août 1998 15,467
9 Septembre 1998 15,505
10 Octobre 1998 15,543
11 Novembre 1998 15,581
12 Décembre 1998 15,618
1 Janvier 1999 15,656
2 Février 1999 15,694
3 Mars 1999 15,732
4 Avril 1999 15,77
5 Mai 1999 15,808
6 Juin 1999 15,846
7 Juillet 1999 15,883
8 Août 1999 15,921
9 Septembre 1999 15,959
10 Octobre 1999 15,997
11 Novembre 1999 16,035
12 Décembre 1999 16,073
1 Janvier 2000 16,11
2 Février 2000 16,148
3 Mars 2000 16,186
4 Avril 2000 16,224
5 Mai 2000 16,262
6 Juin 2000 16,3
7 Juillet 2000 16,338
8 Août 2000 16,375
9 Septembre 2000 16,413
10 Octobre 2000 16,451
11 Novembre 2000 16,489
12 Décembre 2000 16,527
1 Janvier 2001 16,565
2 Février 2001 16,602
3 Mars 2001 16,64
4 Avril 2001 16,678
5 Mai 2001 16,716
6 Juin 2001 16,754
7 Juillet 2001 16,792
8 Août 2001 16,83
9 Septembre 2001 16,867
10 Octobre 2001 16,905
11 Novembre 2001 16,943
12 Décembre 2001 16,981
1 Janvier 2002 17,019
2 Février 2002 17,057
3 Mars 2002 17,094
4 Avril 2002 17,132
5 Mai 2002 17,17
6 Juin 2002 17,208
7 Juillet 2002 17,246
8 Août 2002 17,284
9 Septembre 2002 17,322
10 Octobre 2002 17,359
11 Novembre 2002 17,397
12 Décembre 2002 17,435
1 Janvier 2003 17,473
2 Février 2003 17,511
3 Mars 2003 17,549
4 Avril 2003 17,586
5 Mai 2003 17,624
6 Juin 2003 17,662
7 Juillet 2003 17,7
8 Août 2003 17,738
9 Septembre 2003 17,776
10 Octobre 2003 17,814
11 Novembre 2003 17,851
12 Décembre 2003 17,889
1 Janvier 2004 17,927
2 Février 2004 17,965
3 Mars 2004 18,003
4 Avril 2004 18,041
5 Mai 2004 18,078
6 Juin 2004 18,116
7 Juillet 2004 18,154
8 Août 2004 18,192
9 Septembre 2004 18,23
10 Octobre 2004 18,268
11 Novembre 2004 18,305
12 Décembre 2004 18,343
1 Janvier 2005 18,381
2 Février 2005 18,419
3 Mars 2005 18,457
4 Avril 2005 18,495
5 Mai 2005 18,533
6 Juin 2005 18,57
7 Juillet 2005 18,608
8 Août 2005 18,646
9 Septembre 2005 18,684
10 Octobre 2005 18,722
11 Novembre 2005 18,76
12 Décembre 2005 18,797
1 Janvier 2006 18,835
2 Février 2006 18,873
3 Mars 2006 18,911
4 Avril 2006 18,949
5 Mai 2006 18,987
6 Juin 2006 19,025

Tillard
Messages : 87
Enregistré le : 17 Déc 2004, 10:32

Messagepar Tillard » 13 Mai 2008, 12:33

Code : Tout sélectionner

chronique[1:20,]
> chronique[1:20,]
   mois    mois_l annee     y  x
1     1   Janvier  1984 8.844  0
2     2   Fevrier  1984 8.882  1
3     3      Mars  1984 8.920  2
4     4     Avril  1984 8.958  3
5     5       Mai  1984 8.996  4
6     6      Juin  1984 9.033  5
7     7   Juillet  1984 9.071  6
8     8      Aout  1984 9.109  7
9     9 Septembre  1984 9.147  8
10   10   Octobre  1984 9.185  9
11   11  Novembre  1984 9.223 10
12   12  Decembre  1984 9.260 11
13    1   Janvier  1985 9.298 12
14    2   Fevrier  1985 9.336 13
15    3      Mars  1985 9.374 14
16    4     Avril  1985 9.412 15
17    5       Mai  1985 9.450 16
18    6      Juin  1985 9.488 17
19    7   Juillet  1985 9.525 18
20    8      Aout  1985 9.563 19


plot(y ~ x, chronique)  #c'est une droite

model1 <- lm(y ~ x + I(x^2), chronique)
summary(model1)

Call:
lm(formula = y ~ x + I(x^2), data = chronique)

Residuals:
       Min         1Q     Median         3Q        Max
-5.015e-04 -2.455e-04 -1.188e-06  2.492e-04  4.955e-04

Coefficients:
             Estimate Std. Error   t value Pr(>|t|)   
(Intercept) 8.844e+00  5.250e-05 1.685e+05   <2e-16 ***
x           3.785e-02  9.017e-07 4.197e+04   <2e-16 ***
I(x^2)      3.128e-10  3.245e-09 9.600e-02    0.923   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.0002897 on 267 degrees of freedom
Multiple R-squared:     1,      Adjusted R-squared:     1
F-statistic: 1.4e+10 on 2 and 267 DF,  p-value: < 2.2e-16


le temps a été recodé en mois, de 0 à 269
les coef sont proches de ceux que vous recherchez mais ils ne sont pas identiques
Cordialement
Emmanuel Tillard

UMR ERRC (Elevage des Ruminants en Regions Chaudes)

CIRAD - St PIERRE (La Réunion)

tel: 02 62 49 92 54

Terry Gilles
Messages : 8
Enregistré le : 23 Avr 2008, 12:20

Messagepar Terry Gilles » 13 Mai 2008, 14:10

moi quand j'écris: plot(y ~ x, chronique), ça me sort pas une courbe mais des carrés coloriés en dégradé.

Puis, model1 <- lm(y ~ x + I(x^2), chronique), ça m'affiche un message d'erreur:

Erreur dans storage.mode(y) <- "double" :
la modification du mode de stockage d'un objet 'factor' n'est pas autorisée
De plus : Warning message:
In model.response(mf, "numeric") :
using type="numeric" with a factor response will be ignored


je vois pas bien ce qui colle pas avec moi

Nicolas Péru
Messages : 1408
Enregistré le : 07 Aoû 2006, 08:13

Messagepar Nicolas Péru » 13 Mai 2008, 14:38

bonjour,

peut être un problème de class sur y ou x (plus vraissemblablement sur x) : x devrait être de la class"vector" est ce le cas ?

Mais il faudrait d'abord savoir comment ont été obtenu les valeurs des coefficients que vous recherchez? Comment le temps a t'il été pris en compte par exemple, variable continue ou discrète ?
De ces réponses dépend la manière dont vous aller coder vos variable dans R.

Nicolas.

Terry Gilles
Messages : 8
Enregistré le : 23 Avr 2008, 12:20

Messagepar Terry Gilles » 13 Mai 2008, 14:57

Ok, je verrai de mon côté comment je peux régler ça.

Merci quand même à tous pour votre aide!

je reviendrai sans doute sur le forum demain.

a+


Retourner vers « Questions en cours »

Qui est en ligne

Utilisateurs parcourant ce forum : Google [Bot] et 1 invité