সমাধান — অধ্যায় ৩.৬ · Markov Chains & Introduction to MCMC¶

অধ্যায় ফাইল: part-3-convergence-processes/03-06-markov-chains-mcmc.md (§৭ অনুশীলনী)। সংখ্যাগত উত্তর numpy দিয়ে যাচাইযোগ্য (seed উল্লেখ থাকলে reproducible)। মূল বস্তু — states, transition matrix $P$ (প্রতিটি সারি probability distribution), $\mu_n=\mu_0 P^n$, stationary $\pi=\pi P$, detailed balance $\pi_i P_{ij}=\pi_j P_{ji}$, এবং Metropolis–Hastings। চলমান উদাহরণ: E1 আবহাওয়া-chain $P=\begin{pmatrix}0.8&0.2\\0.4&0.6\end{pmatrix}$ (সারি/কলাম ক্রম Sunny, Rainy), stationary $\pi=(2/3,\,1/3)$; E2 graph random walk; E3 stationary $\pi$; E4 Metropolis MCMC।

ক · ধারণাগত (conceptual)¶

সমাধান ১ (★)¶

Markov property: একটা process $X_0,X_1,X_2,\dots$ Markov হয় যদি ভবিষ্যৎ state-এর distribution শুধু এখনকার state-এর উপর নির্ভর করে, অতীতের পুরো পথের উপর নয়: $$ P(X_{n+1}=j \mid X_n=i,\,X_{n-1}=i_{n-1},\dots,X_0=i_0)=P(X_{n+1}=j\mid X_n=i)=P_{ij}. $$

"অতীত অপ্রাসঙ্গিক" নাকি "অতীত বর্তমানে ধরা আছে"? — দ্বিতীয়টি, এবং এই পার্থক্যটাই সূক্ষ্ম কিন্তু গুরুত্বপূর্ণ। Markov property বলে না অতীতের কোনো মূল্য নেই; বলে অতীত সম্পর্কে যা-কিছু ভবিষ্যতের জন্য দরকারি, তার সবটুকু বর্তমান state-এ সংকুচিত (encoded) আছে। বর্তমান জানলে অতীত আর অতিরিক্ত কিছু বলে না — এটাই "conditional independence of future and past given present"।

Figure 1-এর সাথে: chain যদি আজ Sunny node-এ থাকে, কালকের distribution ঠিক $[P_{SS},P_{SR}]=[0.8,0.2]$ — সেটা গতকাল Sunny ছিল না Rainy, বা গত সপ্তাহে কী ছিল, তার উপর মোটেও নির্ভর করে না। transition graph-এ কোনো node-এ পৌঁছানোর পরে "কোন তীর বেয়ে এলাম" তা ভুলে যাওয়া যায়; শুধু "এখন কোন node-এ আছি" যথেষ্ট।

non-Markov উদাহরণ: এমন কিছু যেখানে প্রবণতা/গতি (momentum) লাগে। যেমন: আজকের আবহাওয়া যদি শুধু গতকালের নয়, গত তিন দিনের ধারা-র উপর নির্ভর করে ("পরপর তিন দিন রোদ থাকলে কাল বৃষ্টির সম্ভাবনা বাড়ে") — তখন একদিনের state ($X_n$) যথেষ্ট নয়, প্রক্রিয়াটা first-order Markov নয়। আরেকটা: কোনো গেমে "যদি আগের দুবার হেরে থাক তবে কৌশল বদলাও" — আগের একটি ফলাফল যথেষ্ট তথ্য নয়। (লক্ষণীয় কৌশল: state-কে বড় করে — যেমন state $=$ "গত $k$ দিনের আবহাওয়া" — প্রায় যেকোনো finite-memory process-কে আবার Markov বানানো যায়; এটাকে বলে state-space augmentation।)

সমাধান ২ (★)¶

stationary distribution $\pi$: এটা states-এর উপর এমন একটা probability distribution যা chain চালালে বদলায় না — একবার $\pi$-তে পৌঁছালে এক ধাপ পরেও distribution $\pi$-ই থাকে। গাণিতিকভাবে $\pi$ একটা row vector যে $$ \pi = \pi P,\qquad \sum_i \pi_i = 1,\quad \pi_i\ge 0. $$

"$\pi=\pi P$" শব্দে: "যদি এই মুহূর্তে state-গুলোর উপর সম্ভাবনা-ভর $\pi$ অনুযায়ী ছড়ানো থাকে, তবে এক step transition চালানোর পরেও ভর ঠিক $\pi$ অনুযায়ীই থাকবে।" অর্থাৎ $\pi$ হলো transition-এর একটা fixed point (ভারসাম্য/equilibrium)। প্রতিটি state থেকে যত ভর বেরিয়ে যায়, ঠিক তত ভর অন্য states থেকে তাতে ঢোকে — তাই net কোনো পরিবর্তন নেই।

stationary-তে পৌঁছে আরেক ধাপ চালালে বদলায় না কেন: সংজ্ঞা অনুসারেই — $\pi_{n}=\pi$ হলে $\pi_{n+1}=\mu_n P=\pi P=\pi$। তাই $\pi$ একটা "শোষক ভারসাম্য": ঢুকে গেলে আর বের হয় না (distribution অর্থে; chain নিজে অবশ্যই state-থেকে-state নড়াচড়া করতেই থাকে, কিন্তু সম্ভাবনাগুলো স্থির)।

"chain শুরু ভুলে যায়" Figure 2-তে কোথায়: তিনটি curve তিন আলাদা শুরু ($[1,0]$, $[0,1]$, $[0.5,0.5]$) থেকে রওনা দিয়েও মাত্র $5$–$6$ ধাপে একই সবুজ রেখা ($\pi_S=2/3$)-তে গিয়ে মিশে স্থির হয়ে যায়। যেহেতু গন্তব্য শুরুর উপর নির্ভর করে না, যথেষ্ট সময় পরে "কোথা থেকে এসেছিলাম" তথ্যটা হারিয়ে যায় — এটাই memory loss / convergence to equilibrium (irreducible + aperiodic chain-এর ergodic ধর্ম)।

সমাধান ৩ (★★)¶

detailed balance ⟹ stationarity, কিন্তু উল্টোটা নয় — স্বজ্ঞা।

detailed balance (DB) একটা local (per-pair) শর্ত: প্রতিটি জোড়া state $(i,j)$-এর জন্য $$ \pi_i P_{ij}=\pi_j P_{ji}. $$ এর ভৌত অর্থ: equilibrium-এ $i\!\to\!j$ দিকে যত probability-প্রবাহ, ঠিক তত $j\!\to\!i$ দিকে — প্রতিটি জোড়ার মধ্যে আলাদা করে প্রবাহ ভারসাম্যপূর্ণ। একে বলে "reversibility": ফিল্ম উল্টো চালালেও পরিসংখ্যান একই দেখায়।

stationarity (global balance) একটা global (per-state) শর্ত: প্রতিটি state $j$-এর জন্য মোট-ঢোকা = মোট-বেরোনো ($\sum_i \pi_i P_{ij}=\pi_j$)।

DB থেকে global আসে কারণ প্রতিটি জোড়া আলাদা ভারসাম্যপূর্ণ হলে, কোনো একটা state-এর চারপাশে সব জোড়ার প্রবাহ যোগ করলেও ভারসাম্য থাকবে (আনুষ্ঠানিক প্রমাণ §Q9)। উল্টোটা সবসময় আসে না: global balance শুধু চায় প্রতিটি state-এ net প্রবাহ শূন্য — কিন্তু এটা চক্রাকার (circulating) প্রবাহ দিয়েও মেটানো যায়। কল্পনা করুন তিনটি state $A\to B\to C\to A$ একমুখী একটা রিং-এ ভর ঘুরছে: প্রতিটি state-এ যত ঢোকে তত বেরোয় (global balance ✓), কিন্তু কোনো জোড়ায় দুদিকে সমান প্রবাহ নেই ($A\!\to\!B$ আছে, $B\!\to\!A$ নেই) — DB ✗। তাই DB কঠোরতর (strictly stronger) শর্ত।

MCMC-তে DB কেন ব্যবহার করি যখন শুধু stationary হলেই চলত: কারণ DB নকশা করা সহজ। আমাদের লক্ষ্য এমন $P$ বানানো যার stationary distribution ঠিক target $\pi$। সরাসরি "$\pi=\pi P$ মেটাও" শর্তে $P$ ডিজাইন করা কঠিন (সব state একসাথে জড়িত)। কিন্তু DB শর্তটা প্রতিটি জোড়ার জন্য আলাদা ও স্থানীয় — Metropolis accept-নিয়মটা ঠিক এমনভাবে বানানো যে প্রতিটি প্রস্তাবিত move স্বয়ংক্রিয়ভাবে DB মানে (§Q11), ফলে target যে stationary হবে তা বিনামূল্যে গ্যারান্টি পাওয়া যায়। অর্থাৎ DB হলো "stationarity পাওয়ার একটা সহজ যথেষ্ট-শর্ত (sufficient condition)" — গাণিতিকভাবে দরকার না হলেও প্রকৌশলগতভাবে অমূল্য।

সমাধান ৪ (★★)¶

(ক) burn-in কেন ফেলি। MCMC chain-কে আমরা একটা নির্বিচার শুরুর state থেকে ছাড়ি, যা সাধারণত target distribution-এর প্রতিনিধি নয় (Figure 3-এ ইচ্ছাকৃতভাবে $x=8$, যেখানে target ঘনত্ব প্রায় শূন্য)। chain target-এ "গড়িয়ে পৌঁছাতে" কিছু সময় নেয় ($\mu_n\to\pi$, ঠিক Figure 2-র মতো)। এই transient পর্বের নমুনাগুলো এখনো stationary distribution থেকে আসছে না — এগুলো রাখলে গড়/histogram-এ শুরুর-state-এর দিকে bias আসে। তাই প্রথম কিছু (যেমন Figure 3-এ $150$, Figure 4-এ $2000$) নমুনা ফেলে দিই (burn-in / warm-up), যাতে বাকিগুলো আনুমানিক target থেকে আসে।

(খ) proposal step-size — mixing। random-walk Metropolis-এ proposal $x'=x+\mathcal N(0,\text{step}^2)$: - খুব ছোট step: প্রতিটি proposal $x$-এর খুব কাছে, তাই প্রায় সবই accept হয় — কিন্তু chain পিঁপড়ের পায়ে নড়ে; এক mode থেকে অন্যে যেতে বহু iteration লাগে, এমনকি আটকে যেতে পারে। ফল: ধীর exploration, পরপর নমুনা প্রবলভাবে সম্পর্কিত (high autocorrelation), কার্যকর নমুনা-সংখ্যা কম। - খুব বড় step: proposal প্রায়ই target-এর কম-ঘনত্ব অঞ্চলে গিয়ে পড়ে, তাই $f(x')/f(x)$ ছোট হয়ে বেশির ভাগ proposal reject হয়; chain একই জায়গায় বহুবার আটকে থাকে (trace-এ লম্বা সমতল ধাপ)। ফলেও high autocorrelation, খারাপ mixing। - মাঝারি step (Figure 3-এর $1.6$): যথেষ্ট দূরে লাফায় (mode-থেকে-mode যেতে পারে) অথচ accept-rate যুক্তিসঙ্গত — এটাই ভালো mixing। (অভিজ্ঞতা-নিয়ম: ১-মাত্রিকে accept-rate $\sim 0.44$, উচ্চ-মাত্রায় $\sim 0.23$ ঘিরে রাখা ভালো।)

(গ) target শুধু অনুপাতে জানলেই চলে কেন। Metropolis accept-probability: $$ \alpha(x,x')=\min!\Big(1,\ \frac{f(x')}{f(x)}\Big). $$ এখানে target ঢোকে কেবল অনুপাত $f(x')/f(x)$ হিসেবে। যদি আসল (normalized) target $\pi(x)=f(x)/Z$ হয় (যেখানে $Z=\int f$), তবে $$ \frac{\pi(x')}{\pi(x)}=\frac{f(x')/Z}{f(x)/Z}=\frac{f(x')}{f(x)}, $$ — $Z$ কাটাকাটি হয়ে যায়। তাই অজানা/দুর্নিরূপণযোগ্য normalizing constant $Z$ লাগেই না; শুধু un-normalized $f$ জানলেই sampler চলে। এই একটা কারণেই MCMC Bayesian posterior-এ এত শক্তিশালী, কারণ সেখানে $Z=\int p(\text{data}\mid\theta)p(\theta)\,d\theta$ (evidence) প্রায়ই হিসাব-অসম্ভব (Part IV, 4.10)।

খ · গণনামূলক (computational)¶

সমাধান ৫ (★)¶

$P=\begin{pmatrix}0.8&0.2\\0.4&0.6\end{pmatrix}$, ক্রম (Sunny, Rainy); আজ Sunny অর্থাৎ state-vector $[1,0]$।

(ক) কাল Rainy: সরাসরি প্রথম সারি থেকে $P_{SR}=0.2$। (অর্থাৎ $[1,0]P=[0.8,0.2]$, Rainy-উপাদান $0.2$।)

(খ) পরশু Rainy $=(P^2)_{SR}$, Chapman–Kolmogorov: $$ (P^2){SR}=\sum. $$ স্বজ্ঞা: দুই দিনে Sunny→Rainy যাওয়ার দুটো পথ — "Sunny→Sunny→Rainy" (}P_{Sk}P_{kR}=P_{SS}P_{SR}+P_{SR}P_{RR}=0.8\cdot0.2+0.2\cdot0.6=0.16+0.12=\boxed{0.28$0.8\times0.2$) আর "Sunny→Rainy→Rainy" ($0.2\times0.6$) — যোগ করলে $0.28$। (যাচাই: পুরো $P^2=\begin{pmatrix}0.72&0.28\\0.56&0.44\end{pmatrix}$, $(P^2)_{SR}=0.28$ ✓।)

সমাধান ৬ (★)¶

$\pi=\pi P$ লিখলে (প্রথম উপাদান): $$ \pi_S=0.8\,\pi_S+0.4\,\pi_R\ \Longrightarrow\ 0.2\,\pi_S=0.4\,\pi_R\ \Longrightarrow\ \pi_S=2\,\pi_R. $$ normalization $\pi_S+\pi_R=1$ বসিয়ে: $2\pi_R+\pi_R=1\Rightarrow \pi_R=\tfrac13,\ \pi_S=\tfrac23$। $$ \boxed{\pi=(\pi_S,\pi_R)=(\tfrac23,\ \tfrac13)\approx(0.667,\ 0.333)}. $$ (দ্বিতীয় সমীকরণ $\pi_R=0.2\pi_S+0.6\pi_R$ একই তথ্য দেয় — $P$-র সারি-যোগ $1$ হওয়ায় সমীকরণ দুটো linearly dependent, তাই normalization লাগে।) Figure 2-র সবুজ ভাঙা-রেখা ঠিক $\pi_S=0.667$-এ — মিলে গেল।

সমাধান ৭ (★★)¶

E2 — path graph A–B–C–D-তে random walk (প্রতি ধাপে এক প্রতিবেশীতে সমান সম্ভাবনায়)।

(ক) transition matrix (ক্রম A,B,C,D)। A-র একমাত্র প্রতিবেশী B (তাই A→B সম্ভাবনা $1$); B-র প্রতিবেশী A,C (প্রতিটি $1/2$); C-র প্রতিবেশী B,D (প্রতিটি $1/2$); D-র একমাত্র প্রতিবেশী C ($1$): $$ P=\begin{pmatrix} 0 & 1 & 0 & 0\ \tfrac12 & 0 & \tfrac12 & 0\ 0 & \tfrac12 & 0 & \tfrac12\ 0 & 0 & 1 & 0 \end{pmatrix}. $$

(খ) stationary distribution। undirected graph-এর simple random walk-এ পরিচিত ফল: $\pi_i\propto \deg(i)$। degrees: $\deg(A)=1,\ \deg(B)=2,\ \deg(C)=2,\ \deg(D)=1$; মোট degree $=6$ ($=2\times$প্রান্ত-সংখ্যা, এখানে ৩টি প্রান্ত)। তাই $$ \boxed{\pi=\Big(\tfrac16,\ \tfrac26,\ \tfrac26,\ \tfrac16\Big)=\Big(\tfrac16,\tfrac13,\tfrac13,\tfrac16\Big)}. $$

detailed balance যাচাই (একটা প্রতিনিধি জোড়া, A–B): $$ \pi_A P_{AB}=\tfrac16\cdot 1=\tfrac16,\qquad \pi_B P_{BA}=\tfrac26\cdot\tfrac12=\tfrac16. $$ সমান ✓। (B–C জোড়া: $\pi_B P_{BC}=\tfrac26\cdot\tfrac12=\tfrac16$ আর $\pi_C P_{CB}=\tfrac26\cdot\tfrac12=\tfrac16$ ✓।) আসলে যেকোনো প্রান্ত $\{i,j\}$-এ $\pi_i P_{ij}=\tfrac{\deg(i)}{2E}\cdot\tfrac{1}{\deg(i)}=\tfrac1{2E}=\pi_j P_{ji}$ — তাই undirected graph-এর random walk সর্বদা reversible, আর $\pi_i\propto\deg(i)$ সবসময় stationary।

সমাধান ৮ (★★)¶

শুরু $\mu_0=[0,1]$ (Rainy)। $$ \pi_1=\mu_0 P=[0,1]\begin{pmatrix}0.8&0.2\0.4&0.6\end{pmatrix}=[0.4,\ 0.6]. $$ $$ \pi_2=\pi_1 P=[0.4,0.6]\begin{pmatrix}0.8&0.2\0.4&0.6\end{pmatrix} =[\,0.4\cdot0.8+0.6\cdot0.4,\ \ 0.4\cdot0.2+0.6\cdot0.6\,]=[\,0.32+0.24,\ 0.08+0.36\,]=[0.56,\ 0.44]. $$ তাই $\pi_2(\text{Sunny})=\boxed{0.56}$। Figure 2-র নীল curve ($[0,1]$ থেকে শুরু) $n=2$-এ ঠিক $0.56$-তে — মিলে গেল। (লক্ষ করুন এটা ইতিমধ্যে stationary $0.667$-এর দিকে উঠছে: $0\to0.4\to0.56\to\cdots\to0.667$।)

গ · প্রমাণভিত্তিক (proof-based)¶

সমাধান ৯ (★★)¶

দাবি: detailed balance $\pi_i P_{ij}=\pi_j P_{ji}$ (সব $i,j$) ⟹ $\pi=\pi P$।

প্রমাণ। যেকোনো নির্দিষ্ট state $j$-এর জন্য $(\pi P)_j$ গণনা করি: $$ (\pi P)j=\sum_i \pi_i P=\pi_j\cdot 1=\pi_j. $$ ধাপগুলো: (১) ম্যাট্রিক্স-গুণফলের সংজ্ঞা; (২) প্রতিটি পদে detailed balance বসিয়ে }\ \overset{\text{(DB)}}{=}\ \sum_i \pi_j P_{ji}=\pi_j\sum_i P_{ji$\pi_i P_{ij}$-কে $\pi_j P_{ji}$ দিয়ে বদলানো; (৩) $\pi_j$ যোগফলের বাইরে নেওয়া (এটা $i$-নিরপেক্ষ); (৪) $\sum_i P_{ji}=1$ কারণ $P$-র প্রতিটি সারি একটা probability distribution (state $j$ থেকে কোথাও-না-কোথাও যেতেই হবে)। যেহেতু এটা প্রতিটি $j$-র জন্য সত্য, $\pi P=\pi$। ∎

(সারমর্ম: detailed balance per-pair ভারসাম্য দেয়; প্রতিটি $j$-এর চারপাশে সব $i$-এর উপর যোগ করলে per-state global ভারসাম্য $=$ stationarity বেরিয়ে আসে।)

সমাধান ১০ (★★)¶

দাবি (Chapman–Kolmogorov): $m,n\ge0$-এর জন্য $(P^{m+n})_{ij}=\sum_k (P^m)_{ik}(P^n)_{kj}$।

প্রমাণ ১ (matrix বীজগণিত): সংজ্ঞা অনুযায়ী $n$-ধাপ transition matrix হলো ম্যাট্রিক্স-ঘাত $P^n$ (আবেশে: $P^1=P$, আর $(P^{n})_{ij}=\sum_k (P^{n-1})_{ik}P_{kj}$ — total probability দিয়ে)। তাই $P^{m+n}=P^m P^n$, আর ম্যাট্রিক্স-গুণফলের $(i,j)$-ভুক্তি ঠিক $\sum_k (P^m)_{ik}(P^n)_{kj}$। ∎

প্রমাণ ২ (probabilistic, যা স্বজ্ঞা দেয়): total probability-র নিয়মে, $(m+n)$ ধাপ পরে $i$ থেকে $j$-তে পৌঁছানোর ঘটনাকে মধ্যবর্তী state $X_m=k$ অনুযায়ী ভাগ করি: $$ (P^{m+n}){ij}=P(X=j,\ X_m=k\mid X_0=i). $$ শর্ত-ভাঙি ও }=j\mid X_0=i)=\sum_k P(X_{m+nMarkov property খাটাই (একবার $X_m=k$ জানলে আগের অংশ ভুলে যাওয়া যায়): $$ =\sum_k \underbrace{P(X_m=k\mid X_0=i)}{(P^m)}}\ \underbrace{P(X_{m+n}=j\mid X_m=k){(P^n)=\sum_k (P^m)}{ik}(P^n). \qquad\blacksquare $$

স্বজ্ঞা: "$i$ থেকে $j$-তে $m+n$ ধাপে যাওয়া" $=$ "প্রথম $m$ ধাপে কোনো এক মধ্যবর্তী থামার-জায়গা $k$-তে যাও, তারপর সেখান থেকে বাকি $n$ ধাপে $j$-তে যাও" — সব সম্ভাব্য মাঝপথ $k$-র উপর যোগফল। দীর্ঘ যাত্রাকে দুই খণ্ডে ভেঙে মাঝবিন্দুর উপর যোগ — এটাই Chapman–Kolmogorov-এর হৃদয়।

সমাধান ১১ (★★★)¶

দাবি: target (un-normalized) $f$, প্রতিসম proposal $q(x'\mid x)=q(x\mid x')$, accept-probability $\alpha(x,x')=\min\!\big(1,\,f(x')/f(x)\big)$ হলে, off-diagonal transition kernel $K(x\to x')=q(x'\mid x)\,\alpha(x,x')$ ($x'\ne x$) target $f$-এর সাপেক্ষে detailed balance মানে: $$ f(x)\,K(x\to x')=f(x')\,K(x'\to x). $$

প্রমাণ। $x=x'$ হলে দুই পাশ অভিন্ন, তাই ধরি $x\ne x'$। সাধারণতা না হারিয়ে ধরি $f(x')\ge f(x)$ (বিপরীত ক্ষেত্রটি $x\leftrightarrow x'$ অদলবদলে একই যুক্তি, কারণ বিবৃতি প্রতিসম)। তখন: $$ \alpha(x,x')=\min!\Big(1,\tfrac{f(x')}{f(x)}\Big)=1\quad(\text{যেহেতু } f(x')/f(x)\ge1), $$ $$ \alpha(x',x)=\min!\Big(1,\tfrac{f(x)}{f(x')}\Big)=\tfrac{f(x)}{f(x')}\quad(\text{যেহেতু } f(x)/f(x')\le1). $$ এবার দুই পাশ: $$ \text{বাঁ}=f(x)\,K(x\to x')=f(x)\,q(x'\mid x)\,\alpha(x,x')=f(x)\,q(x'\mid x)\cdot 1=f(x)\,q(x'\mid x), $$ $$ \text{ডান}=f(x')\,K(x'\to x)=f(x')\,q(x\mid x')\,\alpha(x',x)=f(x')\,q(x\mid x')\cdot\frac{f(x)}{f(x')}=f(x)\,q(x\mid x'). $$ proposal প্রতিসম বলে $q(x'\mid x)=q(x\mid x')$, তাই বাঁ $=$ ডান। ∎

সিদ্ধান্ত: $K$ যেহেতু $f$ (অতএব normalize-করা $\pi=f/Z$)-এর সাপেক্ষে detailed balance মানে, §Q9 অনুযায়ী $\pi$ এই chain-এর stationary distribution। তাই Metropolis chain যথেষ্ট সময় চালালে target $\pi$ থেকে নমুনা দেয়।

দুটি মূল লক্ষণীয় বিষয়: (১) প্রমাণজুড়ে $f$ ঢুকেছে কেবল অনুপাত $f(x')/f(x)$ আকারে — normalizing constant $Z$ কোথাও লাগেনি (§Q4গ-র সাথে মিল)। (২) এখানে proposal প্রতিসম ধরা হয়েছে; সাধারণ Metropolis–Hastings-এ asymmetric $q$-এর জন্য accept-অনুপাতে একটা সংশোধন-গুণক যোগ হয়: $\alpha=\min\!\big(1,\ \tfrac{f(x')\,q(x\mid x')}{f(x)\,q(x'\mid x)}\big)$ — যা ঠিক সেই asymmetry-কে কাটায় যাতে detailed balance বজায় থাকে।

ঘ · কোডিং (coding)¶

সমাধান ১২ (★)¶

import numpy as np
import matplotlib.pyplot as plt

P = np.array([[0.8, 0.2],
              [0.4, 0.6]])          # ক্রম: Sunny, Rainy
pi = np.array([0.0, 1.0])           # [0,1] = Rainy থেকে শুরু

probs = []
for n in range(13):                  # n = 0..12
    probs.append(pi[0])              # P(Sunny) ট্র্যাক করি
    pi = pi @ P                      # এক ধাপ এগোও:  pi_{n+1} = pi_n P

plt.figure(figsize=(8, 4.5))
plt.plot(range(13), probs, "o-", lw=2, label=r"start $[0,1]$ (Rainy)")
plt.axhline(2/3, color="green", ls="--", lw=2, label=r"stationary $\pi_S = 2/3$")
plt.xlabel("step  n"); plt.ylabel(r"$\mu_n(\mathrm{Sunny})$")
plt.ylim(-0.03, 1.03); plt.legend(); plt.title("Convergence to stationary (E1)")
plt.tight_layout(); plt.show()

ফলাফল: ছাপলে দেখা যায় $\mu_n(\text{Sunny})=0,\,0.4,\,0.56,\,0.624,\,0.6496,\dots\to 0.667$ — curve দ্রুত (৫–৬ ধাপে) সবুজ রেখায় পৌঁছায়, ঠিক §Q8-এর হাতে-কষা $\pi_1=0.4,\pi_2=0.56$-এর সাথে মিলে। মূল শিক্ষা: distribution-এর evolution মানেই বারবার $P$ দিয়ে গুণ (pi @ P)।

সমাধান ১৩ (★★)¶

import numpy as np

P = np.array([[0.8, 0.2],
              [0.4, 0.6]])

# (ক) উচ্চ-ঘাত পদ্ধতি: P^50 -এর প্রতিটি সারি একই pi-তে গড়ায়
Pn = np.linalg.matrix_power(P, 50)
print("P^50 =\n", np.round(Pn, 6))
# -> দুই সারিই [0.666667, 0.333333]; শুরু যাই হোক, ৫০ ধাপ পরে distribution একই

# (খ) eigen পদ্ধতি: P^T -এর eigenvalue 1 -এর (বাঁ) eigenvector
w, v = np.linalg.eig(P.T)
idx = np.argmin(np.abs(w - 1.0))      # eigenvalue == 1
pi = np.real(v[:, idx])
pi = pi / pi.sum()                    # normalize -> probability
print("eigen pi =", np.round(pi, 6)) # -> [0.666667, 0.333333]

ফলাফল ও ব্যাখ্যা: দুই পদ্ধতিই $\pi=[0.6667,0.3333]$ দেয় (= হাতে-কষা $2/3,1/3$, §Q6)। কেন কাজ করে: (ক) irreducible+aperiodic chain-এ $P^n$-এর প্রতিটি সারি $\to\pi$ (যেকোনো শুরু একই গন্তব্যে — Figure 2)। (খ) $\pi=\pi P$ মানে $P^{\mathsf T}\pi^{\mathsf T}=\pi^{\mathsf T}$, অর্থাৎ $\pi^{\mathsf T}$ হলো $P^{\mathsf T}$-এর eigenvalue-$1$ eigenvector; stochastic matrix-এর সর্বদা একটি eigenvalue ঠিক $1$ থাকে (Perron–Frobenius), তাই এই eigenvector normalize করলেই stationary distribution।

সমাধান ১৪ (★★★)¶

import numpy as np
import matplotlib.pyplot as plt

# ----- target (un-normalized): 0.6 N(-2, 0.7^2) + 0.4 N(2.2, 1^2) -----
def f(x):
    return (0.6*np.exp(-0.5*((x+2.0)/0.7)**2)/0.7 +
            0.4*np.exp(-0.5*((x-2.2)/1.0)**2)/1.0)

def metropolis(n, step, x0=0.0, seed=0):
    rng = np.random.default_rng(seed)
    x, fx = x0, f(x0)
    chain = np.empty(n)
    for i in range(n):
        xp = x + rng.normal(0.0, step)          # random-walk proposal (symmetric)
        fp = f(xp)
        if rng.random() < fp / fx:               # accept w.p. min(1, f(x')/f(x))
            x, fx = xp, fp
        chain[i] = x
    return chain

# (ক) ভালো step -> 60k iteration, প্রথম 2000 burn-in ফেলি
chain = metropolis(60000, step=1.6, x0=0.0, seed=7)
samples = chain[2000:]

# (খ) trace + histogram, histogram-এ normalize-করা target overlay
xs = np.linspace(-6, 7, 800)
pdf = f(xs) / np.trapezoid(f(xs), xs)            # grid-এ normalize (Z আনুমান)

fig, (a1, a2) = plt.subplots(1, 2, figsize=(12, 4.2))
a1.plot(chain[:1500], lw=0.7, color="#1b6ca8"); a1.set_title("trace (step=1.6)")
a1.set_xlabel("iteration"); a1.set_ylabel("x")
a2.hist(samples, bins=70, density=True, alpha=0.45, color="#1b6ca8",
        label="MCMC samples")
a2.plot(xs, pdf, "r", lw=2.5, label="true target")
a2.legend(); a2.set_title("histogram vs target"); a2.set_xlabel("x")
plt.tight_layout(); plt.show()

# (গ) poor mixing: খুব ছোট ও খুব বড় step
for s in (0.1, 20.0):
    c = metropolis(3000, step=s, x0=0.0, seed=1)
    # accept-rate ≈ পরপর-ভিন্ন মানের ভগ্নাংশ
    moved = np.mean(c[1:] != c[:-1])
    print(f"step={s:>4}:  accept-rate ≈ {moved:.2f}")

ফলাফল ও ব্যাখ্যা: - (ক)–(খ) trace দুই mode ($-2$ ও $2.2$) ঘিরে ওঠানামা করে, আর histogram লাল target curve-এর গায়ে নিখুঁতভাবে বসে (Figure 3, 4 পুনরুৎপাদন) — যদিও আমরা $f$ থেকে সরাসরি নমুনা তুলিনি, শুধু অনুপাত $f(x')/f(x)$ ব্যবহার করেছি। - (গ) ছাপা accept-rate: step=0.1 → accept-rate প্রায় $0.97$ (প্রায় সব accept, কিন্তু ক্ষুদ্র পদক্ষেপে chain এক mode-এই হামাগুড়ি দেয় — $3000$ iteration-এ অন্য mode-এ পৌঁছায়ই না, trace প্রায় সমতল-স্থানীয়); step=20 → accept-rate প্রায় $0.06$ (বেশির ভাগ proposal দূরে কম-ঘনত্বে পড়ে reject হয়, trace লম্বা সমতল ধাপে আটকে)। দুটোই poor mixing-এর দুই বিপরীত চেহারা — মাঝারি step ($\sim1.6$) দুদিকের ভারসাম্য রাখে। (সাধারণ বাস্তব-সতর্কতা: chain যথেষ্ট-দীর্ঘ চালানো ও autocorrelation/effective sample size দেখা জরুরি।)

শেষ কথা। §৭-এর প্রতিটি ফল একে অন্যকে ধরে রাখে: Chapman–Kolmogorov (Q5, Q10) multi-step probability দেয়; stationary সমীকরণ (Q6, Q8, Q13) long-run distribution দেয়; detailed balance (Q3, Q7, Q9) stationary-র সহজ যথেষ্ট-শর্ত দেয়; আর Metropolis (Q4, Q11, Q14) ঠিক সেই detailed balance ব্যবহার করে যেকোনো target থেকে নমুনা তোলে — Markov-তত্ত্ব থেকে MCMC-র সেতু এখানেই সম্পূর্ণ।