Markov chain আর HMM — অবস্থা থেকে অবস্থায়

আগের তিনটি পর্বে (পর্ব ৫০ Bayes-এর উপপাদ্য, ৫১ Naive Bayes, ৫২ Bayesian network) অনিশ্চয়তা ছিল এক মুহূর্তের ছবি: প্রমাণ এল, বিশ্বাস বদলাল। এই পর্বে নতুন একটি মাত্রা যোগ হচ্ছে: সময়. এমন ক্রম নিয়ে যুক্তির সবচেয়ে পুরোনো আর এখনো কাজের হাতিয়ার হলো Markov chain আর তার লুকানো রূপ, hidden Markov model. সম্ভাবনার মৌলিক নিয়ম (conditional probability, total probability) আছে পর্ব ৩৫-এ, আর matrix গুণ ও eigenvalue পর্ব ১৭-এ; এখানে শুধু দরকারটুকু মনে করিয়ে দেওয়া হবে। এই পর্বের সব সম্ভাবনা বাজারঘর (এই কোর্সের কাল্পনিক অনলাইন দোকান) নিয়ে বানানো, কোনো আসল ডেটা নয়।

১ · ক্রম আর Markov property — “ভবিষ্যৎ শুধু বর্তমানের ওপর”

অনেক ডেটা আসে ক্রম (sequence) হয়ে: প্রতি মাসে একজন ক্রেতার অবস্থা, প্রতি ঘণ্টায় রাস্তার অবস্থা, একটি বাক্যের শব্দগুলো। এমন ক্রমের পুরো সম্ভাবনা লিখতে পর্ব ৩৫-এর chain rule লাগে:

P(X₁, X₂, …, X_T) = P(X₁) · P(X₂ | X₁) · P(X₃ | X₁, X₂) ⋯ P(X_T | X₁, …, X_T−1)

এখানে X_t মানে সময় t-তে অবস্থা, আর “|” মানে “জানা থাকলে”। সমস্যা হলো শেষের পদগুলো: ১২ মাসের ইতিহাসে প্রতিটি মাসে ৪টি সম্ভাব্য অবস্থা থাকলে ইতিহাসের রকম 4¹¹, ৪১ লাখের বেশি। প্রতিটির জন্য আলাদা সম্ভাবনা মাপার মতো ডেটা কারও কাছে নেই।

১৯০৬ সালের দিকে রুশ গণিতবিদ Andrey Markov একটি সাহসী সরলীকরণ নিয়ে কাজ শুরু করেন, যা আজ তাঁর নামে পরিচিত। Markov property বলে:

P(X_t+1 | X₁, …, X_t) = P(X_t+1 | X_t)

কথায়: বর্তমান অবস্থা জানলে অতীত আর কোনো বাড়তি খবর দেয় না। তখন পুরো ক্রমের সম্ভাবনা হয়ে যায় P(X₁) আর পাশাপাশি জোড়ার সম্ভাবনার গুণফল, আর মাপতে হয় শুধু “এই অবস্থা থেকে ওই অবস্থায়” যাওয়ার সম্ভাবনা। এমন অনুমান মেনে চলা ক্রমের মডেলকে বলে Markov chain.

কখন অনুমানটি যুক্তিসংগত, কখন নয়

লক্ষ করুন, “Markov কি না” অনেকটা নির্ভর করে আপনি অবস্থা কীভাবে সংজ্ঞায়িত করছেন তার ওপর। অবস্থায় যত বেশি তথ্য রাখবেন (যেমন “গত দুই মাসের অবস্থা” একসঙ্গে), অনুমান তত সত্যের কাছে যাবে, কিন্তু অবস্থার সংখ্যা তত বাড়বে আর মাপতে তত বেশি ডেটা লাগবে। এ ছাড়া এই পর্বে ধরে নেওয়া হবে যে স্থানান্তরের সম্ভাবনা সময়ের সঙ্গে বদলায় না (time-homogeneous chain).

২ · বাজারঘরের ক্রেতার Markov chain — ছবি আর transition matrix

বাজারঘরের ব্যবস্থাপক প্রতি মাসের শেষে প্রত্যেক ক্রেতাকে চারটি দলের একটিতে রাখেন: Active (এ মাসে একাধিকবার কিনেছেন), Occasional (একবার কিনেছেন), Dormant (কেনেননি, কিন্তু অ্যাকাউন্ট আছে) আর Churned (তিন মাস কোনো সাড়া নেই, মানে ছেড়ে গেছেন বলে ধরা)। “ফিরে আসুন” ছাড়ের অফার পেয়ে কিছু Churned ক্রেতা আবার ফেরেন। নিচের সম্ভাবনাগুলো উদাহরণের জন্য বানানো।

একই তথ্য একটি টেবিলে লিখলে তাকে বলে transition matrix P. সারি মানে “এই মাসের অবস্থা”, কলাম মানে “পরের মাসের অবস্থা”, আর ঘর P_ij = P(X_t+1 = j | X_t = i).

এক ধাপ মানে vector × matrix

শুধু একজন ক্রেতা নয়, “কোন অবস্থায় থাকার সম্ভাবনা কত”, এটিকে একটি সারি-vector x_t দিয়ে লেখা যায়। আজ ক্রেতাটি নিশ্চিতভাবে Active, তাই x₀ = [1, 0, 0, 0]. পরের মাসের vector:

x_t+1 = x_t P,   মানে x_t+1[j] = ∑_i x_t[i] · P_ij

এটি পর্ব ৩৫-এর total probability: পরের মাসে j-তে পৌঁছানোর সব পথ (প্রতিটি সম্ভাব্য এখনকার অবস্থা i থেকে) যোগ। প্রথম মাসে x₁ = [1, 0, 0, 0] P = P-এর প্রথম সারি = [0.70, 0.20, 0.08, 0.02]. দ্বিতীয় মাসে কী হয়, ধাপে ধাপে দেখুন। এটি পর্ব ১৭-এর matrix-vector গুণ, শুধু vector বাঁ দিকে।

1x1 = [0.70, 0.20, 0.08, 0.02]2for j in range(4):3    x2[j] = sum(x1[i] * P[i, j])4print(x2, x2.sum())

৩ · n ধাপের সম্ভাবনা — matrix power আর absorbing state

দুই ধাপ পরের vector x₂ = x₁ P = x₀ P P = x₀ P². সাধারণভাবে x_n = x₀ Pⁿ, আর Pⁿ-এর ঘর (i, j) মানে “i থেকে শুরু করে ঠিক n ধাপ পরে j-তে থাকার সম্ভাবনা”। এই সম্পর্কটির নাম Chapman–Kolmogorov সমীকরণ: P^m+n = P^m Pⁿ.

হাতে ২ ধাপ: একটি ঘর

Dormant একজন ক্রেতা ২ মাস পরে Active হবেন, তার সম্ভাবনা (P²)_D,A. মাঝের মাসে তিনি চারটির যেকোনোটিতে থাকতে পারেন:

• D → A → A: 0.10 × 0.70 = 0.070
• D → O → A: 0.20 × 0.30 = 0.060
• D → D → A: 0.50 × 0.10 = 0.050
• D → C → A: 0.20 × 0.10 = 0.020

যোগ 0.200. মানে আজ ঘুমিয়ে থাকা ক্রেতার ২ মাস পরে আবার নিয়মিত হওয়ার সম্ভাবনা ২০%, যদিও সরাসরি এক মাসে মাত্র ১০%. পুরো P² (লেখকের যন্ত্রে NumPy দিয়ে হিসাব করা):

            Active  Occas.  Dormant Churned
Active      0.5600  0.2370  0.1470  0.0560
Occasional  0.3600  0.2725  0.2515  0.1160
Dormant     0.2000  0.2100  0.3180  0.2720
Churned     0.1700  0.0900  0.0855  0.6545

প্রথম সারিটি stepper-এর x₂, আর Dormant সারির প্রথম ঘর হাতের 0.200. P²-এর প্রতিটি সারির যোগও 1: দুটি stochastic matrix-এর গুণফলও stochastic. NumPy-তে Pⁿ লেখা হয় np.linalg.matrix_power(P, n) দিয়ে; P ** n লিখলে প্রতিটি ঘরের আলাদা ঘাত হয়, যা একেবারে ভুল উত্তর দেয় (পর্ব ১৭-এর * বনাম @-এর ফাঁদ)।

Absorbing state — অফার বন্ধ হলে

ধরুন বাজারঘর “ফিরে আসুন” অফার বন্ধ করল, আর Churned ক্রেতা আর ফেরেন না। তখন Churned সারি হয় [0, 0, 0, 1]: একবার ঢুকলে আর বের হওয়া যায় না। এমন অবস্থাকে বলে absorbing state, আর বাকি অবস্থাগুলোকে transient (ক্ষণস্থায়ী)। যেকোনো transient অবস্থা থেকে absorbing-এ পৌঁছানোর পথ থাকলে শেষ পর্যন্ত সব সম্ভাবনা সেখানেই জমে: লেখকের যন্ত্রে হিসাবে এই নতুন chain-এ Active থেকে শুরু করে ২৪ মাস পরে Churned-এ থাকার সম্ভাবনা প্রায় 0.82.

কত মাস লাগবে গড়ে? transient অংশের ৩ × ৩ matrix Q নিয়ে (Active, Occasional, Dormant-এর মধ্যে স্থানান্তর) N = (I − Q)⁻¹ বানাতে হয়, যাকে fundamental matrix বলে; এখানে I হলো identity matrix আর ⁻¹ মানে inverse. N-এর ঘর (i, j) হলো i থেকে শুরু করে শোষিত হওয়ার আগে j-তে গড়ে কত মাস কাটে। প্রতিটি সারির যোগই শোষিত হওয়ার গড় সময়। এই উদাহরণে তা Active থেকে প্রায় ১৫.২৬ মাস, Occasional থেকে ১৩.৬৮, Dormant থেকে ১০.৫৩ (অনুশীলন ৪-এ কোড আর simulation দিয়ে যাচাই)।

৪ · Stationary distribution — তিন পথে একই উত্তর

অফারসহ আসল chain-এ ফিরে আসি। Active থেকে শুরু করলে প্রথম মাসগুলোতে vector দ্রুত বদলায়, কিন্তু অনেক মাস পরে তা একটি নির্দিষ্ট vector-এর কাছে স্থির হয়, আর শুরু কোথা থেকে হয়েছিল তা প্রায় ভুলে যায়। সেই vector π-কে বলে stationary distribution: এমন একটি distribution, যাকে P দিয়ে গুণ করলে বদলায় না।

π P = π,   ∑_i π_i = 1

তিনটি পথ

1. বারবার গুণ (power iteration): যেকোনো x₀ থেকে x ← x P অনেকবার। চার অবস্থার এই ছোট chain-এ ২০০ বার যথেষ্টের চেয়ে বেশি।
2. Eigenvector: π P = π মানে π হলো P-এর বাঁ দিকের eigenvector, যার eigenvalue 1. NumPy-র np.linalg.eig ডান দিকের eigenvector দেয় (A v = λ v), তাই P-এর transpose নিতে হয়: P^T π^T = π^T. তারপর eigenvector-কে তার যোগ দিয়ে ভাগ করে distribution বানাতে হয়, কারণ eigenvector-এর দৈর্ঘ্য বা চিহ্ন নির্দিষ্ট নয়।
3. Simulation: একজন কাল্পনিক ক্রেতাকে অনেক মাস ধরে random ভাবে চালিয়ে দেখা, কত ভাগ মাস তিনি কোন অবস্থায় কাটালেন। এটি পর্ব ৩৬-এর Monte Carlo.

নিচের কোড প্রথমে অংশ ২ আর ৩-এর সংখ্যা মেলায়, তারপর তিন পথে π বের করে, আর শেষে simulation-এর প্রথম ১২ মাসের পথ (অক্ষরের প্রথম অক্ষর: A, O, D, C) দেখায়। প্রথমবার চালালে ব্রাউজারে NumPy নামতে কয়েক সেকেন্ড লাগতে পারে; SciPy থাকা ব্লকে আরও কয়েক সেকেন্ড।

import numpy as np
from bisect import bisect_right
np.set_printoptions(precision=4, suppress=True)

states = ["Active", "Occasional", "Dormant", "Churned"]
P = np.array([[0.70, 0.20, 0.08, 0.02],      # সারি = এই মাস
              [0.30, 0.40, 0.25, 0.05],      # কলাম = পরের মাস
              [0.10, 0.20, 0.50, 0.20],
              [0.10, 0.05, 0.05, 0.80]])
print("সারির যোগ:", P.sum(axis=1))
x0 = np.array([1.0, 0, 0, 0])                 # ক্রেতা এখন Active
print("১ মাস পরে:", x0 @ P)
print("২ মাস পরে:", x0 @ np.linalg.matrix_power(P, 2))
for n in (6, 12, 50):
    print(f"P^{n}-এর Active সারি:", np.linalg.matrix_power(P, n)[0])

x = x0.copy()                                 # পথ ১: বারবার গুণ
for _ in range(200):
    x = x @ P
print("বারবার গুণ :", x)
w, V = np.linalg.eig(P.T)                     # পথ ২: πP = π
v = V[:, np.argmin(np.abs(w - 1))].real
print("eigenvalue  :", np.sort(w.real)[::-1])
print("eigenvector :", v / v.sum())

rng = np.random.default_rng(42)               # পথ ৩: simulation
cum = P.cumsum(axis=1).tolist()
s, counts, path = 0, np.zeros(4), []
for t, u in enumerate(rng.random(20_000)):
    s = min(bisect_right(cum[s], u), 3)
    counts[s] += 1
    if t < 12:
        path.append(states[s][0])
print("প্রথম ১২ মাস:", " ".join(path))
print("simulation  :", counts / counts.sum())

সারির যোগ: [1. 1. 1. 1.]
১ মাস পরে: [0.7  0.2  0.08 0.02]
২ মাস পরে: [0.56  0.237 0.147 0.056]
P^6-এর Active সারি: [0.3909 0.2201 0.1961 0.1929]
P^12-এর Active সারি: [0.3553 0.2031 0.1855 0.2561]
P^50-এর Active সারি: [0.3499 0.1998 0.1827 0.2676]
বারবার গুণ : [0.3499 0.1998 0.1827 0.2676]
eigenvalue  : [1.     0.7251 0.4979 0.177 ]
eigenvector : [0.3499 0.1998 0.1827 0.2676]
প্রথম ১২ মাস: O O D D A D D D O O O D
simulation  : [0.3513 0.2042 0.1834 0.2611]

Ergodicity — কখন π একটিই, আর সবাই সেখানে পৌঁছায়

সব chain-এ এমন সুন্দর ফল আসে না। সসীম অবস্থার একটি chain-কে ergodic বলা হয়, সহজ কথায়, যখন দুটি শর্ত মেলে:

• সব অবস্থা থেকে সব অবস্থায় পৌঁছানো যায় (কয়েক ধাপে হলেও); একে irreducible বলে। Absorbing Churned থাকলে এটি ভাঙে: সেখান থেকে আর কোথাও যাওয়া যায় না।
• কোনো বাঁধা চক্র নেই; একে aperiodic বলে। উদাহরণ: P = [[0, 1], [1, 0]] হলে অবস্থা প্রতিবার উল্টায়: A, B, A, B… π = [0.5, 0.5] আছে, কিন্তু Active থেকে শুরু করা vector কখনো সেখানে স্থির হয় না, শুধু দোলে।

এই দুটি শর্ত মানলে π একটিই, সব শুরু থেকেই x_n তার দিকে যায়, আর একটি লম্বা simulation-এ প্রতিটি অবস্থায় কাটানো সময়ের ভাগও π-এর কাছে যায়। বাজারঘরের chain-এর সব ঘর ধনাত্মক, তাই দুটি শর্তই মেলে।

৫ · PageRank — বিখ্যাত একটি Markov chain

Stationary distribution-এর সবচেয়ে বিখ্যাত প্রয়োগগুলোর একটি এসেছে web search থেকে। ১৯৯৮ সালে Sergey Brin আর Lawrence Page তাঁদের প্রবন্ধ “The Anatomy of a Large-Scale Hypertextual Web Search Engine”-এ (WWW7 সম্মেলন) PageRank বর্ণনা করেন। ধারণাটি হলো একজন random surfer: তিনি একটি পাতায় আছেন, পাতার link-গুলোর একটিতে random ভাবে চাপেন, তারপর আবার, এভাবে চলতে থাকেন। মাঝে মাঝে বিরক্ত হয়ে তিনি link ছেড়ে একেবারে যেকোনো একটি পাতায় চলে যান। অনেকক্ষণ পরে তিনি কোন পাতায় কত ভাগ সময় কাটান, সেটিই সেই পাতার PageRank. মানে এটি একটি random walk-এর stationary distribution.

একটি ছোট বানানো উদাহরণ: বাজারঘরের সাইটের চারটি পাতা। home থেকে offers আর blog-এ link; offers থেকে শুধু home-এ; blog থেকে home, offers আর help-এ; help থেকে offers-এ। link-গুলোর মধ্যে সম্ভাবনা সমান ভাগ হয়: home-এর দুটি link, তাই প্রতিটি ½; blog-এর তিনটি, প্রতিটি ⅓. এভাবে পাওয়া matrix M. তারপর damping factor d যোগ করা হয়:

G = d · M + (1 − d) / n

এখানে n পাতার সংখ্যা (৪), আর d সম্ভাবনায় surfer link ধরে চলেন, 1 − d সম্ভাবনায় যেকোনো পাতায় লাফ দেন। Brin আর Page লিখেছিলেন যে তাঁরা সাধারণত d = 0.85 নিতেন। প্রবন্ধে সূত্রটি একটু ভিন্নভাবে লেখা (n দিয়ে ভাগ ছাড়া), তবে সেখানেই বলা আছে যে PageRank-গুলো পাতার ওপর একটি সম্ভাবনার distribution হিসেবে দেখা যায়; এখানে সেই সম্ভাবনার রূপটি নেওয়া হয়েছে। লাফের অংশটি বদ্ধ পাতার দল আর চক্রের সমস্যা সারায়। G-এর সব ঘর ধনাত্মক, তাই অংশ ৪-এর ergodicity মেলে।

import numpy as np
np.set_printoptions(precision=4, suppress=True)

pages = ["home", "offers", "blog", "help"]
links = {0: [1, 2], 1: [0], 2: [0, 1, 3], 3: [1]}   # কোন পাতা কোথায় link দেয়
n = len(pages)
M = np.zeros((n, n))
for i, outs in links.items():
    M[i, outs] = 1 / len(outs)                # link-গুলোর মধ্যে সমান ভাগ

def pagerank(d, tol=1e-10):
    G = d * M + (1 - d) / n                   # damping: (1-d) সম্ভাবনায় লাফ
    r = np.full(n, 1 / n)
    for k in range(1, 500):
        new = r @ G
        if np.abs(new - r).sum() < tol:
            return new, k, G
        r = new

r, k, G = pagerank(0.85)
print("G (d = 0.85):\n", G)
print("সারির যোগ:", G.sum(axis=1))
print(f"power iteration, {k} ধাপে:", r)
w, V = np.linalg.eig(G.T)
v = V[:, np.argmin(np.abs(w - 1))].real
print("eigenvector দিয়ে:        ", v / v.sum())
print("ক্রম:", [pages[i] for i in np.argsort(-r)])
for d in (0.5, 0.99):
    r2, k2, _ = pagerank(d)
    print(f"d = {d}: {r2} ({k2} ধাপ)")

G (d = 0.85):
 [[0.0375 0.4625 0.4625 0.0375]
 [0.8875 0.0375 0.0375 0.0375]
 [0.3208 0.3208 0.0375 0.3208]
 [0.0375 0.8875 0.0375 0.0375]]
সারির যোগ: [1. 1. 1. 1.]
power iteration, 66 ধাপে: [0.3763 0.3328 0.1974 0.0934]
eigenvector দিয়ে:         [0.3763 0.3328 0.1974 0.0934]
ক্রম: ['home', 'offers', 'blog', 'help']
d = 0.5: [0.3182 0.3182 0.2045 0.1591] (26 ধাপ)
d = 0.99: [0.3985 0.3334 0.1997 0.0684] (120 ধাপ)

৬ · n-gram language model — বাংলা লেখার Markov chain

অবস্থা যদি শব্দ হয়? তাহলে Markov chain একটি language model: এখনকার শব্দ দেখে পরের শব্দের সম্ভাবনা দেয়। শুধু আগের একটি শব্দ দেখলে তাকে bigram model বলে; আগের দুটি দেখলে trigram; সাধারণভাবে আগের n − 1টি শব্দ দেখলে n-gram model. Bigram-এ সম্ভাবনা মাপা হয় সরাসরি গুনে:

P(b | a) = count(a, b) / count(a, যেকোনো পরের শব্দ)

এটি আসলে maximum likelihood-এর অনুমান। Transition matrix-এর প্রতিটি সারি একটি শব্দ, প্রতিটি কলাম তার পরের সম্ভাব্য শব্দ। নতুন লেখা তৈরি করতে একটি শুরুর শব্দ (seed word) দিয়ে, প্রতিটি ধাপে সারির distribution থেকে পরের শব্দ sample করতে হয়, ঠিক অংশ ৪-এর simulation-এর মতো।

নিচের অনুচ্ছেদটি এই পর্বের জন্য লেখা, কোথাও থেকে নেওয়া নয়। শব্দ ভাঙতে একটি সরল regex ব্যবহার হয়েছে: ফাঁকা জায়গা, দাঁড়ি (।) আর কমা দিয়ে ভাঙা, আর দাঁড়ি-কমাকেও আলাদা token রাখা, যাতে model বাক্যের শেষ “শেখে”। বাংলা লেখার tokenization-এর পুরো আলোচনা পর্ব ৬৭-এ। কোডটি গণনা দেখায়, কয়েকটি সম্ভাবনা হিসাব করে, তারপর একই seed word “সুমি” আর তিনটি ভিন্ন random seed দিয়ে নতুন লেখা বানায়।

import re
from collections import Counter, defaultdict
import numpy as np
text = ("সকালে বাজারঘরের গুদামে আলো জ্বলে। সুমি আপা চা নিয়ে অর্ডারের খাতা "
        "খোলেন। খাতা দেখে তিনি প্যাকেট সাজান। প্যাকেট সাজানো হলে কুরিয়ার "
        "ভাই আসেন। কুরিয়ার ভাই প্যাকেট নিয়ে শহরের পথে যান। পথে জ্যাম "
        "থাকলে ডেলিভারি দেরি হয়। দেরি হলে ক্রেতা ফোন করেন। ফোন পেয়ে "
        "সুমি আপা খাতা দেখেন আর হাসিমুখে বলেন, প্যাকেট পথে আছে। বিকেলে "
        "আবার নতুন অর্ডার আসে, আর সুমি আপা আবার চা নিয়ে খাতা খোলেন।")
tokens = re.findall(r"[^।,\s]+|[।,]", text)   # শব্দ আর দাঁড়ি-কমা আলাদা token
nxt = defaultdict(Counter)
for a, b in zip(tokens, tokens[1:]):
    nxt[a][b] += 1                             # bigram গণনা
print("token:", len(tokens), "| আলাদা:", len(set(tokens)))
print("'প্যাকেট'-এর পরে:", dict(nxt["প্যাকেট"]))

def prob(b, a):                                # P(b | a)
    return nxt[a][b] / sum(nxt[a].values())

print("P(সাজান | প্যাকেট) =", round(prob("সাজান", "প্যাকেট"), 3))
print("P(খাতা দেখে তিনি) =", round(prob("দেখে", "খাতা") * prob("তিনি", "দেখে"), 3))
print("P(প্যাকেট দেরি) =", prob("দেরি", "প্যাকেট"))

def generate(word, n, seed):
    rng = np.random.default_rng(seed)
    out = [word]
    for _ in range(n):
        words = list(nxt[word])
        counts = np.array(list(nxt[word].values()), dtype=float)
        word = str(rng.choice(words, p=counts / counts.sum()))
        out.append(word)
    return " ".join(out).replace(" ।", "।").replace(" ,", ",")

for seed in (1, 2, 3):
    print(f"seed {seed}:", generate("সুমি", 22, seed))

token: 77 | আলাদা: 45
'প্যাকেট'-এর পরে: {'সাজান': 1, 'সাজানো': 1, 'নিয়ে': 1, 'পথে': 1}
P(সাজান | প্যাকেট) = 0.25
P(খাতা দেখে তিনি) = 0.25
P(প্যাকেট দেরি) = 0.0
seed 1: সুমি আপা আবার নতুন অর্ডার আসে, আর হাসিমুখে বলেন, আর সুমি আপা আবার নতুন অর্ডার আসে, প্যাকেট সাজানো হলে কুরিয়ার
seed 2: সুমি আপা চা নিয়ে অর্ডারের খাতা দেখে তিনি প্যাকেট সাজানো হলে ক্রেতা ফোন করেন। কুরিয়ার ভাই প্যাকেট নিয়ে শহরের পথে জ্যাম থাকলে
seed 3: সুমি আপা চা নিয়ে শহরের পথে জ্যাম থাকলে ডেলিভারি দেরি হয়। পথে জ্যাম থাকলে ডেলিভারি দেরি হয়। দেরি হয়। বিকেলে

৭ · Hidden Markov model — লুকানো ট্রাফিক, দেখা যায় শুধু ডেলিভারির সময়

বাজারঘরের একজন রাইডার প্রতি ঘণ্টায় একটি করে ডেলিভারি দেন। ডেলিভারির সময় তিন ভাগে লেখা হয়: fast, normal, slow. সময়টি আসলে নির্ভর করে রাস্তার অবস্থার ওপর: Free (ফাঁকা), Heavy (ভারী), Jam (জ্যাম)। কিন্তু অফিসে বসে কেউ রাস্তা দেখতে পান না; দেখা যায় শুধু ডেলিভারির সময়। রাস্তার অবস্থা নিজে একটি Markov chain (জ্যাম থাকলে পরের ঘণ্টাতেও জ্যাম থাকার সম্ভাবনা বেশি), আর প্রতিটি অবস্থা একটি সম্ভাবনায় একটি ফল “নির্গত” করে। এই কাঠামোকে বলে hidden Markov model (HMM).

একটি HMM-এ তিন সেট সম্ভাবনা লাগে (এখানে সব বানানো):

B-এর ঘর হলো emission probability, P(দেখা ফল | লুকানো অবস্থা)। A আর B দুটিরই প্রতিটি সারির যোগ 1. HMM-এর দুটি অনুমান: (১) লুকানো অবস্থা Markov property মানে; (২) এখনকার দেখা ফল শুধু এখনকার লুকানো অবস্থার ওপর নির্ভর করে, আগের ফল বা অবস্থার ওপর নয়। তাহলে একটি নির্দিষ্ট লুকানো পথ Z আর দেখা ক্রম O-এর যৌথ সম্ভাবনা সরল গুণফল:

P(Z, O) = π_z1 B_z1(o₁) · A_{z1 z2} B_z2(o₂) ⋯ A_{z(T−1) zT} B_zT(o_T)

উদাহরণ: তিন ঘণ্টায় দেখা গেল fast, slow, slow. পথ Free, Jam, Jam ধরলে: 0.5 × 0.7 × 0.1 × 0.6 × 0.5 × 0.6 = 0.0063. এটি পরে Viterbi-র উত্তর হিসেবে ফিরে আসবে।

HMM-এর তিনটি মূল প্রশ্ন

Lawrence Rabiner-এর ১৯৮৯ সালের বহুল পঠিত tutorial (“A Tutorial on Hidden Markov Models and Selected Applications in Speech Recognition”, Proceedings of the IEEE) HMM-এর কাজকে তিনটি প্রশ্নে সাজায়, আর এই ভাগটি আজও পাঠ্যবইয়ে প্রচলিত:

৮ · Forward algorithm আর Viterbi — trellis ভরে ভরে

P(O) বের করার সরাসরি পথ হলো সব লুকানো পথের P(Z, O) যোগ করা। তিন ঘণ্টায় 3³ = 27টি পথ, কিন্তু ১০০ ঘণ্টায় 3¹⁰⁰, যা কোনো কম্পিউটারে গোনা সম্ভব নয়। Forward algorithm মূল কৌশলটি নেয় dynamic programming থেকে: প্রতিটি সময় আর অবস্থার জন্য একটি মাত্র সংখ্যা রাখা,

α_t(j) = P(o₁, …, o_t, Z_t = j)

মানে “প্রথম t-টি ফল দেখা গেছে আর সময় t-তে অবস্থা j”. আগের কলাম থেকে পরেরটি:

α₁(j) = π_j B_j(o₁),   α_t+1(j) = [∑_i α_t(i) A_ij] · B_j(o_t+1),   P(O) = ∑_j α_T(j)

বন্ধনীর ভেতরটি অংশ ২-এর vector × matrix, শুধু তারপর দেখা ফলের সম্ভাবনা দিয়ে গুণ। সময় × অবস্থার এই ছকটিকে বলে trellis (জালির মতো দেখতে বলে)। খরচ N অবস্থা আর T সময়ে মোটামুটি N²T, সূচকীয় N^T নয়। দেখা ক্রম fast, slow, slow-এর জন্য trellis ভরুন:

1alpha = pi * B[:, fast]2for o in [slow, slow]:3    alpha = (alpha @ A) * B[:, o]4print(alpha.sum())

Viterbi — যোগের বদলে সর্বোচ্চ, সঙ্গে পেছনের তীর

Decoding প্রশ্নে চাই একটি পথ, যোগফল নয়। Viterbi algorithm forward-এর হুবহু কাঠামো নেয়, শুধু ∑-এর জায়গায় max বসায়:

δ_t+1(j) = max_i [δ_t(i) A_ij] · B_j(o_t+1),   ψ_t+1(j) = argmax_i [δ_t(i) A_ij]

δ_t(j) হলো “সময় t-তে j-তে শেষ হওয়া সেরা একক পথের সম্ভাবনা”, আর ψ হলো backpointer: সেই সেরা পথ আগের ধাপে কোন অবস্থা থেকে এসেছিল। শেষে সবচেয়ে বড় δ থেকে শুরু করে তীর ধরে পেছনে হাঁটলে পুরো পথ পাওয়া যায়। Andrew Viterbi ১৯৬৭ সালে এটি প্রকাশ করেন error-correcting code-এর decoding-এর জন্য; পরে এটি HMM-এর প্রধান হাতিয়ার হয়ে ওঠে।

এবার দুটি algorithm কোডে, প্রতিটি কলামের মান ছাপিয়ে। শেষ অংশে ১০০০ ঘণ্টার একটি random ক্রমে সরাসরি গুণ বনাম log-space: সরাসরি গুণে সংখ্যাগুলো এত ছোট হয় যে float64-এ underflow হয়ে 0.0 হয়ে যায় (পর্ব ৩৫ আর পর্ব ৩৮-এর সমস্যা)। সমাধান হলো log α রাখা; কিন্তু forward-এ যোগ আছে, আর log-এর যোগ সরাসরি করা যায় না। তাই log-sum-exp কৌশল লাগে: log(e^a + e^b) = m + log(e^a−m + e^b−m), যেখানে m = max(a, b); SciPy-তে scipy.special.logsumexp. Viterbi-তে কাজটি আরও সহজ: max আর গুণ log-এ হয়ে যায় max আর যোগ।

import numpy as np
from scipy.special import logsumexp
np.set_printoptions(precision=6, suppress=True)
S = ["Free", "Heavy", "Jam"]                       # লুকানো অবস্থা
O = {"fast": 0, "normal": 1, "slow": 2}           # দেখা যায় এমন ফল
pi = np.array([0.5, 0.3, 0.2])
A = np.array([[0.6, 0.3, 0.1], [0.3, 0.4, 0.3], [0.1, 0.4, 0.5]])
B = np.array([[0.7, 0.2, 0.1], [0.2, 0.5, 0.3], [0.1, 0.3, 0.6]])
obs = [O[o] for o in ["fast", "slow", "slow"]]

alpha = pi * B[:, obs[0]]                          # forward
print("alpha1:", alpha)
for t, o in enumerate(obs[1:], start=2):
    alpha = (alpha @ A) * B[:, o]
    print(f"alpha{t}:", alpha)
print("P(O) =", round(alpha.sum(), 6))

delta, back = pi * B[:, obs[0]], []                # Viterbi
for o in obs[1:]:
    scores = delta[:, None] * A                    # [i, j]: i থেকে j-তে
    back.append(scores.argmax(axis=0))
    delta = scores.max(axis=0) * B[:, o]
    print("delta:", delta, "| পেছনের তীর:", [S[i] for i in back[-1]])
path = [int(delta.argmax())]
for bp in reversed(back):
    path.insert(0, int(bp[path[0]]))
print("Viterbi path:", [S[i] for i in path], "| P =", round(delta.max(), 6))

rng = np.random.default_rng(0)                     # লম্বা ক্রম: underflow
seq = rng.integers(0, 3, size=1000)
a, la = pi * B[:, seq[0]], np.log(pi * B[:, seq[0]])
for o in seq[1:]:
    a = (a @ A) * B[:, o]
    la = logsumexp(la[:, None] + np.log(A), axis=0) + np.log(B[:, o])
print("১০০০ ধাপ, সরাসরি:", a.sum(), "| log-space:", round(logsumexp(la), 2))

alpha1: [0.35 0.06 0.02]
alpha2: [0.023  0.0411 0.0378]
alpha3: [0.002991 0.011538 0.020118]
P(O) = 0.034647
delta: [0.021  0.0315 0.021 ] | পেছনের তীর: ['Free', 'Free', 'Free']
delta: [0.00126 0.00378 0.0063 ] | পেছনের তীর: ['Free', 'Heavy', 'Jam']
Viterbi path: ['Free', 'Jam', 'Jam'] | P = 0.0063
১০০০ ধাপ, সরাসরি: 0.0 | log-space: -1114.76

৯ · শেখা, hmmlearn-এর সত্যিকারের run আর HMM-এর জায়গা

Baum–Welch — শুধু দেখা ক্রম থেকে শেখা (ধারণা)

এতক্ষণ π, A, B জানা ধরে নিয়েছি। বাস্তবে বাজারঘরের কাছে থাকে শুধু হাজারো ডেলিভারির সময়, রাস্তার আসল অবস্থা কোথাও লেখা নেই। লুকানো অবস্থা জানা থাকলে কাজটি সহজ হতো: শুধু গুনে ভাগ করা, অংশ ৬-এর bigram-এর মতো। জানা নেই বলে Baum–Welch algorithm একটি চক্র চালায়, যা expectation–maximization (EM)-এর একটি রূপ:

1. শুরু: π, A, B-এর একটি অনুমান (প্রায়ই random).
2. E-ধাপ: এখনকার model দিয়ে forward আর backward (পেছন থেকে একই ধরনের হিসাব) চালিয়ে প্রতিটি ঘণ্টায় প্রতিটি অবস্থায় থাকার, আর প্রতিটি স্থানান্তর ঘটার প্রত্যাশিত সংখ্যা বের করা। পূর্ণ গণনার বদলে “আংশিক গণনা”।
3. M-ধাপ: সেই প্রত্যাশিত গণনা দিয়ে নতুন π, A, B, ঠিক যেন গুনে ভাগ করা হচ্ছে।
4. আবার: likelihood আর উল্লেখযোগ্যভাবে না বাড়া পর্যন্ত।

প্রতিটি চক্রে likelihood কমে না, এটি প্রমাণিত; কিন্তু চক্রটি একটি local maximum-এ আটকে যেতে পারে (পর্ব ৪৪-এর hill climbing-এর মতো সমস্যা)। পদ্ধতিটি Leonard Baum আর সহকর্মীদের ১৯৬০-এর দশকের শেষ আর ১৯৭০ সালের প্রবন্ধগুলো থেকে এসেছে; EM-এর সাধারণ কাঠামো ১৯৭৭ সালে Dempster, Laird আর Rubin প্রকাশ করেন। নিচের run-এ এই সমস্যাগুলো চোখে দেখা যাবে।

hmmlearn — সত্যিকারের run

HMM-এর একটি পরিচিত Python library হলো hmmlearn, যা scikit-learn-এর মতো fit/predict ধরনের API দেয়। এটি ব্রাউজারের Python-এ নেই, তাই লেখকের যন্ত্রে একটি আলাদা venv-এ চালানো হয়েছে (সত্যিকারের run, Python 3.14.4, uv 0.11.28; লেখার সময়ে PyPI-র সর্বশেষ hmmlearn 0.3.3). দুটি বিষয় লক্ষ করুন। প্রথমত, এর README-তে প্রকল্পটিকে “limited-maintenance mode”-এ বলা আছে। দ্বিতীয়ত, ক্রমিক চিহ্নের (fast, normal, slow-এর মতো) HMM-এর class এখন CategoricalHMM. পুরোনো টিউটোরিয়ালে এই কাজে MultinomialHMM লেখা থাকে, কিন্তু সেই নাম এখন একটি ভিন্ন model-এর (আসল multinomial distribution); বানালেই এই সতর্কবার্তা আসে:

MultinomialHMM has undergone major changes. The previous version was implementing a CategoricalHMM (a special case of MultinomialHMM). This new implementation follows the standard definition for a Multinomial distribution …

প্রথমে প্রকল্প তৈরি:

$ cd ~/bazarghar-ai
$ uv venv --python 3.14
$ uv pip install hmmlearn
 …
 + hmmlearn==0.3.3
 + numpy==2.5.4
 + scikit-learn==1.9.1
 + scipy==1.18.1
 …

import numpy as np
from hmmlearn import hmm
import hmmlearn

np.set_printoptions(precision=4, suppress=True)
print("hmmlearn", hmmlearn.__version__, "| numpy", np.__version__)

model = hmm.CategoricalHMM(n_components=3)          # Free, Heavy, Jam
model.startprob_ = np.array([0.5, 0.3, 0.2])
model.transmat_ = np.array([[0.6, 0.3, 0.1],
                            [0.3, 0.4, 0.3],
                            [0.1, 0.4, 0.5]])
model.emissionprob_ = np.array([[0.7, 0.2, 0.1],     # fast, normal, slow
                                [0.2, 0.5, 0.3],
                                [0.1, 0.3, 0.6]])

X = np.array([[0], [2], [2]])                         # fast, slow, slow
logp = model.score(X)
print("log P(O) =", round(logp, 4), "| P(O) =", round(np.exp(logp), 6))
vlogp, states = model.decode(X, algorithm="viterbi")
names = np.array(["Free", "Heavy", "Jam"])
print("Viterbi:", names[states], "| P(path, O) =", round(np.exp(vlogp), 6))
print("posterior:\n", model.predict_proba(X))

# Baum-Welch: শুধু observation থেকে parameter শেখা
Xs, Zs = model.sample(3000, random_state=7)
best = None
for seed in range(6):
    m = hmm.CategoricalHMM(n_components=3, n_iter=500, tol=1e-4,
                           random_state=seed).fit(Xs)
    ll = m.score(Xs)
    print(f"seed {seed}: log-lik {ll:.1f}, iter {m.monitor_.iter}")
    if best is None or ll > best.score(Xs):
        best = m
print("true model log-lik:", round(model.score(Xs), 1))
order = np.argsort(-best.emissionprob_[:, 0])      # fast-বেশি অবস্থা আগে
print("learned transmat (reordered):\n", best.transmat_[order][:, order])
print("learned emissionprob (reordered):\n", best.emissionprob_[order])

$ python traffic_hmm.py
hmmlearn 0.3.3 | numpy 2.5.4
log P(O) = -3.3625 | P(O) = 0.034647
Viterbi: ['Free' 'Jam' 'Jam'] | P(path, O) = 0.0063
posterior:
 [[0.6879 0.2135 0.0985]
 [0.1394 0.3915 0.4691]
 [0.0863 0.333  0.5807]]
seed 0: log-lik -3289.7, iter 500
seed 1: log-lik -3289.5, iter 500
seed 2: log-lik -3259.1, iter 500
seed 3: log-lik -3258.3, iter 500
seed 4: log-lik -3256.9, iter 500
seed 5: log-lik -3259.3, iter 500
true model log-lik: -3263.0
learned transmat (reordered):
 [[0.5497 0.2916 0.1587]
 [0.0007 0.8612 0.1382]
 [0.9997 0.0002 0.0001]]
learned emissionprob (reordered):
 [[0.8382 0.1024 0.0594]
 [0.1682 0.4178 0.414 ]
 [0.0407 0.5917 0.3676]]

HMM কোথায় ব্যবহার হয়েছে, আর আজ কোথায় দাঁড়িয়ে

• Speech recognition: ১৯৭০-এর দশক থেকে কয়েক দশক ধরে অনেক বড় speech system-এর কেন্দ্রে ছিল HMM: লুকানো অবস্থা শব্দের ধ্বনি-অংশ, দেখা ফল শব্দতরঙ্গের feature.
• Part-of-speech tagging: লুকানো অবস্থা ব্যাকরণের পদ (বিশেষ্য, ক্রিয়া…), দেখা ফল শব্দ।
• Gene finding আর জীববিজ্ঞান: DNA-র ক্রমে লুকানো অবস্থা “gene-এর অংশ কি না”, আর প্রোটিন পরিবার খোঁজায় profile HMM; ১৯৯০-এর দশকের কয়েকটি পরিচিত gene-খোঁজার টুল HMM-ধরনের model ব্যবহার করত বলে উল্লেখ আছে।

২০১০-এর দশক থেকে অনেক কাজে, বিশেষ করে speech আর ভাষায়, neural sequence model (recurrent network, পরে Transformer) বড় অংশে HMM-কে সরিয়ে দিয়েছে, কারণ তারা দীর্ঘ প্রসঙ্গ আর জটিল feature সরাসরি শিখতে পারে। তবু HMM পুরো হারিয়ে যায়নি: অল্প ডেটায়, যেখানে ব্যাখ্যাযোগ্যতা আর সম্ভাবনার পরিষ্কার মানে দরকার, বা জীববিজ্ঞানের কিছু কাজে এটি এখনো ব্যবহৃত হয় বলে দেখা যায়। কোনটি কোথায় “সেরা”, তা কাজ আর ডেটার ওপর নির্ভর করে।

১০ · AI-এর যুগে — AI সহকারীর লেখা Markov আর HMM কোড যাচাই

“একটি HMM বানিয়ে দাও” বা “stationary distribution বের করো” লিখলে AI সহকারী দ্রুত কোড দেয়, কিন্তু ভুলগুলো প্রায়ই নিঃশব্দ। AI সহকারীরা প্রায়ই যা করে:

• সারি আর কলামের রীতি গুলিয়ে ফেলা: P @ x লিখে ধরে নেওয়া যে এটি পরের distribution, অথচ সারি-stochastic P-তে লাগে x @ P; eigenvector-এ P.T ভুলে যাওয়া, বা eigenvector normalize না করা (ঋণাত্মক “সম্ভাবনা”)।
• পুরোনো বা বানানো API: discrete HMM-এ পুরোনো অর্থে MultinomialHMM; np.matrix_power-এর মতো অস্তিত্বহীন নাম (আসলটি np.linalg.matrix_power); P ** n দিয়ে matrix power.
• লম্বা ক্রমে সরাসরি গুণ: log-space বা scaling ছাড়া forward, ফলে 0.0 বা -inf, তারপর warning চাপা দেওয়া।
• শেখা model নিয়ে অতিরিক্ত দাবি: একটি random শুরুর Baum–Welch-কে “আসল” parameter বলা, শেখা অবস্থাকে নিজে থেকে “Jam” নাম দেওয়া, train ডেটার likelihood দিয়েই model-এর মান বিচার।

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

• প্রতিটি transition আর emission matrix-এর সারির যোগ 1, আর কোডে সারি মানে “এখন”, কলাম মানে “পরে”, এটি নিশ্চিত করেছি
• ২-৩ ধাপের একটি ছোট উদাহরণ হাতে হিসাব করে কোডের ফলের সঙ্গে মিলিয়েছি (যেমন P(fast, slow, slow) = 0.034647)
• Stationary distribution অন্তত দুই পথে (বারবার গুণ আর eigenvector) বের করে মিলিয়েছি; eigenvector normalize করা, সব মান ধনাত্মক
• লম্বা ক্রমে log-space বা scaling আছে; কোনো warning বা nan চাপা দেওয়া নেই
• Baum–Welch কয়েকটি random শুরু দিয়ে চালিয়েছি, আলাদা ডেটায় likelihood দেখেছি, আর শেখা অবস্থার নাম নিজে যাচাই করেছি
• নিজের version-এর অফিসিয়াল docs-এ প্রতিটি অচেনা class আর parameter মিলিয়েছি; আসল ক্রেতার ডেটা AI-তে পেস্ট করিনি