Call Center Operations को Probability Distributions के साथ Model करना¶

import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats
import kagglehub

plt.style.use("seaborn-v0_8")
sns.set_context("talk")

# Data download करें
path = kagglehub.dataset_download("akash1vishwakarma/call-center-dataset")

df = pd.read_excel(path+'/01 Call-Center-Dataset.xlsx')
df.columns = [c.strip().replace(" ", "_").replace("(Y/N)", "YN") for c in df.columns]
df

# Date और time को clean करें और timestamp variable बनाएं
df['Date'] = pd.to_datetime(df['Date'], errors='coerce')
df['Time_td'] = pd.to_timedelta(df['Time'].astype(str)) # HH:MM:SS -> timedelta
df['timestamp'] = df['Date'] + df['Time_td']

# Events को timestamp के अनुसार sort करें। (सभी subsequent computations chronological order पर निर्भर करते हैं।)
df = df.dropna(subset=['timestamp']).sort_values('timestamp').reset_index(drop=True)

df['Answered_YN'] = df['Answered_YN'].apply(lambda x: 1 if x=='Y' else 0)
df['Resolved'] = df['Resolved'].apply(lambda x: 1 if x=='Y' else 0)

# Talk duration को seconds में parse करें
def parse_duration(s):
    if pd.isna(s): return np.nan
    try:
        h, m, s2 = 0, 0, 0
        parts = str(s).split(':')
        parts = [int(p) for p in parts]
        if len(parts) == 3: h, m, s2 = parts
        elif len(parts) == 2: m, s2 = parts
        else: return float(parts[0])
        return h*3600 + m*60 + s2
    except Exception:
        return np.nan

df['TalkDuration_seconds'] = df['AvgTalkDuration'].apply(parse_duration)

# Data type conversions
df['Speed_of_answer_in_seconds'] = pd.to_numeric(df['Speed_of_answer_in_seconds'], errors='coerce')

# Calculate: counts/hour
hourly_counts = df.set_index('timestamp').groupby([pd.Grouper(freq='h'), 'Agent']).size().reset_index()
hourly_counts.columns = ['timestamp', 'Agent', 'no_of_calls']

# एक dataset बनाएं जिसमें agent-timestamp (hours) का सभी combination हो
df['dummy_col'] = 1
outer_join_df = df.groupby('Agent').dummy_col.max().reset_index().merge(
    df.groupby('timestamp').dummy_col.max().reset_index(), on='dummy_col')

# datasets को join करें ताकि complete dataset मिल सके
hourly_counts = outer_join_df.merge(hourly_counts, on=['timestamp', 'Agent'], how='outer')
hourly_counts['no_of_calls'] = hourly_counts.no_of_calls.fillna(0)

# Metrics
lambda_hat = hourly_counts.no_of_calls.mean()
var_hat = hourly_counts.no_of_calls.var(ddof=1)

# Sample lambda distribution
max_k = int(hourly_counts.no_of_calls.max())
k_vals = np.arange(0, max_k+1)
pmf = stats.poisson.pmf(k_vals, mu=lambda_hat)

# Distribution को plot करें
fig, ax = plt.subplots(figsize=(10,6))
sns.histplot(hourly_counts.no_of_calls, bins=np.arange(0, max_k+2), stat='probability', color="gold", ax=ax)
ax.plot(k_vals, pmf, 'o-', color="darkgreen", label=f"Poisson(λ={lambda_hat:.2f}) PMF")
ax.set_title('Calls per Hour vs Poisson PMF')
ax.set_xlabel('Calls per hour')
ax.legend()
plt.show();

# Exponential: inter-arrival times
ts = df['timestamp'].sort_values().values
inter_arrival_sec = np.diff(ts) / np.timedelta64(1, 's')
inter_arrival_sec = inter_arrival_sec[inter_arrival_sec > 0]
inter_arrival_min = inter_arrival_sec / 60.0

# Metrics
lambda_exp_hat = 1.0 / np.mean(inter_arrival_min)

fig, ax = plt.subplots(figsize=(10,6))
sns.histplot(inter_arrival_min, bins=40, stat='density', color="gold", ax=ax)
x = np.linspace(0, np.percentile(inter_arrival_min, 99), 200)
ax.plot(x, stats.expon.pdf(x, scale=1/lambda_exp_hat), color="green", lw=2, 
        label=f"Expon PDF (λ={lambda_exp_hat:.3f} PDF)")
ax.set_title('Inter-arrival Times (minutes) with Exponential Fit')
ax.set_xlabel('Minutes')
ax.legend()
plt.show();

# Failures before first success (per agent)

# प्रत्येक agent के chronological sequence of resolved flags के लिए, अगली Y (पहली success) से पहले N (failures) की संख्या गिनें।
geo_failures = []
for agent, grp in df.sort_values('timestamp').groupby('Agent'):
    seq = grp['Resolved'].dropna().values
    fail = 0
    for v in seq:
        if v == 1:
            geo_failures.append(fail)
            fail = 0
        elif v == 0:
            fail += 1

geo_failures = np.array(geo_failures)

# Agents के बीच “success से पहले failures” पर Geometric distribution को fit करें; 
# success probability p (MLE via mean failures) का अनुमान लगाएं।
if len(geo_failures) > 0:
    mean_fail = float(np.mean(geo_failures))
    p_hat = 1.0 / (mean_fail + 1.0)
    k_max_g = int(np.percentile(geo_failures, 99))
    k_vals_g = np.arange(0, max(5, k_max_g)+1)
    pmf_geom = stats.geom.pmf(k_vals_g+1, p_hat)

    fig, ax = plt.subplots(figsize=(10,6))
    sns.histplot(geo_failures, bins=range(0, k_vals_g[-1]+2), stat='probability', 
                 color="gold", ax=ax)
    ax.plot(k_vals_g, pmf_geom, 'o-', color="green", label=f"Geom(p={p_hat:.3f}) PMF (shifted)")
    ax.set_title('Failures before First Success vs Geometric PMF')
    ax.set_xlabel('Failures before success')
    ax.legend()
    plt.show();

# Poisson distribution
chi2_stat_poisson, p_poisson = stats.chisquare(hourly_counts.groupby('no_of_calls').timestamp.count(), 
                                               f_exp=pmf*hourly_counts.no_of_calls.count())
print(f"Poisson Chi-Square p-value: {p_poisson:.4f}")

fig, ax = plt.subplots(figsize=(10,6))
res = stats.probplot(hourly_counts.no_of_calls, dist=stats.poisson, 
                     sparams=(hourly_counts.no_of_calls.mean(),), plot=ax)
ax.set_title("Poisson QQ-Plot for lambda = "+str(round(hourly_counts.no_of_calls.mean(),2)))
plt.show();

# Exponential distribution
s_D, ks_p = stats.kstest(inter_arrival_min, 'expon', args=(0, 1/lambda_exp_hat))
print(f"Exponential Kolmogorov-Smirnov test p-value: {p_poisson:.4f}")

fig, ax = plt.subplots(figsize=(10,6))
probs = np.linspace(0.01, 0.99, 99)
emp_q = np.quantile(inter_arrival_min, probs)
theo_q = stats.expon.ppf(probs, scale=1/lambda_exp_hat)
ax.plot(theo_q, emp_q, 'o', alpha=0.7)
ax.plot(theo_q, theo_q, 'r--', lw=2)
ax.set_title('Exponential QQ-Plot')
ax.set_xlabel('Theoretical quantiles')
ax.set_ylabel('Empirical quantiles')
plt.show();

# Geometric distribution
n_fail = np.mean(geo_failures)
p_hat = 1 / (mean_fail + 1)
k_vals_g = np.arange(0, int(np.percentile(geo_failures, 99))+1)
pmf_geom = stats.geom.pmf(k_vals_g+1, p_hat)

emp_counts = pd.Series(geo_failures).value_counts().sort_index()
obs_g = np.array([emp_counts.get(k, 0) for k in k_vals_g])
exp_g = pmf_geom * obs_g.sum()

chi2_stat_geom, p_geom = stats.chisquare(obs_g, f_exp=exp_g*np.sum(obs_g)/np.sum(exp_g))
print(f"Geometric Chi-Square p-value: {p_geom:.4f}")


fig, ax = plt.subplots(figsize=(10,6))
ax.scatter(obs_g, exp_g, s=40, alpha=0.75, edgecolor='k')
max_q = max(obs_g.max(), exp_g.max())
ax.plot([0, max_q], [0, max_q], 'r--', lw=2)
ax.set_xlabel('Theoretical quantiles (Geometric, failures)')
ax.set_ylabel('Empirical quantiles (sample)')
ax.set_title(f'Geometric Q–Q Plot (p̂={p_hat:.3f})')

Poisson Chi-Square p-value: 0.0000

Exponential Kolmogorov-Smirnov test p-value: 0.0000

Geometric Chi-Square p-value: 0.7642





Text(0.5, 1.0, 'Geometric Q–Q Plot (p̂=0.730)')

जब आप किसी भी population से repeated samples लेते हैं और उनके means की गणना करते हैं, तो उन sample means का distribution एक normal distribution के करीब पहुंचता है जब sample size बढ़ता है—चाहे मूल population का आकार कुछ भी हो।

sample_size = 30
n_samples = 1000

# Sample means generate करें
sample_means = []
for _ in range(n_samples):
    sample = np.random.choice(hourly_counts.no_of_calls, size=sample_size, replace=True)
    sample_means.append(sample.mean())

# Sample means का histogram plot करें
plt.figure(figsize=(10,6))
plt.hist(sample_means, bins=10, color="#86BC25", edgecolor='#2E7C26', density=True)
plt.title(f"CLT in Action: Sample Means का Distribution (n={sample_size})")
plt.xlabel("Sample Mean (calls per hour)")
plt.ylabel("Density")

# Normal curve overlay करें
mean = np.mean(sample_means)
std = np.std(sample_means)
x = np.linspace(min(sample_means), max(sample_means), 200)
plt.plot(x, stats.norm.pdf(x, loc=mean, scale=std), 'r--', lw=2, label="Normal Approximation")
plt.legend()
plt.show();