Modeling Call Center Operations with Probability Distributions¶

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")

# Download the data
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

# Clean date and time and create 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']

# Sort events by timestamp. (All subsequent computations rely on 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)

# Parse talk duration to seconds
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)

# Datatype 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']

# Create a dataset that has all combination of agent-timestamp (hours)
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')

# join the datasets to get 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)

# Plot the distribution
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)

# For each agent’s chronological sequence of resolved flags, count the number of N (failures) before the next Y (first success).
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)

# Fit a Geometric distribution on “failures before success” across agents; 
# estimate 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)')

When you take repeated samples from any population and compute their means, the distribution of those sample means approaches a normal distribution as the sample size grows—regardless of the original population’s shape.

sample_size = 30
n_samples = 1000

# Generate sample means
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())

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

# Overlay normal curve
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();