Distribution আর sampling — normal, binomial আর simulation

পর্ব ১৬-এ np.random.default_rng দিয়ে random সংখ্যা বানিয়েছেন, আর একটি ছোট Monte Carlo দিয়ে দেরিতে ডেলিভারির সম্ভাবনা অনুমান করেছেন। পর্ব ৩৫-এ (AI-এর জন্য সম্ভাবনা) event, conditional probability আর independence-এর নিয়ম আসছে। এই পর্বে সেই দুটি জোড়া লাগবে: একটি অনিশ্চিত সংখ্যার পুরো চেহারা কীভাবে লেখা যায় (distribution), তা থেকে কীভাবে নমুনা তোলা যায় (sampling), আর হাতে হিসাব কঠিন হলে কীভাবে অনেকবার ঘটিয়ে উত্তর পাওয়া যায় (simulation). উদাহরণ বাজারঘরের, এই কোর্সের কাল্পনিক অনলাইন দোকান। পাতার কোড ব্রাউজারেই চলে; প্রথমবার SciPy আর pandas নামাতে কয়েক সেকেন্ড লাগতে পারে।

১ · Random variable — অনিশ্চিত সংখ্যার একটি নাম

বাজারঘরের ম্যানেজার সকাল ১০টায় জানতে চান: “আগামী এক ঘণ্টায় কয়টি অর্ডার আসবে?” উত্তরটি আগে থেকে জানা নেই। হতে পারে ০, ২ বা ৫টি। কিন্তু একেবারে অজানাও নয়: অতীতের খাতা দেখে বলা যায় ২টি আসা বেশ সম্ভব, ১০টি প্রায় অসম্ভব। এমন একটি অনিশ্চিত সংখ্যাকে নাম দিলে, ধরুন X, তাকে বলে random variable: একটি নিয়ম, যা প্রতিটি সম্ভাব্য ফলাফলকে একটি সংখ্যা দেয়, আর সেই সংখ্যাগুলোর প্রত্যেকটির একটি সম্ভাবনা আছে। “Random” মানে এখানে “এলোমেলো” নয়, “আগে থেকে নিশ্চিত নয়”; আর “variable” মানে প্রোগ্রামের variable নয়, গণিতের একটি পরিমাণ।

Random variable দুই রকম, আর দুটির হিসাবের ভাষা আলাদা:

একটি সূক্ষ্মতা: দাম টাকায় আর পয়সায় মাপা হয়, তাই কাগজে-কলমে তা discrete (বাজারঘর দাম রাখে পূর্ণসংখ্যার পয়সায়)। কিন্তু মান এত বেশি আর এত ঘন যে বিশ্লেষণে তাকে continuous ধরা সুবিধার। লেখার রীতি: বড় হাতের X মানে random variable নিজে, ছোট হাতের x মানে তার একটি নির্দিষ্ট মান। P(X = 2) পড়ুন “X-এর মান ২ হওয়ার সম্ভাবনা”।

২ · PMF, PDF আর CDF — সম্ভাবনা কোথায় থাকে

তিনটি সংক্ষিপ্ত নাম এই পর্ব জুড়ে বারবার আসবে। প্রথমে সহজ ভাষায়, তারপর ছবিতে।

• PMF (probability mass function): discrete X-এর জন্য, প্রতিটি মানের সম্ভাবনা। P(X = 2) = 0.35 মানে “ঠিক ২টি অর্ডার আসার সম্ভাবনা ০.৩৫”। সব মানের PMF যোগ করলে 1.
• PDF (probability density function): continuous X-এর জন্য, একটি বক্ররেখা। রেখার উচ্চতা সম্ভাবনা নয়, ঘনত্ব (density): কোন এলাকায় মান বেশি ঘন। সম্ভাবনা হলো রেখার নিচের ক্ষেত্রফল; পুরো রেখার নিচে মোট ক্ষেত্রফল 1.
• CDF (cumulative distribution function): দুই ধরনের জন্যই, F(x) = P(X ≤ x), মানে “x বা তার কম” হওয়ার সম্ভাবনা। এটি শূন্য থেকে ধীরে ধীরে বেড়ে 1-এ পৌঁছায়, কখনো কমে না। পর্ব ৩০-এ ডেটা থেকে আঁকা ECDF এরই নমুনা-সংস্করণ।

PDF-এর উচ্চতা 1-এর বেশিও হতে পারে, কারণ তা সম্ভাবনা নয়: σ = 0.1-এর normal-এর চূড়া প্রায় 3.99, তবু নিচের মোট ক্ষেত্রফল 1. এবার গণনা থেকে একটি PMF আর CDF নিজের হাতে বানান। নিচে বাজারঘরের ঢাকা শাখার ২০ ঘণ্টার খাতা: প্রতিটি সংখ্যা এক ঘণ্টার অর্ডার।

৩ · Expected value আর variance — গড়ে কত, ছড়ানো কত

পুরো distribution একটি টেবিল বা ছবি। প্রায়ই দুটি সংখ্যায় তার সারাংশ লাগে: কেন্দ্র কোথায়, আর মানগুলো কতটা ছড়ানো। Expected value (প্রত্যাশিত মান, E[X] বা μ, গ্রিক অক্ষর “মিউ”) হলো প্রতিটি মানকে তার সম্ভাবনা দিয়ে গুণ করে যোগফল:

E[X] = ∑ x · P(X = x)

এখানে ∑ (গ্রিক “সিগমা”, বড় হাতের) মানে “সব মানের জন্য যোগ করুন”। এটি ওজনসহ গড়: বেশি সম্ভব মান বেশি ওজন পায়। Variance (Var(X) বা σ²) হলো গড় থেকে দূরত্বের বর্গের প্রত্যাশিত মান, আর তার বর্গমূল standard deviation σ:

Var(X) = ∑ (x − μ)² · P(X = x) = E[X²] − μ²

হাতে হিসাব — অংশ ২-এর PMF দিয়ে

তাহলে ঘণ্টায় গড়ে 1.65টি অর্ডার, আর σ = √1.1275 ≈ 1.0618. “1.65টি অর্ডার” কোনো ঘণ্টায় আসে না; expected value একটি দীর্ঘমেয়াদি গড়, একটি সম্ভাব্য মান নয়। শর্টকাট সূত্রেও মেলান: E[X²] = 0 + 0.30 + 1.40 + 1.35 + 0.80 = 3.85, তাই Var = 3.85 − 1.65² = 3.85 − 2.7225 = 1.1275. এবার একই হিসাব কোডে, আর তারপর simulation দিয়ে: এই PMF থেকে অনেকগুলো কাল্পনিক ঘণ্টা তুলে তাদের গড় আর variance দেখা। নিচেরটি লেখকের যন্ত্রে সত্যিকারের run (Python 3.14.4, NumPy 2.5.1, SciPy 1.18.0); পাতার চারটি চালানো যায় এমন ব্লক পরের অংশগুলোর জন্য রাখা, তাই এটি কাগজে দেওয়া। কোডটি যেকোনো ব্লকে বসিয়ে মেলাতে পারেন।

import numpy as np
from scipy import stats

x = np.array([0, 1, 2, 3, 4])                  # ঘণ্টায় অর্ডার
p = np.array([0.15, 0.30, 0.35, 0.15, 0.05])   # PMF
E = np.sum(x * p)
var = np.sum((x - E) ** 2 * p)
print("E[X] =", round(E, 4), "| Var(X) =", round(var, 4), "| sd =", round(var ** 0.5, 4))

X = stats.rv_discrete(values=(x, p))           # একই PMF, scipy object
print("scipy:", X.mean().round(4), X.var().round(4))

rng = np.random.default_rng(42)
for n in (100, 10_000, 1_000_000):
    s = rng.choice(x, size=n, p=p)
    print(f"n={n:>9,}  mean={s.mean():.4f}  var={s.var():.4f}")

E[X] = 1.65 | Var(X) = 1.1275 | sd = 1.0618
scipy: 1.65 1.1275
n=      100  mean=1.5700  var=0.9851
n=   10,000  mean=1.6423  var=1.1250
n=1,000,000  mean=1.6503  var=1.1272

দুটি নিয়ম পরে বারবার কাজে লাগবে। এক, E[aX + b] = a·E[X] + b আর Var(aX + b) = a²·Var(X): প্রতিটি অর্ডারে ৬০ টাকা ডেলিভারি চার্জ হলে ঘণ্টায় চার্জের গড় 60 × 1.65 = 99 টাকা, আর variance 3600 × 1.1275. দুই, যোগফলের গড় সবসময় গড়ের যোগফল, কিন্তু variance যোগ হয় শুধু স্বাধীন (independent, পর্ব ৩৫) random variable-এর ক্ষেত্রে।

৪ · Bernoulli, binomial আর geometric — হ্যাঁ/না-এর গণনা

অংশ ২-এর PMF ছিল ডেটা থেকে মাপা। অনেক সময় পরিস্থিতির গঠন থেকেই PMF-এর সূত্র লেখা যায়, কয়েকটি parameter (মডেলের ঠিক করে দেওয়া সংখ্যা) দিয়ে। এমন পরিচিত ছাঁচগুলোকে নামে চেনা হয়। সবচেয়ে সরলটি দিয়ে শুরু।

Bernoulli: একটি পরীক্ষা, দুটি ফল: হ্যাঁ (1) বা না (0). বাজারঘরের একটি অর্ডার ফেরত আসবে কি না, সম্ভাবনা p. P(X = 1) = p, P(X = 0) = 1 − p; E[X] = p, Var(X) = p(1 − p). p = 0.1 হলে variance 0.09; p = 0.5-এ variance সবচেয়ে বেশি (0.25), কারণ তখনই ফল সবচেয়ে অনিশ্চিত।

Binomial: একই রকম n-টি স্বাধীন Bernoulli পরীক্ষায় মোট কয়টি “হ্যাঁ”। একটি দিনে ১০টি অর্ডার, প্রতিটি আলাদাভাবে ০.১ সম্ভাবনায় ফেরত আসে; X = দিনে ফেরত আসা অর্ডার। সূত্র:

P(X = k) = C(n, k) · p^k · (1 − p)^{n − k}

এখানে C(n, k) (“n choose k”) হলো n-টি থেকে k-টি বাছার উপায়ের সংখ্যা, Python-এ math.comb(n, k). যুক্তিটি সরল: নির্দিষ্ট k-টি ফেরত আর বাকি n − k-টি না ফেরতের একটি নির্দিষ্ট ক্রমের সম্ভাবনা p^k(1 − p)^{n − k}, আর এমন ক্রম আছে C(n, k)-টি। হাতে: P(X = 0) = 0.9¹⁰ ≈ 0.3487; P(X = 1) = 10 × 0.1 × 0.9⁹ ≈ 0.3874; P(X = 2) = 45 × 0.01 × 0.9⁸ ≈ 0.1937. E[X] = np = 1, Var(X) = np(1 − p) = 0.9.

Geometric (সংক্ষেপে): প্রথম “হ্যাঁ” আসতে কয়টি পরীক্ষা লাগে। প্রথম ফেরত আসা অর্ডারটি k-তম হওয়ার সম্ভাবনা (1 − p)^{k − 1}·p, গড় 1/p: p = 0.1 হলে গড়ে দশম অর্ডারে প্রথম ফেরত। scipy.stats.geom এই সংজ্ঞাই নেয় (k শুরু 1 থেকে); কিছু বই “প্রথম হ্যাঁ-এর আগে কয়টি না” গোনে (শুরু 0 থেকে), তাই সূত্র মেলানোর আগে সংজ্ঞা দেখে নিন।

scipy.stats-এর এক ছাঁচ

scipy.stats-এর প্রতিটি distribution একই method দেয়। Parameter দিয়ে একবার বানালে তাকে বলে “frozen” distribution, যেমন X = stats.binom(10, 0.1); তারপর:

ppf (percent point function) নামটি অপরিচিত লাগলেও কাজটি পরিচিত: percentile. random_state-এ একটি পূর্ণসংখ্যা বা np.random.default_rng(42)-এর মতো Generator দেওয়া যায় (পর্ব ১৬); এই কোর্সে সবসময় Generator. নিচের কোডে binomial-এর PMF সূত্র দিয়ে আর scipy দিয়ে বের করে মেলানো হয়েছে, তারপর ১০,০০০টি কাল্পনিক দিন simulate করে তাদের আসল কম্পাঙ্ক (frequency) একই ছবিতে।

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

n, p = 10, 0.1               # দিনে ১০টি অর্ডার, প্রতিটি ফেরতের সম্ভাবনা 0.1
k = np.arange(0, 6)
by_hand = np.array([math.comb(n, i) * p**i * (1 - p)**(n - i) for i in k])
X = stats.binom(n, p)        # frozen distribution object
print("k        :", k)
print("formula  :", by_hand.round(4))
print("scipy pmf:", X.pmf(k).round(4))
print("P(X >= 3):", round(1 - X.cdf(2), 4), "| sf(2):", round(X.sf(2), 4))
print("mean, var:", X.mean(), round(X.var(), 4), "| np, np(1-p):", n * p, round(n * p * (1 - p), 4))
print("ppf(0.95):", X.ppf(0.95))

rng = np.random.default_rng(42)
days = X.rvs(size=10_000, random_state=rng)      # ১০,০০০টি কাল্পনিক দিন
freq = np.bincount(days, minlength=n + 1)[:6] / days.size
print("simulated:", freq.round(4))
print("sim mean, var:", days.mean().round(4), days.var().round(4))

fig, ax = plt.subplots(figsize=(6, 3.2), layout="constrained")
ax.bar(k - 0.2, X.pmf(k), width=0.4, label="PMF (exact)")
ax.bar(k + 0.2, freq, width=0.4, label="simulated, 10,000 days")
ax.set_xlabel("returned orders out of 10")
ax.set_ylabel("probability")
ax.set_title("Binomial(n=10, p=0.1)")
ax.legend()
plt.show()

k        : [0 1 2 3 4 5]
formula  : [0.3487 0.3874 0.1937 0.0574 0.0112 0.0015]
scipy pmf: [0.3487 0.3874 0.1937 0.0574 0.0112 0.0015]
P(X >= 3): 0.0702 | sf(2): 0.0702
mean, var: 1.0 0.9 | np, np(1-p): 1.0 0.9
ppf(0.95): 3.0
simulated: [0.353  0.3884 0.1896 0.0557 0.0112 0.0019]
sim mean, var: 0.9902 0.9015

৫ · Poisson আর exponential — ঘণ্টায় কয়টি অর্ডার, পরেরটি কখন

Binomial-এ n জানা ছিল: দিনে ১০টি অর্ডার। কিন্তু “আগামী এক ঘণ্টায় কয়টি অর্ডার আসবে” প্রশ্নে n কী? বাজারঘরের লাখো সম্ভাব্য ক্রেতা, প্রত্যেকের এই ঘণ্টায় অর্ডার দেওয়ার সম্ভাবনা খুব ছোট। জানা শুধু গড় হার: ঘণ্টায় গড়ে ৩টি। এমন “অনেক সুযোগ, প্রতিটি বিরল, জানা গড় হার” পরিস্থিতির মডেল Poisson distribution, যার একটিমাত্র parameter λ (গ্রিক “ল্যাম্বডা”): একটি নির্দিষ্ট সময় বা জায়গায় ঘটনার গড় সংখ্যা।

P(X = k) = e^−λ · λ^k / k!

e ≈ 2.718 (স্বাভাবিক logarithm-এর ভিত্তি), আর k! = 1 × 2 × … × k (০! = ১)। λ = 3 হলে: P(X = 0) = e⁻³ ≈ 0.0498, মানে প্রায় ৫% ঘণ্টায় একটিও অর্ডার আসে না; P(X = 3) = e⁻³ × 27 / 6 ≈ 0.2240. Poisson-এর একটি বিশেষ বৈশিষ্ট্য: গড় আর variance দুটিই λ। বাস্তব ডেটায় variance গড়ের চেয়ে অনেক বেশি হলে (ঈদের আগে, বিজ্ঞাপনের পরে হঠাৎ ভিড়) Poisson ঠিক মডেল নয়; একে overdispersion বলে, আর তখন negative binomial-এর মতো মডেল লাগে।

Binomial-এর সীমা হিসেবে Poisson

Poisson আসলে binomial থেকেই আসে। ধরুন এক ঘণ্টাকে n-টি ছোট টুকরায় ভাগ করলেন, প্রতিটি টুকরায় একটি অর্ডার আসার সম্ভাবনা p = λ/n. n বড় আর p ছোট হতে থাকলে, কিন্তু গুণফল np = λ স্থির রাখলে, binomial(n, λ/n)-এর PMF Poisson(λ)-এর দিকে যায়। ব্যবহারিক নিয়ম: n বড় আর p ছোট হলে binomial-এর জায়গায় Poisson দিয়ে হিসাব করা যায়, যা সহজ। অংশ ৬-এর কোডে n = 10 থেকে 1,00,000 পর্যন্ত দুটির PMF-এর সবচেয়ে বড় পার্থক্য সত্যিই মাপা হয়েছে।

Exponential (সংক্ষেপে) — পরের অর্ডার পর্যন্ত অপেক্ষা

অর্ডারগুলো যদি Poisson-এর নিয়মে আসে (ঘণ্টায় গড়ে λ = 3টি), তাহলে দুটি পরপর অর্ডারের মাঝের সময় T একটি continuous random variable, যার distribution exponential: গড় অপেক্ষা 1/λ = ⅓ ঘণ্টা = ২০ মিনিট, আর P(T > t) = e^−λt. আধা ঘণ্টায় একটিও অর্ডার না আসার সম্ভাবনা e^−1.5 ≈ 0.2231. একই প্রশ্ন Poisson দিয়েও করা যায় (আধা ঘণ্টায় গড় 1.5টি, শূন্যটি আসার সম্ভাবনা e^−1.5): দুটি মডেল একই ঘটনার দুই দিক।

Exponential-এর একটি অদ্ভুত বৈশিষ্ট্য memorylessness (স্মৃতিহীনতা): ১০ মিনিট ধরে কোনো অর্ডার আসেনি, তাতে পরের ৩০ মিনিটে না আসার সম্ভাবনা বদলায় না; P(T > 40 মিনিট | T > 10 মিনিট) = P(T > 30 মিনিট)। সূত্রে: e^−λ·40 / e^−λ·10 = e^−λ·30. “এতক্ষণ আসেনি, এবার নিশ্চয়ই আসবে” ভাবনাটি এই মডেলে ভুল। পুরোনো যন্ত্রের নষ্ট হওয়ার মতো অপেক্ষা স্মৃতিহীন নয়; সেখানে exponential ভুল মডেল।

৬ · Uniform আর normal — μ, σ, ৬৮–৯৫–৯৯.৭ আর z-score

Uniform: একটি পরিসরের সব মান সমান সম্ভব। Discrete uniform-এর উদাহরণ ছক্কা; continuous uniform-এর উদাহরণ rng.random(), যা [0, 1)-এর যেকোনো জায়গায় সমান ঘনত্বে পড়ে। [a, b]-এ গড় (a + b)/2, variance (b − a)²/12. scipy-তে একটি ফাঁদ: stats.uniform(loc, scale)-এর পরিসর [loc, loc + scale], [loc, scale] নয়; stats.uniform(2, 4) মানে [2, 6]. অংশ ৯-এর π-এর অনুমান পুরোটাই uniform বিন্দুর ওপর দাঁড়িয়ে।

Normal (Gaussian বা “ঘণ্টার মতো” বক্ররেখা) দুটি parameter দিয়ে পুরো ঠিক হয়: μ (কেন্দ্র, যেখানে চূড়া) আর σ (ছড়ানো, চূড়া থেকে বাঁকের জায়গা পর্যন্ত দূরত্ব)। বক্ররেখা μ-এর দুই পাশে সমান (symmetric), আর দুই দিকে লেজ অনন্ত পর্যন্ত যায়, যদিও খুব দ্রুত সরু হয়। অনেক ছোট স্বাধীন প্রভাবের যোগফল প্রায়ই normal-এর মতো হয়; এর গাণিতিক রূপ Central Limit Theorem, পর্ব ৩৭-এর বিষয়।

Normal-এর সবচেয়ে কাজের নিয়ম ৬৮–৯৫–৯৯.৭ নিয়ম: মানগুলোর প্রায় ৬৮% থাকে μ ± 1σ-এর মধ্যে, প্রায় ৯৫% থাকে μ ± 2σ-এর মধ্যে, আর প্রায় ৯৯.৭% থাকে μ ± 3σ-এর মধ্যে। বাজারঘরের ঢাকা শাখায় দৈনিক অর্ডার যদি মোটামুটি normal হয়, μ = 500, σ = 40, তাহলে প্রায় ৯৫% দিনে অর্ডার 420 থেকে 580-এর মধ্যে।

z-score — যেকোনো normal-কে একটি মাপকাঠিতে আনা

৫৬০টি অর্ডারের দিন কতটা অস্বাভাবিক? উত্তর নির্ভর করে σ-এর ওপর। z-score বলে একটি মান গড় থেকে কত σ দূরে:

z = (x − μ) / σ

z = 1.5 মানে “গড়ের দেড় σ ওপরে”। X normal হলে z-ও normal, μ = 0 আর σ = 1, যাকে standard normal বলে। তাই যেকোনো normal-এর সম্ভাবনা একটি টেবিল বা একটি function (stats.norm.cdf) দিয়েই পাওয়া যায়। ধাপে ধাপে দেখুন।

এবার তিনটি দাবি একসঙ্গে যাচাই: Poisson সত্যিই binomial-এর সীমা কি না, ৬৮–৯৫–৯৯.৭ নিয়মের আসল সংখ্যা, আর exponential-এর স্মৃতিহীনতা।

import numpy as np
from scipy import stats

lam = 3                                   # ঘণ্টায় গড়ে ৩টি অর্ডার
k = np.arange(0, 9)
P = stats.poisson(lam)
print("Poisson(3) pmf:", P.pmf(k).round(4))
print("P(X=0):", round(P.pmf(0), 4), "| e^-3:", round(np.exp(-3), 4), "| mean, var:", P.mean(), P.var())

for n in (10, 100, 1000, 100_000):        # n*p = 3 স্থির রেখে n বাড়ানো
    B = stats.binom(n, lam / n)
    gap = np.abs(B.pmf(k) - P.pmf(k)).max()
    print(f"n={n:>7,}  p={lam / n:.5f}  max|binom - poisson| = {gap:.6f}")

Z = stats.norm(0, 1)
for z in (1, 2, 3):
    print(f"P(|Z| < {z}) = {Z.cdf(z) - Z.cdf(-z):.4f}")
print("one-sided 95%:", round(Z.ppf(0.95), 4), "| two-sided 95%:", round(Z.ppf(0.975), 4))

rng = np.random.default_rng(42)
x = stats.norm(500, 40).rvs(size=100_000, random_state=rng)   # দৈনিক অর্ডার
print("simulated within 1, 2, 3 sd:",
      [round(float(np.mean(np.abs(x - 500) < z * 40)), 4) for z in (1, 2, 3)])

T = stats.expon(scale=1 / lam)            # অর্ডারের মাঝের সময়, ঘণ্টায়
print("P(T > 30 min):", round(T.sf(0.5), 4), "| mean wait (min):", round(T.mean() * 60, 1))
print("P(T > 40 | T > 10):", round(T.sf(40 / 60) / T.sf(10 / 60), 4))

Poisson(3) pmf: [0.0498 0.1494 0.224  0.224  0.168  0.1008 0.0504 0.0216 0.0081]
P(X=0): 0.0498 | e^-3: 0.0498 | mean, var: 3.0 3.0
n=     10  p=0.30000  max|binom - poisson| = 0.042786
n=    100  p=0.03000  max|binom - poisson| = 0.003432
n=  1,000  p=0.00300  max|binom - poisson| = 0.000337
n=100,000  p=0.00003  max|binom - poisson| = 0.000003
P(|Z| < 1) = 0.6827
P(|Z| < 2) = 0.9545
P(|Z| < 3) = 0.9973
one-sided 95%: 1.6449 | two-sided 95%: 1.96
simulated within 1, 2, 3 sd: [0.6817, 0.9536, 0.997]
P(T > 30 min): 0.2231 | mean wait (min): 20.0
P(T > 40 | T > 10): 0.2231

৭ · যখন ডেটা normal নয় — skew, log-normal আর ভারী লেজ

Normal এত পরিচিত যে অনেকে ধরে নেন সব ডেটাই normal. বাজারঘরের অর্ডারের মূল্য (order value) দেখলেই ভুলটি ধরা পড়ে। বেশির ভাগ অর্ডার ৩০০ থেকে ১,৫০০ টাকার, কিন্তু মাঝে মাঝে কেউ ফ্রিজ বা ল্যাপটপ কেনেন, ৬০,০০০ টাকার। মূল্য কখনো শূন্যের নিচে যায় না, অথচ ওপরের দিকে লম্বা লেজ। এমন ডেটাকে বলে ডান দিকে skewed (right-skewed). এর লক্ষণ: গড় median-এর চেয়ে অনেক বেশি, কারণ কয়েকটি বড় মান গড়কে টেনে নেয় (পর্ব ১৬-এর outlier-এর আলোচনা)।

Skewed ডেটায় normal-এর মডেল বসালে দুটি ক্ষতি হয়। এক, মডেল ঋণাত্মক মূল্য বা ঋণাত্মক সময়কে সম্ভব বলে। দুই, ওপরের লেজকে মডেল অনেক ছোট ধরে, তাই “৫০,০০০ টাকার বেশি অর্ডার প্রায় অসম্ভব” বলে, অথচ সেগুলো নিয়মিত আসে।

ভারী লেজ — যখন চরম মানই আসল খবর

কিছু distribution-এর লেজ normal তো বটেই, log-normal-এর চেয়েও মোটা: চরম মান অনেক বেশি ঘন ঘন আসে। এদের বলে heavy-tailed. এমন ডেটায় কয়েকটি বড় মান মোট যোগফলের বড় অংশ নিতে পারে, আর নমুনার গড় খুব ধীরে স্থির হয়; কিছু চরম ক্ষেত্রে (তাত্ত্বিক Cauchy distribution) গড়ই সংজ্ঞায়িত নয়, তাই law of large numbers কাজ করে না। কোন বাস্তব ডেটা কতটা ভারী লেজের, তা নিয়ে বিতর্ক আছে (ভাবনার প্রশ্ন ৫)। ব্যবহারিক শিক্ষাটি হলো: শুধু গড় আর σ দিয়ে ঝুঁকি মাপবেন না; লেজ আলাদা করে দেখুন, যেমন সর্বোচ্চ ১%-এর মান, আর পর্ব ৩০-এর histogram বা ECDF.

৮ · Sampling — population, নমুনা, stratified আর পক্ষপাত

এ পর্যন্ত distribution জানা ছিল, আর আমরা তা থেকে সংখ্যা তুলেছি। বাস্তবে উল্টো: distribution অজানা, হাতে শুধু কিছু ডেটা। Population হলো যাদের নিয়ে প্রশ্ন, তাদের সবাই: বাজারঘরের এ বছরের সব অর্ডার। নমুনা (sample) হলো তাদের যে অংশটি আসলে দেখা হয়: ২০০টি অর্ডারের ক্রেতাকে ফোন করে সন্তুষ্টির জরিপ। নমুনা থেকে population সম্পর্কে বলতে হলে নমুনাটি population-এর ছোট কিন্তু সৎ প্রতিচ্ছবি হতে হবে।

• Simple random sampling: population-এর প্রতিটি সদস্যের বাছাই হওয়ার সম্ভাবনা সমান, আর বাছাই random. pandas-এ orders.sample(n=200, random_state=rng). সহজ আর পক্ষপাতহীন, কিন্তু ছোট নমুনায় ছোট দলগুলো দৈবাৎ কম বা বেশি আসতে পারে।
• Stratified sampling: population-কে আগে স্তরে (stratum, বহুবচনে strata) ভাগ করে, যেমন জেলা, তারপর প্রতিটি স্তর থেকে আলাদাভাবে random নমুনা, সাধারণত স্তরের আকারের অনুপাতে। তাহলে বরিশালের ৩% অর্ডার নমুনাতেও প্রায় ৩%. pandas-এ orders.groupby("district").sample(frac=0.02, random_state=...): পর্ব ২৪-এর split-apply-combine, যেখানে “apply” ধাপটি নমুনা তোলা। মডেলের train/test ভাগে class imbalance থাকলে একই ভাবনা (stratified split) লাগে।

নিচের কোডে ১০,০০০টি অর্ডারের একটি কাল্পনিক population, ছয়টি জেলায়। দুই পদ্ধতিতে ২০০টি করে নমুনা ১০০ বার তোলা হয়েছে, আর প্রতিবার মাপা হয়েছে: নমুনায় কোনো জেলার ভাগ population-এর ভাগ থেকে সবচেয়ে বেশি কতটা সরে গেছে।

import numpy as np
import pandas as pd

rng = np.random.default_rng(42)
share = {"ঢাকা": 0.50, "চট্টগ্রাম": 0.20, "খুলনা": 0.12,
         "রাজশাহী": 0.10, "সিলেট": 0.05, "বরিশাল": 0.03}
orders = pd.DataFrame({"district": rng.choice(list(share), size=10_000,
                                              p=list(share.values()))})
true = orders["district"].value_counts(normalize=True)
print(true.round(4).to_string())

def worst_share_error(sample):
    s = sample["district"].value_counts(normalize=True)
    return (s.reindex(true.index, fill_value=0) - true).abs().max()

simple, strat = [], []
for i in range(100):                             # ১০০ বার নমুনা তোলা
    a = orders.sample(n=200, random_state=i)     # simple random
    b = orders.groupby("district").sample(frac=0.02, random_state=i)
    simple.append(worst_share_error(a))
    strat.append(worst_share_error(b))
print("simple random: mean", round(np.mean(simple), 4), "| worst", round(np.max(simple), 4))
print("stratified   : mean", round(np.mean(strat), 4), "| worst", round(np.max(strat), 4))
print("stratified sample size:", len(b))
print(b["district"].value_counts().to_dict())

district
ঢাকা         0.5015
চট্টগ্রাম    0.2043
খুলনা        0.1180
রাজশাহী      0.0977
সিলেট        0.0480
বরিশাল       0.0305
simple random: mean 0.0395 | worst 0.1035
stratified   : mean 0.004 | worst 0.004
stratified sample size: 201
{'ঢাকা': 100, 'চট্টগ্রাম': 41, 'খুলনা': 24, 'রাজশাহী': 20, 'সিলেট': 10, 'বরিশাল': 6}

Sampling bias — বড় নমুনাও বাঁচায় না

Random না হলে নমুনা যত বড়ই হোক, ভুল দিকে নিখুঁত হয়। একে বলে sampling bias: নমুনায় ঢোকার সম্ভাবনা নিজেই এমন কিছুর ওপর নির্ভর করে, যা আমরা মাপতে চাই। একটি হাতে হিসাবের উদাহরণ। ধরুন বাজারঘরের সব ক্রেতার সত্যিকারের রেটিং এমন: ১★ ১০%, ২★ ১০%, ৩★ ২০%, ৪★ ৩০%, ৫★ ৩০%; সত্যিকারের গড় 0.1 + 0.2 + 0.6 + 1.2 + 1.5 = 3.6. কিন্তু review লেখেন বেশি খুশিরা: ১★ আর ২★ ক্রেতার ৫%, ৩★-এর ১০%, ৪★-এর ২০%, ৫★-এর ৪০%. তাহলে review-এর মধ্যে ভাগ হয় 0.005, 0.005, 0.02, 0.06, 0.12 (মোট 0.21), আর review-এর গড় (0.005 + 0.01 + 0.06 + 0.24 + 0.60) / 0.21 = 0.915 / 0.21 ≈ 4.36. দশ লাখ review হলেও গড় 4.36-এর কাছেই থাকবে, 3.6-এর দিকে যাবে না। বাস্তবে অসন্তুষ্টরাও অনেক সময় বেশি লেখেন; মূল কথা, কারা নমুনায় ঢুকছে তা random নয়।

পক্ষপাতের পরিচিত কয়েকটি রূপ:

• Survivorship bias: শুধু যারা “টিকে আছে” তাদের দেখা। “বাজারঘরে ৫ বছর ধরে যে বিক্রেতারা আছেন, তাঁদের সবার রেটিং ভালো, তাই আমাদের প্ল্যাটফর্মে সবাই সফল” — অথচ যাঁরা টেকেননি, তাঁরা ডেটাতেই নেই।
• সুবিধার নমুনা (convenience sample): শুধু ঢাকার ক্রেতাকে ফোন, কারণ তাঁদের নম্বর হাতের কাছে।
• উত্তর না দেওয়া (non-response): ১,০০০ জনকে জরিপ পাঠানো হলো, ১৫০ জন উত্তর দিলেন; যাঁরা দিলেন না, তাঁরা হয়তো সবচেয়ে ব্যস্ত বা সবচেয়ে বিরক্ত।

AI-এর জন্য এটি সরাসরি প্রাসঙ্গিক: একটি মডেল তার training ডেটার distribution শেখে। ডেটা যদি একটি পক্ষপাতদুষ্ট নমুনা হয় (শুধু শহরের ক্রেতা, শুধু ইংরেজি review), মডেল বাকিদের জন্য খারাপ কাজ করতে পারে, আর একই রকম test ডেটায় মাপলে কেউ টের পায় না।

৯ · Monte Carlo আর বড় সংখ্যার নিয়ম — π, ক্ষমতা ছাড়ানো আর 1/√n

পর্ব ১৬-এ Monte Carlo দিয়ে দেরিতে ডেলিভারির সম্ভাবনা অনুমান করেছেন: অনেকবার random ঘটনা ঘটিয়ে “হ্যাঁ”-এর ভাগ গোনা। এখানে দুটি নতুন প্রয়োগ আর তার পেছনের নিয়ম। প্রথমটি বিখ্যাত খেলনা-সমস্যা: π-এর অনুমান। এক বাহুর একটি বর্গের ভেতরে uniform random বিন্দু ফেলুন। বর্গের এক কোণকে কেন্দ্র করে ব্যাসার্ধ ১-এর একটি পোয়া বৃত্ত আঁকা; তার ক্ষেত্রফল π/4, বর্গের ক্ষেত্রফল ১। তাই বিন্দুর মোটামুটি π/4 ভাগ বৃত্তের ভেতরে পড়বে, আর

π ≈ 4 × (বৃত্তের ভেতরের বিন্দু) / (মোট বিন্দু)

একটি বিন্দু (x, y) ভেতরে কি না, তা পরীক্ষা x² + y² ≤ 1. নিচে লেখকের যন্ত্রে default_rng(42) দিয়ে একটি run, ধাপে ধাপে বাড়ানো (NumPy 2.5.1).

দ্বিতীয় প্রয়োগটি ব্যবসার: বাজারঘরের গুদাম দিনে ৫২০টি অর্ডার প্যাক করতে পারে, আর দৈনিক অর্ডার মোটামুটি Poisson, λ = 480. কোনো দিন ক্ষমতা ছাড়ানোর সম্ভাবনা কত? এখানে সঠিক উত্তর stats.poisson(480).sf(520) দিয়েই পাওয়া যায়, আর ঠিক সেই কারণেই উদাহরণটি ভালো: Monte Carlo-র ভুল সত্যিই মাপা যায়। মডেল জটিল হলে সূত্র থাকে না, কিন্তু simulation একই রকম লেখা যায়।

নিচের কোড প্রতিটি n-এ অনুমানটি ৩০ বার চালিয়ে ভুলের RMS (বর্গের গড়ের বর্গমূল, মানে একটি সাধারণ run কতটা দূরে) মেপেছে, তারপর log-log ছবিতে এঁকেছে: দুই অক্ষেই log স্কেল (log scale), যেখানে ভুল = c/√n একটি ঢাল −½-এর সরলরেখা। ডানের ছবিতে law of large numbers: একটি run-এর চলতি অনুমান (running estimate) কীভাবে π-এর দিকে স্থির হয়।

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

rng = np.random.default_rng(42)
ns = np.array([100, 1_000, 10_000, 100_000])
cap_exact = stats.poisson(480).sf(520)          # P(দিনের অর্ডার > 520)
pi_err, cap_err = [], []
for n in ns:
    e1, e2 = [], []
    for _ in range(30):                         # প্রতিটি n-এ ৩০টি run
        xy = rng.random((n, 2))
        e1.append(4 * np.mean((xy**2).sum(axis=1) <= 1) - np.pi)
        e2.append(np.mean(rng.poisson(480, size=n) > 520) - cap_exact)
    pi_err.append(np.sqrt(np.mean(np.square(e1))))   # RMS ভুল
    cap_err.append(np.sqrt(np.mean(np.square(e2))))
    print(f"n={n:>7,}  pi RMS err={pi_err[-1]:.4f}  overflow RMS err={cap_err[-1]:.5f}")
print("exact P(orders > 520) =", round(cap_exact, 4))
print("log-log slope:", np.polyfit(np.log10(ns), np.log10(pi_err), 1)[0].round(3),
      np.polyfit(np.log10(ns), np.log10(cap_err), 1)[0].round(3))

xy = rng.random((5_000, 2))
running = 4 * np.cumsum((xy**2).sum(axis=1) <= 1) / np.arange(1, 5_001)
print("running estimate at n = 10, 100, 1000, 5000:", running[[9, 99, 999, 4999]].round(4))
fig, (a1, a2) = plt.subplots(1, 2, figsize=(7, 3.2), layout="constrained")
a1.loglog(ns, pi_err, "o-", label="pi")
a1.loglog(ns, cap_err, "s-", label="P(overflow)")
a1.loglog(ns, 1.64 / np.sqrt(ns), "k--", label="1.64/sqrt(n)")
a1.set_xlabel("n (samples)"); a1.set_ylabel("RMS error"); a1.legend(fontsize=8)
a2.plot(np.arange(1, 5_001), running); a2.axhline(np.pi, color="k", ls="--")
a2.set_xscale("log"); a2.set_xlabel("n"); a2.set_ylabel("running estimate of pi")
a2.set_title("Law of large numbers")
plt.show()

n=    100  pi RMS err=0.1539  overflow RMS err=0.02363
n=  1,000  pi RMS err=0.0485  overflow RMS err=0.00522
n= 10,000  pi RMS err=0.0157  overflow RMS err=0.00191
n=100,000  pi RMS err=0.0064  overflow RMS err=0.00045
exact P(orders > 520) = 0.0336
log-log slope: -0.464 -0.559
running estimate at n = 10, 100, 1000, 5000: [3.2    3.08   3.168  3.1432]

Law of large numbers সহজ ভাষায়: একই distribution থেকে স্বাধীন নমুনা যত বাড়ে, নমুনার গড় তত distribution-এর expected value-এর কাছে যায় (যদি expected value থাকে; অংশ ৭-এর ভারী লেজের সতর্কতা মনে রাখুন)। এটি বলে না যে ভুলগুলো “শোধ” হয়: ১০টি বিন্দুর পর অনুমান কম ছিল বলে পরের বিন্দুগুলো ভেতরে পড়ার সম্ভাবনা বাড়ে না, যা জুয়াড়ির ভুল ধারণা (gambler's fallacy). আগের ভুল শুধু অনেক নতুন নমুনায় “পাতলা” হয়ে যায়। আর নমুনার গড় কতটা ছড়ায় ও কেন তার distribution প্রায় normal হয়, সেই Central Limit Theorem আসবে পর্ব ৩৭-এ।

পুনরুৎপাদনের নিয়ম পর্ব ১৬-এর মতোই: একটি rng = np.random.default_rng(seed) বানিয়ে সব জায়গায় সেটিই পাঠান, scipy-র rvs-এ random_state=rng, আর pandas-এর sample-এও একটি পূর্ণসংখ্যা বা Generator. Seed-এর পাশে NumPy আর SciPy-র version লিখে রাখুন (পর্ব ১২), আর একটি seed-এর ফল নয়, কয়েকটি seed-এ ফলের ছড়ানো রিপোর্ট করুন।

১০ · AI-এর যুগে — distribution আর simulation-এর কোড যাচাই

“অর্ডারের ডেটা থেকে একটি distribution fit করো” বা “ক্ষমতা ছাড়ানোর সম্ভাবনা simulate করো” লিখলে AI সহকারী দ্রুত কোড আর আত্মবিশ্বাসী সংখ্যা দেবে। এই বিষয়ে তাদের সাধারণ ভুলগুলো:

• পুরোনো বা বানানো API: np.random.seed আর np.random.poisson-এর global পথ; NumPy 2-এ সরানো np.float; pandas-এর সরানো append; আর এমন parameter যা নেই, যেমন stats.norm(mean=500, sd=40) (আসল নাম loc আর scale). stats.uniform(a, b)-কে [a, b] ভাবা আরেকটি পরিচিত ভুল।
• Parameter-এর সংজ্ঞা গুলিয়ে ফেলা: exponential-এ scale=λ লেখা, অথচ scale = 1/λ; geometric শুরু 0 না 1 থেকে, তা না দেখা।
• সব কিছুকে normal ধরা: skewed মূল্যের ডেটায় μ ± 2σ দিয়ে “সীমা” বলা, ঋণাত্মক সময় বা দামকে সম্ভব বানানো।
• Plausible কিন্তু ভুল পরিসংখ্যান: continuous-এ pdf(x)-কে সম্ভাবনা বলা; এক দিকের আর দুই দিকের ৯৫% (1.645 বনাম 1.96) গুলিয়ে ফেলা; অল্প trial-এর Monte Carlo ফল ভুলের মাপ ছাড়া দেওয়া।
• Sampling-এ ফাঁক: পুরো ডেটা দেখে distribution বা scaling ঠিক করে তারপর train/test ভাগ (data leakage), train ডেটাতেই মাপা ফল দেখানো, আর পক্ষপাতদুষ্ট নমুনা থেকে population-এর দাবি।

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

• প্রতিটি scipy.stats distribution-এর loc/scale-এর অর্থ নিজের SciPy version-এর অফিসিয়াল docs-এ মিলিয়েছি, আর X.mean() ছেপে দেখেছি তা আমার ধারণার সঙ্গে মেলে
• অন্তত একটি সম্ভাবনা হাতে বা সূত্রে হিসাব করে কোডের সঙ্গে মিলিয়েছি (যেমন P(X = 0))
• ডেটার histogram বা ECDF দেখে তবেই normal বা অন্য মডেল ধরেছি; skew আর লেজ আলাদা করে দেখেছি
• Simulation-এর ফলের সঙ্গে n আর standard error আছে, আর যেখানে সঠিক উত্তর জানা, সেখানে simulation তার সঙ্গে মিলেছে
• একটি default_rng(seed) সব জায়গায় পাঠানো, seed আর library-র version লেখা; কয়েকটি seed-এ ফলের ছড়ানো দেখেছি
• নমুনা কীভাবে তোলা হয়েছে, কারা বাদ পড়েছে, তা লিখেছি; train/test ভাগ সবকিছুর আগে, আর আসল ক্রেতার ডেটা AI-তে পেস্ট করিনি