# Laboratorio 02 de Inferencia Estadística
# En myCompiler, selecciona R, pega el archivo en main.r y pulsa Ejecutar.
alpha <- 0.05

# P1(a). S es la suma de las n observaciones.
# La parte de la log-verosimilitud que depende de lambda es S*log(lambda) - n*lambda.
cat("P1(a) Derivada de la log-verosimilitud\n")
print(D(expression(S * log(lambda) - n * lambda), "lambda"))
cat("Al igualar a cero: lambda estimado = S/n = media muestral.\n")

# P1(b). Calculamos la media ponderada antes de agrupar la cola.
colonias <- c(0:11, 19)
frec <- c(56, 104, 80, 62, 42, 27, 9, 9, 5, 3, 2, 0, 1)
n <- sum(frec)
total_colonias <- sum(colonias * frec)
lambda <- weighted.mean(colonias, frec)
cat(sprintf("\nP1(b) n = %d, total de colonias = %d, lambda = %.3f\n",
            n, total_colonias, lambda))

# P1(c). Categorías 0, 1, ..., 6 y 7 o más.
# La cola incluye todos los enteros desde 7, incluido el caso de 19 colonias.
categorias <- c(as.character(0:6), "7+")
observadas <- c(frec[colonias < 7], sum(frec[colonias >= 7]))
prob <- c(dpois(0:6, lambda), ppois(6, lambda, lower.tail = FALSE))
esperadas <- n * prob
aportes <- (observadas - esperadas)^2 / esperadas
cat("\nP1(c) Frecuencias por categoría y aportes a chi-cuadrado\n")
print(data.frame(Categoria = categorias, O = observadas, p = round(prob, 8),
                 E = round(esperadas, 6), Aporte = round(aportes, 6)),
      row.names = FALSE)

# P1(d). Salida inicial de R. El argumento p contiene las probabilidades.
prueba1 <- chisq.test(x = observadas, p = prob)
cat("\nP1(d) Bondad de ajuste, salida inicial de R\n")
print(prueba1)

# P1(e). Restamos un grado de libertad por estimar lambda de la muestra.
gl1 <- length(observadas) - 1 - 1
x2_1 <- unname(prueba1$statistic)
p1 <- pchisq(x2_1, df = gl1, lower.tail = FALSE)
critico1 <- qchisq(1 - alpha, df = gl1)
cat(sprintf("\nP1(e) X2 = %.6f, gl = %d, crítico = %.6f, p corregido = %.8e\n",
            x2_1, gl1, critico1, p1))
if (p1 < alpha) cat("Se rechaza H0. El modelo Poisson no se ajusta a los conteos.\n")

# P2. Las filas son partidos y las columnas son respuestas de la encuesta.
tabla <- matrix(c(42, 309, 31, 64, 246, 46, 20, 115, 27),
               nrow = 3, byrow = TRUE,
               dimnames = list(Partido = c("Democrata", "Republicano", "Otro"),
                               Opinion = c("Se sabe todo", "Se ocultan hechos", "No sabe")))
cat("\nP2 Frecuencias observadas con sus totales\n")
print(addmargins(tabla))
prueba2 <- chisq.test(tabla)
cat("\nFrecuencias esperadas bajo independencia\n")
print(round(prueba2$expected, 4))
cat("\nAportes (O - E)^2 / E\n")
print(round((tabla - prueba2$expected)^2 / prueba2$expected, 6))
cat("\nPrueba de independencia de Pearson\n")
print(prueba2)
critico2 <- qchisq(1 - alpha, df = prueba2$parameter)
cat(sprintf("Valor crítico = %.6f, valor-p = %.9f\n", critico2, prueba2$p.value))
if (prueba2$p.value < alpha) cat("Se rechaza H0. La afiliación y la opinión están asociadas.\n")
cat("\nPorcentaje que responde que se ocultan hechos, dentro de cada partido\n")
print(round(100 * prop.table(tabla, margin = 1)[, "Se ocultan hechos"], 2))

cat("\nResumen al 5 %\n")
cat(sprintf("P1  lambda = %.3f\n", lambda))
cat(sprintf("    X2 = %.6f, gl = %d\n", x2_1, gl1))
cat(sprintf("    p corregido = %.8e\n", p1))
cat(sprintf("P2  X2 = %.6f, gl = %d\n", unname(prueba2$statistic), prueba2$parameter))
cat(sprintf("    p = %.9f\n", prueba2$p.value))