1. Grâce à ce formulaire vous pouvez proposer une astuce ou un script sur R.
  2. Votre script doit pouvoir être lancé en l'état. Veuillez penser à :
    • inclure le chargement des "library" nécessaires
    • construire un petit jeu de données si besoin est.
    • commenter les lignes de codes pour en faciliter la compréhension.
  3. Le titre de votre script ou astuce doit être clair et explicite.
  4. Pensez à mettre votre code entre les balises [R] et [/R]; Pour cela, vous pouvez utiliser le bouton

Vous pouvez utiliser vos comptes Facebook, twitter ou google pour vous identifer (google est compatible yahoo, openID...)
L'ideal étant de vous connecter si vous avez un compte utilisateur, ou faire une demande de compte utilisateur si vous n'en avez pas encore.
Créer un compte va vous permettre de pouvoir éditer vos codes et de mettre en avant votre site internet.
Sinon vous pouvez soumettre anonymement en remplissant les champs ci-après.







Choisissez les catégories correspondantes à votre Code:

  • algorithmique
  • Analyse de survie
  • base indispensable
  • bayésien
  • configuration de R
  • exportation de données
  • fonctions utiles
  • graphique
  • importation de données
  • inférence
  • manipulation de données
  • message d'erreur
  • modélisation
  • Non classé
  • optimisation
  • planification
  • programmer avec R
  • regression linéaire
  • Test
  • tidyverse
  • Transformation de données


Comment ajuster des fonctions avec les moindres carrés ? nls

7
dans inférence, modélisation
La manière la plus simple d’ajuster une fonction à des données est la méthode « géométrique » des moindres carrés (minimiser la somme des carrés des écarts correspond à maximiser la vraisemblance avec une loi normale).

La fonction nls de R permet de réaliser ceci de manière simple.

Voyons deux exemples :

 
#exemple modèle de croissance exponentiel
 
#on crée une fonction qui correspond à un modèle de Malthus
 
malthus<-function(t,N0,r){N0*exp(r*t)}
 
t<-seq(0:100)
NO<-2
r<-0.05
 
plot(malthus(t,NO,r)~t,type='l',col='green',lwd=2)
 
#on crée des données en ajoutant du bruit
sim<-malthus(t,NO,r)+rnorm(t,sd=0.3*malthus(t,NO,r))
plot(sim~t,pch=20)
 
#on ajuste la fonction sur les données simulées en utilisant les moindres carrés
fitmalthus<-nls(sim~malthus(t,a,b),start=list(a=1,b=0.01))
fitmalthus
summary(fitmalthus)
 
#on vérifie
plot(sim~t,pch=20)
lines(malthus(t,NO,r)~t,type='l',col='green',lwd=2)
lines(malthus(t,coef(fitmalthus)[1],coef(fitmalthus)[2])~t,type='l',col='red',lwd=2)
 
#exemple 2 avec une fonction logistique
 
logist<-function(t,v,w,z){v/(1+w*exp(-z*t))}
 
v<-10;w<-8;z<-0.09;x<-seq(from=-100,to=100)
 
plot(logist(x,v,w,z)~x,type='l',col='blue',lwd=2)
 
#on crée des données avec du bruit
noisy<-rnorm(x,mean=logist(x,v,w,z),sd=0.9)
 
plot(noisy,pch=20)#notez qu'il y a des valeurs négatives pas réalistes...
 
fitlog<-nls(noisy~logist(x,a,b,c),start=list(a=5,b=5,c=0.02))
fitlog
summary(fitlog)
 
#on vérifie
plot(noisy~x,pch=20,col="grey")
lines(logist(x,v,w,z)~x,type='l',col='blue',lwd=2)
lines(logist(x,coef(fitlog)[1],coef(fitlog)[2],coef(fitlog)[3])~x,type='l',col='red',lwd=2) 
 


Ces deux exemples cachent en fait de nombreuses difficultés :

Pour que l’algorithme converge, vous devez donner des valeurs de départ des paramètres assez proches de la réalité –> essayez d’ajuster vos courbes à l’œil afin d’obtenir les meilleurs paramètres de départ possibles. Sinon vous pouvez utiliser la fonction nls2 et l’algorithme brute force, ou passer à la vraisemblance voire au merveilleux monde des mcmc…

La fonction nls présente, comme beaucoup de fonctions R, de nombreux arguments qui vous permettront peut être d’ajuster vos fonctions comme vous le souhaitez : control, trace, weights etc.

Enfin vous pouvez récupérer de nombreuses informations à partir de l’objet créé lors de l’ajustement : résidus, valeur des paramètres estimés etc…

Proposé par Melen.

Ce script vous à rendu service? remerciez l'auteur en votant ici:
- ça ne sert à rien -- c\'est interessant - (score de +11 sur 13 votes)
Loading ... Loading ...

7 Comments

  1. Stan
    Posté le 3 mars 2012 a 15 h 53 min | Permalink

    Juste une remarque, qui vaut aussi pour question : on peut aussi retrouver les coefficients du modèle ajusté via un modèle linéaire tout simple non ?

     
    #On crée le jeu de données (repris de la fonction de Malthus précédente avec le bruit)
    donnees=rnorm(t,mean=2*exp(0.05*t),sd=0.3*2*exp(0.05*t))
     
    #On y applique le modèle linéaire avec notre hypothèse de loi exponentielle
    LMdon=lm(log(donnees)~t)
    summary(LMdon)
     
    #Graphe de vérification
    plot(donnees~t,pch=20)
    lines(exp(coef(LMdon)[2]*t+coef(LMdon)[1])~t,type="l",col="darkred",lwd=2)
    lines(2*exp(0.05*t)~t,type="l",col="darkgreen",lwd=2) 
     

    Dans ce cas, qu’apporte le NLS par rapport au LM ?

  2. Posté le 7 mars 2012 a 12 h 40 min | Permalink

    Dans le cas où la fonction est linéarisable, utiliser lm ne doit en théorie pas changer beaucoup de chose. Par contre, c’est quand même très féquent qu’il ne soit pas possible de linéariser la fonction à ajuster…dans ce cas là si on veut utiliser les moindres carrés…Après en terme d’inférence statistique c’est souvent mieux d’utiliser une vraissemblance adaptée…

  3. Stan
    Posté le 8 mars 2012 a 18 h 34 min | Permalink

    Ah oui, avec par exemple une fonction du type y = a*x^b+c non linéarisable…

     
     #On crée une fonction 'non linéarisable'
    t=seq(0,100,1)
    donnees=rnorm(t,mean=2*t^0.2+5,sd=0.5)
    plot(donnees~t,pch=20)
     
    #On y applique le NLS avec notre hypothèse de loi...
    NLSdon=nls(donnees~a*t^b+c,start=list(a=1,b=0.1,c=1))
    summary(NLSdon)
    lines(coef(NLSdon)[1]*t^coef(NLSdon)[2]+coef(NLSdon)[3]~t,type="l",lwd=3,col="darkred")
     
    #On vérifie si le modèle trouvé par la fonction NLS est proche de la loi de départ...
    lines(2*t^0.2+5~t,type="l",col="darkgreen",lwd=3) 
     
  4. AC
    Posté le 13 juin 2012 a 15 h 55 min | Permalink

    Justement, comment récupérer les fameux paramètres calculés ? J’arrive à les afficher avec summary ou print, mais j’aimerais pouvoir les réutiliser sans avoir à les recopier à la main. Quelqu’un a une idée de la manière de les récupérer ?

  5. Vincent
    Posté le 26 juin 2012 a 17 h 29 min | Permalink

    @AC

    le resultat doit etre stocké dans un objet (une liste ou un tableau) ensuite vos pourrez y avoir accés.

  6. Melen
    Posté le 30 juillet 2012 a 9 h 24 min | Permalink

    c’est exactement ce qui est fait en utilisant la fonction coef(modèle); Il faut regarder dans quel ordre les coefficients sont stocké dans l’objet modèle et ensuite il suffit d’utiliser coef()[i] pour utiliser ces valeurs!

  7. Posté le 23 août 2012 a 9 h 04 min | Permalink

    @melen

    merci pour ce post, il me sert de pense bete :)

Poster un commentaire


Votre email ne seras jamais communiqué champs requis désigné par une *

*
*


× 8 = vingt quatre

Formation logiciel R