Gaussian mixture আর EM algorithm — নরম clustering

AI-এর ভিত্তি কোর্সের পর্ব ৬৫-এ k-means হাতে আর NumPy-তে, k-means++, elbow আর silhouette এসেছে। সেখানে প্রতিটি বিন্দু ঠিক একটি দলে পড়ে। এই কোর্সের পর্ব ১৮ শিখিয়েছে likelihood সবচেয়ে বড় করে parameter বের করা, আর পর্ব ১৯ শিখিয়েছে multivariate Gaussian-এর density, উপবৃত্তাকার contour আর covariance-এর MLE. এই পর্ব তিনটিকে জোড়া দেয়: কয়েকটি Gaussian মিশিয়ে ডেটার একটি model বানানো, আর লুকানো দল-পরিচয় না জেনেও তার parameter শেখা। চালানো যায় এমন চারটি ব্লক পাতাতেই চলে (ব্রাউজারের Python: NumPy, SciPy, matplotlib আর scikit-learn সহ)। প্রথমবার প্যাকেজ নামতে কয়েক সেকেন্ড লাগে, scikit-learn-এ প্রায় ১০ সেকেন্ড।

১ · শক্ত ভাগ না নরম ভাগ — k-means থেকে Gaussian mixture

Unsupervised learning-এ label নেই, তাই দলগুলো ডেটা থেকেই খুঁজতে হয়। দুই ধরনের উত্তর হতে পারে। Hard clustering প্রতিটি বিন্দুকে ঠিক একটি দলে বসায়, যেমন k-means. Soft clustering প্রতিটি বিন্দুর জন্য সব দলের সম্ভাবনা দেয়, যার যোগফল ১. দুই দলের মাঝখানের বিন্দু তখন নিজেই জানিয়ে দেয় যে সে অনিশ্চিত।

k-means চুপচাপ ধরে নেয় যে প্রতিটি দল গোল আর সব দলের ছড়ানো সমান। বাজারঘরে এটি প্রায়ই খাটে না: ছোট বাস্কেটের ছাত্রছাত্রীরা আঁটসাঁট দল, পাইকারি ক্রেতারা ছড়ানো, হেলানো আর সংখ্যায় কম। GMM প্রতিটি দলকে নিজের গড় μ_k, নিজের covariance matrix Σ_k আর নিজের অনুপাত π_k দেয়। তাই দলগুলো হেলানো উপবৃত্ত, ভিন্ন আকার আর ভিন্ন আকৃতির হতে পারে।

k-means আসলে GMM-এর একটি সীমা। ধরুন সব component-এর covariance একই আর গোল, Σ_k = εI, আর সব π_k সমান। তাহলে বিন্দু x_i-এর জন্য component k-এর responsibility (অংশ ২-এ এর সংজ্ঞা আসবে):

r_ik = exp(−‖x_i − μ_k‖² / 2ε) ÷ ∑_j exp(−‖x_i − μ_j‖² / 2ε)

এটি দূরত্বের বর্গের ওপর একটি softmax, আর ε তার “তাপমাত্রা”। ε বড় হলে সব r প্রায় সমান, ভাগ খুব নরম। ε → ০ হলে নিকটতম কেন্দ্রের r → ১ আর বাকিগুলো → ০. Σ আর π স্থির রেখে শুধু গড় হালনাগাদ করলে M-step-এর ওজন দেওয়া গড় হয়ে যায় দলের সাধারণ গড়, আর EM হুবহু k-means-এর “বসানো, তারপর গড় নেওয়া” loop. তাই k-means-কে বলা যায় শক্ত, গোল আর সমান মাপের GMM.

২ · Mixture density আর লুকানো চলক — ডেটা তৈরির গল্প

GMM-এর density হলো K-টি Gaussian-এর ওজন দেওয়া যোগফল:

p(x) = π₁N(x; μ₁, Σ₁) + … + π_KN(x; μ_K, Σ_K), π_k ≥ ০, ∑ π_k = ১

এখানে N(x; μ, Σ) হলো পর্ব ১৯-এর multivariate Gaussian density, আর π_k-কে বলে mixing weight. ওজনগুলোর যোগফল ১ বলে p(x)-এর মোট integral-ও ১, অর্থাৎ এটি একটি বৈধ density. সূত্রটিকে ডেটা তৈরির একটি গল্প (generative story) হিসেবে পড়া যায়: প্রতিটি ক্রেতার জন্য প্রথমে একটি পাশা চালা হয়, যার k নম্বর দিক পড়ার সম্ভাবনা π_k; ফল z হলো ক্রেতার দল। তারপর সেই দলের Gaussian থেকে x আসে। আমরা x দেখি, z দেখি না। এমন চলককে বলে latent variable, বা লুকানো চলক। GMM তাই একটি generative model, যা থেকে নতুন নমুনাও বানানো যায়।

Responsibility আসে Bayes-এর নিয়ম থেকে। z-এর prior P(z = k) = π_k, আর likelihood p(x | z = k) = N(x; μ_k, Σ_k). তাই x দেখার পর z-এর posterior (AI-এর ভিত্তি কোর্সের পর্ব ৫০):

r_ik = P(z_i = k | x_i) = π_kN(x_i; μ_k, Σ_k) ÷ ∑_j π_jN(x_i; μ_j, Σ_j)

প্রতিটি সারির r-এর যোগফল ১. লবটি “দল k-এর ক্রেতা হওয়া আর ঠিক এই x দেখা”-র ঘনত্ব, হরটি সব দল মিলিয়ে x দেখার ঘনত্ব, অর্থাৎ p(x). বাজারঘরের উদাহরণে (অংশ ৫-এর fit থেকে): ১,৫০০ টাকার অর্ডারে সাধারণ দিনের দলের ঘনত্ব খুব কম, কারণ সেটি গড় ৮৮১ থেকে প্রায় তিন standard deviation দূরে। উৎসবের দলের ঘনত্ব বেশি, যদিও ওই দলের ওজন মাত্র ০.৩২. ফল r ≈ (০.১৬, ০.৮৪).

Parameter কতগুলো? d মাত্রায় K-টি full-covariance component-এ গড়ের জন্য Kd-টি, covariance-এর জন্য K·d(d + ১)/২টি, আর ওজনের জন্য K − ১টি সংখ্যা (শেষ ওজনটি বাকিগুলো থেকে আসে)। d = ২, K = ৩ হলে ৬ + ৯ + ২ = ১৭টি। d = ২০ হলে প্রতিটি covariance একাই ২১০টি সংখ্যা চায়। অংশ ৭-এর covariance type এই সংখ্যা কমানোর উপায়।

৩ · সরাসরি MLE কেন কঠিন — যোগফলের log

parameter θ = (π, μ, Σ) শেখার স্বাভাবিক পথ maximum likelihood: n-টি বিন্দুর log-likelihood সবচেয়ে বড় করা।

ℓ(θ) = ∑_{i = 1..n} log ( ∑_{k = 1..K} π_kN(x_i; μ_k, Σ_k) )

একটি Gaussian হলে (পর্ব ১৮) log সরাসরি exp-কে কেটে দেয়, ℓ হয় μ-এর একটি দ্বিঘাত রাশি, আর gradient শূন্য করলেই নমুনার গড়। Mixture-এ log বসে আছে একটি যোগফলের বাইরে। নিচের stepper দেখায় এই একটি পার্থক্য কেন সব বদলে দেয়।

আরও দুটি সমস্যা আছে। প্রথমত ℓ concave নয় (পর্ব ১২-এর ভাষায় non-convex সমস্যা)। label switching-এর জন্যই অন্তত K!-টি সমান উঁচু চূড়া থাকে, আর তার বাইরে নিচু চূড়া আর saddle point-ও থাকে। দ্বিতীয়ত, ℓ-এর কোনো সর্বোচ্চ মানই নেই: একটি component-এর গড় একটি বিন্দুর ওপর বসিয়ে তার σ শূন্যের দিকে নিলে ℓ → ∞ (অংশ ৮)। তাই “GMM-এর MLE” বলতে বাস্তবে বোঝায় একটি যুক্তিসংগত local maximum.

৪ · EM-এর দুই ধাপ — E-step আর M-step

Expectation–maximization (EM) stepper-এর শেষ ধাপের কৌশলটিকে একটি নিয়মিত loop বানায়। বর্তমান parameter θ^(t) থেকে শুরু করে প্রতি পাকে:

M-step-এর সূত্রগুলো পর্ব ১৯-এর MLE-রই ওজন দেওয়া রূপ: গড় হলো ওজন দেওয়া গড়, covariance হলো ওজন দেওয়া বর্গ-বিচ্যুতির গড় (n − ১ নয়, N_k দিয়ে ভাগ)। π_k-এর সূত্রে একটি শর্ত আছে, ∑ π_k = ১. পর্ব ১৬-এর Lagrange multiplier দিয়ে ∑_k N_k log π_k + λ(১ − ∑ π_k)-এর derivative শূন্য করলে N_k/π_k = λ, আর যোগফলের শর্ত থেকে λ = n. তাই π_k = N_k/n, অর্থাৎ দল k-এর ওপর মোট responsibility-র অংশ। M-step-এ শর্তগুলো আপনা-আপনি রক্ষা পায়: π-এর যোগফল ১, আর Σ_k ওজন দেওয়া outer product-এর যোগফল বলে PSD (পর্ব ১০)।

পরীক্ষায় প্রায়ই এক মাত্রায় একটি পাক হাতে করতে বলা হয়। নিচের stepper চারটি বিন্দুতে তা-ই করে। সংখ্যাগুলো একটি scratch script-এ যাচাই করা।

কখন থামবেন? সাধারণ নিয়ম: পরপর দুই পাকে ℓ-এর পরিবর্তন একটি ছোট সীমা tol-এর নিচে নামলে। scikit-learn-এর GaussianMixture প্রতি নমুনার গড় log-likelihood দেখে, ডিফল্ট tol=1e-3 আর max_iter=100. এই tol বেশ ঢিলা; সূক্ষ্ম তুলনার সময় (অংশ ৬ আর ৭-এর মতো) একে ছোট করা ভালো।

৫ · NumPy-তে EM — আর log-likelihood কেন কখনো কমে না

এবার বাজারঘরের ৩০০টি কাল্পনিক অর্ডার: ২১০টি সাধারণ দিনের (গড় ৯০০ টাকা, sd ২৫০), ৯০টি উৎসবের (গড় ২,৪০০, sd ৬০০)। Model জানে না কোনটি কোন দলের। শুরুর অনুমান ইচ্ছা করে দূরে রাখা হয়েছে: গড় ৫০০ আর ৩,০০০, দুটিরই sd ৮০০. প্রতিটি লাইনে আগে সেই পাকের শুরুর parameter-এ ℓ, তীরের পরে M-step-এর নতুন parameter.

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
rng = np.random.default_rng(42)            # কাল্পনিক বাজারঘর: অর্ডারের মূল্য (টাকা)
x = np.concatenate([rng.normal(900, 250, 210),     # সাধারণ দিনের ক্রেতা
                    rng.normal(2400, 600, 90)])    # উৎসবের ক্রেতা
pi = np.array([0.5, 0.5])                  # শুরুর অনুমান, ইচ্ছা করে দূরে
mu = np.array([500.0, 3000.0])
sd = np.array([800.0, 800.0])
lls = []
for it in range(30):
    w = pi * norm.pdf(x[:, None], mu, sd)  # n × 2: πₖ·N(xᵢ; μₖ, σₖ)
    lls.append(np.log(w.sum(axis=1)).sum())          # এই θ-তে log-likelihood
    r = w / w.sum(axis=1, keepdims=True)   # E-step: responsibility
    Nk = r.sum(axis=0)                     # M-step: ওজন দেওয়া গণনা
    pi = Nk / len(x)
    mu = (r * x[:, None]).sum(axis=0) / Nk
    sd = np.sqrt((r * (x[:, None] - mu) ** 2).sum(axis=0) / Nk)
    if it in (0, 1, 2, 4, 9, 29):
        print(f"iter {it + 1:2d}: loglik {lls[-1]:9.2f} -> pi {pi.round(3)}"
              f" mu {mu.round()} sd {sd.round()}")
print("loglik never fell:", bool(np.all(np.diff(lls) >= -1e-9)))
for v in (1200, 1500, 1800):               # মাঝামাঝি অর্ডার: নরম ভাগ
    p = pi * norm.pdf(v, mu, sd)
    print(f"order {v} tk: r = {(p / p.sum()).round(3)}")
g = np.linspace(0, 5000, 400)
fig, ax = plt.subplots(figsize=(6, 3.6), layout="constrained")
ax.hist(x, bins=30, density=True, alpha=0.4, label="orders")
for k in (0, 1):
    ax.plot(g, pi[k] * norm.pdf(g, mu[k], sd[k]), "--", label=f"component {k + 1}")
ax.plot(g, (pi * norm.pdf(g[:, None], mu, sd)).sum(axis=1), "k", label="mixture")
ax.set_xlabel("order value (taka)"); ax.set_ylabel("density")
ax.set_title("Two-component GMM fitted by EM"); ax.legend(fontsize=8)
plt.show()

iter  1: loglik  -2536.12 -> pi [0.727 0.273] mu [ 950. 2363.] sd [336. 702.]
iter  2: loglik  -2325.16 -> pi [0.716 0.284] mu [ 914. 2398.] sd [257. 642.]
iter  3: loglik  -2303.21 -> pi [0.705 0.295] mu [ 895. 2388.] sd [228. 618.]
iter  5: loglik  -2297.90 -> pi [0.694 0.306] mu [ 886. 2355.] sd [216. 630.]
iter 10: loglik  -2297.51 -> pi [0.684 0.316] mu [ 882. 2318.] sd [213. 655.]
iter 30: loglik  -2297.50 -> pi [0.681 0.319] mu [ 881. 2306.] sd [212. 663.]
loglik never fell: True
order 1200 tk: r = [0.896 0.104]
order 1500 tk: r = [0.162 0.838]
order 1800 tk: r = [0.001 0.999]

কেন ℓ কখনো কমে না: একটি নিম্নসীমা। লুকানো z-এর ওপর যেকোনো বণ্টন q(z) নিন। log একটি concave function, তাই Jensen's inequality (পর্ব ১২) বলে “গড়ের log ≥ log-এর গড়”। এ থেকে প্রতিটি বিন্দুর জন্য:

log p(x | θ) = log ∑_z q(z)·[p(x, z | θ) / q(z)] ≥ ∑_z q(z) log [p(x, z | θ) / q(z)] = ELBO(q, θ)

ডান দিকের রাশির নাম ELBO, অর্থাৎ evidence lower bound. দুই দিকের ফাঁকটুকু ঠিক KL divergence KL(q ‖ p(z | x, θ)), যা কখনো ঋণাত্মক নয় আর শূন্য হয় শুধু q যখন আসল posterior (AI-এর ভিত্তি কোর্সের পর্ব ৩৮)। এখন EM-কে এভাবে পড়ুন। E-step-এ q-কে posterior বানানো হয়, অর্থাৎ q(z_i = k) = r_ik. KL শূন্য হয়, আর নিম্নসীমা বর্তমান θ-তে ℓ-কে ছুঁয়ে ফেলে। M-step-এ q স্থির রেখে θ বদলে ELBO বাড়ানো হয়। ELBO-র θ-নির্ভর অংশই complete-data log-likelihood-এর প্রত্যাশা, তাই বন্ধ সূত্রগুলো আসে। নতুন θ-তে ℓ ≥ নতুন ELBO ≥ পুরোনো ELBO = পুরোনো ℓ. তাই প্রতি পাকে ℓ বাড়ে বা একই থাকে।

তিনটি ফল মনে রাখুন। (১) EM-এ learning rate নেই, আর ℓ কখনো কমে না। (২) নিশ্চয়তা শুধু stationary বিন্দুর, সাধারণত একটি local maximum (অংশ ৭)। (৩) শেষের দিকে গতি ধীর, বিশেষ করে component-গুলো মিশে থাকলে। ব্লকের ১০ থেকে ৩০ নম্বর পাকে ℓ মাত্র ০.০১ বাড়ে, অথচ μ₂ তখনও ২,৩১৮ থেকে ২,৩০৬-এ সরছে। সমতল চূড়ায় ℓ প্রায় স্থির থাকলেও parameter নড়তে পারে, তাই খুব ঢিলা tol বিপজ্জনক। ELBO-র এই দৃষ্টিভঙ্গি পর্ব ৪৮ (Variational inference-এর ধারণা — ELBO আর approximation)-এর ভিত্তি। সেখানে q-কে আর নিখুঁত posterior বানানো যায় না, একটি সরল পরিবার থেকে বেছে নিতে হয়।

৬ · দুই মাত্রায় full covariance — নিজের EM বনাম GaussianMixture

বাজারঘরের ৩৬০ জন কাল্পনিক ক্রেতার দুটি feature: গড় বাস্কেটের মূল্য (শত টাকায়) আর মাসে অর্ডারের সংখ্যা। তিনটি লুকানো দল: ছোট বাস্কেটের ১৫০ জন, পরিবারের ১৫০ জন, আর ৬০ জন পাইকারি ক্রেতা, যাঁদের ছড়ানো বেশি আর হেলানো উল্টো দিকে। নিচের ব্লক EM নিজে লেখে (covariance-এর M-step একটি einsum-এ), তারপর scikit-learn-কে একই শুরুর মান দেয়: weights_init, means_init আর precisions_init (covariance-এর inverse দিয়ে)। একই শুরু, একই tol হলে দুটির ফল মেলা উচিত।

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal as mvn
from sklearn.mixture import GaussianMixture
rng = np.random.default_rng(42)            # (গড় বাস্কেট শত টাকা, মাসে অর্ডার)
X = np.vstack([rng.multivariate_normal([6, 3], [[4, 1], [1, 1]], 150),
               rng.multivariate_normal([15, 7], [[9, 4], [4, 4]], 150),
               rng.multivariate_normal([28, 4.5], [[30, -6], [-6, 3]], 60)])
K, n = 3, len(X); mu = X[rng.choice(n, K, replace=False)]  # শুরু: এলোমেলো ৩ সারি
cov = np.array([np.cov(X.T)] * K); pi = np.full(K, 1 / K); old = -np.inf
init = dict(weights_init=pi.copy(), means_init=mu.copy(), precisions_init=np.linalg.inv(cov))
for it in range(1, 501):
    w = np.column_stack([pi[k] * mvn.pdf(X, mu[k], cov[k]) for k in range(K)])
    ll = np.log(w.sum(axis=1)).mean()      # গড় log-likelihood
    r = w / w.sum(axis=1, keepdims=True)   # E-step
    Nk = r.sum(axis=0); pi = Nk / n        # M-step
    mu = (r.T @ X) / Nk[:, None]
    D = X[:, None, :] - mu                 # n × K × 2
    cov = np.einsum("nk,nki,nkj->kij", r, D, D) / Nk[:, None, None] + 1e-6 * np.eye(2)
    if ll - old < 1e-6: break
    old = ll
gm = GaussianMixture(K, init_params="random_from_data", tol=1e-6, **init).fit(X)
print("iterations:", it, gm.n_iter_, "| loglik/n:", round(ll, 6), round(gm.lower_bound_, 6))
print("weights:", pi.round(3), gm.weights_.round(3))
print("means:\n", mu.round(2), "\n", gm.means_.round(2))
print("max |cov diff|:", f"{np.abs(cov - gm.covariances_).max():.1e}")
fig, ax = plt.subplots(figsize=(6, 4), layout="constrained")
ax.scatter(X[:, 0], X[:, 1], c=gm.predict(X), s=8, cmap="viridis")
t = np.linspace(0, 2 * np.pi, 100)
for k in range(K):                         # Mahalanobis 2-এর উপবৃত্ত (পর্ব ১৯)
    lam, V = np.linalg.eigh(gm.covariances_[k])
    e = gm.means_[k][:, None] + V @ (2 * np.sqrt(lam)[:, None] * [np.cos(t), np.sin(t)])
    ax.plot(e[0], e[1], "r")
ax.set_xlabel("avg basket (100 taka)"); ax.set_ylabel("orders per month")
ax.set_title("GMM, full covariance: 2-sd ellipses"); plt.show()

iterations: 33 33 | loglik/n: -4.890552 -4.890552
weights: [0.405 0.432 0.163] [0.405 0.432 0.163]
means:
 [[ 6.02  2.95]
 [14.88  6.96]
 [28.45  4.07]]
 [[ 6.02  2.95]
 [14.88  6.96]
 [28.45  4.07]]
max |cov diff|: 2.8e-14

উপবৃত্তগুলো কোথা থেকে। প্রতিটি component-এর Σ_k-কে eigh দিয়ে ভাঙলে অক্ষের দিক আসে eigenvector থেকে আর অর্ধ-অক্ষের দৈর্ঘ্য c√λ (পর্ব ১৯-এর অংশ ২)। এখানে c = ২, অর্থাৎ Mahalanobis দূরত্ব ২-এর contour. দুই মাত্রায় এর ভেতরে প্রায় ৮৬.৫% থাকার কথা। ছবিতে পাইকারি দলের উপবৃত্ত সবচেয়ে বড় আর অন্য দুটির উল্টো দিকে হেলানো: বাস্কেট যত বড়, অর্ডার তত কম। k-means এমন আকৃতি ধরতে পারে না।

৭ · Covariance type, BIC/AIC দিয়ে K বাছাই, আর restart

Full covariance সবচেয়ে নমনীয়, কিন্তু parameter-ও সবচেয়ে বেশি। covariance_type দিয়ে Σ_k-এর ওপর গঠন বসানো যায়, যাকে বলে covariance type:

শেষ কলামের “+ 6 + 2” হলো গড়ের Kd = ৬ আর ওজনের K − ১ = ২. এবার K কত? ℓ দিয়ে বাছাই করা যায় না, কারণ component বাড়ালে training ডেটায় ℓ প্রায় সবসময় বাড়ে। তাই parameter-এর জন্য শাস্তি যোগ করা হয়:

BIC = −2ℓ + p·ln n, AIC = −2ℓ + 2p

এখানে ℓ fit-এর সর্বোচ্চ log-likelihood, p parameter-এর সংখ্যা, n সারির সংখ্যা; দুটিই যত কম, তত ভালো। AIC (Hirotugu Akaike, ১৯৭৩–৭৪) ভবিষ্যৎ ডেটায় অনুমানের মান লক্ষ্য করে, BIC (Gideon Schwarz, ১৯৭৮) আসে Bayesian যুক্তি থেকে। n ≥ ৮ হলেই ln n > ২, তাই BIC কম component বেছে নেয়। scikit-learn-এ gm.bic(X) আর gm.aic(X) ঠিক এই সূত্র। নিচে অংশ ৫-এর অর্ডার-ডেটায় K = ১ বনাম K = ২ হাতে তুলনা।

নিচের ব্লক অংশ ৬-এর ক্রেতা-ডেটায় চারটি covariance type আর K = ১ থেকে ৬ পর্যন্ত সব জোড়া fit করে BIC আর AIC তুলনা করে। তারপর K = ৩-এ শুধু এক শুরু থেকে (n_init=1) ছয়টি ভিন্ন seed চালায়, যাতে local maximum চোখে পড়ে।

import numpy as np
import matplotlib.pyplot as plt
from sklearn.mixture import GaussianMixture
rng = np.random.default_rng(42)            # একই কাল্পনিক ক্রেতা-ডেটা
X = np.vstack([rng.multivariate_normal([6, 3], [[4, 1], [1, 1]], 150),
               rng.multivariate_normal([15, 7], [[9, 4], [4, 4]], 150),
               rng.multivariate_normal([28, 4.5], [[30, -6], [-6, 3]], 60)])
Ks, types = range(1, 7), ["full", "tied", "diag", "spherical"]
opt = dict(tol=1e-5, max_iter=500, random_state=0)
full = [GaussianMixture(k, init_params="random_from_data", n_init=3, **opt).fit(X) for k in Ks]
bic, aic = {}, {}
for t in types:                            # শুরু: full fit-এর গড় (নিচের সতর্কতা)
    fits = [GaussianMixture(k, covariance_type=t, means_init=f.means_, init_params="random",
                            **opt).fit(X)
            for k, f in zip(Ks, full)]
    bic[t] = [g.bic(X) for g in fits]; aic[t] = [g.aic(X) for g in fits]
    print(f"{t:9s} best K: BIC {np.argmin(bic[t]) + 1}, AIC {np.argmin(aic[t]) + 1}"
          f" | BIC at K=3 {bic[t][2]:6.0f}")
print("full BIC K=1..6:", np.round(bic["full"]).astype(int))
print("full AIC K=1..6:", np.round(aic["full"]).astype(int))
for seed in range(6):                      # n_init=1: একেক শুরু, একেক শেষ
    g = GaussianMixture(3, init_params="random_from_data", tol=1e-6,
                        max_iter=500, random_state=seed).fit(X)
    print(f"seed {seed}: loglik/n {g.lower_bound_:.4f}  weights {np.sort(g.weights_).round(3)}")
fig, ax = plt.subplots(figsize=(6, 3.6), layout="constrained")
for t in types:
    ax.plot(Ks, bic[t], "o-", label=t)
ax.set_xlabel("number of components K"); ax.set_ylabel("BIC (lower is better)")
ax.set_title("BIC by covariance type"); ax.legend(fontsize=8); plt.show()

full      best K: BIC 3, AIC 6 | BIC at K=3   3621
tied      best K: BIC 6, AIC 6 | BIC at K=3   3754
diag      best K: BIC 5, AIC 6 | BIC at K=3   3716
spherical best K: BIC 6, AIC 6 | BIC at K=3   3901
full BIC K=1..6: [4204 3719 3621 3648 3669 3685]
full AIC K=1..6: [4184 3676 3555 3558 3557 3549]
seed 0: loglik/n -4.8906  weights [0.163 0.405 0.432]
seed 1: loglik/n -4.8906  weights [0.163 0.405 0.432]
seed 2: loglik/n -4.8906  weights [0.163 0.405 0.432]
seed 3: loglik/n -5.0455  weights [0.006 0.208 0.787]
seed 4: loglik/n -4.8906  weights [0.163 0.406 0.431]
seed 5: loglik/n -5.0467  weights [0.037 0.179 0.784]

~/bazarghar-ml$ python tied_init.py
python 3.14.4 | numpy 2.5.4 | scikit-learn 1.9.1
k-means++        means [[13.5, 4.9], [13.6, 4.9], [13.5, 4.9]]  loglik/n -5.798
random_from_data means [[13.5, 4.9], [13.5, 4.9], [13.5, 4.9]]  loglik/n -5.798
kmeans           means [[16.3, 7.3], [30.7, 3.5], [6.7, 3.2]]  loglik/n -5.124
K=1 for reference means [[13.5, 4.9]]  loglik/n -5.798

Local maximum থেকে বাঁচার সাধারণ পথ restart: n_init-টি ভিন্ন শুরু থেকে EM চালিয়ে যে fit-এর ℓ সবচেয়ে বেশি, সেটি রাখা। scikit-learn-এ ডিফল্ট n_init=1, অর্থাৎ restart নেই। BIC বা AIC হিসাব করার আগে প্রতিটি K-তে যথেষ্ট restart দরকার। নইলে একটি খারাপ local maximum একটি ভালো K-কে খারাপ দেখাতে পারে।

৮ · নরম সিদ্ধান্ত, ভেঙে পড়া component আর reg_covar

predict_proba প্রতিটি নতুন ক্রেতার responsibility দেয়, predict শুধু সবচেয়ে বড়টির নম্বর। দুই দলের মাঝের ক্রেতার বেলায় এই পার্থক্য সিদ্ধান্ত বদলায়। ব্লকের দ্বিতীয় অংশ উল্টো দিকের একটি বিপদ দেখায়। বাজারঘরের অর্ডার-তালিকায় ৫টি হুবহু এক দামের (৯,৯৯৯ টাকা) “বান্ডেল” অর্ডার আছে। একটি component সেগুলোর ওপর বসে পড়লে কী হয়?

import numpy as np
from scipy.stats import norm
from sklearn.mixture import GaussianMixture
rng = np.random.default_rng(42)            # একই কাল্পনিক ক্রেতা-ডেটা
X = np.vstack([rng.multivariate_normal([6, 3], [[4, 1], [1, 1]], 150),
               rng.multivariate_normal([15, 7], [[9, 4], [4, 4]], 150),
               rng.multivariate_normal([28, 4.5], [[30, -6], [-6, 3]], 60)])
gm = GaussianMixture(3, init_params="random_from_data", n_init=3, tol=1e-6,
                     max_iter=500, random_state=0).fit(X)
order = np.argsort(gm.means_[:, 0])        # দল: ছোট, মাঝারি, বড় বাস্কেট
new = np.array([[6, 3], [10, 4.5], [11, 5.5], [21, 5.5], [24, 5]])
P = gm.predict_proba(new)[:, order]
for x, p in zip(new, P):
    print(f"customer {x}: soft {p.round(3)}  hard -> group {p.argmax() + 1}")

x = np.r_[rng.normal(900, 250, 40), [9999.0] * 5]   # ৫টি হুবহু এক দামের বান্ডেল
pi, mu, sd = np.array([0.9, 0.1]), np.array([900.0, 8000.0]), np.array([250.0, 1500.0])
for it in range(1, 7):                     # reg ছাড়া নিজের EM
    w = pi * norm.pdf(x[:, None], mu, sd)
    print(f"iter {it}: loglik {np.log(w.sum(axis=1)).sum():9.2f}  sd2 {sd[1]:.3g}")
    r = w / w.sum(axis=1, keepdims=True); Nk = r.sum(axis=0)
    pi, mu = Nk / len(x), (r * x[:, None]).sum(axis=0) / Nk
    sd = np.sqrt((r * (x[:, None] - mu) ** 2).sum(axis=0) / Nk)
    if sd[1] < 1e-6: print(f"sd2 -> {sd[1]:.3g}: density at 9999 -> infinity"); break
Xd = x[:, None]
for reg in (0.0, 1e-6, 100.0):
    try:
        g = GaussianMixture(2, reg_covar=reg, means_init=[[900], [8000]],
                            init_params="random_from_data", random_state=0).fit(Xd)
        sd = np.sqrt(g.covariances_.ravel())
        print(f"reg_covar={reg}: sd {sd[0]:.1f}, {sd[1]:.3g}  loglik {g.score(Xd) * len(x):.1f}")
    except ValueError as e:
        print(f"reg_covar={reg}: ValueError:", str(e)[:60], "...")

customer [6. 3.]: soft [0.99 0.01 0.  ]  hard -> group 1
customer [10.   4.5]: soft [0.45 0.55 0.  ]  hard -> group 2
customer [11.   5.5]: soft [0.039 0.96  0.001]  hard -> group 2
customer [21.   5.5]: soft [0.    0.044 0.956]  hard -> group 3
customer [24.  5.]: soft [0. 0. 1.]  hard -> group 3
iter 1: loglik   -343.30  sd2 1.5e+03
iter 2: loglik   -317.75  sd2 24.2
sd2 -> 0: density at 9999 -> infinity
reg_covar=0.0: sd 275.7, 3.64e-12  loglik -172.6
reg_covar=1e-06: sd 275.7, 0.001  loglik -267.3
reg_covar=100.0: sd 275.9, 10  loglik -313.3

কেন এমন হয়। এক মাত্রায় একটি বিন্দু x-এর ওপর বসা component-এর density N(x; x, σ²) = ১/(σ√(২π)), যা σ → ০ হলে অসীম। বাকি বিন্দুগুলোকে অন্য component সামলায়, তাই মোট ℓ সীমাহীনভাবে বাড়ে। এটি likelihood singularity, আর এই কারণেই GMM-এর “global MLE” অর্থহীন। হুবহু এক মান বারবার থাকলে (বান্ডেলের দাম, এন্ট্রির ডিফল্ট মান) সমস্যা সহজে আসে। d মাত্রায় d-টির কম বিন্দু একটি component-এ পড়লেই তার Σ singular হয়।

প্রতিকার। (১) Covariance regularization: reg_covar প্রতিটি Σ_k-এর কর্ণে একটি ধ্রুবক যোগ করে, ফলে variance তার নিচে নামতে পারে না। ডিফল্ট 1e-6 মাপা হয় ডেটার এককে, টাকা²-এ। ৯,৯৯৯ টাকার উদাহরণে তা কার্যত কোনো বাধা নয়। তাই আগে feature scale করুন (যেমন StandardScaler), তারপর reg_covar বাছুন। (২) হুবহু এক মানগুলো সত্যিই আলাদা দল, না ভুল, দেখুন। (৩) খুব ছোট ওজনের component (π < ০.০১) সন্দেহ করুন। (৪) Covariance-এর ওপর prior বসালে MAP অনুমান (পর্ব ১৮) singularity এড়ায়।

Bayesian GMM, একটি নোট। scikit-learn-এর BayesianGaussianMixture ওজন, গড় আর covariance-এর ওপর prior বসায় আর fit করে variational inference দিয়ে, EM দিয়ে নয়। ডিফল্ট weight_concentration_prior_type="dirichlet_process". Dirichlet process prior-এর সুবিধা: n_components একটু বেশি দিলে অপ্রয়োজনীয় component-এর ওজন প্রায় শূন্যে নেমে যেতে পারে। তবে ফল prior-এর মানের ওপর নির্ভর করে। এই Bayesian Gaussian mixture-এর পেছনের পদ্ধতি পর্ব ৪৮-এ (Variational inference-এর ধারণা) আসবে।

EM শুধু GMM-এর জন্য নয়। লুকানো চলক জানলে যেখানে fit সহজ, সেখানেই EM খাটে: HMM-এর Baum–Welch (AI-এর ভিত্তি কোর্সের পর্ব ৫৩-এর HMM), missing value-সহ Gaussian, factor analysis আর অনেক topic model. পরের পর্বগুলোতে GMM আবার আসবে: পর্ব ৪০ (Density estimation — histogram, KDE আর mixture model)-এ parametric density হিসেবে, আর পর্ব ৪৩ (Anomaly detection — isolation forest, LOF আর one-class SVM)-এ score_samples দিয়ে কম-ঘনত্বের বিন্দু খোঁজায়।

৯ · AI-এর যুগে — GMM আর EM-এর কোড যাচাই

AI assistant-এর GMM কোড দেখতে বিশ্বাসযোগ্য, কিন্তু ভুল লুকিয়ে থাকে parameter-এর নাম, ডিফল্ট মান আর K বাছাইয়ের যুক্তিতে:

AI-এর দেওয়া GMM কোড হাতে পেলে — যাচাইয়ের চেকলিস্ট

• প্রতিটি parameter-এর নাম আর ডিফল্ট মান নিজের scikit-learn version-এর docs-এ মিলিয়েছি (n_init, init_params, tol, reg_covar)
• feature scale করেছি, আর K বেছেছি BIC/AIC বা held-out log-likelihood দিয়ে, training ℓ দিয়ে নয়
• কয়েকটি seed-এ চালিয়ে lower_bound_, ওজন আর গড় তুলনা করেছি; কোনো component-এর ওজন প্রায় শূন্য বা গড় মিশে গেছে কি না দেখেছি
• নিজে লেখা EM-কে ছোট seed দেওয়া উদাহরণে GaussianMixture-এর সঙ্গে একই শুরু থেকে মিলিয়েছি, আর ℓ কখনো কমেনি তা পরীক্ষা করেছি
• segment-গুলোর গড় প্রোফাইল দেখে ব্যবসায়িক মানে যাচাই করেছি, আর একটি সরল baseline (যেমন k-means বা নিয়মভিত্তিক ভাগ) এর সঙ্গে তুলনা করেছি
• কোনো আসল গ্রাহকের ব্যক্তিগত ডেটা AI-তে পেস্ট করিনি