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

df = pd.read_csv('superstore_sales.csv')
df['Order Date'] = pd.to_datetime(df['Order Date'])
print(df.shape)

p_tech = (df['Category'] == 'Technology').mean()
n = 20
print(f'Real P(Technology) per order line: {p_tech:.4f}')
print(f'Modeling batches of n={n} order lines')

rng = np.random.default_rng(seed=7)
n_batches = 10000
simulated_counts = rng.binomial(n=n, p=p_tech, size=n_batches)
print(f'Simulated mean Technology count per batch: {simulated_counts.mean():.3f}')
print(f'Theoretical mean (n*p): {n * p_tech:.3f}')

k_values = np.arange(0, n + 1)
pmf = stats.binom.pmf(k_values, n, p_tech)
for k, prob in zip(k_values[:6], pmf[:6]):
    print(f'P(exactly {k} Technology orders in a batch of {n}) = {prob:.4f}')

plt.figure(figsize=(9, 4))
plt.bar(k_values, pmf, color='steelblue', alpha=0.7, label='Theoretical PMF')
sim_probs = pd.Series(simulated_counts).value_counts(normalize=True).sort_index()
plt.scatter(sim_probs.index, sim_probs.values, color='darkred', zorder=3, label='Simulated frequencies')
plt.title(f'Binomial(n={n}, p={p_tech:.3f}): Theoretical vs. Simulated')
plt.xlabel('Number of Technology Orders in Batch')
plt.ylabel('Probability')
plt.legend()
plt.show()

mean_binom, var_binom = stats.binom.stats(n, p_tech, moments='mv')
print(f'Theoretical mean: {float(mean_binom):.3f}')
print(f'Theoretical variance: {float(var_binom):.3f}')
print(f'Simulated mean: {simulated_counts.mean():.3f}, simulated variance: {simulated_counts.var():.3f}')

p_at_most_2 = stats.binom.cdf(2, n, p_tech)
p_at_least_5 = 1 - stats.binom.cdf(4, n, p_tech)
print(f'P(at most 2 Technology orders) = {p_at_most_2:.4f}')
print(f'P(at least 5 Technology orders) = {p_at_least_5:.4f}')

daily_counts = df.groupby(df['Order Date'].dt.date).size()
lam = daily_counts.mean()
print(f'Number of real days with at least one order: {len(daily_counts)}')
print(f'Real average daily order-line count (lambda): {lam:.3f}')

k_values2 = np.arange(0, daily_counts.max() + 1)
poisson_pmf = stats.poisson.pmf(k_values2, lam)
observed_freq = daily_counts.value_counts(normalize=True).sort_index()
plt.figure(figsize=(9, 4))
plt.bar(k_values2, poisson_pmf, color='steelblue', alpha=0.7, label=f'Poisson(lambda={lam:.2f}) PMF')
plt.scatter(observed_freq.index, observed_freq.values, color='darkred', zorder=3, label='Real observed frequencies')
plt.title('Real Daily Order Counts vs. a Fitted Poisson')
plt.xlabel('Order Lines per Day')
plt.ylabel('Probability')
plt.legend()
plt.show()

var_real = daily_counts.var()
print(f'Real daily count mean: {lam:.3f}')
print(f'Real daily count variance: {var_real:.3f}')
print(f'Poisson assumes mean == variance; real ratio (variance/mean) = {var_real / lam:.2f}')

p_zero_or_one = stats.poisson.cdf(1, lam)
p_busy_day = 1 - stats.poisson.cdf(15, lam)
print(f'Modeled P(0 or 1 order lines in a day) = {p_zero_or_one:.4f}')
print(f'Modeled P(more than 15 order lines in a day) = {p_busy_day:.4f}')
real_p_busy = (daily_counts > 15).mean()
print(f'Real observed P(more than 15 order lines) = {real_p_busy:.4f}')

big_n, small_p = 2000, p_tech / 100
approx_lambda = big_n * small_p
binom_pmf_tail = stats.binom.pmf(np.arange(0, 8), big_n, small_p)
poisson_pmf_tail = stats.poisson.pmf(np.arange(0, 8), approx_lambda)
comparison = pd.DataFrame({'k': np.arange(0, 8), 'Binomial': binom_pmf_tail, 'Poisson_approx': poisson_pmf_tail})
print(comparison.round(5))