-
-
Notifications
You must be signed in to change notification settings - Fork 266
RStan Getting Started (Français)
(Translation: Virgile Baudrot)
RStan est une interface R à Stan. Pour plus d'informations sur Stan, et son langage de modélisation, vous pouvez visiter le site web :
La quasi totalité des instructions ci-dessous sont pour la version de RStan mentionnée ci-dessus.
Pour installer RStan, suivez le lien correspondant à votre plateforme :
Si vous suivez les instructions d'installation mais que l'installation ne fonctionne pas, vous pouvez obtenir de l'aide sur la mailing liste (liste de distribution) des utilisateurs de Stan (à noter que le forum est en anglais):
Le reste de ce document suppose que vous avez intallé RStan en suivant les instructions d'un des liens ci-dessus.
Le nom du package étant rstan (en minuscule), vous pouvez le charger en utilisant la fonction library("rstan")
.
library("rstan") # observez les messages de démarrage
Comme l'indique le message de démarrage, si vouzutilisez rstan locallement sur une machine multi-coeur et que vous avez beaucoup de RAM pour estimer votre modèle en parallèle, alors exécutez :
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
Ces options vous permettent (i) de sauvegarder automatiquement une version d'un programme Stan compilé sur le disque dur afin de ne pas avoir besoin de recompiler ainsi que (ii) d'exécuter plusieurs chaînes de Markov en parallèle.
C'est un exemple de la Section 5.5. de Gelman et al (2003), qui a étudié les effets de l'entraînement au sein de huit écoles. Pour simplifier, nous appellerons cet exemple "eight schools".
Nous commençons par écrire un programme Stan pour le model et le sauvons dans
un nouveau fichier 8schools.stan
.
// sauvegarder dans 8schools.stan
data {
int<lower=0> J; // number of schools
real y[J]; // estimated treatment effects
real<lower=0> sigma[J]; // s.e. of effect estimates
}
parameters {
real mu;
real<lower=0> tau;
real eta[J];
}
transformed parameters {
real theta[J];
for (j in 1:J)
theta[j] = mu + tau * eta[j];
}
model {
target += normal_lpdf(eta | 0, 1);
target += normal_lpdf(y | theta, sigma);
}
Dans ce modèle, nous déclarons theta
dans le bloc transformed parameters
puisqu'on le construit à partir des paramètres mu
et eta
plutôt que directement
comme un paramètre (qui serait alors dans le bloc parameters).
En paramètrant ainsi, l'échantillonneur sera plus efficace (voir l'explication détaillée)
Nous pouvons préparer les données dans R avec :
schools_dat <- list(J = 8,
y = c(28, 8, -3, 7, -1, 1, 18, 12),
sigma = c(15, 10, 16, 11, 9, 11, 10, 18))
Et nous pouvons obtenir un ajustement avec les commandes R suivantes. Notez que l'argument
pour file =
doit indiquer l'emplacement du fichier dans votre système, sauf si vous l'avez
placé dans l'espace de travail (le working directory) de R auquel cas les commandes ci-dessous fonctionnerons
fit <- stan(file = '8schools.stan', data = schools_dat,
iter = 1000, chains = 4)
Il est également possible de spécifier un modèle Stan en utilisant une chaîne de
caractère que l'on indiquera dans l'argument model_code
de la fonction Stan.
Cependant, cette manière de faire n'est pas recommandée puisque l'utilisation d'un fichier
permet d'afficher les déclarations et de déréférencer les numéros de lignes dans les messages d'erreurs.
L'objet fit
, qui est retourné par la fonction stan
est un objet S4 de classe
stanfit
. Les méthodes telles que print
, plot
, et pairs
sont associées aux
au résultat ajusté ce qui permet d'utiliser le code suivant pour faire une vérification des éléments inclus dans fit
.
print
renvoie un résumé des paramètres du modèle ainsi que le log-posterior dénoté lp__
(voir les sorties de l'exemple suivant).
Pour plus de méthodes et de détails de la classe stanfit
, utilisez l'aide de cette classe.
La fonction extract
appliquée à un objet stanfit
renvoie les données d'échantillonnage.
extract
extrait les échantillonnages d'un objet stanfit
sous forme d'une list de tableaux des
paramètres d'intérêt, ou juste un tableau.
De plus, les fonctions S3 as.array
et as.matrix
sont définies pour les objets stanfit
(utilisez la commande help("as.array.stanfit")
pour afficher l'aide dans R).
print(fit)
plot(fit)
pairs(fit, pars = c("mu", "tau", "lp__"))
la <- extract(fit, permuted = TRUE) # return a list of arrays
mu <- la$mu
### retourner un tableau de trois dimensions: iterations, chains, parameters
a <- extract(fit, permuted = FALSE)
### utiliser les fonctions S3 as.array (ou as.matrix) sur les objets stanfit
a2 <- as.array(fit)
m <- as.matrix(fit)
> print(fit, digits = 1)
Inference for Stan model: 8schools.
4 chains, each with iter=1000; warmup=500; thin=1;
post-warmup draws per chain=500, total post-warmup draws=2000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
mu 8.2 0.2 5.4 -1.9 4.8 8.1 11.3 19.3 480 1
tau 6.8 0.3 6.2 0.3 2.5 5.2 9.2 21.7 425 1
eta[1] 0.4 0.0 1.0 -1.5 -0.3 0.4 1.0 2.2 2000 1
eta[2] 0.0 0.0 0.8 -1.7 -0.6 0.0 0.5 1.7 2000 1
eta[3] -0.2 0.0 1.0 -2.1 -0.9 -0.2 0.4 1.7 2000 1
eta[4] -0.1 0.0 0.9 -1.8 -0.7 -0.1 0.5 1.7 2000 1
eta[5] -0.4 0.0 0.9 -2.1 -1.0 -0.4 0.2 1.4 2000 1
eta[6] -0.2 0.0 0.9 -1.9 -0.8 -0.2 0.4 1.5 1731 1
eta[7] 0.3 0.0 0.9 -1.4 -0.2 0.4 0.9 2.0 1507 1
eta[8] 0.0 0.0 0.9 -1.9 -0.6 0.0 0.7 1.8 1988 1
theta[1] 11.5 0.3 8.8 -2.4 5.9 10.1 15.6 32.9 977 1
theta[2] 7.8 0.1 6.2 -4.7 4.1 7.9 11.6 20.3 2000 1
theta[3] 6.1 0.2 7.7 -11.2 2.1 6.4 10.5 20.2 2000 1
theta[4] 7.6 0.1 6.5 -4.9 3.8 7.8 11.4 21.3 2000 1
theta[5] 5.0 0.1 6.6 -9.3 1.2 5.6 9.3 16.7 2000 1
theta[6] 6.2 0.2 6.7 -8.2 2.2 6.5 10.5 18.5 2000 1
theta[7] 10.8 0.2 7.0 -1.3 6.1 10.1 15.1 26.8 2000 1
theta[8] 8.7 0.2 8.2 -7.3 3.9 8.4 12.8 27.2 1446 1
lp__ -39.5 0.1 2.6 -45.1 -41.2 -39.4 -37.7 -35.1 590 1
Samples were drawn using NUTS(diag_e) at Fri May 5 10:41:43 2017.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).
De plus, comme dans BUGS (ou JAGS), CmdStan (l'interface en ligne de commande de Stan) requiert que
toutes les données soit dans un fichier dump de R. Dans le cas ou ce fichier existe, rstan fournit
la fonction read_rdump
permettant de lire toutes les données dans une liste de R. Par exemple, si
le fichier nommé 8scholls.rdump contenant les lignes suivantes est dans l'espace de travail (le working directory).
J <- 8
y <- c(28, 8, -3, 7, -1, 1, 18, 12)
sigma <- c(15, 10, 16, 11, 9, 11, 10, 18)
Alors, nous pouvons lire le fichier 8schools.rdump avec :
schools_dat <- read_rdump('8schools.rdump')
Le fichier R dump peut également être sourcé dans l'environnement global en utilisant la fonction source
de R.
Dans ce cas, nous pouvons omettre l'argument data
et stan
cherchera dans l'environnement pour les objets qui ont les mêmes noms
que ceux indiqués dans le bloc data de 8schools.stan. Par exemple:
source('8schools.rdump')
fit <- stan(file = '8schools.stan', iter = 1000, chains = 4)
L'exemple Rats est aussi un exemple assez classique. On trouve par exemple une version OpenBUGS à ce lien. Il a pour origine la publication de Gelfand et al (1990).
Les données font référence à la croissance hebdomadaire de 30 rats durant 5 semaines.
Dans le tableau suivant, nous avons listé les données, dans lesquelles x
correspond
aux dates de collecte des données. Vous pouvez essayer cet exemple en téléchargeant les données rats.txt
et le code de modèle rats.stan.
Rat | x=8 | x=15 | x=22 | x=29 | x=36 | Rat | x=8 | x=15 | x=22 | x=29 | x=36 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|
1 | 151 | 199 | 246 | 283 | 320 | 16 | 160 | 207 | 248 | 288 | 324 | |
2 | 145 | 199 | 249 | 293 | 354 | 17 | 142 | 187 | 234 | 280 | 316 | |
3 | 147 | 214 | 263 | 312 | 328 | 18 | 156 | 203 | 243 | 283 | 317 | |
4 | 155 | 200 | 237 | 272 | 297 | 19 | 157 | 212 | 259 | 307 | 336 | |
5 | 135 | 188 | 230 | 280 | 323 | 20 | 152 | 203 | 246 | 286 | 321 | |
6 | 159 | 210 | 252 | 298 | 331 | 21 | 154 | 205 | 253 | 298 | 334 | |
7 | 141 | 189 | 231 | 275 | 305 | 22 | 139 | 190 | 225 | 267 | 302 | |
8 | 159 | 201 | 248 | 297 | 338 | 23 | 146 | 191 | 229 | 272 | 302 | |
9 | 177 | 236 | 285 | 350 | 376 | 24 | 157 | 211 | 250 | 285 | 323 | |
10 | 134 | 182 | 220 | 260 | 296 | 25 | 132 | 185 | 237 | 286 | 331 | |
11 | 160 | 208 | 261 | 313 | 352 | 26 | 160 | 207 | 257 | 303 | 345 | |
12 | 143 | 188 | 220 | 273 | 314 | 27 | 169 | 216 | 261 | 295 | 333 | |
13 | 154 | 200 | 244 | 289 | 325 | 28 | 157 | 205 | 248 | 289 | 316 | |
14 | 171 | 221 | 270 | 326 | 358 | 29 | 137 | 180 | 219 | 258 | 291 | |
15 | 163 | 216 | 242 | 281 | 312 | 30 | 153 | 200 | 244 | 286 | 324 |
y <- as.matrix(read.table('https://raw.github.com/wiki/stan-dev/rstan/rats.txt', header = TRUE))
x <- c(8, 15, 22, 29, 36)
xbar <- mean(x)
N <- nrow(y)
T <- ncol(y)
rats_fit <- stan(file = 'https://raw.githubusercontent.com/stan-dev/example-models/master/bugs_examples/vol1/rats/rats.stan')
Vous pouvez essayer beaucoup des exemples de BUGS et bien d'autres qui ont été créé dans Stan en éxécutant
model <- stan_demo()
et choisir l'un des exemples de modèle à partir de la liste qui s'affiche. La première fois que vous utiliserez stan_demo()
, on vous demande
si vous souhaitez télécharger ces exemples. Vous devez choisir l'option 1 pour les mettre dans le dossier où rstan est installé afin qu'ils puissent
être utilisé dans le future sans avoir à les re-télécharger. L'objet model
ci-dessus est de la classe stanfit
, donc vous pouvez utiliser
les fonctions print
, plot
, pairs
, extract
, etc. pour l'analyser.
Plus de détails à propos de RStan peuvent être dans la documentation qui inclut la vignette du package rstan.
Par exemple, en utilisant help(stan)
et help("stanfit-class")
pour obtenir l'aide sur les fonctions de stan
et la classe S4 stanfit
.
Aussi, vous pouvez voir Stan's modeling language manual pour les détails sur l'échantillonneur Stan, et le langage de modélisation Stan.
Enfin, la Stan User's Mailing list peut être utilisée afin de discuter de l'utilisation de Stan, poster des exemples ou poser des questions sur (R)Stan. Si vous avez besoin d'aide, il est important de fournir assez d'information telles que:
- le code du modèle dans le langage Stan
- les données
- le code R si nécessaire
- les messages d'erreurs en utilisant
verbose=TRUE
etcores=1
lors de l'appel de la fonctionstan
- la version du compileur C++, par exemple en utilisant
g++ -v
(dans un terminal) sigcc
est utilisé - les informations sur R en utilisant
sessionInfo()
dans R
- Gelman, A., Carlin, J. B., Stern, H. S., and Rubin, D. B. (2003). Bayesian Data Analysis, CRC Press, London, 2nd Edition.
- Stan Development Team. Stan Modeling Language User's Guide and Reference Manual.
- Gelfand, A. E., Hills S. E., Racine-Poon, A., and Smith A. F. M. (1990). "Illustration of Bayesian Inference in Normal Data Models Using Gibbs Sampling", Journal of the American Statistical Association, 85, 972-985.
- Stan
- R
- BUGS
- OpenBUGS
- JAGS
- Rcpp