#=========================================
# Laboratorio 02 - Inferencia Estadistica
# Integrantes:
# Valentina Herrera
# Cristian Araneda
# Matias Dimontti
# Davidson Pierre
# ========================================

# DATOS (se definen una sola vez)
colonias <- c(0:11, 19)
frec     <- c(56, 104, 80, 62, 42, 27, 9, 9, 5, 3, 2, 0, 1)
n        <- sum(frec)
alpha    <- 0.05

tabla_p2 <- matrix(c(42, 309, 31,
                     64, 246, 46,
                     20, 115, 27),
                   nrow = 3, byrow = TRUE)
rownames(tabla_p2) <- c("Democrata", "Republicano", "Otro")
colnames(tabla_p2) <- c("Sabemos todo", "Se ocultan hechos", "No sabe")

# P1(a): Estimador de maxima verosimilitud 
# L(lambda)  = prod( lambda^xi * exp(-lambda) / xi! )
# l(lambda)  = (sum xi) * log(lambda) - n*lambda - sum(log(xi!))
# l'(lambda) = (sum xi)/lambda - n = 0  =>  lambda_hat = sum(xi)/n = media muestral
# l''(lambda) = -(sum xi)/lambda^2 < 0  =>  es un maximo
cat("---- P1(a) ----\n")
cat("EMV: lambda_hat = media muestral (x barra)\n\n")

# P1(b): Calculo de lambda 
cat("---- P1(b) ----\n")
lambda_hat <- weighted.mean(colonias, frec)
cat("n =", n, "\n")
cat("lambda_hat =", format(round(lambda_hat, 3), nsmall = 3), "\n\n")

# P1(c): Categorias, probabilidades y esperados 
cat("---- P1(c) ----\n")
categorias <- c("0", "1", "2", "3", "4", "5", "6", ">=7")
obs        <- c(frec[colonias <= 6], sum(frec[colonias >= 7]))
probs      <- dpois(0:6, lambda_hat)
probs      <- c(probs, 1 - sum(probs))   # ultima categoria por complemento
esperado   <- n * probs

tabla_p1 <- data.frame(Categoria    = categorias,
                       Observado    = obs,
                       Probabilidad = round(probs, 4),
                       Esperado     = round(esperado, 2))
print(tabla_p1)

if (any(esperado < 5)) {
  cat("AVISO: hay categorias con esperado < 5\n")
} else {
  cat("OK: todos los esperados son >= 5\n")
}
cat("\n")

# P1(d): Test de bondad de ajuste (gl de R) 
cat("---- P1(d) ----\n")
cat("H0: el numero de colonias distribuye Poisson\n")
cat("H1: el numero de colonias NO distribuye Poisson\n")
test_gof    <- chisq.test(obs, p = probs)
estadistico <- unname(test_gof$statistic)
cat("Estadistico =", round(estadistico, 3), "\n")
cat("gl (segun R) =", unname(test_gof$parameter), "\n")
cat("valor-p (segun R) =", format(test_gof$p.value, scientific = TRUE, digits = 4), "\n")
if (test_gof$p.value < alpha) {
  cat("Conclusion: se rechaza H0 con alpha = 0.05; los datos no siguen una Poisson.\n\n")
} else {
  cat("Conclusion: no se rechaza H0 con alpha = 0.05.\n\n")
}
# En simple: el modelo Poisson no describe bien los datos
# (subestima los ceros y la cola larga; hay sobredispersion).

# P1(e): Grados de libertad correctos 
cat("---- P1(e) ----\n")
k     <- length(obs)
m     <- 1                      # parametros estimados (lambda)
gl_ok <- k - 1 - m
p_ok  <- pchisq(estadistico, df = gl_ok, lower.tail = FALSE)
crit  <- qchisq(0.95, gl_ok)
cat("gl correctos = k - 1 - m =", k, "- 1 -", m, "=", gl_ok, "\n")
cat("valor-p correcto =", format(p_ok, scientific = TRUE, digits = 4), "\n")
cat("valor-p de R     =", format(test_gof$p.value, scientific = TRUE, digits = 4), "\n")
cat("Valor critico (alpha = 0.05, gl =", gl_ok, ") =", round(crit, 3), "\n")
cat("Estadistico > valor critico:", estadistico > crit, "\n")
cat("Por que R falla: chisq.test usa gl = k - 1 y no sabe que lambda se estimo\n")
cat("con los mismos datos; cada parametro estimado quita un grado de libertad.\n")
cat("Conclusion: la decision no cambia, se rechaza H0.\n\n")

# P2: Test de independencia 
cat("---- P2 ----\n")
cat("H0: afiliacion politica y opinion son independientes\n")
cat("H1: afiliacion politica y opinion son dependientes\n")
test_ind <- chisq.test(tabla_p2)
cat("Estadistico =", round(unname(test_ind$statistic), 3), "\n")
cat("gl =", unname(test_ind$parameter), "\n")
cat("valor-p =", format(test_ind$p.value, scientific = TRUE, digits = 4), "\n\n")
cat("Tabla de esperados:\n")
print(round(test_ind$expected, 2))
cat("\nResiduos estandarizados (|valor| > 2 indica celdas influyentes):\n")
print(round(test_ind$stdres, 2))
if (all(test_ind$expected >= 5)) cat("\nOK: todos los esperados son >= 5\n")
if (test_ind$p.value < alpha) {
  cat("Conclusion: se rechaza H0; hay evidencia de dependencia entre partido y opinion.\n")
}
cat("Nota: es asociacion, no causalidad.\n")
# En simple: la opinion sobre el encubrimiento cambia segun el partido.