On cherche souvent \u00e0 mod\u00e9liser un \u00e9chantillon par une loi de probabilit\u00e9.
A partir d’un jeu de donn\u00e9es, comment peut-on trouver les param\u00e8tres d’une loi pr\u00e9alablement fix\u00e9e?<\/p>\n
Plusieurs m\u00e9thodes peuvent \u00eatre utilis\u00e9es.<\/p>\n
On prend l’exemple ici du d\u00e9lai entre l’infection d’un individu et la d\u00e9tection de cet individu comme malade.
On mod\u00e9lise ici ce d\u00e9lai (en jours) par une loi de Weibull (on peut aussi essayer les lois gamma et lognormale par exemple)<\/p>\n
La m\u00e9thode 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 \u00e7a donne pour une loi Weibull.<\/p>\n
#on charge le package Mass
library(MASS)
# on appelle z le d\u00e9lai entre infection et d\u00e9tection, z est un vecteur contenant les donn\u00e9es
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\u00e9es, ici un histogramme est le plus adapt\u00e9!
hist(z)
# on utilise la fonction fitdistr pour une loi Weibull, regarder ?fitdistr
paraw <- fitdistr(z,densfun=\"weibull\")
logLik(paraw)\u00a0 # on peut avoir le loglikelihood
# on visualise les r\u00e9sultats sur un graphique : histogramme+ loi
hist(z,freq=FALSE) # penser \u00e0 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\u00e8tres estim\u00e9s
#regarder ?dweibull,
#on peut simuler de nombreuses lois sur R par exemple, gamma avec dgamma(), normale avec dnorm()...
<\/code><\/pre>\nOn peut \u00e9galement \u00e9crire directement la vraisemblance et trouver le maximum de vraisemblance.
Dans cet exemple on utilise le package bbmle<\/p>\n
library(bbmle)
#on \u00e9crit la vraisemblance pour une weibull (loglikelihood ici)
weiblikfun<-function(shape,scale){
\u00a0\u00a0 \u00a0-sum((dweibull(z,shape=shape,scale=scale,log=TRUE)))}
# on cherche \u00e0 minimiser le loglikelihood avec mle2, regarder ?mle2!!!!
w<-mle2(weiblikfun,start=list(shape= 10,scale=20))
w #donne les param\u00e8tres estim\u00e9s
confint(w) #donne des intervalles de confiance
<\/code><\/pre>\nEnfin on peut estimer les param\u00e8tres dans un cadre bay\u00e9sien en lan\u00e7ant 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\u00a0par exemple.
Enfin on peut analyser le resultat des MCMC via le package coda.<\/p>\n
library(R2WinBUGS)
library(coda)
#on met en forme les donn\u00e9es pour Winbugs (en liste)
data<-list(T=z,N=length(z))
<\/code><\/pre>\nIl faut \u00e9crire un code Winbugs dans un fichier txt et le mettre dans le r\u00e9pertoire courant.
Voil\u00e0 ce que \u00e7a peut donner pour une weibull (ceci correspond \u00e0 modelweibull.txt plac\u00e9 dans le r\u00e9pertoire courant):<\/p>\n
#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)
}
}
<\/code><\/pre>\nEnfin on lance le tout et on regarde le comportement des MCMC!<\/p>\n
rweibull<-bugs(data,inits=NULL,parameters.to.save=c(\"v\",\"s\"),
\u00a0\u00a0 \u00a0\u00a0\u00a0 \u00a0model.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)
<\/code><\/pre>\nCe long exemple peut \u00eatre utilis\u00e9 pour d’autres lois. Je vous conseille vivement de bien regarder les expressions analytiques utilis\u00e9es pour chaque fonction. Par exemple la loi de Weibull peut \u00eatre \u00e9crite de diff\u00e9rentes fa\u00e7ons. Vous aurez peut \u00eatre remarqu\u00e9 que les fonctions utilis\u00e9es dans R et winbugs n’utilisent pas la m\u00eame expression: ATTENTION!!<\/p>\n
Enfin, il est souvent possible de mod\u00e9liser des donn\u00e9es par diff\u00e9rentes lois, il s’agit ensuite de trouver celle qui correspond le mieux….<\/p>\n","protected":false},"excerpt":{"rendered":"
On cherche souvent \u00e0 mod\u00e9liser un \u00e9chantillon par une loi de probabilit\u00e9. A partir d’un jeu de donn\u00e9es, comment peut-on trouver les param\u00e8tres d’une loi pr\u00e9alablement fix\u00e9e? Plusieurs m\u00e9thodes peuvent \u00eatre utilis\u00e9es. On prend l’exemple ici du d\u00e9lai entre l’infection d’un individu et la d\u00e9tection de cet individu comme malade.On mod\u00e9lise ici ce d\u00e9lai (en jours) par une loi de Weibull (on peut aussi essayer les lois gamma et lognormale par exemple) La m\u00e9thode 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 \u00e7a donne pour une loi Weibull. #onRead More →<\/a><\/p>\n","protected":false},"author":6,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"_acf_changed":false,"content-type":"","rop_custom_images_group":[],"rop_custom_messages_group":[],"rop_publish_now":"initial","rop_publish_now_accounts":{"twitter_399453572_399453572":""},"rop_publish_now_history":[],"rop_publish_now_status":"pending","jetpack_post_was_ever_published":false,"_jetpack_newsletter_access":"","_jetpack_dont_email_post_to_subs":false,"_jetpack_newsletter_tier_id":0,"_jetpack_memberships_contains_paywalled_content":false,"_jetpack_memberships_contains_paid_content":false,"footnotes":""},"categories":[5,11],"tags":[],"class_list":{"0":"entry","1":"post","2":"publish","3":"author-melen","4":"post-553","6":"format-standard","7":"category-bayesien","8":"category-inference"},"acf":[],"jetpack_featured_media_url":"","jetpack_shortlink":"https:\/\/wp.me\/p9O7Sx-8V","jetpack_sharing_enabled":true,"_links":{"self":[{"href":"https:\/\/thinkr.fr\/abcdr\/wp-json\/wp\/v2\/posts\/553","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/thinkr.fr\/abcdr\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/thinkr.fr\/abcdr\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/thinkr.fr\/abcdr\/wp-json\/wp\/v2\/users\/6"}],"replies":[{"embeddable":true,"href":"https:\/\/thinkr.fr\/abcdr\/wp-json\/wp\/v2\/comments?post=553"}],"version-history":[{"count":2,"href":"https:\/\/thinkr.fr\/abcdr\/wp-json\/wp\/v2\/posts\/553\/revisions"}],"predecessor-version":[{"id":4085,"href":"https:\/\/thinkr.fr\/abcdr\/wp-json\/wp\/v2\/posts\/553\/revisions\/4085"}],"wp:attachment":[{"href":"https:\/\/thinkr.fr\/abcdr\/wp-json\/wp\/v2\/media?parent=553"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/thinkr.fr\/abcdr\/wp-json\/wp\/v2\/categories?post=553"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/thinkr.fr\/abcdr\/wp-json\/wp\/v2\/tags?post=553"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}