Density estimation — histogram, KDE আর mixture model

AI-এর ভিত্তি কোর্সের পর্ব ৩০-এ histogram আর seaborn-এর kdeplot ছবি আঁকার যন্ত্র হিসেবে এসেছে, আর পর্ব ৩৬-এ PDF মানে কী, তা দেখেছেন। এই কোর্সের পর্ব ১৮ একটি জানা আকারের বণ্টনের parameter MLE দিয়ে বের করেছে, আর পর্ব ১৯ অনেক মাত্রার Gaussian density দেখিয়েছে। এই পর্বে প্রশ্নটি উল্টো: বণ্টনের আকার আগে থেকে জানা না থাকলে, শুধু নমুনা থেকে পুরো density কীভাবে অনুমান করবেন, আর সেই অনুমান কতটা ভালো, তা কীভাবে মাপবেন? চালানো যায় এমন চারটি ব্লক পাতাতেই চলে (ব্রাউজারের Python: NumPy, SciPy, matplotlib আর scikit-learn); প্রথমবার প্যাকেজ নামতে কয়েক সেকেন্ড লাগে, scikit-learn-এ প্রায় ১০ সেকেন্ড।

১ · Density estimation কী, আর কোথায় কাজে লাগে

আনুষ্ঠানিকভাবে: একটি অজানা বণ্টন p থেকে স্বাধীনভাবে আসা নমুনা x₁, …, x_n দেওয়া আছে। লক্ষ্য একটি function p̂, যা সব জায়গায় অ-ঋণাত্মক, মোট integral ১, আর p-এর যতটা সম্ভব কাছে। একে বলে density estimation. এটি unsupervised learning-এর সবচেয়ে মৌলিক কাজগুলোর একটি: কোনো label নেই, শুধু “ডেটা কোথায় ঘন, কোথায় পাতলা” শেখা। একটি ভালো p̂ হাতে থাকলে অনেক কাজ হয়:

দুটি বড় পথ আছে। Parametric: আকার আগে ধরে নিন (Gaussian, gamma, log-normal), তারপর কয়েকটি parameter MLE দিয়ে বের করুন (পর্ব ১৮)। আকার ঠিক হলে অল্প ডেটাতেই নির্ভুল; ভুল হলে যত ডেটাই দিন, ভুল থেকে যায়। Nonparametric: আকার ধরে না নিয়ে ডেটাকেই কথা বলতে দিন: histogram আর KDE. এদের “parameter” আসলে পুরো ডেটাসেট, তাই এরা nonparametric model. মাঝামাঝি পথ mixture model: কয়েকটি সরল আকারের ওজন-দেওয়া যোগফল, যা পর্ব ৩৯ (Gaussian mixture আর EM algorithm — নরম clustering)-এর বিষয়।

এই পর্বের ডেটা। বাজারঘরের ৪০০টি কাল্পনিক ডেলিভারি সময় (ঘণ্টায়), seed ৪২ দিয়ে বানানো। তিন ধরনের অর্ডার: ২৫% ঢাকার express (exponential, গড় ৪ ঘণ্টা), ৪৫% ঢাকার সাধারণ (gamma, গড় ২৪ ঘণ্টা), ৩০% ঢাকার বাইরে (gamma, গড় ৬৩ ঘণ্টা)। ডেটা আমরাই বানিয়েছি, তাই আসল density জানা; বাস্তবে এই সুবিধা থাকে না, কিন্তু শেখার সময় অনুমানকে সত্যের সঙ্গে মেলানো যায়।

২ · Histogram — সবচেয়ে পুরোনো অনুমান আর bin-এর সমস্যা

সংখ্যারেখাকে w চওড়া খোপে (bin) ভাগ করুন। x যে খোপে পড়ে, সেই খোপে n_kটি নমুনা থাকলে histogram density হলো:

p̂(x) = n_k / (n · w)

n দিয়ে ভাগ করলে অনুপাত, w দিয়ে ভাগ করলে “প্রতি ঘণ্টায়” ঘনত্ব, আর সব খোপের ক্ষেত্রফলের যোগ ১. matplotlib-এর hist(..., density=True) আর NumPy-র np.histogram(..., density=True) ঠিক এটিই করে। দুটি সিদ্ধান্ত লাগে: খোপের চওড়া w আর প্রথম খোপ কোথা থেকে শুরু (origin). চওড়া w বড় হলে প্রতিটি খোপে অনেক নমুনা, গণনা স্থির, কিন্তু খোপের ভেতরের আকার মুছে যায়: bias বেশি। w ছোট হলে আকার ধরা পড়ে, কিন্তু প্রতিটি খোপে দু-চারটি নমুনা, তাই উচ্চতা লাফায়: variance বেশি। এটি bias–variance tradeoff-এর সবচেয়ে চোখে-দেখা রূপ।

হিসাব করলে দেখা যায় bias² প্রায় w²-এর সমানুপাতিক আর variance প্রায় 1/(nw)-এর। দুটির যোগ ন্যূনতম হয় যখন w ∝ n^−1/3, আর তখন গড় বর্গ-ভুল (MISE) কমে n^−2/3 হারে। সমানুপাতিক ধ্রুবকটি আসল density-র ওপর নির্ভর করে, যা অজানা। তাই তিনটি প্রচলিত নিয়ম:

এখানে σ̂ নমুনার standard deviation, আর IQR = Q3 − Q1, মাঝের ৫০% মানের বিস্তার। একটি outlier σ̂-কে অনেক টানে, IQR-কে প্রায় নড়ায় না। নিচের ব্লক দুটি নিয়ম হাতে হিসাব করে, NumPy-র np.histogram_bin_edges-এর সঙ্গে মেলায়, আর চার রকম bin-এ histogram আঁকে।

import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(42)            # বাজারঘরের ৪০০টি ডেলিভারি (ঘণ্টা)
n = 400
g = rng.choice(3, size=n, p=[0.25, 0.45, 0.30])   # express, ঢাকা, ঢাকার বাইরে
x = np.where(g == 0, rng.exponential(4, n),
     np.where(g == 1, rng.gamma(4, 6, n), rng.gamma(9, 7, n)))
q1, q3 = np.percentile(x, [25, 75])
print(f"n={n}  min={x.min():.2f}  median={np.median(x):.1f}  max={x.max():.1f}")
print(f"sd={x.std(ddof=1):.2f}  IQR={q3 - q1:.2f}")
fd = 2 * (q3 - q1) * n ** (-1 / 3)          # Freedman–Diaconis
scott = 3.49 * x.std(ddof=1) * n ** (-1 / 3)
print(f"FD width by hand={fd:.2f}  Scott width by hand={scott:.2f}")
for rule in ("fd", "scott", "sturges"):
    e = np.histogram_bin_edges(x, bins=rule)
    print(f"numpy {rule:8s}: {len(e) - 1:3d} bins, width {e[1] - e[0]:.2f}")

fig, axes = plt.subplots(2, 2, figsize=(7, 5), sharex=True)
for ax, b, name in zip(axes.flat, [4, "fd", "scott", 120],
                       ["4 bins (too wide)", "Freedman-Diaconis",
                        "Scott", "120 bins (too narrow)"]):
    ax.hist(x, bins=b, density=True, color="#93c5fd", edgecolor="#1e3a8a")
    ax.set_title(name, fontsize=10)
for ax in axes[1]:
    ax.set_xlabel("delivery time (hours)")
for ax in axes[:, 0]:
    ax.set_ylabel("density")
fig.subplots_adjust(hspace=0.35, wspace=0.25)
plt.show()

n=400  min=0.04  median=22.9  max=129.2
sd=26.70  IQR=36.26
FD width by hand=9.84  Scott width by hand=12.65
numpy fd      :  14 bins, width 9.23
numpy scott   :  11 bins, width 11.74
numpy sturges :  10 bins, width 12.92

Histogram-এর আরও দুটি দুর্বলতা আছে। প্রথমত, origin অর্ধেক bin সরালেই ছবি বদলে যেতে পারে, অথচ ডেটা একই। দ্বিতীয়ত, p̂ ধাপের মতো: খোপের সীমানায় হঠাৎ লাফ, ভেতরে সমতল। অনেকগুলো সরানো origin-এর histogram-এর গড় নিলে (averaged shifted histogram) দুটি সমস্যাই কমে, আর সরানোর সংখ্যা অসীমের দিকে নিলে যা পাওয়া যায়, তা-ই পরের অংশের KDE (একটি ত্রিভুজ-আকারের kernel-সহ)।

৩ · KDE — প্রতিটি বিন্দুতে একটি ছোট ঢিবি

Kernel density estimation (KDE)-এর সূত্র:

p̂_h(x) = (1 / (n h)) · ∑_i=1ⁿ K((x − x_i) / h)

এখানে K একটি kernel: শূন্যের চারপাশে প্রতিসম, অ-ঋণাত্মক, মোট integral ১. সবচেয়ে প্রচলিত Gaussian kernel φ(u) = e^−u²/2/√(2π). h হলো bandwidth: প্রতিটি ঢিবির চওড়া। K((x − x_i)/h)/h হলো x_i-কে কেন্দ্র করে h-মাপের একটি ঢিবি, যার ক্ষেত্রফল ১; n দিয়ে ভাগ করলে প্রতিটির ক্ষেত্রফল ১/n, আর যোগফলের ক্ষেত্রফল ঠিক ১. তাই p̂ নিজেই একটি বৈধ density. এটি পর্ব ২৯-এর RBF kernel-এর আত্মীয়: Gaussian KDE-র মান আসলে x আর সব ডেটা বিন্দুর RBF similarity-র গড়, একটি ধ্রুবক দিয়ে মাপা। তবে সেখানে kernel ছিল দুটি বিন্দুর মিল, এখানে মসৃণ করার ছাঁকনি।

Kernel-এর আকার কতটা জরুরি? Gaussian ছাড়াও আছে Epanechnikov K(u) = ¾(1 − u²) যখন |u| ≤ 1, tophat (সমতল বাক্স) ইত্যাদি। তাত্ত্বিকভাবে গড় বর্গ-ভুলের মাপে Epanechnikov সেরা, কিন্তু Gaussian-এর দক্ষতা তার প্রায় ৯৫%, আর বাকিগুলোও কাছাকাছি। তাই বাস্তবে kernel-এর চেয়ে h বাছাই অনেক বেশি গুরুত্বপূর্ণ। দুটি লাইব্রেরির আচরণ মনে রাখুন। SciPy-র gaussian_kde(data, bw_method=c)-এ c একটি গুণক: আসল চওড়া = c × ডেটার sd. scikit-learn-এর KernelDensity(bandwidth=h)-এ h সরাসরি চওড়া, আর score_samples density নয়, log-density ফেরত দেয়; np.exp লাগে। ভেতরে এটি KD-tree বা ball tree দিয়ে দূরের বিন্দু বাদ দিয়ে হিসাব দ্রুত করে।

৪ · Bandwidth বাছাই — Silverman-এর নিয়ম বনাম held-out likelihood

Histogram-এর মতোই bias আর variance টানাটানি করে, শুধু হার ভালো। ছোট h-এ bias প্রায় (h²/2)·p″(x) (বাঁক যত তীক্ষ্ণ, ভুল তত বেশি), আর variance প্রায় p(x)·R(K)/(nh), যেখানে R(K) = ∫K². যোগফল ন্যূনতম হয় h ∝ n^−1/5-এ, আর তখন MISE কমে n^−4/5 হারে, histogram-এর n^−2/3-এর চেয়ে দ্রুত। কিন্তু ধ্রুবকে আছে অজানা p″. দুটি সমাধান প্রচলিত।

(১) Normal-reference নিয়ম। ধরে নিন আসল density normal, তাহলে সূত্রটি বন্ধ রূপে আসে: h = ১.০৬ σ n^−1/5. Bernard Silverman তাঁর ১৯৮৬ সালের বইয়ে একটি একটু সাবধানী রূপ সুপারিশ করেন, যা বাঁকা বা দুই-চূড়ার ডেটায় কম মসৃণ করে:

h = 0.9 · min(σ̂, IQR / 1.34) · n^−1/5

১.৩৪ কেন? Normal বণ্টনে IQR ≈ ১.৩৪৯σ, তাই IQR/১.৩৪ হলো σ-এর একটি robust অনুমান; দুটির ছোটটি নিলে একটি outlier σ̂ ফুলিয়ে h-কে বড় করতে পারে না।

(২) Cross-validated likelihood. আকার ধরে না নিয়ে সরাসরি জিজ্ঞেস করুন: কোন h দিলে নতুন ডেটা সবচেয়ে সম্ভাব্য দেখায়? ডেটাকে কয়েক ভাগে (fold) ভাগ করে, প্রতিবার বাকি ভাগে KDE বানিয়ে বাদ রাখা ভাগের log-likelihood ∑ log p̂(x_j) যোগ করুন। এটিই held-out log-likelihood. Training ডেটার নিজের likelihood কেন চলবে না, তা নিচের stepper দেখায়।

scikit-learn-এ কাজটি সহজ: KernelDensity.score(X) ফেরত দেয় X-এর মোট log-likelihood, আর GridSearchCV estimator-এর score-কেই ডিফল্ট মাপ ধরে। তাই GridSearchCV(KernelDensity(), {"bandwidth": ...}) কোনো label ছাড়াই cross-validation দিয়ে held-out log-likelihood সর্বোচ্চ করে। (AI-এর ভিত্তি কোর্সের পর্ব ৬৬-এর grid search, শুধু মাপটি এখানে likelihood.) নিচের ব্লক নিজের KDE তিনটি বিন্দুতে SciPy আর scikit-learn-এর সঙ্গে মেলায়, Silverman আর cross-validation-এর h তুলনা করে, আর চারটি h-এ বক্ররেখা আঁকে।

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import gaussian_kde, iqr
from sklearn.neighbors import KernelDensity
from sklearn.model_selection import GridSearchCV

rng, n = np.random.default_rng(42), 400    # আগের ব্লকের একই ৪০০টি ডেলিভারি
g = rng.choice(3, size=n, p=[0.25, 0.45, 0.30])
x = np.where(g == 0, rng.exponential(4, n),
     np.where(g == 1, rng.gamma(4, 6, n), rng.gamma(9, 7, n)))
def kde(t, data, h):                       # নিজের KDE: প্রতিটি বিন্দুতে একটি Gaussian ঢিবি
    u = (t[:, None] - data[None, :]) / h
    return np.exp(-0.5 * u**2).sum(axis=1) / (len(data) * h * np.sqrt(2 * np.pi))

s = x.std(ddof=1)
h_silv = 0.9 * min(s, iqr(x) / 1.34) * n ** (-1 / 5)   # Silverman-এর নিয়ম
t = np.array([2.0, 24.0, 60.0])
sp = gaussian_kde(x, bw_method=h_silv / s)(t)          # scipy নেয় h/sd
sk = np.exp(KernelDensity(bandwidth=h_silv).fit(x[:, None]).score_samples(t[:, None]))
print(f"Silverman h = {h_silv:.2f} hours")
for row in zip(t, kde(t, x, h_silv), sp, sk):
    print("t={:4.0f}  mine={:.5f}  scipy={:.5f}  sklearn={:.5f}".format(*row))
grid = GridSearchCV(KernelDensity(), {"bandwidth": np.linspace(1, 15, 29)}, cv=5)
grid.fit(x[:, None])                       # score = held-out fold-এর মোট log-likelihood
h_cv = grid.best_params_["bandwidth"]
print(f"GridSearchCV h = {h_cv:.1f}")
sk_rule = KernelDensity(bandwidth="silverman").fit(x[:, None]).bandwidth_
print(f"sklearn bandwidth='silverman' gives h = {sk_rule:.3f}  (not times sd!)")
tt = np.linspace(-10, 140, 400)
fig, ax = plt.subplots(figsize=(7, 3.6), layout="constrained")
ax.hist(x, bins="fd", density=True, color="#e5e7eb")
for h, lab in [(1, "too small"), (h_cv, "CV"), (h_silv, "Silverman"), (20, "too large")]:
    ax.plot(tt, kde(tt, x, h), label=f"h = {h:.1f} ({lab})")
ax.set_xlabel("delivery time (hours)"); ax.set_ylabel("density")
ax.set_title("Gaussian KDE with four bandwidths"); ax.legend(); plt.show()

Silverman h = 7.25 hours
t=   2  mine=0.01575  scipy=0.01575  sklearn=0.01575
t=  24  mine=0.01609  scipy=0.01609  sklearn=0.01609
t=  60  mine=0.00539  scipy=0.00539  sklearn=0.00539
GridSearchCV h = 2.5
sklearn bandwidth='silverman' gives h = 0.320  (not times sd!)

৫ · ধনাত্মক ডেটার সীমানা — boundary bias, log রূপান্তর আর নতুন নমুনা

ডেলিভারি সময় কখনো ঋণাত্মক নয়। কিন্তু Gaussian ঢিবি দুই দিকে ছড়ায়: ০.৩ ঘণ্টার একটি express অর্ডারের ঢিবির (h = ২.৫) প্রায় অর্ধেক পড়ে শূন্যের নিচে, যেখানে কোনো ডেলিভারি সম্ভব নয়। ফলে দুটি ক্ষতি: শূন্যের নিচে অর্থহীন সম্ভাবনা, আর শূন্যের ঠিক ওপরে ঘনত্ব কম। আসল density সীমানায় মসৃণভাবে ধনাত্মক হলে, সীমানার ঠিক ওপরে সাধারণ KDE প্রায় অর্ধেক মান দেয়, কারণ সেখানে ঢিবি শুধু এক দিক থেকে আসে। এটিই boundary bias.

সমাধান: log জগতে কাজ। y = log x নিন। y-এর কোনো সীমানা নেই, তাই y-তে সাধারণ KDE বানান, তারপর ফিরে আসুন। কিন্তু শুধু p̂_Y(log x) বসালে ভুল হবে: density একটি “প্রতি একক” মাপ, আর log একককে টানে-চাপে। চলক বদলের নিয়ম (পর্ব ১১-এর Jacobian-এর এক মাত্রার রূপ):

p̂_X(x) = p̂_Y(log x) · |dy/dx| = p̂_Y(log x) / x

১/x গুণকটি না দিলে বক্ররেখার ক্ষেত্রফল ১ থাকে না। এর আরও কয়েকটি বিকল্প আছে: reflection (প্রতিটি বিন্দুর −x_i-এও ঢিবি বসিয়ে শূন্যের নিচের অংশ ভাঁজ করে ওপরে আনা), বিশেষ boundary kernel, বা gamma-আকারের kernel. Log পথটি সবচেয়ে সরল, আর ডানে লম্বা লেজওয়ালা ডেটার জন্য বাড়তি সুবিধা দেয়: log জগতে লেজ ছোট হয়, তাই একটি h সব জায়গায় ভালো খাটে।

নতুন নমুনা। KDE থেকে নমুনা নেওয়া খুব সহজ: এলোমেলোভাবে একটি ডেটা বিন্দু x_I বাছুন, তারপর তার সঙ্গে kernel-এর একটি noise যোগ করুন: x* = x_I + h·ε, যেখানে Gaussian kernel-এ ε ~ N(০, ১). এটি আসলে পর্ব ২১-এর bootstrap, একটু “মসৃণ” করে (smoothed bootstrap). scikit-learn-এ KernelDensity.sample() এটি করে, তবে শুধু Gaussian আর tophat kernel-এ। Log জগতে নমুনা নিয়ে exp করলে সব নমুনা নিশ্চিতভাবে ধনাত্মক।

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import expon, gamma, norm, iqr
from sklearn.neighbors import KernelDensity
rng, n = np.random.default_rng(42), 400    # একই ৪০০টি ডেলিভারি
g = rng.choice(3, size=n, p=[0.25, 0.45, 0.30])
x = np.where(g == 0, rng.exponential(4, n),
     np.where(g == 1, rng.gamma(4, 6, n), rng.gamma(9, 7, n)))
def true_pdf(t):                           # ডেটা আমরাই বানিয়েছি, তাই আসল density জানা
    return (0.25 * expon.pdf(t, scale=4) + 0.45 * gamma.pdf(t, 4, scale=6)
            + 0.30 * gamma.pdf(t, 9, scale=7))
h = 2.5                                    # আগের ব্লকের cross-validated h
plain = KernelDensity(bandwidth=h).fit(x[:, None])
print(f"plain KDE mass below 0 hours: {norm.cdf(-x / h).mean():.3f}")
y = np.log(x)                              # log জগতে KDE, তারপর ফেরত
hy = 0.9 * min(y.std(ddof=1), iqr(y) / 1.34) * n ** (-1 / 5)
logk = KernelDensity(bandwidth=hy).fit(y[:, None])
def log_kde_pdf(t):                        # f_X(t) = f_Y(log t) / t   (Jacobian)
    return np.exp(logk.score_samples(np.log(t)[:, None])) / t
t = np.array([0.5, 2.0, 24.0])
print(f"Silverman h on log scale = {hy:.3f}")
for row in zip(t, true_pdf(t), np.exp(plain.score_samples(t[:, None])), log_kde_pdf(t)):
    print("t={:4.1f}  true={:.4f}  plain={:.4f}  log-KDE={:.4f}".format(*row))
grid = np.geomspace(1e-4, 2000, 5000)      # log-ফাঁকে বিন্দু
print(f"log-KDE integrates to {np.trapezoid(log_kde_pdf(grid), grid):.3f}")
s1 = plain.sample(1000, random_state=42).ravel()          # নতুন নমুনা
s2 = np.exp(logk.sample(1000, random_state=42).ravel())
print(f"new samples < 0: plain {(s1 < 0).sum()}, log-KDE {(s2 < 0).sum()}")
tt = np.linspace(0.05, 40, 400)
fig, ax = plt.subplots(figsize=(7, 3.6), layout="constrained")
ax.plot(tt, true_pdf(tt), "k--", label="true density")
ax.plot(tt, np.exp(plain.score_samples(tt[:, None])), label=f"plain KDE, h={h}")
ax.plot(tt, log_kde_pdf(tt), label="KDE on log scale, back-transformed")
ax.set_xlabel("delivery time (hours)"); ax.set_ylabel("density")
ax.set_title("Boundary bias near zero"); ax.legend(); plt.show()

plain KDE mass below 0 hours: 0.046
Silverman h on log scale = 0.325
t= 0.5  true=0.0552  plain=0.0228  log-KDE=0.0526
t= 2.0  true=0.0382  plain=0.0264  log-KDE=0.0372
t=24.0  true=0.0155  plain=0.0189  log-KDE=0.0152
log-KDE integrates to 1.000
new samples < 0: plain 38, log-KDE 0

৬ · অনেক মাত্রায় — multivariate KDE, mixture model আর curse of dimensionality

KDE সহজেই d মাত্রায় যায়: প্রতিটি বিন্দুতে একটি d মাত্রার ঢিবি, p̂(x) = (1/(n h^d))∑K(‖x − x_i‖/h). scikit-learn সব মাত্রায় একই h নেয়, তাই আগে feature-গুলো একই মাপে আনুন (standardize বা log); নইলে km আর মিনিটের মতো আলাদা এককে একটি h অর্থহীন। SciPy-র gaussian_kde ডেটার covariance matrix দিয়ে ঢিবির আকার নেয়, তাই হেলানো মেঘের সঙ্গে নিজে থেকেই খাপ খায়।

Parametric বিকল্প: mixture. একটি Gaussian mixture model density-কে K-টি Gaussian-এর ওজন-দেওয়া যোগ ধরে: p(x) = ∑_k π_k N(x | μ_k, Σ_k), ∑π_k = ১. KDE-কে চরম mixture ভাবা যায়: n-টি উপাদান, প্রতিটির ওজন ১/n, গড় একেকটি ডেটা বিন্দু, covariance h²I. GMM সেখানে কয়েকটি উপাদানে পুরো মেঘ ধরে: parameter-এর সংখ্যা ডেটার সঙ্গে বাড়ে না, মূল্যায়ন দ্রুত (n-টি ঢিবির বদলে K-টি), আর আকার ঠিক হলে অল্প ডেটায় ভালো। EM দিয়ে fit করা আর BIC দিয়ে K বাছাই পর্ব ৩৯-এর বিষয়; এখানে আমরা GMM-কে শুধু density estimate হিসেবে দেখব।

Curse of dimensionality. KDE-র সেরা MISE হার d মাত্রায় n^−4/(4+d). d = ১-এ n^−0.8, d = ১০-এ n^−0.29: একই নির্ভুলতা পেতে ডেটা লাগে বহু গুণ বেশি। কারণটি সহজ: একটি বিন্দুর “কাছে” যথেষ্ট প্রতিবেশী না থাকলে স্থানীয় গড় সম্ভব নয়, আর মাত্রা বাড়লে প্রতিবেশী উধাও হয় (AI-এর ভিত্তি কোর্সের পর্ব ৩৩-এ এর প্রথম দেখা)। নিচের ব্লক প্রথমে বাজারঘরের দুই মাত্রার ডেটায় (log সময়, log দূরত্ব) KDE আর বিভিন্ন K-এর GMM-কে আলাদা test ডেটার গড় log-likelihood দিয়ে তুলনা করে, তারপর দেখায় মাত্রা বাড়লে ব্যাসার্ধ ১-এর ভেতরে কত অংশ বিন্দু থাকে।

import numpy as np
from scipy.stats import chi2
from sklearn.neighbors import KernelDensity
from sklearn.mixture import GaussianMixture
from sklearn.model_selection import GridSearchCV

def orders(n, seed):                       # (log সময়, log দূরত্ব), তিন ধরনের ডেলিভারি
    rng = np.random.default_rng(seed)
    g = rng.choice(3, size=n, p=[0.25, 0.45, 0.30])
    mt, md = np.array([1.0, 3.1, 4.1])[g], np.array([1.2, 2.1, 5.0])[g]
    t = mt + rng.normal(0, 0.45, n)
    return np.column_stack([t, md + 0.6 * (t - mt) + rng.normal(0, 0.35, n)])
train, test = orders(600, 1), orders(2000, 2)

kde = GridSearchCV(KernelDensity(), {"bandwidth": np.linspace(0.05, 0.6, 23)}, cv=5)
kde.fit(train)
h = kde.best_params_["bandwidth"]
print(f"2-D KDE  h={h:.3f}  held-out mean log-lik = {kde.score(test) / len(test):.3f}")
for k in (1, 2, 3, 4, 6):
    gm = GaussianMixture(k, init_params="k-means++", n_init=3, random_state=0).fit(train)
    print(f"GMM k={k}  held-out mean log-lik = {gm.score(test):.3f}  BIC = {gm.bic(train):7.1f}")

rng = np.random.default_rng(0)             # মাত্রা বাড়লে "কাছের" বিন্দু কত?
print(" d  share within radius 1 (sim)  theory")
for d in (1, 2, 5, 10, 20):
    z = rng.standard_normal((20000, d))
    inside = (np.linalg.norm(z, axis=1) <= 1).mean()
    print(f"{d:2d}  {inside:12.5f}  {chi2.cdf(1, d):12.2e}")

2-D KDE  h=0.175  held-out mean log-lik = -2.102
GMM k=1  held-out mean log-lik = -2.981  BIC =  3594.0
GMM k=2  held-out mean log-lik = -2.236  BIC =  2821.6
GMM k=3  held-out mean log-lik = -2.056  BIC =  2626.0
GMM k=4  held-out mean log-lik = -2.059  BIC =  2660.9
GMM k=6  held-out mean log-lik = -2.064  BIC =  2732.6
 d  share within radius 1 (sim)  theory
 1       0.68425      6.83e-01
 2       0.39570      3.93e-01
 5       0.03730      3.74e-02
10       0.00015      1.72e-04
20       0.00000      1.71e-10

বেশি মাত্রায় তাই সাধারণ কৌশল: কম মাত্রায় নামানো (PCA, পর্ব ০৮), feature-গুলো স্বাধীন ধরে এক মাত্রার density-র গুণফল (Naive Bayes, AI-এর ভিত্তি কোর্সের পর্ব ৫১), বা mixture-এর মতো কাঠামোবদ্ধ model. ছবি, শব্দ বা লেখার মতো হাজার মাত্রার ডেটায় density শেখে normalizing flow, VAE বা diffusion-এর মতো deep generative model, যা পরে deep learning আর generative AI কোর্সে আসবে।

৭ · অনুমান কতটা ভালো — মূল্যায়ন আর পদ্ধতি বাছাই

Classification-এ accuracy আছে; density-তে মূল মাপ held-out গড় log-likelihood: (১/m)∑ log p̂(x_j), যেখানে x_j শেখার সময় দেখা হয়নি। এটি পর্ব ১৮-এর MLE-র নীতিকে নতুন ডেটায় খাটায়, আর AI-এর ভিত্তি কোর্সের পর্ব ৩৮-এর ভাষায় এটি আসল বণ্টন আর p̂-এর cross-entropy-র ঋণাত্মক অনুমান: যত বড়, KL divergence তত কম। তিনটি নিয়ম মানুন:

• একই জগতে তুলনা। Log সময়ের ওপর বানানো density-র log-likelihood আর ঘণ্টার ওপর বানানোটির log-likelihood সরাসরি তুলনীয় নয়। ঘণ্টার জগতে আনতে প্রতিটি বিন্দুতে −log x যোগ করুন (Jacobian).
• একটি শূন্য সব নষ্ট করে। কোনো test বিন্দুতে p̂ = ০ হলে log = −∞. Tophat বা Epanechnikov kernel-এর সীমিত ঢিবি, বা ফাঁকা bin-ওয়ালা histogram-এ এটি সহজেই হয়; Gaussian kernel-এ হয় না।
• সংখ্যার সঙ্গে চোখ। Log-likelihood-এর ছোট পার্থক্য seed বদলালে উল্টে যেতে পারে; কয়েকটি seed বা fold-এ গড় আর ছড়ানো দেখুন (পর্ব ২১)। আর ছবি আঁকুন: শূন্যের নিচে ভর, লেজে ঢেউ, অর্থহীন চূড়া সংখ্যায় ধরা পড়ে না।

৮ · AI-এর যুগে — density estimation-এর কোড যাচাই

AI assistant KDE-র কোড দ্রুত লেখে, কিন্তু ভুলগুলো চুপচাপ থাকে: কোড চলে, ছবি সুন্দর, অথচ সংখ্যা ভুল। নিচের ভুলগুলো বারবার দেখা যায়:

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

• একটি ছোট উদাহরণে (৪টি বিন্দু) হাতের হিসাব, নিজের NumPy সূত্র আর লাইব্রেরি তিনটি মিলিয়েছি
• Bandwidth বা K বেছেছি held-out log-likelihood দিয়ে, training likelihood দিয়ে নয়, আর grid-এর কিনারায় সেরা মান পড়েনি কি না দেখেছি
• ধনাত্মক বা সীমাবদ্ধ ডেটায় সীমানার বাইরে কতটা ভর পড়ছে মেপেছি, আর রূপান্তর করলে Jacobian দিয়েছি
• Density-র integral ১ আর নমুনার পরিসর সংখ্যাগতভাবে যাচাই করেছি
• KDE-কে একটি সরল parametric baseline আর একটি GMM-এর সঙ্গে কয়েকটি seed-এ তুলনা করেছি
• কোনো আসল গ্রাহকের ব্যক্তিগত ডেটা AI-তে পেস্ট করিনি; দরকার হলে বানানো নমুনা দিয়ে প্রশ্ন করেছি