data("iris")
View(iris)
# Cargar librería
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.4.3
# Cargar el conjunto de datos iris
data(iris)
# Gráfico bivariado
p<-ggplot(iris, aes(x = Sepal.Width, y = Petal.Width, color = Species)) +
geom_point(size = 2, alpha = 0.8) + # puntos
geom_smooth(method = "lm", se = FALSE, linetype = "dashed", size = 1) + # línea punteada
labs(
title = "Relación entre Sepal.Width y Petal.Width por Especie",
x = "Ancho del Sépalo (cm)",
y = "Ancho del Pétalo (cm)",
color = "Especie"
) +
theme_minimal(base_size = 14)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
ggsave("grafico_iris.png", plot = p, width = 7, height = 5, dpi = 300)
## `geom_smooth()` using formula = 'y ~ x'
#B
## 0) Datos y preparación -------------------
data(iris)
X <- as.matrix(iris[, c("Sepal.Width", "Petal.Width")]) # p = 2
g <- iris$Species
groups <- levels(g)
N <- nrow(X) # total
p <- ncol(X) # número de variables (2)
k <- length(groups) # número de grupos (3)
## 1) Medias por grupo y global -------------
n_i <- tapply(g, g, length)
mu_i <- lapply(groups, function(cl) colMeans(X[g==cl, , drop=FALSE]))
mu_i <- do.call(cbind, mu_i) # p x k
colnames(mu_i) <- groups
mu <- colMeans(X) # media global (p)
mu_i
## setosa versicolor virginica
## Sepal.Width 3.428 2.770 2.974
## Petal.Width 0.246 1.326 2.026
## 2) Matrices de covarianza por grupo ------
# S_i con divisor (n_i - 1)
S_i <- lapply(groups, function(cl) {
Xi <- X[g==cl, , drop=FALSE]
sweep(Xi, 2, colMeans(Xi)) # centra (no guarda)
cov(Xi) # por defecto usa (n_i - 1)
})
names(S_i) <- groups
## 3) SSCP Within (W), Between (B), Total (T)
# W = sum_i (n_i - 1) S_i
W <- Reduce("+", Map(function(S, ni) (ni - 1) * S, S_i, as.list(n_i)))
# B = sum_i n_i (mu_i - mu)(mu_i - mu)' (usando medias)
B <- matrix(0, p, p)
for (j in seq_len(k)) {
d <- mu_i[, j] - mu
B <- B + n_i[j] * (d %*% t(d))
}
# T = W + B (y también T = sum sobre todo (x - mu)(x - mu)')
Tmat <- W + B
## 4) Wilks' Lambda y prueba (Bartlett) -----
Lambda <- det(W) / det(Tmat)
# Estadístico aproximado Chi-cuadrado:
# -[ (N - 1) - (p + k)/2 ] * ln(Lambda) ~ Chi2_{p(k-1)}
m <- (N - 1) - (p + k)/2
X2 <- - m * log(Lambda)
df <- p * (k - 1)
pval_wilks <- 1 - pchisq(X2, df = df)
## 5) Covarianzas iguales: Box's M "a mano" --
# Sp (covarianza agrupada) con divisor (N - k)
Sp <- W / (N - k)
# M de Box:
# M = (N - k) * ln|Sp| - sum_i (n_i - 1) * ln|S_i|
term1 <- (N - k) * log(det(Sp))
term2 <- sum(mapply(function(S, ni) (ni - 1) * log(det(S)), S_i, as.list(n_i)))
Mbox <- term1 - term2
# Corrección c (Bartlett) y aproximación Chi-cuadrado:
# c = [ (2p^2 + 3p - 1) * ( sum 1/(n_i - 1) - 1/(N - k) ) ] / [ 6(p+1)(k-1) ]
sum_inv <- sum(1 / (n_i - 1)) - 1 / (N - k)
c_corr <- ((2 * p^2 + 3 * p - 1) * sum_inv) / (6 * (p + 1) * (k - 1))
X2_box <- (1 - c_corr) * Mbox
df_box <- (k - 1) * p * (p + 1) / 2
pval_box <- 1 - pchisq(X2_box, df = df_box)
## 6) Salidas legibles -----------------------
cat("\n=== Tamaños por grupo ===\n"); print(n_i)
##
## === Tamaños por grupo ===
## setosa versicolor virginica
## 50 50 50
cat("\nMedia global mu:\n"); print(mu)
##
## Media global mu:
## Sepal.Width Petal.Width
## 3.057333 1.199333
cat("\nMedias por grupo (columnas):\n"); print(mu_i)
##
## Medias por grupo (columnas):
## setosa versicolor virginica
## Sepal.Width 3.428 2.770 2.974
## Petal.Width 0.246 1.326 2.026
cat("\nMatriz W (Within SSCP):\n"); print(W)
##
## Matriz W (Within SSCP):
## Sepal.Width Petal.Width
## Sepal.Width 16.9620 4.8084
## Petal.Width 4.8084 6.1566
cat("\nMatriz B (Between SSCP):\n"); print(B)
##
## Matriz B (Between SSCP):
## Sepal.Width Petal.Width
## [1,] 11.34493 -22.93267
## [2,] -22.93267 80.41333
cat("\nMatriz T = W + B (Total SSCP):\n"); print(Tmat)
##
## Matriz T = W + B (Total SSCP):
## Sepal.Width Petal.Width
## Sepal.Width 28.30693 -18.12427
## Petal.Width -18.12427 86.56993
cat("\n--- Wilks' Lambda ---\n")
##
## --- Wilks' Lambda ---
cat(sprintf("Lambda = %.6f\nChi2 approx = %.4f, df = %d, p-valor = %.6g\n",
Lambda, X2, df, pval_wilks))
## Lambda = 0.038316
## Chi2 approx = 477.8676, df = 4, p-valor = 0
cat("\n--- Box's M (igualdad de covarianzas) ---\n")
##
## --- Box's M (igualdad de covarianzas) ---
cat(sprintf("M = %.6f, corrección c = %.6f\nChi2 approx = %.4f, df = %d, p-valor = %.6g\n",
Mbox, c_corr, X2_box, df_box, pval_box))
## M = 52.832419, corrección c = 0.019652
## Chi2 approx = 51.7941, df = 6, p-valor = 2.05128e-09
## ===========================
## QDA "a mano" en R base
## iris: Sepal.Width, Petal.Width
## ===========================
## 1) Datos
data(iris)
X <- as.matrix(iris[, c("Sepal.Width", "Petal.Width")]) # p = 2
y <- iris$Species
groups <- levels(y)
k <- length(groups)
p <- ncol(X)
## 2) Particiones por grupo, tamaños, medias y covarianzas (sin paquetes)
split_idx <- split(seq_len(nrow(X)), y)
n_i <- sapply(split_idx, length)
mu_i <- sapply(groups, function(g) colMeans(X[ y==g, , drop = FALSE ]))
colnames(mu_i) <- groups
S_i <- lapply(groups, function(g){
Xi <- X[ y==g, , drop = FALSE ]
cov(Xi) # divisor (n_i - 1)
})
names(S_i) <- groups
## 3) Densidad normal bivariada y puntuación QDA
dmvnorm2 <- function(x, mu, Sigma) {
# x y mu vectores de longitud p; Sigma matriz p x p
x <- as.numeric(x); mu <- as.numeric(mu)
p <- length(mu)
S_inv <- solve(Sigma)
quad <- (x - mu)
quad <- as.numeric(t(quad) %*% S_inv %*% quad)
cte <- 1 / ((2*pi)^(p/2) * sqrt(det(Sigma)))
cte * exp(-0.5 * quad)
}
dQ <- function(x, mu, Sigma, logpi = 0) {
# d_i^Q(x) = log p_i - 0.5 log|Sigma| - 0.5 (x-mu)'Sigma^{-1}(x-mu)
x <- as.numeric(x); mu <- as.numeric(mu)
-0.5 * log(det(Sigma)) - 0.5 * t(x - mu) %*% solve(Sigma) %*% (x - mu) + logpi
}
## 4) Priors iguales y punto a clasificar
priors <- rep(1/k, k)
logpi <- log(priors)
x0 <- c(3.5, 1.75)
## 5) Puntuaciones QDA, posteriors y clase predicha
scores <- sapply(seq_len(k), function(j) dQ(x0, mu_i[, j], S_i[[j]], logpi[j]))
names(scores) <- groups
# Posteriores (Bayes): p_i f_i(x0) / sum_j p_j f_j(x0)
num <- sapply(seq_len(k), function(j) priors[j] * dmvnorm2(x0, mu_i[, j], S_i[[j]]))
post <- num / sum(num); names(post) <- groups
pred <- names(which.max(scores))
## 6) Salidas
cat("\nMedias por grupo (columnas):\n"); print(round(mu_i, 3))
##
## Medias por grupo (columnas):
## setosa versicolor virginica
## Sepal.Width 3.428 2.770 2.974
## Petal.Width 0.246 1.326 2.026
cat("\nCovarianzas por grupo:\n")
##
## Covarianzas por grupo:
for (g in groups) { cat("\n", g, ":\n", sep=""); print(round(S_i[[g]], 4)) }
##
## setosa:
## Sepal.Width Petal.Width
## Sepal.Width 0.1437 0.0093
## Petal.Width 0.0093 0.0111
##
## versicolor:
## Sepal.Width Petal.Width
## Sepal.Width 0.0985 0.0412
## Petal.Width 0.0412 0.0391
##
## virginica:
## Sepal.Width Petal.Width
## Sepal.Width 0.1040 0.0476
## Petal.Width 0.0476 0.0754
cat("\nPuntuaciones QDA d_i^Q(x0):\n"); print(round(scores, 6))
##
## Puntuaciones QDA d_i^Q(x0):
## setosa versicolor virginica
## -104.871856 -1.055387 -2.325402
cat("\nPosteriores Bayes P(pi_i | x0):\n"); print(round(post, 6))
##
## Posteriores Bayes P(pi_i | x0):
## setosa versicolor virginica
## 0.000000 0.780745 0.219255
cat("\nClasificación de x0 = [3.5, 1.75]: ", pred, "\n", sep = "")
##
## Clasificación de x0 = [3.5, 1.75]: versicolor
## LDA "a mano" con covarianza común (iris: Sepal.Width, Petal.Width)
data(iris)
X <- as.matrix(iris[, c("Sepal.Width", "Petal.Width")])
g <- iris$Species
groups <- levels(g); k <- length(groups); p <- ncol(X)
## Medias y covarianzas por grupo
mu_i <- sapply(groups, function(cl) colMeans(X[g==cl, , drop=FALSE]))
S_i <- lapply(groups, function(cl) cov(X[g==cl, , drop=FALSE]))
n_i <- sapply(groups, function(cl) sum(g==cl))
N <- nrow(X)
## Covarianza agrupada Sp = W/(N-k), con W = sum (n_i-1) S_i
W <- Reduce("+", Map(function(S,ni) (ni-1)*S, S_i, as.list(n_i)))
Sp <- W/(N - k)
Sinv <- solve(Sp)
## Priors iguales
priors <- rep(1/k, k); logpi <- log(priors)
## Puntuación lineal de LDA
delta <- function(x, mu) { as.numeric(t(x) %*% Sinv %*% mu - 0.5 * t(mu) %*% Sinv %*% mu) }
x0 <- c(3.5, 1.75)
scores_lda <- sapply(1:k, function(j) delta(x0, mu_i[, j]) + logpi[j])
names(scores_lda) <- groups
## Posteriores con covarianza común (densidad normal con Sp)
dmvnorm_common <- function(x, mu, S) {
d <- x - mu
detS <- det(S); invS <- solve(S)
(1/((2*pi)^(p/2)*sqrt(detS))) * exp(-0.5 * t(d) %*% invS %*% d)
}
dens <- sapply(1:k, function(j) dmvnorm_common(x0, mu_i[, j], Sp))
post_lda <- (priors * dens) / sum(priors * dens); names(post_lda) <- groups
pred_lda <- names(which.max(scores_lda))
## Mostrar
cat("\nSp (covarianza común):\n"); print(round(Sp, 5))
##
## Sp (covarianza común):
## Sepal.Width Petal.Width
## Sepal.Width 0.11539 0.03271
## Petal.Width 0.03271 0.04188
cat("\nPuntuaciones LDA delta_i(x0):\n"); print(round(scores_lda, 6))
##
## Puntuaciones LDA delta_i(x0):
## setosa versicolor virginica
## 27.01748 57.75737 56.81905
cat("\nPosteriores LDA P(pi_i | x0):\n"); print(round(post_lda, 6))
##
## Posteriores LDA P(pi_i | x0):
## setosa versicolor virginica
## 0.000000 0.718759 0.281241
cat("\nClasificación LDA de x0 = [3.5, 1.75]: ", pred_lda, "\n", sep = "")
##
## Clasificación LDA de x0 = [3.5, 1.75]: versicolor
data(iris)
X <- as.matrix(iris[, c("Sepal.Width", "Petal.Width")])
g <- iris$Species
groups <- levels(g); k <- length(groups); p <- ncol(X)
# 2) Medias por grupo y covarianzas
mu_i <- sapply(groups, function(cl) colMeans(X[g==cl, , drop=FALSE]))
colnames(mu_i) <- groups
n_i <- sapply(groups, function(cl) sum(g==cl))
S_i <- lapply(groups, function(cl) cov(X[g==cl, , drop=FALSE]))
# 3) Covarianza común agrupada Sp = W/(N-k), con W = sum (n_i-1) S_i
N <- nrow(X)
W <- Reduce("+", Map(function(S,ni) (ni-1)*S, S_i, as.list(n_i)))
Sp <- W/(N - k)
Sinv <- solve(Sp)
# 4) Priors iguales y función lineal d_i(x) (Regla 11.56)
priors <- rep(1/k, k)
logpi <- log(priors)
d_lin <- function(x, mu, Sinv, logpi_i=0) {
# d_i(x) = x' Sinv mu_i - 0.5 mu_i' Sinv mu_i + log p_i
as.numeric(t(x) %*% Sinv %*% mu - 0.5 * t(mu) %*% Sinv %*% mu) + logpi_i
}
# 5) Clasificación de x0
x0 <- c(3.5, 1.75)
scores_lda <- sapply(1:k, function(j) d_lin(x0, mu_i[, j], Sinv, logpi[j]))
names(scores_lda) <- groups
post_unnorm <- sapply(1:k, function(j) {
# Posteriores bajo LDA usando densidad con Sp (opcional, para reporte)
d <- x0 - mu_i[, j]
detSp <- det(Sp)
(priors[j] / ((2*pi)^(p/2) * sqrt(detSp))) * exp(-0.5 * t(d) %*% Sinv %*% d)
})
post_lda <- as.numeric(post_unnorm / sum(post_unnorm)); names(post_lda) <- groups
pred_lda <- names(which.max(scores_lda))
cat("\n=== Regla (11.56): LDA con covarianza común ===\n")
##
## === Regla (11.56): LDA con covarianza común ===
cat("Sp (covarianza común):\n"); print(round(Sp, 5))
## Sp (covarianza común):
## Sepal.Width Petal.Width
## Sepal.Width 0.11539 0.03271
## Petal.Width 0.03271 0.04188
cat("\nPuntuaciones lineales d_i(x0):\n"); print(round(scores_lda, 6))
##
## Puntuaciones lineales d_i(x0):
## setosa versicolor virginica
## 27.01748 57.75737 56.81905
cat("\nPosteriores (LDA) P(pi_i | x0):\n"); print(round(post_lda, 6))
##
## Posteriores (LDA) P(pi_i | x0):
## setosa versicolor virginica
## 0.000000 0.718759 0.281241
cat(sprintf("\nClasificación de x0 = [%.2f, %.2f]: %s\n", x0[1], x0[2], pred_lda))
##
## Clasificación de x0 = [3.50, 1.75]: versicolor
x1_rng <- range(X[,1]); x2_rng <- range(X[,2])
pad1 <- diff(x1_rng)*0.1; pad2 <- diff(x2_rng)*0.1
x1_seq <- seq(x1_rng[1]-pad1, x1_rng[2]+pad1, length.out = 400)
x2_seq <- seq(x2_rng[1]-pad2, x2_rng[2]+pad2, length.out = 400)
grid <- expand.grid(Sepal.Width = x1_seq, Petal.Width = x2_seq)
## --- REEMPLAZO COMPLETO DEL BLOQUE "desde aquí no corre" ---
## 1) Grilla numérica
x1_rng <- range(X[,1]); x2_rng <- range(X[,2])
pad1 <- diff(x1_rng)*0.10; pad2 <- diff(x2_rng)*0.10
x1_seq <- seq(x1_rng[1]-pad1, x1_rng[2]+pad1, length.out = 400)
x2_seq <- seq(x2_rng[1]-pad2, x2_rng[2]+pad2, length.out = 400)
grid_df <- expand.grid(Sepal.Width = x1_seq, Petal.Width = x2_seq)
grid_mat <- as.matrix(grid_df) # 160000 x 2 (numérico)
## 2) Pre-cálculos de LDA (si no los tienes aún)
# V = Sinv %*% mu_i (2 x 3); c_i = -0.5 * mu_i' Sinv mu_i + log p_i
V <- Sinv %*% mu_i
quad <- colSums(mu_i * (Sinv %*% mu_i))
cvec <- as.numeric(-0.5 * quad + logpi) # longitud 3, una por clase
names(cvec) <- groups
## 3) Puntuaciones d_i(x) para TODA la grilla (sin apply)
# Dmat = grid_mat %*% V -> (Ngrid x 3); luego sumamos constantes por columna
Dmat <- grid_mat %*% V
Dmat <- sweep(Dmat, 2, cvec, FUN = "+")
## 4) Guardar d_i en columnas y obtener predicción por argmax
colnames(Dmat) <- paste0("d_", groups)
grid_df <- cbind(grid_df, as.data.frame(Dmat))
imax <- max.col(Dmat, ties.method = "first")
grid_df$pred <- factor(groups[imax], levels = groups)
## 5) Fronteras: d_i - d_j = 0
grid_df$diff_12 <- grid_df[[paste0("d_", groups[1])]] - grid_df[[paste0("d_", groups[2])]]
grid_df$diff_13 <- grid_df[[paste0("d_", groups[1])]] - grid_df[[paste0("d_", groups[3])]]
grid_df$diff_23 <- grid_df[[paste0("d_", groups[2])]] - grid_df[[paste0("d_", groups[3])]]
## 6) (Opcional) Graficar regiones + fronteras
library(ggplot2)
p2 <- ggplot() +
geom_raster(data = grid_df, aes(Sepal.Width, Petal.Width, fill = pred), alpha = 0.28) +
geom_point(data = iris, aes(Sepal.Width, Petal.Width, color = Species), size = 1.8, alpha = 0.9) +
geom_contour(data = grid_df, aes(Sepal.Width, Petal.Width, z = diff_12), breaks = 0, linewidth = 0.6) +
geom_contour(data = grid_df, aes(Sepal.Width, Petal.Width, z = diff_13), breaks = 0, linewidth = 0.6) +
geom_contour(data = grid_df, aes(Sepal.Width, Petal.Width, z = diff_23), breaks = 0, linewidth = 0.6) +
labs(title = "Regiones de clasificación LDA (Regla 11.56)",
x = "Sepal.Width", y = "Petal.Width", fill = "Región", color = "Especie") +
theme_minimal(base_size = 13)
ggsave("grafico_iris.png", plot = p2, width = 7, height = 5, dpi = 300)
## 0) Datos
data(iris)
X <- as.matrix(iris[, c("Sepal.Width", "Petal.Width")]) # p = 2
g <- iris$Species
groups <- levels(g)
k <- length(groups); p <- ncol(X); N <- nrow(X)
## 1) Utilidades: estimadores y puntuaciones LDA -----------------------------
# Medias por grupo (2 x k) y tamaños por grupo
get_mu_n <- function(X, g, groups) {
mu <- sapply(groups, function(cl) colMeans(X[g==cl, , drop=FALSE]))
n <- sapply(groups, function(cl) sum(g==cl))
list(mu=mu, n=n)
}
# W (SSCP within) = sum_i (n_i - 1) S_i
get_W <- function(X, g, groups) {
Reduce("+", lapply(groups, function(cl){
Xi <- X[g==cl, , drop=FALSE]
(nrow(Xi) - 1) * cov(Xi)
}))
}
# Covarianza común agrupada: Sp = W / (N - k)
get_Sp <- function(W, N, k) W / (N - k)
# Puntuación lineal LDA: d_i(x) = x' Sinv mu_i - 0.5 mu_i' Sinv mu_i + log p_i
# Versión vectorizada para TODAS las filas de X (N x p)
lda_scores <- function(X, mu, Sinv, logpi) {
# X: N x p; mu: p x k; Sinv: p x p; logpi: long k
V <- Sinv %*% mu # p x k
base <- -0.5 * colSums(mu * (Sinv %*% mu)) + logpi # long k
D <- X %*% V # N x k
sweep(D, 2, base, "+") # N x k
}
# Clasifica devolviendo factor de etiquetas predichas
predict_classes <- function(score_mat, groups){
groups[max.col(score_mat, ties.method="first")]
}
## 2) Priors iguales y objetos de Parte (d) ---------------------------------
priors <- rep(1/k, k); logpi <- log(priors)
mn <- get_mu_n(X, g, groups); mu <- mn$mu; n_i <- mn$n
W <- get_W(X, g, groups)
Sp <- get_Sp(W, N, k)
Sinv <- solve(Sp)
## 3) Clasificación por re-sustitución (APER) -------------------------------
scores_resub <- lda_scores(X, mu, Sinv, logpi) # N x k
yhat_resub <- predict_classes(scores_resub, groups)
# Matriz de confusión y APER
tab_resub <- table(Real=g, Pred=yhat_resub)
APER <- 1 - sum(diag(tab_resub)) / N
cat("\n=== Re-sustitución (APER) ===\n")
##
## === Re-sustitución (APER) ===
print(tab_resub)
## Pred
## Real setosa versicolor virginica
## setosa 50 0 0
## versicolor 0 49 1
## virginica 0 4 46
cat(sprintf("APER = %.6f\n", APER))
## APER = 0.033333
## 4) Lachenbruch (LOOCV) ---------------------------------------------------
# Para cada i: eliminar i, recalcular mu, W, Sp, Sinv y clasificar x_i
yhat_loo <- character(N)
for (i in 1:N) {
# conjunto de entrenamiento sin i
idx <- (1:N) != i
Xi <- X[idx, , drop=FALSE]
gi <- g[idx]
Ni <- nrow(Xi)
# Si alguna clase quedara vacía se debería saltar, pero en iris no ocurre
mn_i <- get_mu_n(Xi, gi, groups); mu_i <- mn_i$mu; n_i2 <- mn_i$n
Wi <- get_W(Xi, gi, groups)
Spi <- get_Sp(Wi, Ni, k)
Sinvi <- solve(Spi)
scores_i <- lda_scores(matrix(X[i,], nrow=1), mu_i, Sinvi, logpi) # 1 x k
yhat_loo[i] <- predict_classes(scores_i, groups)
}
yhat_loo <- factor(yhat_loo, levels=groups)
tab_loo <- table(Real=g, Pred=yhat_loo)
Ehat_AER <- 1 - sum(diag(tab_loo)) / N
cat("\n=== Lachenbruch (LOOCV) ===\n")
##
## === Lachenbruch (LOOCV) ===
print(tab_loo)
## Pred
## Real setosa versicolor virginica
## setosa 50 0 0
## versicolor 0 48 2
## virginica 0 4 46
cat(sprintf("Ê(AER) por Lachenbruch = %.6f\n", Ehat_AER))
## Ê(AER) por Lachenbruch = 0.040000
## 5) Resumen breve ----------------------------------------------------------
cat("\nResumen:\n")
##
## Resumen:
cat(sprintf("APER (resub) = %.6f\n", APER))
## APER (resub) = 0.033333
cat(sprintf("Ê(AER) Lachenbruch (LOO)= %.6f\n", Ehat_AER))
## Ê(AER) Lachenbruch (LOO)= 0.040000
p
## [1] 2
p2
