Va de si hay una o dos lambdas

r

Un año, el 2016, mueren 1160 personas en accidentes de tráfico. El anterior, 1131, i.e., 29 menos. Ruido estadístico aparte, ¿aumentan?

Comenzamos a optar. Primera elección subjetiva: son muestras de una Poisson de parámetro desconocido. La pregunta: ¿el mismo?

Una manera de estudiar lo anterior es plantear

1160 ~ poisson(lambda * (1 + incr))
1131 ~ poisson(lambda)

y estudiar la distribución de incr. Que a saber qué distribución tendrá (teóricamente). Pero, ¿importa?

Mejor que rebuscar a ver qué distribución podría tener la cosa, basta con envolverlo en un poco de seudo-C++,

stancode <- '
  parameters{
    real<lower = 0> lambda;
    real<lower = -1, upper=1> incr;
  }

  model{
    1160 ~ poisson(lambda * (1 + incr));
    1131 ~ poisson(lambda);
  }
'

y luego, dejar que la magia surta efecto:

library(rstan)
fit <- stan(model_code = stancode, chains = 1,
            iter = 12000, warmup = 2000, thin = 10)

summary(fit)$summary[1:2, c("mean", "2.5%", "97.5%")]
#               mean         2.5%       97.5%
#lambda 1.131041e+03 1065.3705780 1200.341736
#incr   2.745496e-02   -0.0531976    0.116086

Y se podrían pintar los histogramas, etc. También podemos de alguna manera estimar la probabilidad de que la accidentalidad haya subido, si es que la pregunta tiene sentido:

tmp <- as.data.frame(fit)
mean(tmp$incr > 0)
#[1] 0.732

Es decir, podíamos decir que hay un 25% de probabilidades de que los neocríticos de la DGT estén equivocados.