# =============================================================================
# 1. FUNGSI
# =============================================================================
# --- kernel dan KDE ---
kernel_gaussian <- function(r2, d) (2 * pi)^(-d / 2) * exp(-0.5 * r2)
bobot_kernel <- function(r2) exp(-0.5 * r2)
jarak2 <- function(A, B) {
pmax(outer(rowSums(A^2), rowSums(B^2), "+") - 2 * A %*% t(B), 0)
}
bandwidth_kuantil <- function(X, q = 0.3) as.numeric(quantile(dist(X), q))
kde <- function(x, X, h) {
n <- nrow(X); d <- ncol(X)
r2 <- colSums((t(X) - x)^2) / h^2
sum(kernel_gaussian(r2, d)) / (n * h^d)
}
kde_grid <- function(X, gx, gy, h) {
n <- nrow(X); d <- ncol(X)
sisa <- if (d > 2) rowSums(X[, -(1:2), drop = FALSE]^2) else rep(0, n)
Kx <- exp(-0.5 * outer(gx, X[, 1], "-")^2 / h^2)
Ky <- exp(-0.5 * outer(gy, X[, 2], "-")^2 / h^2)
Kx <- sweep(Kx, 2, exp(-0.5 * sisa / h^2), "*")
(Kx %*% t(Ky)) * (2 * pi)^(-d / 2) / (n * h^d)
}
# --- mean shift ---
toleransi <- function(X) 0.05 * mean(apply(X, 2, sd))
vektor_ms <- function(P, X, h) {
W <- bobot_kernel(jarak2(P, X) / h^2)
(W %*% X) / rowSums(W) - P
}
mean_shift <- function(X, h, eps = 1e-4, maxit = 500, tol = toleransi(X),
simpan = FALSE, verbose = FALSE) {
P <- X; norma <- numeric(0)
riwayat <- if (simpan) list(P) else NULL
for (it in seq_len(maxit)) {
m <- vektor_ms(P, X, h)
P0 <- P
P <- P0 + m
norma[it] <- max(sqrt(rowSums((P - P0)^2)))
if (simpan) riwayat[[it + 1]] <- P
if (verbose)
cat("iterasi", it, "| maks ||x - x0|| =", format(norma[it], digits = 3), "\n")
if (norma[it] < eps) break
}
list(P = P, kl = cutree(hclust(dist(P), "single"), h = tol),
iter = it, norma = norma, riwayat = riwayat)
}
# --- evaluasi ---
n_kluster <- function(kl) length(unique(kl))
pusat_mode <- function(P, kl) {
do.call(rbind, lapply(split(as.data.frame(P), kl), colMeans))
}
purity <- function(kl, asal) sum(apply(table(kl, asal), 1, max)) / length(asal)
galat_mode <- function(r, pusat) {
M <- pusat_mode(r$P, r$kl)
mean(apply(pusat, 1, function(p) min(sqrt(colSums((t(M) - p)^2)))))
}
# --- bantuan tampilan ---
koma <- function(x) sub("\\.", ",", as.character(x))
atur_plot <- function(baris = 1, kolom = 1, mar = c(4.5, 4.5, 3.5, 1)) {
par(mfrow = c(baris, kolom), mar = mar)
}
palet_kde <- colorRampPalette(c("#ffffff", "#c7e9c0", "#74c476", "#31a354", "#006d2c", "#00441b"))(80)
peta_kde <- function(gx, gy, f, ...) {
image(gx, gy, f^0.6, col = palet_kde, ...)
contour(gx, gy, f^0.6, add = TRUE, drawlabels = FALSE, nlevels = 8,
col = adjustcolor("#00441b", 0.5))
}
lingkaran <- function(pusat, r, ...) {
theta <- seq(0, 2 * pi, length.out = 200)
lines(pusat[1] + r * cos(theta), pusat[2] + r * sin(theta), ...)
}
lingkaran_ganda <- function(pusat, r) {
lingkaran(pusat, r, col = "black", lwd = 4)
lingkaran(pusat, r, col = "white", lwd = 2, lty = 2)
}
label_bg <- function(x, y, teks, cex = 0.72) {
w <- strwidth(teks, cex = cex); hh <- strheight(teks, cex = cex)
rect(x - w / 2 - 0.05 * u, y - hh / 2 - 0.05 * u,
x + w / 2 + 0.05 * u, y + hh / 2 + 0.05 * u,
col = adjustcolor("white", 0.9), border = "grey40")
text(x, y, teks, cex = cex, font = 2)
}
# peta KDE + data + mode + lingkaran bandwidth
plot_peta <- function(Xd, r, h2, judul, gx, gy, normalisasi = FALSE) {
f <- kde_grid(Xd, gx, gy, h2)
if (normalisasi) f <- f / max(f)
peta_kde(gx, gy, f, asp = 1, xlab = "x1", ylab = "x2", main = judul, cex.main = 0.95)
points(Xd[, 1], Xd[, 2], pch = 21, bg = "white", col = "#08519c", cex = 0.4, lwd = 0.5)
points(pusat_mode(r$P[, 1:2], r$kl), pch = 24, bg = "#ffd400", cex = 1.6)
lingkaran_ganda(c(0, 0), h2)
}
# =============================================================================
# 2. DATA
# =============================================================================
set.seed(42)
asal <- rep(1:3, times = c(167, 167, 166))
pusat_asli <- rbind(c(0, 0), c(10, 0), c(5, 9))
# E = sebaran standar tiap blob
E <- rbind(cbind(rnorm(167), rnorm(167)),
cbind(rnorm(167), rnorm(167)),
cbind(rnorm(166), rnorm(166)))
X <- E + pusat_asli[asal, ]
n <- nrow(X); d <- ncol(X)
u <- mean(apply(X, 2, sd))
eps <- 1e-4
tol <- toleransi(X)
warna <- c("#e05a5a", "#5fd35f", "#7b7bd8")
tepi <- c("#a33636", "#2a8a2a", "#44449c")
print(table(Kelompok = asal))
## Kelompok
## 1 2 3
## 167 167 166
atur_plot()
plot(X, pch = 21, cex = 0.9, lwd = 0.6, bg = adjustcolor(warna[asal], 0.8), col = tepi[asal],
xlab = "x1", ylab = "x2", main = "Scatter plot Gaussian Blobs")
legend("topright", legend = paste("Kelompok", 1:3), pch = 21, pt.bg = warna, col = tepi, bty = "n")
# =============================================================================
# 3. EKSPLORASI DATA
# =============================================================================
cat("\n=== 1. VALIDASI DATA ===\n")
##
## === 1. VALIDASI DATA ===
cat("Struktur : matriks", nrow(X), "x", ncol(X), "| tipe:", typeof(X), "\n")
## Struktur : matriks 500 x 2 | tipe: double
cat("Missing :", sum(is.na(X)), "\n")
## Missing : 0
cat("Duplikat :", sum(duplicated(X)), "\n")
## Duplikat : 0
cat("Tak hingga :", sum(is.infinite(X)), "\n")
## Tak hingga : 0
cat("\n=== 2. PEMERIKSAAN PENCILAN (1,5 x IQR) ===\n")
##
## === 2. PEMERIKSAAN PENCILAN (1,5 x IQR) ===
pencilan <- sapply(1:ncol(X), function(j) {
q <- quantile(X[, j], c(0.25, 0.75)); batas <- 1.5 * (q[2] - q[1])
sum(X[, j] < q[1] - batas | X[, j] > q[2] + batas)
})
names(pencilan) <- c("x1", "x2")
cat("Jumlah pencilan per fitur:\n"); print(pencilan)
## Jumlah pencilan per fitur:
## x1 x2
## 0 0
atur_plot()
boxplot(X, names = c("x1", "x2"), col = c("#e05a5a", "#7b7bd8"),
main = "Boxplot (aturan 1,5 x IQR)", ylab = "Nilai")
cat("\n=== 3. SKALA FITUR ===\n")
##
## === 3. SKALA FITUR ===
sdv <- apply(X, 2, sd)
cat("SD x1 =", round(sdv[1], 3), "| SD x2 =", round(sdv[2], 3),
"| Rasio SD =", round(min(sdv) / max(sdv), 3), "(tanpa standardisasi)\n")
## SD x1 = 4.206 | SD x2 = 4.406 | Rasio SD = 0.955 (tanpa standardisasi)
# =============================================================================
# 4. MEAN SHIFT PADA DATA UTAMA
# =============================================================================
h1 <- bandwidth_kuantil(X, 0.3)
cat("\nBandwidth h =", round(h1, 4), "\n")
##
## Bandwidth h = 3.0726
ms <- mean_shift(X, h1, eps = eps, tol = tol, maxit = 200, simpan = TRUE, verbose = TRUE)
## iterasi 1 | maks ||x - x0|| = 3.22
## iterasi 2 | maks ||x - x0|| = 1.26
## iterasi 3 | maks ||x - x0|| = 0.475
## iterasi 4 | maks ||x - x0|| = 0.114
## iterasi 5 | maks ||x - x0|| = 0.0244
## iterasi 6 | maks ||x - x0|| = 0.00507
## iterasi 7 | maks ||x - x0|| = 0.00105
## iterasi 8 | maks ||x - x0|| = 0.000218
## iterasi 9 | maks ||x - x0|| = 4.5e-05
P <- ms$P
kl <- ms$kl
riw <- ms$riwayat
norma <- ms$norma
Tit <- ms$iter
cat("Konvergen:", Tit, "iterasi | kluster:", n_kluster(kl), "\n")
## Konvergen: 9 iterasi | kluster: 3
print(table(kl))
## kl
## 1 2 3
## 167 167 166
sel <- sort(unique(pmin(c(0, 1, 2, 3, 4, Tit), Tit)))
atur_plot(2, 3, mar = c(4, 4, 3, 1))
for (it in sel) {
plot(X, pch = 16, cex = 0.4, col = "grey85", xlab = "x1", ylab = "x2",
main = if (it == 0) "Iterasi 0: data awal" else
paste0("Iterasi ", it, " (maks ||x - x0|| = ", format(norma[it], digits = 2), ")"))
if (it > 0)
segments(riw[[it]][, 1], riw[[it]][, 2], riw[[it + 1]][, 1], riw[[it + 1]][, 2],
col = adjustcolor("red3", 0.2))
points(riw[[it + 1]], pch = 16, cex = 0.5, col = adjustcolor(kl + 1, 0.7))
}
par(mfrow = c(1, 2))
plot(norma, type = "b", pch = 19, col = "blue", log = "y", xlab = "Iterasi",
ylab = "maks ||x - x0|| (skala log)", main = "Cek konvergensi")
abline(h = eps, col = "red", lty = 2, lwd = 2)
plot(X, col = kl + 1, pch = 16, cex = 0.5, xlab = "x1", ylab = "x2",
main = paste0("Hasil: ", n_kluster(kl), " kluster"))
points(pusat_mode(P, kl), pch = 24, bg = "gold", cex = 2)
gx <- seq(min(X[, 1]) - 0.5 * u, max(X[, 1]) + 0.5 * u, length.out = 120)
gy <- seq(min(X[, 2]) - 0.5 * u, max(X[, 2]) + 0.5 * u, length.out = 120)
f <- kde_grid(X, gx, gy, h1)
mode_k <- pusat_mode(P, kl)
f_mode <- apply(mode_k, 1, function(p) kde(p, X, h1))
jauh <- sqrt(rowSums((X - P)^2))
rep_i <- which.max(jauh)
rep_k <- sapply(sort(unique(kl)), function(k) { idx <- which(kl == k); idx[which.max(jauh[idx])] })
lintasan <- function(i) t(sapply(riw, function(M) M[i, ]))
jejak <- lintasan(rep_i)
p0 <- jejak[1, ]; p1 <- jejak[2, ]; pT <- jejak[Tit + 1, ]
m1 <- p1 - p0
cat("\n--- NOMINAL UNTUK ANOTASI ---\n")
##
## --- NOMINAL UNTUK ANOTASI ---
cat("1. KDE f_K(mode) :", sprintf("%.4f", f_mode), "\n")
## 1. KDE f_K(mode) : 0.0052 0.0052 0.0051
cat("2. Kernel Gaussian K(0) :", sprintf("%.4f", kernel_gaussian(0, d)), "\n")
## 2. Kernel Gaussian K(0) : 0.1592
cat("3. m(x) titik contoh :", sprintf("(%.4f, %.4f), ||m|| = %.4f", m1[1], m1[2], sqrt(sum(m1^2))), "\n")
## 3. m(x) titik contoh : (-3.0464, -1.0552), ||m|| = 3.2239
cat("4. x_baru = x + m(x) :", Tit, "iterasi | dari", sprintf("(%.3f, %.3f)", p0[1], p0[2]),
"ke", sprintf("(%.3f, %.3f)", pT[1], pT[2]), "\n")
## 4. x_baru = x + m(x) : 9 iterasi | dari (8.211, 10.245) ke (4.806, 8.977)
cat("5. Bandwidth h = Q_0,3 :", sprintf("%.4f", h1), "\n")
## 5. Bandwidth h = Q_0,3 : 3.0726
cat("6. Kriteria henti :", sprintf("maks ||x - x0|| = %.2e < eps = %g", norma[Tit], eps), "\n")
## 6. Kriteria henti : maks ||x - x0|| = 4.50e-05 < eps = 0.0001
col_titik <- "#08519c"; col_jejak <- "#e31a1c"; col_eps <- "#7b2cbf"; col_box <- "#08519c"
xl <- c(min(gx) - 3.4 * u, max(gx) + 3.4 * u)
yl <- c(min(gy) - 1 * u, max(gy) + 0.3 * u)
atur_plot(mar = c(4.5, 4.5, 4, 1))
peta_kde(gx, gy, f, xlim = xl, ylim = yl, xlab = "x1", ylab = "x2",
main = sprintf("Mean Shift | h = %.3f | %d kluster", h1, n_kluster(kl)))
points(X, pch = 21, bg = "white", col = col_titik, lwd = 0.6, cex = 0.5)
for (i in rep_k) {
jj <- lintasan(i)
q0 <- jj[1, ]; q1 <- jj[2, ]; qT <- jj[Tit + 1, ]
lines(jj, col = col_jejak, lwd = 3)
points(jj, pch = 21, bg = col_jejak, col = "white", cex = 1.1)
lingkaran(q0, h1, lty = 2, lwd = 2.5, col = "black")
arrows(q0[1], q0[2], q1[1], q1[2], col = col_jejak, lwd = 4, length = 0.12)
points(q0[1], q0[2], pch = 21, bg = "white", col = "black", cex = 1.8, lwd = 2.5)
lingkaran(qT, 0.28 * u, col = col_eps, lwd = 3.5)
}
points(mode_k, pch = 24, bg = "#ffd400", col = "black", cex = 2.6, lwd = 1.5)
for (i in seq_len(nrow(mode_k)))
label_bg(mode_k[i, 1], mode_k[i, 2] - 0.55 * u,
sprintf("(%.2f, %.2f) f = %.3f", mode_k[i, 1], mode_k[i, 2], f_mode[i]))
kotak <- function(x, y, no, teks, target, kanan = FALSE) {
w <- strwidth(teks, cex = 0.75); hh <- strheight(teks, cex = 0.75); pad <- 0.1 * u
x0 <- if (kanan) x - w else x
xc <- x0 - pad - 0.18 * u
segments(if (kanan) xc else x0 + w + pad, y, target[1], target[2], col = "grey20", lwd = 1.3)
rect(x0 - pad, y - hh / 2 - pad, x0 + w + pad, y + hh / 2 + pad,
col = "#fffbe6", border = col_box, lwd = 2)
text(x0, y, teks, adj = c(0, 0.5), cex = 0.75)
points(xc, y, pch = 21, bg = col_box, col = "white", cex = 2.6)
text(xc, y, no, col = "white", font = 2, cex = 0.85)
}
ys <- seq(max(gy) - 0.2 * u, min(gy) + 1.0 * u, length.out = 3)
xk <- xl[1] + 0.45 * u; xr <- xl[2] - 0.35 * u
kotak(xk, ys[1], 1, sprintf("Kepadatan (KDE)\nf_K(mode) = %s",
paste(sprintf("%.3f", f_mode), collapse = ", ")),
mode_k[kl[rep_i], ])
kotak(xk, ys[2], 2, sprintf("Kernel Gaussian (radial)\nd = %d, K(0) = %.4f",
d, kernel_gaussian(0, d)), p0)
kotak(xk, ys[3], 3, sprintf("Vektor mean shift m(x)\nm = (%.3f, %.3f)\n||m|| = %.3f",
m1[1], m1[2], sqrt(sum(m1^2))), (p0 + p1) / 2)
kotak(xr, ys[1], 4, sprintf("x_baru = x + m(x)\n%d iterasi\n(%.2f, %.2f) -> (%.2f, %.2f)",
Tit, p0[1], p0[2], pT[1], pT[2]),
jejak[min(3, Tit + 1), ], kanan = TRUE)
kotak(xr, ys[2], 5, sprintf("Bandwidth h = Q0,3(dist(X'))\nh = %.4f (radius)", h1),
p0 + h1 * c(cos(pi / 4), sin(pi / 4)), kanan = TRUE)
kotak(xr, ys[3], 6, sprintf("Kriteria henti\n||x - x0|| = %.1e\n< eps = %g",
norma[Tit], eps), pT, kanan = TRUE)
legend("bottom", ncol = 3, bty = "o", bg = "white", cex = 0.8,
legend = c("titik data", "mode (puncak)", "jejak iterasi", "radius h", "vektor m(x)", "lingkaran eps (diperbesar)"),
pch = c(21, 24, 21, NA, NA, NA), lty = c(NA, NA, 1, 2, 1, 1), lwd = c(NA, NA, 3, 2.5, 4, 3.5),
col = c(col_titik, "black", col_jejak, "black", col_jejak, col_eps),
pt.bg = c("white", "#ffd400", col_jejak, NA, NA, NA), pt.cex = c(1, 1.6, 1, NA, NA, NA))
# =============================================================================
# 5. PENGARUH BANDWIDTH
# =============================================================================
dst <- dist(X)
qs <- c(0.05, 0.1, 0.2, 0.3, 0.5, 0.7)
cand <- setNames(sapply(qs, function(q) as.numeric(quantile(dst, q))), paste0("Q", koma(qs)))
print(round(cand, 3))
## Q0,05 Q0,1 Q0,2 Q0,3 Q0,5 Q0,7
## 0.808 1.198 1.914 3.073 9.361 10.514
ms_cand <- lapply(cand, function(h2) mean_shift(X, h2, eps = eps, tol = tol))
ringkas <- data.frame(
sumber = names(cand), h = round(cand, 3),
iterasi = sapply(ms_cand, `[[`, "iter"),
kluster = sapply(ms_cand, function(r) n_kluster(r$kl)),
purity = sapply(ms_cand, function(r) round(purity(r$kl, asal), 3)),
galat_mode = sapply(ms_cand, function(r) round(galat_mode(r, pusat_asli), 3)),
silhouette = sapply(ms_cand, function(r)
if (n_kluster(r$kl) < 2) NA else
round(mean(cluster::silhouette(r$kl, dst)[, 3]), 3)),
row.names = NULL
)
print(ringkas)
## sumber h iterasi kluster purity galat_mode silhouette
## 1 Q0,05 0.808 34 3 1.000 0.155 0.82
## 2 Q0,1 1.198 15 3 1.000 0.111 0.82
## 3 Q0,2 1.914 9 3 1.000 0.099 0.82
## 4 Q0,3 3.073 9 3 1.000 0.167 0.82
## 5 Q0,5 9.361 9 1 0.334 5.888 NA
## 6 Q0,7 10.514 8 1 0.334 5.888 NA
atur_plot(2, 3, mar = c(4, 4, 4, 1))
for (i in seq_along(cand))
plot_peta(X, ms_cand[[i]], cand[[i]],
sprintf("%s\nh = %.2f | %d kluster", names(cand)[i], cand[[i]], n_kluster(ms_cand[[i]]$kl)),
gx, gy)
# =============================================================================
# 6. PENGARUH DIMENSI DATA
# =============================================================================
dims <- c(2, 3, 5, 10, 20, 50)
jmin <- 10 # jarak antar-pusat terdekat: (0,0) ke (10,0)
set.seed(123)
Xd_list <- list(); ms_dim <- list(); hasil_d <- list(); h_dim <- numeric(length(dims))
for (i in seq_along(dims)) {
dd <- dims[i]
Eb <- if (dd == 2) E else cbind(E, matrix(rnorm(nrow(E) * (dd - 2)), nrow(E)))
pus <- cbind(pusat_asli, matrix(0, 3, dd - 2))
Xd <- Eb + pus[asal, ]
h2 <- bandwidth_kuantil(Xd, 0.3)
wkt <- system.time(r <- mean_shift(Xd, h2, eps = eps))[["elapsed"]]
Xd_list[[i]] <- Xd; ms_dim[[i]] <- r; h_dim[i] <- h2
hasil_d[[i]] <- data.frame(d = dd, h = round(h2, 2), h_per_jarak = round(h2 / jmin, 2),
iterasi = r$iter, kluster = n_kluster(r$kl),
purity = round(purity(r$kl, asal), 3), detik = round(wkt, 2))
}
dat_d <- do.call(rbind, hasil_d)
print(dat_d)
## d h h_per_jarak iterasi kluster purity detik
## 1 2 3.07 0.31 9 3 1.000 0.63
## 2 3 3.54 0.35 13 3 1.000 0.94
## 3 5 4.35 0.43 164 1 0.334 4.94
## 4 10 5.65 0.56 20 1 0.334 0.43
## 5 20 7.55 0.76 11 1 0.334 0.33
## 6 50 11.18 1.12 8 1 0.334 0.31
gx7 <- seq(-6, 16, length.out = 100)
gy7 <- seq(-6, 15, length.out = 100)
atur_plot(2, 3, mar = c(4, 4, 4, 1))
for (i in seq_along(dims))
plot_peta(Xd_list[[i]], ms_dim[[i]], h_dim[i],
sprintf("d = %d | h = %.2f\nh/jarak = %.2f | %d kluster",
dims[i], h_dim[i], h_dim[i] / jmin, n_kluster(ms_dim[[i]]$kl)),
gx7, gy7, normalisasi = TRUE)
atur_plot(1, 2)
plot(dat_d$d, dat_d$h_per_jarak, type = "b", pch = 21, cex = 1.6, log = "x",
bg = c("#d73027", "#1a9850")[(dat_d$kluster == 3) + 1],
xlab = "Dimensi (d)", ylab = "h / jarak antar-pusat", main = "Rasio h/jarak vs dimensi")
abline(h = 0.4, lty = 2)
legend("topleft", legend = c("3 kluster", "bukan 3"), pt.bg = c("#1a9850", "#d73027"),
pch = 21, bty = "n")
plot(dat_d$d, dat_d$iterasi, type = "b", pch = 19, col = "blue", log = "x",
xlab = "Dimensi (d)", ylab = "Iterasi sampai konvergen", main = "Iterasi vs dimensi")
