Fig_08

# ¿BINOMIAL NEGATIVA?

#---- VARIABLES
espacio <- 2500             # espacio total
m <- 8                      # número de grupos
n <- 50                     # número de individuos/grupo
radio <- 40                 # radio del grupo
nH <- 1000                  # número de hospedadores
veces <- 1000               # veces que se repite la simulación
picax <- 0:4                # nº posible de picaduras/hospedador
clases <- length(picax)
picam <- numeric(clases)         # tabla media picaduras

#---- BUCLE DE LAS VECES
for (v in 1:veces){
  S <- integer(espacio)       # posiciones hospedadores
  Si <- S                     # posiciones hospedadores infestados

  #---- REPARTO DE HOSPEDADORES
  S[sample(espacio, nH, replace = TRUE)] <- 1
  nHreal <- sum(S)

  #---- REPARTO DE PARÁSITOS
  mx <- sample(espacio, m, replace = TRUE)    # centro gravedad de cada grupo
  for (i in 1:m) {            # poblar los m grupos
    for (j in 1:n) {
      nx <- sample((mx[i] - radio):(mx[i] + radio), 1)
      nx <- nx %% espacio
      if (nx == 0) {nx <- 1}
      if (S[nx] == 1) { # ha pillado a un hospedador
        Si[nx] <- Si[nx] + 1     # lo marcamos como infestado
      }
    }
  }

  #---- DISTRIBUCIÓN DE PICADURAS
  pica <- as.numeric(table(factor(Si, levels = picax)))
  pica[1] <- nHreal - sum(pica[2:clases])
  picam <- picam + pica

} # fin bucle veces

picam <- picam / veces

#---- RESULTADO GRÁFICO
pdf("fig_08.pdf", width = 6, height = 6)
par(bty = "n", las = 1)
plot(picax, log10(picam), ylim = c(0, log10(nH)), xlim = c(-0.3, clases - 0.5), type = "h", xlab = "Nº ataques / hospedador", ylab = "Nº hospedadores", col = "black", lwd = 4, yaxt = "n")
etiquetas <- c(1, 10, 25, 50, 100, 250, 500, 1000)
axis(2, at = log10(etiquetas), labels = etiquetas)
text(picax, log10(picam), round(picam), cex = 0.8, pos = 3, col = "black")
nHreal <- sum(picam)
media <- sum(picax*picam)/nHreal
varia <- sum(picam*(picax-media)^2)/(nHreal - 1)
poi <- nHreal*dpois(picax, media)
points(picax - 0.1, log10(poi), type = "h", col = "cadetblue1", lwd = 4)
text(picax - .1, log10(poi), round(poi), cex = 0.8, pos = 2, col = "darkgrey")

#---- AJUSTE BINOMIAL
SDC2 <- function(par) { # mínimos cuadrados
  ft <- nHreal*dnbinom(picax, size = par, mu = media)
  SDCn = sum((picam - ft)^2)/(var(picam)*clases)
  return(SDCn)
}
SDC <- function(par) { # mínimo chi2
  ft <- nHreal*dnbinom(picax, size = par, mu = media)
  SDCn = sum(((picam - ft)^2)/ft)
  return(SDCn)
}
ajusteBN <- optim(par = c(0.18), SDC)
k <- ajusteBN$par
ft <- (1 + media/k)^(-k)
pbn <- nHreal*dnbinom(picax, size = k, mu = media)
points(picax + .1, log10(pbn), type = "h", col = "royalblue", lwd = 4)
text(picax + .1, log10(pbn), round(pbn), cex = 0.8, pos = 4, col = "royalblue")
legend("topright", c("Experimento", "Binomial negativa", "Poisson"), lty = 1, lwd = 2, col = c("black", "royalblue", "cadetblue1"), cex = 0.8, bty = "n")
dev.off()