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
  • sig - cartographie
  • Test
  • tidyverse
  • Transformation de données


inférence

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...
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 estimer les paramètres d’une loi de probabilité avec R ?

1
dans bayésien, inférence
- ça ne sert à rien -- c\'est interessant - (score de +13 sur 13 votes)
Loading...
On cherche souvent à modéliser un échantillon par une loi de probabilité.
A partir d’un jeu de données, comment peut-on trouver les paramètres d’une loi préalablement fixée?

Plusieurs méthodes peuvent être utilisées.

On prend l’exemple ici du délai entre l’infection d’un individu et la détection de cet individu comme malade.
On modélise ici ce délai (en jours) par une loi de Weibull (on peut aussi essayer les lois gamma et lognormale par exemple)

La méthode la plus simple est d’utiliser la fonction fitdistr du package MASS.
Cette fonction permet d’ajuster de nombreuses lois par maximum de vraisemblance. Regardons ce que ça donne pour une loi Weibull.

 
 
#on charge le package Mass
 library(MASS) 
 
# on appelle z le délai entre infection et détection, z est un vecteur contenant les données
 
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)
 
#toujours visualiser ses données, ici un histogramme est le plus adapté!
 
 hist(z) 
 
# on utilise la fonction fitdistr pour une loi Weibull, regarder ?fitdistr
 
 paraw <- fitdistr(z,densfun="weibull") 
 logLik(paraw)  # on peut avoir le loglikelihood
 
# on visualise les résultats sur un graphique : histogramme+ loi 
 
 hist(z,freq=FALSE) # penser à frep=FALSE!!
 lines(dweibull(0:max(z),shape=paraw$estimate[1],scale=paraw$estimate[2]),type='l',col='green',lwd=2)
 
#on simule une loi de weibull avec dweibull en utilisant les paramètres estimés
#regarder ?dweibull, 
#on peut simuler de nombreuses lois sur R par exemple, gamma avec dgamma(), normale avec dnorm()...
 
 


On peut également écrire directement la vraisemblance et trouver le maximum de vraisemblance.
Dans cet exemple on utilise le package bbmle

 
 
library(bbmle) 
 
#on écrit la vraisemblance pour une weibull (loglikelihood ici)
 
 weiblikfun<-function(shape,scale){
    -sum((dweibull(z,shape=shape,scale=scale,log=TRUE)))}
 
# on cherche à minimiser le loglikelihood avec mle2, regarder ?mle2!!!!
 
 w<-mle2(weiblikfun,start=list(shape= 10,scale=20)) 
 
w #donne les paramètres estimés
confint(w) #donne des intervalles de confiance
 
 


Enfin on peut estimer les paramètres dans un cadre bayésien en lançant des MCMC.
Le plus simple est d’utiliser le Gibbs sampler via Winbugs ou OpenBugs.
On peut lancer ces derniers via R et le package R2WinBUGS par exemple.
Enfin on peut analyser le resultat des MCMC via le package coda.

 
 
 library(R2WinBUGS)
library(coda) 
 
#on met en forme les données pour Winbugs (en liste)
 data<-list(T=z,N=length(z))
 
 


Il faut écrire un code Winbugs dans un fichier txt et le mettre dans le répertoire courant.
Voilà ce que ça peut donner pour une weibull (ceci correspond à modelweibull.txt placé dans le répertoire courant):

 
 
 #model
 
model{
 
#priors
 
v~dgamma(0.01,0.01)
 
#v~dlnorm(0,0.0001)
lambda~dbeta(1,1)
 
s<-pow((1/lambda),(1/v))
 
#likelihood
for(i in 1:N) {
 
T[i]~dweib(v,lambda)
 
}
 
} 
 


Enfin on lance le tout et on regarde le comportement des MCMC!
 
 
 rweibull<-bugs(data,inits=NULL,parameters.to.save=c("v","s"),
        model.file="modelweibull.txt",n.chains=4,n.iter=100000,codaPkg=TRUE)
 
rweibull.coda<-read.bugs(rweibull)
summary(rweibull.coda)
xyplot(rweibull.coda)
acfplot(rweibull.coda)
densityplot(rweibull.coda,col="blue")
rejectionRate(rweibull.coda)
 
 


Ce long exemple peut être utilisé pour d’autres lois. Je vous conseille vivement de bien regarder les expressions analytiques utilisées pour chaque fonction. Par exemple la loi de Weibull peut être écrite de différentes façons. Vous aurez peut être remarqué que les fonctions utilisées dans R et winbugs n’utilisent pas la même expression: ATTENTION!!

Enfin, il est souvent possible de modéliser des données par différentes lois, il s’agit ensuite de trouver celle qui correspond le mieux….

Proposé par Melen.

1 commentaire. Cliquez ici pour réagir.
Formation logiciel R