Multivariate Gaussian — covariance, Mahalanobis দূরত্ব আর শর্তসাপেক্ষ বণ্টন

AI-এর ভিত্তি কোর্সের পর্ব ৩৬-এ এক মাত্রার normal distribution আর Monte Carlo, পর্ব ৩৭-এ covariance, correlation আর z-score এসেছে। এই কোর্সের পর্ব ০৮ প্রতিসম matrix-এর eigen-decomposition দেখিয়েছে, আর পর্ব ১০ দেখিয়েছে covariance কেন PSD, Cholesky দিয়ে x = μ + Lz নমুনা কীভাবে বানাতে হয়, আর Mahalanobis দূরত্বের প্রথম পরিচয়। এই পর্ব সেই টুকরোগুলো জুড়ে পুরো ছবি আঁকে: একসঙ্গে অনেকগুলো সংখ্যার অনিশ্চয়তা কীভাবে একটি বণ্টনে ধরা যায়, আর তা দিয়ে outlier ধরা, অনুমান আর classification কীভাবে হয়। চালানো যায় এমন চারটি ব্লক পাতাতেই চলে (ব্রাউজারের Python: NumPy, SciPy, matplotlib আর scikit-learn); প্রথমবার প্যাকেজ নামতে কয়েক সেকেন্ড লাগে, scikit-learn-এ প্রায় ১০ সেকেন্ড।

১ · এক মাত্রার ঘণ্টা থেকে d মাত্রার পাহাড় — density সূত্র পদে পদে

এক মাত্রায় normal distribution-এর density p(x) = (1/(σ√(2π))) · exp(−(x − μ)²/(2σ²)). এখানে μ গড়, σ standard deviation, আর (x − μ)/σ হলো z-score. d মাত্রায় x একটি vector (দৈর্ঘ্য d), μ গড় vector, আর σ²-এর জায়গায় বসে d × d covariance matrix Σ. সূত্রটি:

p(x) = (2π)^−d/2 · (det Σ)^−1/2 · exp(−½ (x − μ)ᵀ Σ⁻¹ (x − μ))

তিনটি অংশ আলাদা করে পড়ুন। (১) সূচকের ভেতরের রাশি d² = (x − μ)ᵀΣ⁻¹(x − μ) একটি quadratic form, আর এটিই Mahalanobis দূরত্ব-এর বর্গ। এক মাত্রায় এটি ঠিক z² = (x − μ)²/σ². Σ⁻¹-কে বলে precision matrix: covariance যত বড়, ঘনত্ব তত ধীরে কমে। (২) det Σ হলো মেঘের “আয়তনের” মাপ। মেঘ যত ছড়ানো, মোট সম্ভাবনা ১ রাখতে চূড়া তত নিচু হতে হয়, তাই det Σ-এর বর্গমূল দিয়ে ভাগ। (৩) (2π)^d/2 শুধু একটি ধ্রুবক, যাতে পুরো জায়গায় ঘনত্বের integral ঠিক ১ হয়। সূত্রটি অর্থবহ হতে Σ-কে positive definite হতে হবে (পর্ব ১০-এর ভাষায়), নইলে Σ⁻¹ আর √det Σ সংজ্ঞায়িত নয়।

পর্ব ১০-এর বাজারঘরের ডেলিভারি আবার নিই: x = (সময়, দূরত্ব), μ = (৩৫ মিনিট, ৬ km), সময়ের standard deviation ১০, দূরত্বের ৩, correlation ০.৭, তাই Σ = [[100, 21], [21, 9]]. একটি অর্ডার (৪৫ মিনিট, ৮ km)-এ density কত?

দুটি সতর্কতা। প্রথমত, density সম্ভাবনা নয়: এর একক “প্রতি মিনিট প্রতি km”, আর ছোট Σ-এ এটি ১-এর বেশিও হতে পারে; সম্ভাবনা আসে কোনো এলাকার ওপর integral থেকে (AI-এর ভিত্তি কোর্সের পর্ব ৩৬-এর PDF-এর মতো)। দ্বিতীয়ত, বড় d-তে det Σ খুব ছোট বা বড় হয়ে underflow বা overflow করতে পারে, তাই বাস্তব কোডে log-density হিসাব করা হয়: log p = −½(d log 2π + log det Σ + d²), আর log det Σ = ২∑ log L_ii (Cholesky-র কর্ণ থেকে, পর্ব ১০)।

Gaussian-এর একটি বিশেষ গুণ মনে রাখুন: যৌথভাবে Gaussian হলে correlation শূন্য মানেই স্বাধীন. Σ কর্ণ-matrix হলে exp-এর ভেতরের যোগফল আলাদা আলাদা পদে ভাগ হয়, আর density এক মাত্রার density-গুলোর গুণফল। অন্য বণ্টনে এটি সত্য নয়: x আর x² correlation-হীন হতে পারে, অথচ একটি জানলে অন্যটি পুরো জানা যায়।

২ · Contour মানে উপবৃত্ত — অক্ষ আসে Σ-এর eigenvector থেকে

Density একই থাকে যেখানে d² একই। তাই সমান-উচ্চতার রেখা (contour) হলো d² = c² সমীকরণের সমাধান। পর্ব ০৮-এর spectral theorem অনুযায়ী Σ = VΛVᵀ, যেখানে V-এর কলামগুলো পরস্পর লম্ব একক eigenvector আর Λ-এর কর্ণে eigenvalue λ₁, λ₂, …. তাহলে Σ⁻¹ = VΛ⁻¹Vᵀ, আর y = Vᵀ(x − μ) লিখলে:

d² = y₁²/λ₁ + y₂²/λ₂ + … + y_d²/λ_d

এটি অক্ষ-বরাবর একটি উপবৃত্তের সমীকরণ। y হলো eigenvector-গুলোর দিকে মাপা স্থানাঙ্ক, তাই উপবৃত্তের অক্ষগুলো eigenvector বরাবর, আর d = c রেখায় i-তম অর্ধ-অক্ষের দৈর্ঘ্য c√λ_i. বড় eigenvalue মানে সেই দিকে ডেটা বেশি ছড়ানো, অক্ষ লম্বা। আরও এক ধাপ: z_i = y_i/√λ_i নিলে d² = ‖z‖², সাধারণ দৈর্ঘ্যের বর্গ। এই রূপান্তর z = Λ^−1/2Vᵀ(x − μ)-কে বলে whitening: হেলানো উপবৃত্তকে ঘুরিয়ে, প্রতিটি অক্ষ বরাবর চেপে বা টেনে, একটি গোল বৃত্ত বানানো। Mahalanobis দূরত্ব আসলে whitened জগতের সাধারণ দূরত্ব। (পর্ব ১০-এ Cholesky দিয়ে Ly = x − μ আরেক রকম whitening, যার ঘোরানোর ধাপ আলাদা কিন্তু দৈর্ঘ্য একই।)

উপবৃত্তের ভেতরে কত সম্ভাবনা? z-এর উপাদানগুলো স্বাধীন standard normal, তাই d² = z₁² + … + z_d², যার বণ্টনকে বলে d degree of freedom-এর chi-square, χ²_d. ২ মাত্রায় এর একটি সহজ বন্ধ রূপ আছে: P(d² ≤ c²) = 1 − e^−c²/2. নিচের ব্লক density সূত্র নিজে লিখে SciPy-র multivariate_normal.pdf-এর সঙ্গে মেলায়, তারপর eigenvector থেকে তিনটি উপবৃত্ত আঁকে।

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal

mu = np.array([35.0, 6.0])                 # গড়: ৩৫ মিনিট, ৬ km
Sigma = np.array([[100.0, 21.0], [21.0, 9.0]])
def gauss_pdf(x, mu, Sigma):               # সূত্রটি হুবহু, পদে পদে
    d = len(mu)
    diff = x - mu
    m2 = diff @ np.linalg.solve(Sigma, diff)          # Mahalanobis²
    norm = (2 * np.pi) ** (d / 2) * np.sqrt(np.linalg.det(Sigma))
    return np.exp(-0.5 * m2) / norm, m2

for x in ([45.0, 8.0], [35.0, 6.0], [55.0, 3.0]):
    p, m2 = gauss_pdf(np.array(x), mu, Sigma)
    q = multivariate_normal.pdf(x, mean=mu, cov=Sigma)
    print(f"x={x}: d²={m2:6.3f}  formula={p:.6f}  scipy={q:.6f}")
lam, V = np.linalg.eigh(Sigma)             # ছোট থেকে বড় eigenvalue
print("eigenvalues:", lam.round(3), "| sqrt:", np.sqrt(lam).round(3))
print("long axis direction:", V[:, 1].round(3))

t = np.linspace(0, 2 * np.pi, 200)
circle = np.stack([np.cos(t), np.sin(t)])  # একক বৃত্ত
fig, ax = plt.subplots(figsize=(6, 4), layout="constrained")
for c in (1, 2, 3):                        # d = c রেখা: অক্ষ c·√λ
    ell = mu[:, None] + V @ (c * np.sqrt(lam)[:, None] * circle)
    inside = 1 - np.exp(-c * c / 2)        # ২ মাত্রায় ভেতরের সম্ভাবনা
    ax.plot(ell[1], ell[0], lw=1.5, label=f"Mahalanobis {c}: {inside:.1%} inside")
for k in (0, 1):                           # অক্ষ: ±2√λ · eigenvector
    a = 2 * np.sqrt(lam[k]) * V[:, k]
    ax.plot([mu[1] - a[1], mu[1] + a[1]], [mu[0] - a[0], mu[0] + a[0]], "k--", lw=1)
ax.plot([8, 6, 3], [45, 35, 55], "rx")
ax.set_xlabel("distance (km)"); ax.set_ylabel("delivery time (min)")
ax.set_title("Contours of N(mu, Sigma)"); ax.legend(fontsize=8, loc="lower right")
plt.show()

x=[45.0, 8.0]: d²= 1.002  formula=0.004501  scipy=0.004501
x=[35.0, 6.0]: d²= 0.000  formula=0.007429  scipy=0.007429
x=[55.0, 3.0]: d²=15.294  formula=0.000004  scipy=0.000004
eigenvalues: [  4.388 104.612] | sqrt: [ 2.095 10.228]
long axis direction: [-0.977 -0.215]

৩ · Mahalanobis দূরত্ব আর χ² — কোন অর্ডার সত্যিই অস্বাভাবিক

বাজারঘরের অপারেশন দল প্রতিদিন কয়েকশো অর্ডারের মধ্যে সন্দেহজনকগুলো আলাদা করতে চায়: হয়তো রাইডার ভুল ঠিকানায় গিয়েছিলেন, হয়তো ডেটা এন্ট্রিতে ভুল। প্রচলিত সরল নিয়ম: প্রতিটি feature-এর z-score দেখুন, কোনোটি ৩-এর বেশি হলে চিহ্নিত করুন। এই নিয়ম প্রতিটি feature আলাদা করে দেখে, তাই “কাছে অথচ অনেক দেরি” ধরনের outlier ধরতে পারে না, কারণ সময় আর দূরত্ব দুটিই আলাদাভাবে সাধারণ পরিসরে।

Mahalanobis দূরত্বের সুবিধা: d² = ‖z‖² (অংশ ২), আর ডেটা Gaussian হলে d² ~ χ²_d, যার গড় d আর variance 2d. তাই একটি নির্দিষ্ট সম্ভাবনার সীমা টানা যায়, যেমন “সাধারণ অর্ডারের মাত্র ১% এর বাইরে যাবে”। কয়েকটি সীমা (scipy.stats.chi2.ppf থেকে):

লক্ষ করুন, d বাড়লে সীমাও বাড়ে: তিন মাত্রায় গড় থেকে ৩.৩৭ “Mahalanobis একক” দূরে থাকাও ৯৯% ক্ষেত্রে স্বাভাবিক। প্রতিটি feature-এ ৩ σ নিয়ম এই বৃদ্ধি মানে না। নিচের ব্লকে ৪০০টি কাল্পনিক সাধারণ অর্ডার (সময়, দূরত্ব, ওজন) আর চারটি বিশেষ অর্ডার।

import numpy as np
from scipy.stats import chi2

rng = np.random.default_rng(42)           # কাল্পনিক বাজারঘর: ৪০০টি সাধারণ অর্ডার
sd = np.array([10.0, 3.0, 1.5])           # সময় (মিনিট), দূরত্ব (km), ওজন (kg)
R = np.array([[1.0, 0.7, 0.3], [0.7, 1.0, 0.1], [0.3, 0.1, 1.0]])
X = rng.multivariate_normal([35, 6, 3], R * np.outer(sd, sd), size=400)
odd = np.array([[55, 2, 2.5],             # A: কাছে, অথচ অনেক দেরি
                [20, 12, 3],              # B: দূরে, অথচ খুব দ্রুত
                [68, 15, 4],              # C: দূরে আর ধীর (সামঞ্জস্যপূর্ণ)
                [40, 6, 8]])              # D: খুব ভারী পার্সেল
X = np.vstack([X, odd])

mu = X.mean(axis=0)
S = np.cov(X, rowvar=False)               # অনুমান করা Σ (n − 1 দিয়ে ভাগ)
D = X - mu
d2 = np.einsum("ij,ij->i", D @ np.linalg.inv(S), D)   # প্রতিটি সারির d²
z = np.abs(D / X.std(axis=0, ddof=1))     # feature ধরে আলাদা z-score

cut = chi2.ppf(0.99, df=3)                # χ²(3)-এর ৯৯% সীমা
print(f"chi2 99% cutoff for d² (df=3): {cut:.2f}")
print("order  d²     max|z|  Mahalanobis  z>3")
for name, i in zip("ABCD", range(400, 404)):
    print(f"  {name}  {d2[i]:6.2f}  {z[i].max():5.2f}   "
          f"{'flag' if d2[i] > cut else '-':^11}  {'flag' if z[i].max() > 3 else '-'}")

normal = slice(0, 400)
print("normal orders flagged:", (d2[normal] > cut).sum(), "by d²,",
      (z[normal].max(axis=1) > 3).sum(), "by z")
for q in (0.5, 0.9, 0.99):                # d²-এর quantile বনাম χ²(3)
    print(f"quantile {q}: data {np.quantile(d2[normal], q):5.2f}"
          f" | chi2 {chi2.ppf(q, 3):5.2f}")

chi2 99% cutoff for d² (df=3): 11.34
order  d²     max|z|  Mahalanobis  z>3
  A   22.12   1.92      flag      -
  B   22.22   1.88      flag      -
  C   10.63   3.16        -       flag
  D   12.05   3.42      flag      flag
normal orders flagged: 5 by d², 4 by z
quantile 0.5: data  2.12 | chi2  2.37
quantile 0.9: data  6.09 | chi2  6.25
quantile 0.99: data 12.02 | chi2 11.34

৪ · Marginal আর conditional — একটি অংশ জানলে বাকিটা কী

x-কে দুই ভাগে ভাগ করুন: x_a (যা জানতে চাই) আর x_b (যা জানা)। μ আর Σ-ও একইভাবে ভাগ হয়: μ = (μ_a, μ_b), আর Σ-এর চারটি ব্লক Σ_aa, Σ_ab, Σ_ba, Σ_bb. Gaussian-এর দুটি চমৎকার গুণ:

μ_a|b = μ_a + Σ_ab Σ_bb⁻¹ (x_b − μ_b)
Σ_a|b = Σ_aa − Σ_ab Σ_bb⁻¹ Σ_ba

দ্বিতীয় রাশিটির নাম Schur complement. যা বিয়োগ হচ্ছে, Σ_abΣ_bb⁻¹Σ_ba, তা PSD, তাই শর্ত দিলে variance কখনো বাড়ে না: তথ্য পেলে অনিশ্চয়তা কমে বা একই থাকে। দুই মাত্রায় (a = সময় t, b = দূরত্ব) সূত্র দুটি হয় পরিচিত রূপে: গড় = μ_t + ρ(σ_t/σ_d)(x_d − μ_d), আর variance = σ_t²(1 − ρ²)।

সূত্রটি কি সত্যিই ঠিক? নিচের ব্লক Cholesky দিয়ে (পর্ব ১০-এর x = μ + Lz) দুই লাখ জোড়া বানায়, তারপর দূরত্ব d₀-এর ±০.১ km-এর একটি সরু ফালি কেটে শুধু সেই অর্ডারগুলোর সময়ের গড় আর variance মাপে: “শর্ত দেওয়া”-র আক্ষরিক মানে।

import numpy as np
import matplotlib.pyplot as plt

mu = np.array([35.0, 6.0])                # (সময় মিনিট, দূরত্ব km)
Sigma = np.array([[100.0, 21.0], [21.0, 9.0]])
L = np.linalg.cholesky(Sigma)
rng = np.random.default_rng(42)
X = mu + rng.standard_normal((200_000, 2)) @ L.T      # x = μ + Lz

def cond_time(d0):                        # সূত্র: সময় | দূরত্ব = d0
    m = mu[0] + Sigma[0, 1] / Sigma[1, 1] * (d0 - mu[1])
    v = Sigma[0, 0] - Sigma[0, 1] ** 2 / Sigma[1, 1]
    return m, v

for d0 in (3.0, 6.0, 10.0):
    m, v = cond_time(d0)
    band = X[np.abs(X[:, 1] - d0) < 0.1, 0]           # সরু ফালি
    print(f"d0={d0:4.1f} formula: mean {m:5.2f} var {v:4.1f} | "
          f"slice n={band.size}: mean {band.mean():5.2f} var {band.var(ddof=1):4.1f}")

slope, icpt = np.polyfit(X[:, 1], X[:, 0], 1)         # least-squares রেখা
print(f"regression line: time = {icpt:.2f} + {slope:.3f} * distance")
print(f"formula line:    time = {mu[0] - 21 / 9 * mu[1]:.2f} + {21 / 9:.3f} * distance")

fig, ax = plt.subplots(figsize=(6, 4), layout="constrained")
ax.scatter(X[:3000, 1], X[:3000, 0], s=3, alpha=0.3, label="samples")
g = np.linspace(-2, 14, 50)
m, v = cond_time(g)
ax.plot(g, m, "r", label="E[time | distance]")
ax.plot(g, m + 2 * np.sqrt(v), "r--", lw=1, label="mean ± 2 sd")
ax.plot(g, m - 2 * np.sqrt(v), "r--", lw=1)
ax.axvspan(9.9, 10.1, color="orange", alpha=0.5, label="slice at 10 km")
ax.set_xlabel("distance (km)"); ax.set_ylabel("delivery time (min)")
ax.legend(fontsize=8); plt.show()

d0= 3.0 formula: mean 28.00 var 51.0 | slice n=3304: mean 28.26 var 49.4
d0= 6.0 formula: mean 35.00 var 51.0 | slice n=5090: mean 35.09 var 51.5
d0=10.0 formula: mean 44.33 var 51.0 | slice n=2212: mean 44.55 var 48.6
regression line: time = 21.08 + 2.318 * distance
formula line:    time = 21.00 + 2.333 * distance

Precision matrix-এর একটি গোপন অর্থ। Σ⁻¹-এর বাইরের কোনো ঘর শূন্য হলে সেই দুটি feature বাকি সবকিছু জানা থাকলে স্বাধীন (conditional independence). Σ-এর ঘর শূন্য হওয়া আর Σ⁻¹-এর ঘর শূন্য হওয়া আলাদা কথা: প্রথমটি “একা দেখলে সম্পর্কহীন”, দ্বিতীয়টি “বাকিদের প্রভাব সরালে সম্পর্কহীন”। Gaussian graphical model আর graphical lasso এই ধারণার ওপর দাঁড়ানো। অংশ ১-এর stepper-এ Σ⁻¹-এর বাইরের ঘর ঋণাত্মক ছিল, কারণ ২ মাত্রায় সেটি −ρ-এর সমানুপাতিক: correlation যত বেশি, যৌথ সরণ তত “ক্ষমা” পায়।

৫ · Conditional গড় মানেই regression রেখা — আর “গড়ের দিকে ফেরা”

অংশ ৪-এর শর্তসাপেক্ষ গড় E[সময় | দূরত্ব] = ২১ + ২.৩৩৩ × দূরত্ব, একটি সরলরেখা। এর ঢাল Cov(t, d)/Var(d), যা হুবহু least squares-এর ঢাল (AI-এর ভিত্তি কোর্সের পর্ব ৫৮)। তাই বলা যায়: ডেটা যৌথভাবে Gaussian হলে linear regression কোনো অনুমান নয়, সঠিক model; শর্তসাপেক্ষ গড় সত্যিই রৈখিক, আর ভুল (residual) Gaussian, যার variance σ_t²(1 − ρ²) সব x-এ এক। R² এখানে ρ² = ০.৪৯: দূরত্ব সময়ের variance-এর ৪৯% ব্যাখ্যা করে, বাকি ৫১ (মিনিট²) থেকে যায়। একাধিক feature হলে Σ_abΣ_bb⁻¹ হলো সব coefficient-এর vector, যা normal equations-এর (XᵀX)⁻¹Xᵀy-এর জনসংখ্যা-রূপ।

উল্টো দিকও খেয়াল করুন: বাস্তব ডেটা Gaussian না হলে E[y | x] বাঁকা হতে পারে, আর তখন সরলরেখা শুধু “সবচেয়ে ভালো রৈখিক অনুমান”, সত্যিকারের শর্তসাপেক্ষ গড় নয়। আর variance x-এর সঙ্গে বদলালে (দূরের অর্ডারে অনিশ্চয়তা বেশি) সব x-এ এক পরিসর দেওয়া ভুল হবে।

Regression to the mean: standard এককে অনুমান সবসময় ρ গুণ, তাই |ρ| < ১ হলে অনুমান গড়ের দিকে টানা থাকে। দূরত্ব ১.৩৩ sd বেশি হলে সময়ের অনুমান মাত্র ০.৯৩ sd বেশি। এটি কোনো কারণগত শক্তি নয়, শুধু অসম্পূর্ণ correlation-এর গাণিতিক ফল। একই কারণে খুব খারাপ মাসের পর বাজারঘরের যেকোনো “সংস্কার”-এর পরের মাস ভালো দেখাতে পারে, সংস্কার কাজ না করলেও। ভাবনার প্রশ্ন ২-এ এর ইতিহাস আর ফাঁদ।

ML-এ এই সূত্র অনেক জায়গায় ফিরে আসে। পর্ব ৪৫ (Bayesian linear regression)-এ posterior আর predictive বণ্টন, পর্ব ৪৬ (Gaussian process)-এ কিছু বিন্দুর মান জানা থাকলে নতুন বিন্দুর অনুমান, দুটিই আসলে এই conditional Gaussian সূত্র। Kalman filter (১৯৬০) প্রতিটি নতুন মাপে একই হিসাব বারবার চালায়, আর Gaussian mixture-এর EM (পর্ব ৩৯) প্রতিটি উপাদানের জন্য এই পর্বের density ব্যবহার করে।

৬ · ডেটা থেকে μ আর Σ — MLE, n বনাম n − ১, আর parameter-এর হিসাব

বাস্তবে μ আর Σ জানা থাকে না; n-টি সারি x₁, …, x_n থেকে অনুমান করতে হয়। Maximum likelihood (পর্ব ১৮, Maximum likelihood আর MAP-এ বিস্তারিত) মানে log-likelihood সবচেয়ে বড় করা:

ℓ(μ, Σ) = −(n/2)·log det Σ − ½ ∑_i (x_i − μ)ᵀΣ⁻¹(x_i − μ) − (nd/2)·log 2π

μ-এর সাপেক্ষে gradient শূন্য করলে Σ⁻¹∑(x_i − μ) = ০, তাই μ̂ = নমুনার গড়. Σ-এর সাপেক্ষে (পর্ব ১১-এর দুটি সূত্র লাগে: ∂ log det Σ/∂Σ = Σ⁻¹, আর trace-এর কৌশল) পাওয়া যায় Σ̂ = (1/n)∑(x_i − μ̂)(x_i − μ̂)ᵀ. এখানে n দিয়ে ভাগ, n − ১ নয়: MLE একটু biased, কারণ μ̂ নিজেই ডেটা থেকে এসেছে আর ডেটার দিকে টানা। n − ১ দিয়ে ভাগ করলে unbiased. NumPy-র np.cov ডিফল্টভাবে n − ১ নেয় (bias=True দিলে n), scikit-learn-এর EmpiricalCovariance আর LDA-র ভেতরের হিসাব n নেয়। n বড় হলে পার্থক্য তুচ্ছ, কিন্তু ছোট n-এ আর লাইব্রেরির ফল মেলাতে গেলে গুরুত্বপূর্ণ (অংশ ৭-এর ব্লকে দুটিই আছে)।

Parameter কত দ্রুত বাড়ে। μ-তে d-টি সংখ্যা, আর প্রতিসম Σ-এ d(d + ১)/২টি। d = ৩-এ মোট ৯, d = ৪০-এ ৪০ + ৮২০ = ৮৬০, d = ১,০০০-এ পাঁচ লাখের বেশি। তার ওপর Σ̂-এর rank সর্বোচ্চ n − ১, তাই n ≤ d হলে Σ̂ singular: inverse নেই, density নেই, Mahalanobis নেই। অংশ ৮ এই সমস্যার একটি নীতিগত সমাধান দেবে।

উল্টো দিক, নমুনা বানানো। পর্ব ১০-এ দেখেছি Σ = LLᵀ হলে x = μ + Lz ~ N(μ, Σ). NumPy-র rng.multivariate_normal(mean, cov) একই কাজ করে, তবে ডিফল্টভাবে SVD দিয়ে (method="svd"), যা singular Σ-ও সামলায়; method="cholesky" দ্রুত। আর SciPy-র multivariate_normal(mean, cov).rvs(size, random_state=...) একটি বণ্টন-object দেয়, যার pdf, logpdf আর cdf-ও আছে। লক্ষ করুন, একই seed হলেও এই তিনটি পথ আলাদা সংখ্যা দেয়, কারণ ভেতরের বর্গমূল-matrix আলাদা।

৭ · Gaussian দিয়ে classification — LDA আর QDA

ধরুন বাজারঘর জানতে চায় কোন অর্ডার ফেরত আসতে পারে। দুটি feature: ডেলিভারির দেরি (মিনিট) আর ছাড় (%). একটি পথ হলো logistic regression, যা সরাসরি P(ফেরত | x) শেখে, একটি discriminative model (AI-এর ভিত্তি কোর্সের পর্ব ৫৯)। অন্য পথ generative model: প্রতিটি শ্রেণির ডেটা কেমন দেখতে, তা শেখা। শ্রেণি k-এর জন্য ধরা হয় x | y = k ~ N(μ_k, Σ_k), আর শ্রেণির অনুপাত π_k (prior). তারপর Bayes-এর নিয়ম (AI-এর ভিত্তি কোর্সের পর্ব ৫০): P(y = k | x) ∝ π_k·N(x; μ_k, Σ_k). log নিলে প্রতিটি শ্রেণির discriminant function:

δ_k(x) = −½ log det Σ_k − ½ (x − μ_k)ᵀΣ_k⁻¹(x − μ_k) + log π_k

যে শ্রেণির δ বড়, সেটিই অনুমান। প্রতিটি শ্রেণির নিজের Σ_k থাকলে এর নাম QDA (Quadratic Discriminant Analysis), কারণ দুই শ্রেণির সীমারেখা δ₁ = δ₀ একটি দ্বিঘাত বক্ররেখা (উপবৃত্ত, পরাবৃত্ত বা অধিবৃত্ত)। সব শ্রেণির জন্য একটিই Σ (pooled covariance) ধরলে xᵀΣ⁻¹x পদটি দুই দিকে সমান হয়ে কেটে যায়, থাকে শুধু রৈখিক অংশ। সেটিই LDA (Linear Discriminant Analysis):

log [P(1 | x) / P(0 | x)] = wᵀ(x − m) + log(π₁/π₀), w = Σ⁻¹(μ₁ − μ₀), m = (μ₀ + μ₁)/2

লক্ষ করুন, এটি হুবহু logistic regression-এর আকার: P(1 | x) = σ(wᵀx + b). পার্থক্য শুধু w আর b কীভাবে আসে: LDA গড় আর covariance থেকে বন্ধ সূত্রে, logistic regression সরাসরি log loss কমিয়ে। আর সমান prior হলে LDA মানে “Mahalanobis দূরত্বে যে শ্রেণির গড় কাছে, সেটি”। Σ_k-কে কর্ণ-matrix ধরলে পাওয়া যায় Gaussian Naive Bayes (AI-এর ভিত্তি কোর্সের পর্ব ৫১)। নিচের stepper একটি অর্ডারের LDA সিদ্ধান্ত হাতে করে, প্যারামিটারগুলো পরের ব্লকের অনুমান থেকে গোল করে নেওয়া।

এবার পুরো ডেটায়: ৪৫০টি রাখা আর ১৫০টি ফেরত অর্ডার (কাল্পনিক), ফেরতগুলোর দেরিতে ছড়ানো বেশি, তাই দুই শ্রেণির Σ আলাদা। LDA আর QDA নিজে লিখে scikit-learn-এর LinearDiscriminantAnalysis আর QuadraticDiscriminantAnalysis-এর সঙ্গে মেলানো হলো, ৭০/৩০ ভাগে, test set আলাদা রেখে।

import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.discriminant_analysis import (LinearDiscriminantAnalysis,
                                           QuadraticDiscriminantAnalysis)
rng = np.random.default_rng(42)           # (দেরি মিনিট, ছাড় %); 0 = রাখা, 1 = ফেরত
X0 = rng.multivariate_normal([35, 10], [[100, 15], [15, 25]], size=450)
X1 = rng.multivariate_normal([50, 18], [[225, -20], [-20, 49]], size=150)
X = np.vstack([X0, X1]); y = np.repeat([0, 1], [450, 150])
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.3, stratify=y, random_state=42)

pri = np.array([np.mean(ytr == k) for k in (0, 1)])            # π_k
mus = np.array([Xtr[ytr == k].mean(axis=0) for k in (0, 1)])   # μ_k
covs = [np.cov(Xtr[ytr == k], rowvar=False) for k in (0, 1)]   # Σ_k (n_k − 1)
pooled = sum(p * np.cov(Xtr[ytr == k], rowvar=False, bias=True)
             for k, p in enumerate(pri))                       # এক Σ সবার
def lda_score(x, k):                      # রেখা: x-এ রৈখিক
    a = np.linalg.solve(pooled, mus[k])
    return x @ a - 0.5 * mus[k] @ a + np.log(pri[k])

def qda_score(x, k):                      # বক্ররেখা: x-এ দ্বিঘাত
    D = x - mus[k]
    m2 = np.einsum("ij,ij->i", D @ np.linalg.inv(covs[k]), D)
    return -0.5 * np.linalg.slogdet(covs[k])[1] - 0.5 * m2 + np.log(pri[k])
w = np.linalg.solve(pooled, mus[1] - mus[0])
b = lda_score(np.zeros(2), 1) - lda_score(np.zeros(2), 0)
lda = LinearDiscriminantAnalysis(solver="lsqr").fit(Xtr, ytr)
qda = QuadraticDiscriminantAnalysis().fit(Xtr, ytr)
print("w, b scratch:", w.round(4), round(b, 3))
print("w, b sklearn:", lda.coef_[0].round(4), lda.intercept_.round(3))
for name, f, model in (("LDA", lda_score, lda), ("QDA", qda_score, qda)):
    mine = (f(Xte, 1) > f(Xte, 0)).astype(int)
    print(f"{name}: agrees with sklearn on {np.mean(mine == model.predict(Xte)):.0%}"
          f" | test accuracy {np.mean(mine == yte):.3f}")
print("always 'kept' baseline accuracy:", round(np.mean(yte == 0), 3))

w, b scratch: [0.1269 0.314 ] -11.075
w, b sklearn: [0.1269 0.314 ] [-11.075]
LDA: agrees with sklearn on 100% | test accuracy 0.856
QDA: agrees with sklearn on 100% | test accuracy 0.878
always 'kept' baseline accuracy: 0.75

৮ · অল্প ডেটায় বড় Σ — Ledoit–Wolf shrinkage

অংশ ৬-এর সমস্যায় ফিরি। বাজারঘর ৪০টি feature (প্রতিটি পণ্য-বিভাগে খরচ ইত্যাদি) দিয়ে LDA চালাতে চায়, কিন্তু একটি নতুন জেলায় প্রতি শ্রেণিতে মাত্র ২৫টি অর্ডার। Σ-এ ৮২০টি অজানা, ডেটায় ৫০টি সারি। নমুনার covariance S এখানে শুধু অনিশ্চিত নয়, পদ্ধতিগতভাবে বিকৃত: তার বড় eigenvalue-গুলো আসলের চেয়ে বড়, ছোটগুলো আসলের চেয়ে ছোট (n ≤ d হলে একেবারে শূন্য)। আর Σ⁻¹ ঠিক এই ছোট eigenvalue-গুলোকে উল্টে বিশাল বানায়, তাই Mahalanobis দূরত্ব, density আর LDA-র w সবই ভুল দিকে ফুলে ওঠে।

Shrinkage-এর ধারণা: S-কে একটি সরল, স্থিতিশীল লক্ষ্যের দিকে একটু টানুন।

Σ̂ = (1 − α)·S + α·(tr S / d)·I, ০ ≤ α ≤ ১

লক্ষ্য (tr S/d)·I হলো “সব feature-এর variance সমান, কোনো correlation নেই”। α = ০ মানে শুধু ডেটা, α = ১ মানে শুধু লক্ষ্য। ফল: সব eigenvalue গড়ের দিকে সরে, ক্ষুদ্রতমটি আর শূন্য থাকে না, Σ̂ সবসময় invertible. এটি পর্ব ১৭-এর ridge-এর আত্মীয়: একটু bias মেনে অনেক variance কমানো। কিন্তু α কত? Olivier Ledoit আর Michael Wolf (২০০৪) একটি সূত্র দেন, যা ডেটা থেকেই α হিসাব করে, Frobenius দূরত্বে (সব ঘরের বর্গ-ভুলের যোগ) প্রত্যাশিত ভুল প্রায় সবচেয়ে কম রাখতে; কোনো cross-validation লাগে না। scikit-learn-এ sklearn.covariance.LedoitWolf, আর LDA-তে solver="lsqr" বা "eigen"-এর সঙ্গে shrinkage="auto".

নিচের পরীক্ষাটি একটু বড় (৪০ feature, ২০টি seed, কয়েক হাজার test সারি), তাই এটি ডেস্কটপে চালানো হয়েছে। সত্যিকারের run: Python 3.14.4, NumPy 2.5.4, scikit-learn 1.9.1 (লেখার সময়ের সর্বশেষ stable). আসল Σ জানা: Σ_ij = ০.৬^{|i − j|}. প্রতিটি n-এ ২০টি seed-এর গড় আপেক্ষিক ভুল ‖Σ̂ − Σ‖/‖Σ‖, গড় α, আর শেষ seed-এর condition number; তারপর দুই শ্রেণির LDA (প্রথম ৫টি feature-এ গড় ১ করে সরানো), প্রতি শ্রেণিতে ২৫টি training সারি, ১০টি seed.

~/bazarghar-ml$ python ledoit_wolf.py
numpy 2.5.4 | scikit-learn 1.9.1
 n   err(sample)  err(LW)  shrinkage  cond(sample)  cond(LW)
  20       0.975    0.591      0.609      1.17e+18      5.31
  50       0.629    0.479      0.423           875      7.39
 200       0.320    0.294      0.162          40.3      14.6
1000       0.140    0.137      0.037          19.3      16.3
LDA shrinkage=None    test accuracy mean 0.599 sd 0.033 over 10 seeds
LDA shrinkage='auto'  test accuracy mean 0.710 sd 0.020 over 10 seeds

পড়ুন সারি ধরে। n = ২০ (feature-এর অর্ধেক): নমুনার S-এর ভুল প্রায় ১০০% (০.৯৭৫), আর condition number ১০^১৮: কার্যত singular, কারণ rank সর্বোচ্চ ১৯. Ledoit–Wolf গড়ে α = ০.৬১ বেছে নেয়, ভুল নামে ০.৫৯-এ, condition number ৫.৩. n = ৫০-এ α = ০.৪২, n = ২০০-এ ০.১৬, n = ১,০০০-এ ০.০৪: ডেটা যত বাড়ে, সূত্র তত কম টানে, আর দুই অনুমান প্রায় এক হয়ে যায়। এটিই ভালো shrinkage-এর আচরণ। LDA: ৫০টি training সারিতে ৪০ মাত্রার Σ, shrinkage ছাড়া test accuracy গড়ে ০.৫৯৯, shrinkage-সহ ০.৭১০, আর seed-ভেদে ছড়ানোও কম (sd ০.০৩৩ বনাম ০.০২০)। এই কাল্পনিক ডেটায় আসল Σ জানা থাকলে সর্বোচ্চ সম্ভব accuracy প্রায় ০.৭৯ (হিসাব করা Bayes accuracy), তাই shrinkage ফাঁকের অনেকটা পূরণ করেছে। এগুলো একটি local run-এর সংখ্যা, সাধারণ ফল নয়।

Ledoit–Wolf একমাত্র পথ নয়: OAS (Oracle Approximating Shrinkage) Gaussian ডেটায় প্রায়ই একটু ভালো, GraphicalLasso Σ⁻¹-এ শূন্য খোঁজে (অংশ ৪-এর precision matrix), আর লেখার সময়ের scikit-learn-এ (1.8 থেকে) QuadraticDiscriminantAnalysis-এও solver আর shrinkage parameter আছে; পুরোনো version-এ শুধু reg_param ছিল। আপনার version-এর docs দেখে নিন।

৯ · AI-এর যুগে — Gaussian-এর সূত্র আর কোড যাচাই

AI assistant এই বিষয়ের সূত্র প্রায়ই নির্ভুল ছাপে, কিন্তু ভুল হয় মাত্রায়, ভাগে আর অনুমানে। নিচের ভুলগুলো বারবার দেখা যায়:

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

• Density বা d² নিজের হাতের ২ × ২ উদাহরণ আর SciPy, দুটির সঙ্গেই মিলিয়েছি
• Outlier-এর সীমা feature সংখ্যা অনুযায়ী χ² থেকে নিয়েছি, আর সাধারণ ডেটায় কত শতাংশ চিহ্নিত হয় গুনে দেখেছি
• Gaussian অনুমান যাচাই করেছি: histogram, বাঁকা feature-এ log, আর d²-এর quantile বনাম χ²
• ডেটা আগে ভাগ করেছি, test set তালা দেওয়া, আর LDA/QDA-কে একটি সরল baseline আর logistic regression-এর সঙ্গে কয়েকটি seed-এ তুলনা করেছি
• Feature অনেক আর সারি কম হলে shrinkage চেষ্টা করেছি, আর α-এর মান লিখে রেখেছি
• কোনো আসল গ্রাহকের ব্যক্তিগত ডেটা AI-তে পেস্ট করিনি