RStan Getting Started (Français)

RStan est une interface R à Stan. Pour plus d'informations sur Stan, et son langage de modélisation, vous pouvez visister le site web :

mc-stan.org

Dernière Version : 2.16.2 (19 Juillet 2017)

La quasi totalité des instructions ci-dessoute sont pour la version de RStan mentionnée ci-dessus.

Installation

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):

Stan users mailing list

Comment utiliser RStan

Le reste de ce document suppose que vous avez intallé RStan en suivant les instructions d'un des liens ci-dessus.

Chargement du package (ou librairie)

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.

Exemple 1: Eight Schools

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 appelerons 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)

Exemple 2: Rats

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')

Exemple 3: autre

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 d'aide

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 et cores=1 lors de l'appel de la fonction stan
la version du compileur C++, par exemple en utilisant g++ -v (dans un terminal) si gcc est utilisé
les informations sur R en utilisant sessionInfo() dans R

Références

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

Provide feedback

Saved searches

Use saved searches to filter your results more quickly