Gestión de la mendacidad encuestoelectoral: los números

probs <- c(4, 3, 2, 1)
probs <- probs / sum(probs)
partidos <- letters[1:length(probs)]

library(plyr)
library(rstan)
library(ggplot2)
library(reshape2)

n <- 3000
resultados <- data.frame(
  tarjeta = sample(1:nrow(tarjetas), n, replace = T),
  partido = sample(partidos, n, prob = probs, replace = T))
resultados <- data.frame(
  tarjetas[resultados$tarjeta,],
  partido = resultados$partido)
resultados$coincide <- resultados$partido == resultados$partido1 |
  resultados$partido == resultados$partido2

# proporciones reales en la muestra
props.muestra <- table(resultados$partido) / nrow(resultados)

# resultados agregados (por tarjeta)
resultados.agg <- ddply(
    resultados, .(partido1, partido2),
    summarize,
    total = length(partido1),
    coincidencias = sum(coincide))

codigo <- '
data {
  int<lower=1> N;
  int partido1[N];
  int partido2[N];
  int total[N];
  int coincidencias[N];
  int <lower = 1> n_partidos;
}

parameters {
  simplex[n_partidos] pes; //probabilidades
}

model {
  // dejamos la priori a "la cocina"

  // verosimilitud
  for (i in 1:N){
    coincidencias[i] ~ binomial(total[i],
      pes[partido1[i]] + pes[partido2[i]]);
  }
}
'

fit <- stan(model_code = codigo,
    data = list(n_partidos = length(partidos),
      N = nrow(resultados.agg),
      partido1 = match(resultados.agg$partido1, partidos),
      partido2 = match(resultados.agg$partido2, partidos),
      total = resultados.agg$total,
      coincidencias = resultados$coincidencias),
    iter=12000, warmup=2000,
    chains=2, thin=10)

res <- as.data.frame(fit)[,1:length(partidos)]
colnames(res) <- partidos
res <- melt(res)
colnames(res) <- c("partido", "pct")

ggplot(res, aes(x = pct, fill = partido)) +
  geom_density(alpha = 0.3) + ylab("")