پیاده‌سازی Mean Shift در پایتون؛ آموزش عملی از صفر تا scikit-learn

1. مقدمه

هدف این بخش انتقال از «دانستن» به «توانستن» است. فصل اصلی رابطه Mean Shift با برآورد چگالی هسته‌ای، بردار انتقال میانگین، fixed-point iteration، bandwidth، حوزه جذب، mode merging و محدودیت‌های همگرایی را تثبیت کرده است. در اینجا همان قراردادها به یک workflow قابل اجرا تبدیل می‌شوند، بدون آنکه مبانی نظری دوباره تدریس شوند.

سه سطح عملی در پیوست دنبال می‌شود:

1. پیاده‌سازی مرجع Gaussian Mean Shift از صفر برای مشاهده مستقیم update فصل؛
2. مطالعه موردی Iris با flat-kernel Mean Shift در scikit-learn برای نشان‌دادن یک workflow استاندارد Load → Split → Scale → Tune → Fit → Predict → Evaluate؛
3. سناریوی پیشرفته با چگالی نامساوی، عدم‌توازن و نقاط پرت برای نمایش اثر bandwidth، scaling، cluster_all و seed strategy.

در هر مطالعه، هر جا برچسب ground truth موجود است، فقط پس از آموزش برای سنجش پسینی استفاده می‌شود. انتخاب bandwidth و ساختار خوشه‌بندی بدون استفاده از برچسب انجام می‌شود.

.

2. پیش‌نیازهای عملی و ساختار داده ورودی

2.1 محیط اجرایی آزموده‌شده

این بسته روی محیط زیر اجرا و آزمون شده است:

مؤلفه	نسخه آزموده‌شده
Python	3.13.5
NumPy	2.3.5
SciPy	1.17.0
pandas	2.2.3
scikit-learn	1.8.0
Matplotlib	3.10.8
pytest	8.x یا بالاتر

نسخه‌ها در requirements.txt و مشخصات محیط در environment_manifest.json ثبت شده‌اند. مستندات stable فعلی scikit-learn نیز MeanShift را با flat kernel معرفی می‌کنند و برای estimate_bandwidth هشدار می‌دهند که هزینه محاسباتی دست‌کم نسبت به تعداد نمونه‌های مورد استفاده درجه دوم است.

.

2.2 نصب و اجرای تکرارپذیر

python -m venv .venv
source .venv/bin/activate        # Linux/macOS
# .venv\Scripts\activate       # Windows
pip install -r requirements.txt

python examples/small_example.py
python case_studies/iris_case_study.py
python case_studies/noisy_imbalanced_case_study.py
pytest

برای اجرای دقیق این بسته، seedهای تصادفی در فایل manifest ثبت شده‌اند. فایل checksums.sha256 نیز امکان کنترل یکسان‌بودن تمام فایل‌های ورودی و خروجی را فراهم می‌کند.

Input Schema

نسخه پایه این پیوست انتظار دارد:

• داده عددی و دوبعدی با شکل (n_samples, n_features)؛
• مقادیر finite؛
• نبود NaN و Inf؛
• ویژگی‌هایی که بعد از preprocessing مقیاس قابل مقایسه داشته باشند؛
• فاصله اقلیدسی معنادار باشد.

اگر داده categorical، mixed-type، sparse بسیار پُربُعد یا manifold باشد، این پیاده‌سازی مرجع نباید بدون اصلاح metric و update rule استفاده شود.

.

2.3 پیش‌پردازش ضروری

Scaling

Mean Shift مستقیماً به فاصله حساس است. اگر یک ویژگی در بازه [0,1] و دیگری در [0,10000] باشد، ویژگی دوم هندسه فضای ویژگی و در نتیجه KDE یا flat neighborhood را غالب می‌کند. در مطالعه Iris از StandardScaler و در مطالعه دارای outlier از RobustScaler استفاده می‌شود.

.

Missing values

نسخه پایه missing-value aware نیست. پیش از اجرا باید با روش مناسب دامنه مسئله، داده‌های مفقود حذف یا impute شوند. افزودن یک imputer عمومی بدون تحلیل سازوکار مفقودی توصیه نمی‌شود.

.

Outliers

از آنجا که bandwidth کوچک می‌تواند ساختارهای محلی بسیار ریز بسازد، outlier یا گروه کوچک outlierها ممکن است modeهای مستقل تولید کنند. مطالعه دوم این وضعیت را عملاً نشان می‌دهد.

.

3. پیاده‌سازی پایه از صفر: Gaussian Mean Shift

3.1 قرارداد ریاضی پیاده‌سازی

به‌روزرسانی مرجع دقیقاً از رابطه Gaussian فصل استفاده می‌کند:

و شرط توقف اجرایی:

merge_radius در این کد تصمیم اجرایی است، نه قضیه Mean Shift. مقدار پیش‌فرض آن bandwidth / 2 انتخاب شده و در API آشکار نگه داشته شده است تا دانشجو آن را با بخش mode merging فصل اشتباه نگیرد.

.

3.2 ویژگی‌های مهندسی پیاده‌سازی

• اعتبارسنجی صریح ورودی؛
• جلوگیری از NaN/Inf؛
• seed مستقل برای هر trajectory؛
• max_iter و flag همگرایی برای هر seed؛
• ادغام deterministic مُدها بر اساس امتیاز چگالی؛
• fit, predict, fit_predict؛
• type hints و docstring؛
• ۱۰ آزمون خودکار؛
• نبود وابستگی به کتابخانه clustering سطح بالا.

.

3.3 کد کامل پیاده‌سازی مرجع

"""Reference Gaussian Mean Shift implementation aligned with the audited chapter.

The implementation is intentionally explicit: Gaussian weights are used for the
fixed-point update, each seed is iterated independently, and converged modes are
merged in a deterministic post-processing stage.
"""
from __future__ import annotations

from dataclasses import dataclass
from typing import Iterable

import numpy as np
from numpy.typing import ArrayLike, NDArray


FloatArray = NDArray[np.float64]


@dataclass(frozen=True)
class MeanShiftConfig:
    bandwidth: float
    tol: float = 1e-3
    max_iter: int = 300
    merge_radius: float | None = None

    def validate(self) -> None:
        if not np.isfinite(self.bandwidth) or self.bandwidth <= 0:
            raise ValueError("bandwidth must be a finite positive number")
        if not np.isfinite(self.tol) or self.tol <= 0:
            raise ValueError("tol must be a finite positive number")
        if self.max_iter < 1:
            raise ValueError("max_iter must be >= 1")
        if self.merge_radius is not None and self.merge_radius <= 0:
            raise ValueError("merge_radius must be positive when provided")


class GaussianMeanShift:
    """Gaussian Mean Shift for educational and reproducible experiments.

    Parameters
    ----------
    bandwidth:
        Gaussian bandwidth ``h`` from the chapter.
    tol:
        Stop a seed when ``||y_{t+1} - y_t||_2 < tol``.
    max_iter:
        Maximum number of fixed-point updates per seed.
    merge_radius:
        Distance threshold for merging numerically close converged modes.
        Defaults to ``bandwidth / 2``. This is an implementation decision, not
        a theorem of Mean Shift, and is therefore exposed explicitly.
    """

    def __init__(
        self,
        bandwidth: float,
        *,
        tol: float = 1e-3,
        max_iter: int = 300,
        merge_radius: float | None = None,
    ) -> None:
        self.config = MeanShiftConfig(bandwidth, tol, max_iter, merge_radius)
        self.config.validate()

    @staticmethod
    def _as_2d_float(X: ArrayLike) -> FloatArray:
        arr = np.asarray(X, dtype=np.float64)
        if arr.ndim != 2:
            raise ValueError("X must be a 2D array of shape (n_samples, n_features)")
        if arr.shape[0] == 0 or arr.shape[1] == 0:
            raise ValueError("X must be non-empty")
        if not np.isfinite(arr).all():
            raise ValueError("X must contain only finite values")
        return arr

    def _weights(self, y: FloatArray, X: FloatArray) -> FloatArray:
        h = self.config.bandwidth
        squared_distance = np.sum((X - y) ** 2, axis=1)
        return np.exp(-squared_distance / (2.0 * h * h))

    def _step(self, y: FloatArray, X: FloatArray) -> FloatArray:
        weights = self._weights(y, X)
        total = float(weights.sum())
        if total <= np.finfo(np.float64).tiny:
            # Gaussian weights are positive analytically, but extreme floating
            # point underflow can make all of them numerically zero.
            return y.copy()
        return np.sum(weights[:, None] * X, axis=0) / total

    def _iterate_seed(self, seed: FloatArray, X: FloatArray) -> tuple[FloatArray, int, bool]:
        y = seed.astype(np.float64, copy=True)
        for iteration in range(1, self.config.max_iter + 1):
            y_next = self._step(y, X)
            shift = float(np.linalg.norm(y_next - y))
            y = y_next
            if shift < self.config.tol:
                return y, iteration, True
        return y, self.config.max_iter, False

    def _density_score(self, mode: FloatArray, X: FloatArray) -> float:
        # Normalization constants are irrelevant for ranking modes.
        return float(self._weights(mode, X).sum())

    def _merge_modes(self, modes: FloatArray, X: FloatArray) -> FloatArray:
        radius = self.config.merge_radius
        if radius is None:
            radius = self.config.bandwidth / 2.0

        scores = np.asarray([self._density_score(m, X) for m in modes])
        # Stable deterministic ordering: high estimated density first, then index.
        order = np.argsort(-scores, kind="stable")
        accepted: list[FloatArray] = []
        for idx in order:
            candidate = modes[idx]
            if not accepted:
                accepted.append(candidate.copy())
                continue
            distances = np.linalg.norm(np.vstack(accepted) - candidate, axis=1)
            if np.all(distances > radius):
                accepted.append(candidate.copy())
        return np.vstack(accepted)

    def fit(self, X: ArrayLike, *, seeds: ArrayLike | None = None) -> "GaussianMeanShift":
        X_arr = self._as_2d_float(X)
        seeds_arr = X_arr.copy() if seeds is None else self._as_2d_float(seeds)
        if seeds_arr.shape[1] != X_arr.shape[1]:
            raise ValueError("seeds and X must have the same number of features")

        raw_modes = []
        iterations = []
        converged = []
        for seed in seeds_arr:
            mode, n_iter, ok = self._iterate_seed(seed, X_arr)
            raw_modes.append(mode)
            iterations.append(n_iter)
            converged.append(ok)

        raw_modes_arr = np.vstack(raw_modes)
        centers = self._merge_modes(raw_modes_arr, X_arr)

        self.X_ = X_arr
        self.raw_modes_ = raw_modes_arr
        self.cluster_centers_ = centers
        self.n_iter_per_seed_ = np.asarray(iterations, dtype=int)
        self.converged_per_seed_ = np.asarray(converged, dtype=bool)
        self.labels_ = self.predict(X_arr)
        return self

    def predict(self, X: ArrayLike) -> NDArray[np.int64]:
        if not hasattr(self, "cluster_centers_"):
            raise RuntimeError("fit must be called before predict")
        X_arr = self._as_2d_float(X)
        if X_arr.shape[1] != self.cluster_centers_.shape[1]:
            raise ValueError("X has an unexpected number of features")
        diff = X_arr[:, None, :] - self.cluster_centers_[None, :, :]
        distances = np.sum(diff * diff, axis=2)
        return np.argmin(distances, axis=1).astype(np.int64)

    def fit_predict(self, X: ArrayLike, *, seeds: ArrayLike | None = None) -> NDArray[np.int64]:
        return self.fit(X, seeds=seeds).labels_


def gaussian_mean_shift_step(X: ArrayLike, y: ArrayLike, bandwidth: float) -> FloatArray:
    """One transparent Gaussian Mean Shift update for hand-checking examples."""
    model = GaussianMeanShift(bandwidth)
    X_arr = model._as_2d_float(X)
    y_arr = np.asarray(y, dtype=np.float64)
    if y_arr.shape != (X_arr.shape[1],):
        raise ValueError("y must be one vector with n_features elements")
    return model._step(y_arr, X_arr)

3.4 نکته مهم درباره predict

در خوشه‌بندی، «پیش‌بینی» به معنای پیش‌بینی supervised نیست. متد predict در این پیوست یک مشاهده جدید را به نزدیک‌ترین mode نهایی نسبت می‌دهد. اگر نیاز مسئله این است که نمونه‌های دور از همه modeها «ناشناخته/نویز» باقی بمانند، باید علاوه بر نزدیک‌ترین مرکز یک آستانه فاصله یا چگالی تعریف شود.

.

4. مثال آموزشی کوچک

4.1 داده و هدف

همان مثال عددی فصل با داده زیر اجرا می‌شود:

وزن‌های Gaussian محاسبه‌شده توسط کد:

نقطه	وزن
0	0.882497
1	1.000000
4	0.324652

اولین update:

و در نتیجه:

نتیجه با محاسبه دستی فصل هم‌خوان است. اجرای iterative از seed برابر 1 با tolerance بسیار کوچک پس از 20 تکرار به mode عددی حدود 1.073768 می‌رسد.

.

4.2 کد کامل مثال

from __future__ import annotations
import json
from pathlib import Path
import numpy as np
from src.mean_shift_gaussian import GaussianMeanShift, gaussian_mean_shift_step

ROOT = Path(__file__).resolve().parents[1]
OUT = ROOT / "outputs"
OUT.mkdir(exist_ok=True)

X = np.array([[0.0], [1.0], [4.0]])
y0 = np.array([1.0])
h = 2.0

y1 = gaussian_mean_shift_step(X, y0, h)
weights = np.exp(-((X[:, 0] - y0[0]) ** 2) / (2 * h**2))

model = GaussianMeanShift(h, tol=1e-8, max_iter=300, merge_radius=0.5)
model.fit(X, seeds=np.array([[1.0]]))

result = {
    "X": X[:, 0].tolist(),
    "bandwidth": h,
    "seed_y0": float(y0[0]),
    "weights": weights.tolist(),
    "y1": float(y1[0]),
    "mean_shift_vector_first_step": float(y1[0] - y0[0]),
    "final_mode_from_seed_1": float(model.cluster_centers_[0, 0]),
    "iterations": int(model.n_iter_per_seed_[0]),
}
(OUT / "small_example.json").write_text(json.dumps(result, indent=2), encoding="utf-8")
print(json.dumps(result, indent=2))

.

4.3 تفسیر اجرایی

حرکت اول فقط حدود 0.0414 است. این مشاهده عملی یک سوءبرداشت متداول را برطرف می‌کند: Mean Shift الزاماً در هر تکرار «پرش بزرگ» به مرکز خوشه ندارد. اندازه حرکت از توزیع وزن‌های kernel در همسایگی نقطه جاری نتیجه می‌شود و می‌تواند بسیار کوچک باشد.

.

5. مطالعه موردی اول: خوشه‌بندی بدون‌نظارت داده Iris

5.1 تعریف مسئله

مجموعه Iris شامل ۱۵۰ نمونه و ۴ اندازه‌گیری ریخت‌شناختی است. سه گونه شناخته‌شده در داده وجود دارد، اما Mean Shift در این مطالعه این برچسب‌ها را هنگام fitting یا انتخاب bandwidth نمی‌بیند. هدف عملی این است که workflow صحیح خوشه‌بندی بدون‌نظارت با یک داده استاندارد و قابل بازتولید نشان داده شود.

.

5.2.pipeline

1. بارگذاری داده با load_iris؛
2. split تصادفی ۷۵/۲۵ بدون stratification مبتنی بر برچسب؛
3. fit کردن StandardScaler فقط روی train؛
4. برآورد چند bandwidth کاندید با quantileهای مختلف؛
5. fit کردن MeanShift با flat kernel روی train؛
6. انتخاب bandwidth فقط با training silhouette؛
7. assignment نمونه‌های test با predict؛
8. محاسبه Silhouette روی test؛
9. محاسبه ARI/NMI فقط به‌عنوان ارزیابی پسینی.

.

5.3 خروجی واقعی اجرا

پیکربندی انتخاب‌شده:

• quantile = 0.25
• bandwidth = 1.435890
• تعداد mode روی train = 2
• Silhouette آموزش = 0.5723
• ARI آموزش، فقط پسینی = 0.5440

روی test:

• تعداد clusterهای حاضر = 2
• Silhouette = 0.6123
• ARI پسینی = 0.6269
• NMI پسینی = 0.7620

نکته مهم این است که داده دارای سه گونه است، اما criterion بدون‌نظارت در این تنظیم دو mode را ترجیح داده است. این «خطای نرم‌افزار» نیست؛ نشان می‌دهد modeهای KDE/flat-density الزاماً برابر کلاس‌های معنایی نیستند.

Silhouette در sweep پهنای‌باند Iris

تعداد modeها در برابر پهنای‌باند Iris

5.4 اثر seed strategy

با همان bandwidth:

bin_seeding	تعداد mode	زمان اجرا (ثانیه)
False	2	0.1255
True	2	0.0240

در این داده، bin_seeding=True تعداد mode را تغییر نداد و زمان اجرا را کاهش داد. این نتیجه داده‌وابسته است و نباید به همه مجموعه‌ها تعمیم داده شود؛ مطالعه دوم نمونه‌ای را نشان می‌دهد که seed reduction ساختار modeها را نیز تغییر می‌دهد.

.

5.5 کد کامل مطالعه موردی اول

"""Case study 1: reproducible flat-kernel Mean Shift on the UCI Iris data.

The clustering model and bandwidth are selected without using class labels.
Labels are consulted only after fitting to report post-hoc ARI/NMI diagnostics.
"""
from __future__ import annotations
import json
from pathlib import Path
import time
import numpy as np
import pandas as pd
from sklearn.cluster import MeanShift, estimate_bandwidth
from sklearn.datasets import load_iris
from sklearn.metrics import adjusted_rand_score, normalized_mutual_info_score, silhouette_score
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler

ROOT = Path(__file__).resolve().parents[1]
OUT = ROOT / "outputs"
OUT.mkdir(exist_ok=True)
RANDOM_STATE = 42

iris = load_iris(as_frame=True)
X = iris.data.to_numpy(dtype=float)
y = iris.target.to_numpy()

# Random split only; labels are NOT used to stratify or tune the model.
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.25, random_state=RANDOM_STATE, shuffle=True, stratify=None
)
scaler = StandardScaler().fit(X_train)
X_train_s = scaler.transform(X_train)
X_test_s = scaler.transform(X_test)

rows = []
for quantile in [0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, 0.50]:
    bandwidth = float(
        estimate_bandwidth(
            X_train_s,
            quantile=quantile,
            n_samples=len(X_train_s),
            random_state=RANDOM_STATE,
        )
    )
    if bandwidth <= 0:
        continue
    t0 = time.perf_counter()
    model = MeanShift(
        bandwidth=bandwidth,
        bin_seeding=False,
        cluster_all=True,
        max_iter=300,
        n_jobs=1,
    ).fit(X_train_s)
    elapsed = time.perf_counter() - t0
    train_labels = model.labels_
    n_clusters = len(model.cluster_centers_)
    train_silhouette = (
        float(silhouette_score(X_train_s, train_labels))
        if 1 < n_clusters < len(X_train_s)
        else float("nan")
    )
    rows.append(
        {
            "quantile": quantile,
            "bandwidth": bandwidth,
            "n_clusters": n_clusters,
            "train_silhouette": train_silhouette,
            "runtime_seconds": elapsed,
            # y_train is used only AFTER fitting; not for model selection.
            "train_ARI_posthoc": float(adjusted_rand_score(y_train, train_labels)),
        }
    )

sweep = pd.DataFrame(rows)
valid = sweep[np.isfinite(sweep["train_silhouette"])].copy()
# Purely unsupervised selection.
best_idx = valid["train_silhouette"].idxmax()
best = sweep.loc[best_idx].to_dict()

best_model = MeanShift(
    bandwidth=float(best["bandwidth"]),
    bin_seeding=False,
    cluster_all=True,
    max_iter=300,
    n_jobs=1,
).fit(X_train_s)

test_labels = best_model.predict(X_test_s)
test_n_clusters = len(np.unique(test_labels))
test_silhouette = (
    float(silhouette_score(X_test_s, test_labels))
    if 1 < test_n_clusters < len(X_test_s)
    else float("nan")
)

runtime_rows = []
for bin_seeding in [False, True]:
    t0 = time.perf_counter()
    m = MeanShift(
        bandwidth=float(best["bandwidth"]),
        bin_seeding=bin_seeding,
        min_bin_freq=2 if bin_seeding else 1,
        cluster_all=True,
        max_iter=300,
        n_jobs=1,
    ).fit(X_train_s)
    runtime_rows.append(
        {
            "bin_seeding": bin_seeding,
            "n_clusters": int(len(m.cluster_centers_)),
            "runtime_seconds": time.perf_counter() - t0,
        }
    )

sweep.to_csv(OUT / "iris_bandwidth_sweep.csv", index=False)
pd.DataFrame(runtime_rows).to_csv(OUT / "iris_seed_strategy.csv", index=False)
pd.DataFrame(best_model.cluster_centers_).to_csv(OUT / "iris_centers_scaled.csv", index=False)
pd.DataFrame(
    {
        "true_label_posthoc": y_test,
        "predicted_cluster": test_labels,
    }
).to_csv(OUT / "iris_test_assignments.csv", index=False)

summary = {
    "dataset": "sklearn.datasets.load_iris (UCI Iris data)",
    "n_samples": int(X.shape[0]),
    "n_features": int(X.shape[1]),
    "train_size": int(len(X_train)),
    "test_size": int(len(X_test)),
    "selection_rule": "maximum training silhouette; labels excluded from fitting and selection",
    "selected": {
        k: (
            float(v)
            if isinstance(v, (np.floating, float))
            else int(v)
            if isinstance(v, (np.integer, int))
            else v
        )
        for k, v in best.items()
    },
    "test": {
        "n_clusters_present": int(test_n_clusters),
        "silhouette": test_silhouette,
        "ARI_posthoc": float(adjusted_rand_score(y_test, test_labels)),
        "NMI_posthoc": float(normalized_mutual_info_score(y_test, test_labels)),
    },
    "seed_strategy_comparison": runtime_rows,
}
(OUT / "iris_summary.json").write_text(json.dumps(summary, indent=2), encoding="utf-8")
print(json.dumps(summary, indent=2))

.

5.6 Troubleshooting مطالعه اول

تعداد mode بسیار زیاد است

اقدامات به‌ترتیب اولویت:

1. افزایش bandwidth؛
2. بررسی scaling؛
3. بررسی outlier؛
4. در صورت استفاده از bin_seeding، کنترل min_bin_freq.

فقط یک mode پیدا می‌شود

• bandwidth احتمالاً بیش از حد بزرگ است؛
• quantile مورد استفاده در estimate_bandwidth را کاهش دهید؛
• بررسی کنید scaling ساختار فاصله را بیش از حد فشرده نکرده باشد.

.

ARI پایین ولی Silhouette مناسب است

این وضعیت الزاماً خطا نیست. Mean Shift modeهای چگالی را پیدا می‌کند، نه کلاس‌های برچسب‌خورده انسانی را. اول بررسی کنید هدف مسئله واقعاً density-mode clustering است یا classification.

.

نتیجه کد با فرمول Gaussian فصل متفاوت است

این انتظار طبیعی است: sklearn.cluster.MeanShift از flat kernel استفاده می‌کند، در حالی‌که پیاده‌سازی آموزشی بخش ۳ Gaussian است.

.

6. مطالعه موردی دوم: چگالی نامساوی، عدم‌توازن و نقاط پرت

6.1 سناریو

داده مصنوعی کنترل‌شده با سه ساختار اصلی ساخته شده است:

• خوشه اول: ۶۰۰ نمونه؛
• خوشه دوم: ۱۸۰ نمونه؛
• خوشه سوم: ۷۰ نمونه؛
• ۳۰ outlier یکنواخت در گستره‌ای بزرگ‌تر.

این طراحی سه محدودیت فصل را هم‌زمان فعال می‌کند: چگالی متفاوت، عدم‌توازن و outlier.

.

6.2 مقایسه baseline و تنظیم اصلاح‌شده

سناریو	Scaling		مُدها	بدون انتساب	Silhouette	ARI پسینی
Baseline	Standard	0.18	5	62.7%	0.656	0.211
Tuned	Robust	0.36	3	4.4%	0.825	0.954

Baseline با bandwidth کوچک 5 mode ساخت و 62.7% داده را با cluster_all=False بدون انتساب گذاشت. با Robust Scaling و ، سه mode اصلی بازیابی شدند و نرخ بدون انتساب به حدود 4.4% کاهش یافت.

مقایسه baseline و تنظیم اصلاح‌شده

تعداد مُدها در برابر bandwidth

6.3 چرا RobustScaler؟

StandardScaler از میانگین و انحراف معیار استفاده می‌کند و outlierها می‌توانند scale آن را جابه‌جا کنند. RobustScaler بر quantileها متکی است و در این سناریوی خاص هندسه سه ساختار اصلی را بهتر حفظ کرد. این نتیجه نباید به معنای برتری عمومی RobustScaler باشد؛ انتخاب scaler باید با سازوکار داده و نوع outlier سازگار باشد.

.

6.4 نقش cluster_all=False

هنگام fit، نمونه‌هایی که خارج از kernelهای نهایی قرار می‌گیرند می‌توانند label برابر -1 داشته باشند. این قابلیت برای آشکارکردن نمونه‌های دور در آزمایش مفید است. اما یک نکته اجرایی مهم وجود دارد: متد predict برای نمونه جدید آن را به نزدیک‌ترین centroid/mode نهایی نسبت می‌دهد. بنابراین اگر rejection واقعی در inference لازم است، باید threshold فاصله/چگالی جداگانه اضافه شود.

چهار probe زیر در خروجی این مطالعه ذخیره شده‌اند:

نقطه جدید	cluster نسبت‌داده‌شده
(-3.1, -1.1)	0
(1.5, 2.1)	1
(4.1, -0.5)	2
(8.5, 6.5)	1

آخرین نقطه عمداً بسیار دور است، اما predict همچنان نزدیک‌ترین مرکز را برمی‌گرداند. این رفتار باید در کاربردهای anomaly/noise detection صریحاً مدیریت شود.

.

6.5 seed strategy: فقط بهینه‌سازی سرعت نیست

در bandwidth تنظیم‌شده:

bin_seeding	مُدها	بدون انتساب	زمان (s)	ARI پسینی
False	12	2.8%	0.6311	0.954
True	3	4.4%	0.0110	0.954

در این نمونه، کاهش seedها زمان اجرا را به‌شدت پایین آورد، اما همچنین تعداد modeهای گزارش‌شده را از ۱۲ به ۳ کاهش داد. ARI روی داده clean تقریباً ثابت ماند، زیرا modeهای اضافه عمدتاً ساختارهای فرعی/دور را نمایندگی می‌کردند. نتیجه مهم: bin_seeding یک تصمیم محاسباتی است که می‌تواند اثر مدل‌سازی هم داشته باشد.

.

6.6 کد کامل مطالعه موردی دوم

"""Case study 2: unequal densities, imbalance, and outliers with flat-kernel Mean Shift.

Labels are preserved only for post-hoc diagnostics. The tuning rule uses only
cluster count, unassigned fraction, and silhouette score.
"""
from __future__ import annotations
import json
from pathlib import Path
import time
import numpy as np
import pandas as pd
from sklearn.cluster import MeanShift
from sklearn.datasets import make_blobs
from sklearn.metrics import adjusted_rand_score, silhouette_score
from sklearn.preprocessing import RobustScaler, StandardScaler

ROOT = Path(__file__).resolve().parents[1]
OUT = ROOT / "outputs"
OUT.mkdir(exist_ok=True)
RNG = np.random.default_rng(20260808)

X_clean, y_clean = make_blobs(
    n_samples=[600, 180, 70],
    centers=[(-3.0, -1.0), (1.4, 2.0), (4.0, -0.4)],
    cluster_std=[0.55, 0.32, 0.20],
    random_state=7,
)
outliers = RNG.uniform(low=[-8.0, -6.0], high=[9.0, 7.0], size=(30, 2))
X = np.vstack([X_clean, outliers])
y_eval = np.concatenate([y_clean, np.full(len(outliers), -1)])

bandwidths = [0.18, 0.24, 0.30, 0.36, 0.42, 0.50, 0.60, 0.72, 0.85]
rows = []
for scaler_name, scaler in {
    "standard": StandardScaler(),
    "robust": RobustScaler(quantile_range=(10, 90)),
}.items():
    Xs = scaler.fit_transform(X)
    for bandwidth in bandwidths:
        t0 = time.perf_counter()
        # bin_seeding is used during tuning to make the sweep inexpensive and
        # reproducible. A later comparison shows that this is a real modeling
        # decision, not merely an optimization switch.
        model = MeanShift(
            bandwidth=bandwidth,
            bin_seeding=True,
            min_bin_freq=4,
            cluster_all=False,
            max_iter=200,
            n_jobs=1,
        ).fit(Xs)
        labels = model.labels_
        elapsed = time.perf_counter() - t0
        assigned = labels >= 0
        unique_assigned = np.unique(labels[assigned])
        silhouette = float("nan")
        if assigned.sum() > 2 and 1 < len(unique_assigned) < assigned.sum():
            silhouette = float(silhouette_score(Xs[assigned], labels[assigned]))
        rows.append({
            "scaler": scaler_name,
            "bandwidth": bandwidth,
            "n_clusters": int(len(model.cluster_centers_)),
            "unassigned_fraction": float(np.mean(labels < 0)),
            "silhouette_assigned": silhouette,
            "runtime_seconds": elapsed,
            "ARI_clean_posthoc": float(adjusted_rand_score(y_clean, labels[: len(y_clean)])),
        })

sweep = pd.DataFrame(rows)
sweep.to_csv(OUT / "noisy_imbalanced_sweep.csv", index=False)

# Deliberately problematic baseline: small h + ordinary scaling.
baseline = sweep[(sweep.scaler == "standard") & (sweep.bandwidth == 0.18)].iloc[0]

# Tuning uses no true labels: require few orphans and plausible non-degenerate
# cluster count, then choose the highest silhouette.
candidates = sweep[
    (sweep.scaler == "robust")
    & (sweep.n_clusters.between(2, 6))
    & (sweep.unassigned_fraction <= 0.05)
    & np.isfinite(sweep.silhouette_assigned)
].copy()
tuned = candidates.sort_values(["silhouette_assigned", "runtime_seconds"], ascending=[False, True]).iloc[0]

# Compare seed strategies at the tuned bandwidth. This demonstrates the chapter's
# warning that seed reduction affects both scalability and possibly the modes found.
Xs_tuned = RobustScaler(quantile_range=(10, 90)).fit_transform(X)
seed_rows = []
for bin_seeding in [False, True]:
    t0 = time.perf_counter()
    model = MeanShift(
        bandwidth=float(tuned.bandwidth),
        bin_seeding=bin_seeding,
        min_bin_freq=4 if bin_seeding else 1,
        cluster_all=False,
        max_iter=200,
        n_jobs=1,
    ).fit(Xs_tuned)
    seed_rows.append({
        "bin_seeding": bin_seeding,
        "n_clusters": int(len(model.cluster_centers_)),
        "unassigned_fraction": float(np.mean(model.labels_ < 0)),
        "runtime_seconds": time.perf_counter() - t0,
        "ARI_clean_posthoc": float(adjusted_rand_score(y_clean, model.labels_[: len(y_clean)])),
    })

# Fit the selected configuration and demonstrate assignment of new observations.
robust_scaler = RobustScaler(quantile_range=(10, 90)).fit(X)
X_tuned = robust_scaler.transform(X)
tuned_model = MeanShift(
    bandwidth=float(tuned.bandwidth),
    bin_seeding=True,
    min_bin_freq=4,
    cluster_all=False,
    max_iter=200,
    n_jobs=1,
).fit(X_tuned)
probe_points = np.array([[-3.1, -1.1], [1.5, 2.1], [4.1, -0.5], [8.5, 6.5]])
probe_labels = tuned_model.predict(robust_scaler.transform(probe_points))
probe_table = pd.DataFrame({
    "x1": probe_points[:, 0],
    "x2": probe_points[:, 1],
    "assigned_cluster": probe_labels,
})

selected = pd.DataFrame([baseline, tuned], index=["baseline", "robust_tuned"]).reset_index(names="scenario")
selected.to_csv(OUT / "noisy_imbalanced_selected.csv", index=False)
pd.DataFrame(seed_rows).to_csv(OUT / "noisy_seed_strategy.csv", index=False)
probe_table.to_csv(OUT / "noisy_probe_predictions.csv", index=False)
np.savez_compressed(OUT / "noisy_imbalanced_data.npz", X=X, y_eval=y_eval, X_clean=X_clean, y_clean=y_clean)

summary = {
    "n_clean": int(len(X_clean)),
    "n_outliers": int(len(outliers)),
    "imbalance_counts": [600, 180, 70],
    "selection_policy": "robust scaling; 2-6 modes; <=5% unassigned; maximize silhouette; labels excluded from tuning",
    "baseline": baseline.to_dict(),
    "tuned": tuned.to_dict(),
    "seed_strategy_comparison_at_tuned_bandwidth": seed_rows,
    "probe_predictions": probe_table.to_dict(orient="records"),
}
def clean(v):
    if isinstance(v, dict): return {k: clean(x) for k, x in v.items()}
    if isinstance(v, list): return [clean(x) for x in v]
    if isinstance(v, np.integer): return int(v)
    if isinstance(v, np.floating): return float(v)
    if isinstance(v, np.bool_): return bool(v)
    return v
summary = clean(summary)
(OUT / "noisy_imbalanced_summary.json").write_text(json.dumps(summary, indent=2), encoding="utf-8")
print(json.dumps(summary, indent=2))

.

6.7 Troubleshooting پیشرفته

نشانه: ده‌ها mode کوچک یا modeهای تک‌نمونه‌ای

• bandwidth را افزایش دهید؛
• outlierها و scaling را بررسی کنید؛
• mode merging را کنترل کنید؛
• در flat Mean Shift، bin_seeding و min_bin_freq را به‌عنوان ابزار کنترل seedها آزمایش کنید.

.

نشانه: خوشه کوچک ناپدید می‌شود

• bandwidth را کاهش دهید؛
• از bandwidth محلی/تطبیقی در توسعه پیشرفته استفاده کنید؛
• بررسی کنید preprocessing خوشه کوچک را به ساختار بزرگ فشرده نکرده باشد.

.

نشانه: زمان اجرا زیاد است

• تعداد seedها را کاهش دهید؛
• bin_seeding=True را آزمایش کنید؛
• برای estimate_bandwidth از n_samples استفاده کنید؛
• در داده بزرگ، approximate neighbor search یا واریانت‌های scalable را در نظر بگیرید؛
• فراموش نکنید کاهش seedها ممکن است modeهای نهایی را نیز تغییر دهد.

.

نشانه: اجرای from-scratch با scikit-learn یکسان نیست

ابتدا kernel را کنترل کنید. پیاده‌سازی مرجع این پیوست Gaussian و scikit-learn flat است. سپس tolerance، seed set و mode-merging rule را مقایسه کنید.

.

نشانه: الگوریتم به سقف iteration می‌رسد

• مسیرهای converged_per_seed_ را بررسی کنید؛
• tolerance را اندکی بزرگ‌تر کنید؛
• bandwidth بسیار کوچک را بررسی کنید؛
• scaling را کنترل کنید؛
• این وضعیت را با «اثبات عدم همگرایی نظری» یکی ندانید؛ ممکن است صرفاً ناشی از cutoff اجرایی max_iter باشد.

.

7. جمع‌بندی نهایی

نکات اجرایی کلیدی

• kernel را همیشه ثبت کنید. Gaussian Mean Shift و flat Mean Shift ممکن است با bandwidth مشابه جواب متفاوت بدهند.
• bandwidth مهم‌ترین کنترل granularity است. sweep عملی آن باید بخشی از workflow باشد.
• feature scaling اختیاریِ تزئینی نیست. هر metric-based clustering به هندسه scale وابسته است.
• mode merging را صریح کنید. مقدار threshold می‌تواند نتیجه clustering را تغییر دهد.
• seed strategy را هم تصمیم محاسباتی و هم تصمیم مدل‌سازی بدانید.
• برچسب‌های واقعی را وارد tuning بدون‌نظارت نکنید. ARI/NMI در این پیوست فقط پس از fitting گزارش می‌شوند.
• predict در clustering همان classification نیست. در sklearn، نمونه جدید به نزدیک‌ترین مرکز نسبت داده می‌شود؛ rejection جداگانه نیاز به قاعده اضافی دارد.
• estimate_bandwidth می‌تواند bottleneck باشد. برای داده بزرگ subsampling باید کنترل‌شده و با seed ثبت‌شده انجام شود.

.

هشدارهای پیاده‌سازی

1. ادعای «تعداد خوشه کاملاً خودکار» نکنید؛ bandwidth تعداد modeها را کنترل می‌کند.
2. با flat kernel، همسایه درون شعاع و بیرون شعاع رفتار گسسته دارد؛ Gaussian weighting پیوسته است.
3. cluster_all=False را با قابلیت rejection در predict اشتباه نگیرید.
4. تعداد mode بیشتر همیشه به معنی مدل بهتر نیست.
5. Silhouette بالا در داده دارای orphan زیاد می‌تواند گمراه‌کننده باشد؛ نرخ بدون انتساب را همراه آن گزارش کنید.
6. ARI/NMI پایین به‌تنهایی اثبات نمی‌کند clustering غلط است؛ ممکن است modal structure با کلاس معنایی متفاوت باشد.

.

توسعه‌های عملیاتی بعدی

• adaptive/variable bandwidth؛
• approximate nearest-neighbor Mean Shift؛
• GPU implementation؛
• streaming Mean Shift؛
• density-threshold rejection برای predict؛
• مقایسه Gaussian و flat kernel روی یک dataset واحد؛
• benchmark زمان اجرا در برابر , ,  و ؛
• ثبت trajectoryهای seed برای visualization و debugging.

.

راهنمای فایل‌های بسته بازتولیدپذیری

مسیر	کاربرد
src/mean_shift_gaussian.py	پیاده‌سازی Gaussian Mean Shift از صفر
examples/small_example.py	مثال عددی فصل
case_studies/iris_case_study.py	مطالعه موردی استاندارد
case_studies/noisy_imbalanced_case_study.py	مطالعه پیشرفته
tests/test_mean_shift_gaussian.py	۱۰ آزمون خودکار
outputs/*.csv	نتایج عددی قابل تحلیل
outputs/*.json	خلاصه‌های ماشینی نتایج
outputs/*.png	نمودارهای بازتولیدشده
requirements.txt	وابستگی‌های pin شده
environment_manifest.json	محیط و seedها
checksums.sha256	کنترل یکپارچگی فایل‌ها

منابع فنی این پیوست

scikit-learn Developers. MeanShift: https://scikit-learn.org/stable/modules/generated/sklearn.cluster.MeanShift.html

scikit-learn Developers. estimate_bandwidth: https://scikit-learn.org/stable/modules/generated/sklearn.cluster.estimate_bandwidth.html

scikit-learn Developers. Iris dataset API: https://scikit-learn.org/stable/modules/generated/sklearn.datasets.load_iris.html

Comaniciu, D., & Meer, P. (2002). Mean shift: A robust approach toward feature space analysis. IEEE TPAMI, 24(5), 603–619.