Fig_09

#---- DATOS DE M. HASSELL Y ANÁLISIS DE R. MAY
#---- con Cyzenis albicans
picax <- c(0, 1, 2, 3, 4)
picam <- c(1066, 176, 48, 8, 5)

# Gráfico
pdf("fig_09.pdf", width = 6, height = 6)
par(bty = "n", las = 1)
plot(picax, log10(picam), ylim = c(0, 3.5), xlim = c(-0.4, 4.5), type = "h", xlab = "Nº larvas / polilla", ylab = "Nº polillas", 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 de la binomial (mínimo chi2) y gráfico final
SDC <- function(par) {
  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("Datos de Hassell", "Binomial negativa", "Poisson"), lty = 1, lwd = 2, col = c("black", "royalblue", "cadetblue1"), cex = 0.8, bty = "n")

dev.off()