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


modélisation

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

8
dans inférence, modélisation
- ça ne sert à rien -- c\'est interessant - (score de +12 sur 14 votes)
Loading ... Loading ...
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.

Déjà 8 commentaires. Cliquez ici pour réagir.

Comment obtenir le critère d’Akaike avec R ? : AIC

1
dans modélisation
- ça ne sert à rien -- c\'est interessant - (score de +7 sur 7 votes)
Loading ... Loading ...
Le critère d’Akaike (AIC) est un critère utilisé pour la sélection de modèles. Ce critère représente un compromis entre le biais diminuant avec le nombre de paramètres libres et la parcimonie, volonté de décrire les données avec le plus petit nombre de paramètres possible. Il se calcule de la façon suivante -2*log-likelihood + k*npar. Par défaut on a souvent k=2. Le meilleur modèle est celui qui possède l’AIC le plus faible.

On obtient ce critère en utilisant la fonction AIC(objet,k=?), k=2 par défaut. Prenons un exemple

 
 
 
 library(MASS) #pour la fonction fitdistr
 
#  z est un vecteur contenant les données, on essaie de modéliser ces données par une distribution de weibull ou gamma
 
z<-c(14,14,14,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,23)
 
# on utilise la fonction fitdistr pour une loi Weibull, regarder ?fitdistr
 
paraw <- fitdistr(z,densfun="weibull") 
loglikw<-logLik(paraw)  # on peut avoir le loglikelihood
loglikw
 
#on fait la même chose pour une loi gamma
 
parag<-fitdistr(z,densfun="gamma")
loglikg<-logLik(parag)
loglikg
 
#AIC avec k=2 par défaut
 
AIC(paraw)#environ 204
AIC(parag)#environ 209
 
#d'après ce critère on choisirait la loi de weibull
 
#on vérifie en recalculant les AIC directement
 
akaike<-function(npar,loglik,k){-2*loglik+k*npar}
 
akaike(2,loglikw[1],2)
akaike(2,loglikg[1],2)
 
 
 


L’AIC est un critère comme les autres, faites en bon usage et n’en n’abusez pas trop!


Proposé par Melen.

1 commentaire. Cliquez ici pour réagir.

Comment estimer les paramètres d’un modèle d’équations différentielles ordinaires par les moindres carrés avec R ?

0
dans modélisation
- ça ne sert à rien -- c\'est interessant - (score de +4 sur 4 votes)
Loading ... Loading ...
De nombreux systèmes sont modélisés par des équations différentielles ordinaires. R peut permettre de résoudre certains de ces systèmes et aussi d’estimer leurs paramètres.

On prend ici l’exemple d’un modèle épidémiologique temporel SI.

 
 
 #edo SIX
 
library(deSolve) # on utilise deSolve
 
#on définit le système dans une fonction six ici
 
six<-function(t,x,parms){
    with( as.list(c(parms,x)),{
 
    rp<-ap*exp(-bp*t)
    rs<-as*exp(-0.5*(log(t/gs)/bs)^2)
    
    dI<- (rp*X*S+rs*I*S)
    dS<- -(rp*X*S+rs*I*S)
 
    res<-c(dI,dS)
    list(res)})
}
 
#on définit les paramètres pour la simulation du système
 
parms<-c(ap=0.002, bp=0.0084, as=5.9e-7, bs=0.25, gs=1396, X=1, N=1010)
 
#on crée un vecteur pour le temps
times<-seq(0:3000)
 
#valeurs initiales des variables (ici tous les individus sont sains au début)
 
y<- xstart <-(c(I = 0, S = 1010))
 
#on résout le système avec la fonction lsoda
 
out<-as.data.frame(lsoda(xstart, times, six, parms))
 
# toujours regarder le résultat pour détecter les incohérences
 
plot(out$S~out$time, type="l",col='blue')
plot(out$I~out$time, type="l",col='green',lwd=3)
 
#on crée un jeu de données en ajoutant un bruit blanc (gaussien)
 
noisy<-out$I+rnorm(nrow(out),sd=0.15*out$I)
 
plot(noisy)
 
#on ajuste le modèle sur ces données par les moindres carrés
# pour cela on utilise la fonction nls, il faut faire attention aux paramètres initiaux
 
fitnoisy<-nls(noisy~lsoda(xstart,times,six,c(ap=ap, bp=bp, as=as, bs=bs, gs=gs, X=1, N=1010))[,2]
     ,start=list(ap=0.002, bp=0.0084, as=5.9e-7, bs=0.25, gs=1000)
    , control=list(minFactor=1e-20,maxiter=150))
 
summary(fitnoisy)
 
# un beau graphique pour visualiser tout ça!
 
plot(noisy,pch=20,col='grey',main="SIX ode Model Fit",xlab="time",ylab="disease progress")
lines(predict(fitnoisy,times),type='l',col='green',lwd=4.5)
 
 
 


Vous pouvez définir différents systèmes d’équations différentielles en conservant la même structure.
Enfin pour estimer les paramètres il est possible d’utiliser d’autres fonctions et/ou d’autres méthodes (vraisemblance par exemple).
Pour que nls converge il faut donner des valeurs initiales le plus proche possibles des « vraies valeurs ». Je vous conseille donc d’essayer de fitter votre système graphiquement. Si ça ne marche pas vous pouvez utiliser nls2 avec l’algorithme bruteforce…

Proposé par Melen.

Pas encore de commentaire, cliquez ici pour réagir.
Formation logiciel R