import seaborn as sns
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
import hashlib
from scipy.stats import ttest_ind, t

rng = np.random.default_rng(42)

n_samples = 20000          # общий размер выборки
p_zero = 0.4               # доля нулевых значений

# Параметры для положительных значений (смесь двух распределений)
# 1-я компонента: Гамма-распределение (основная масса плательщиков)
gamma_shape = 2.0
gamma_scale = 10.0         # среднее = shape * scale = 20

# 2-я компонента: Логнормальное распределение (киты)
lognorm_mu = 8
lognorm_sigma = 1       # среднее = exp(mu + sigma^2/2) ≈ 4.5, но из-за хвоста будут большие значения

# Вероятности компонент среди плательщиков
prob_gamma = 0.8           # 80% плательщиков – из Гаммы
prob_lognorm = 0.2         # 20% – из Логнормального

# 1. Генерируем индикатор платящий/неплатящий
is_payer = rng.binomial(1, 1 - p_zero, size=n_samples)

# 2. Количество плательщиков
n_payers = np.sum(is_payer)

# 3. Для каждого плательщика выбираем компоненту смеси
components = rng.choice(['gamma', 'lognorm'], size=n_payers, p=[prob_gamma, prob_lognorm])

# 4. Генерируем значения для каждой компоненты
positive_values = np.zeros(n_payers)

# Гамма-часть
mask_gamma = components == 'gamma'
n_gamma = np.sum(mask_gamma)
positive_values[mask_gamma] = np.random.gamma(shape=gamma_shape, scale=gamma_scale, size=n_gamma)

# Логнормальная часть
mask_lognorm = components == 'lognorm'
n_lognorm = np.sum(mask_lognorm)
positive_values[mask_lognorm] = rng.lognormal(mean=lognorm_mu, sigma=lognorm_sigma, size=n_lognorm)

# 5. Формируем итоговый массив: нули + положительные значения
arpu_samples = np.zeros(n_samples)
arpu_samples[is_payer == 1] = positive_values

# Выводим основные статистики
print(f"Всего значений: {len(arpu_samples)}")
print(f"Доля нулей: {np.mean(arpu_samples == 0):.3f}")
print(f"Среднее (включая нули): {np.mean(arpu_samples):.2f}")
print(f"Среднее среди плательщиков: {np.mean(positive_values):.2f}")
print(f"Максимум: {np.max(arpu_samples):.2f}")
print(f"95-й процентиль: {np.percentile(arpu_samples, 95):.2f}")


# Настройка стиля для научной публикации
sns.set_theme(style="white")          # стиль с белым фоном и сеткой
plt.rcParams['font.family'] = 'serif'      # шрифт с засечками (часто используется в статьях)
plt.rcParams['font.size'] = 12
plt.rcParams['figure.figsize'] = (10, 6)

# Положительные значения
positive = arpu_samples[arpu_samples > 0]

# Визуализация (гистограмма в логарифмической шкале по оси y)
plt.figure(figsize=(10, 5))
plt.hist(arpu_samples[arpu_samples > 0], bins=50, edgecolor='black', alpha=0.7, log=True)
plt.xlabel('ARPU')
plt.ylabel('Частота')
plt.title('Рапределение ARPU')
#plt.grid(axis='y', linestyle='--', alpha=0.5)
ax.grid(False)
ax = plt.gca()
ax.yaxis.set_major_formatter(ticker.LogFormatter(labelOnlyBase=False))
plt.ticklabel_format(style='plain', axis='x')
plt.show()

n = len(arpu_samples)
user_ids = np.arange(n).astype(str)  # строковые идентификаторы

# Соль эксперимента
salt = "experiment_11"

# Функция для назначения группы на основе хэша
def assign_group(user_id, salt, split_point=50):
    """split_point - процент пользователей в тесте (0-100)"""
    # Создаём строку для хэширования
    input_str = str(user_id) + "_" + salt
    # Берём MD5-хэш (можно любой другой)
    hash_obj = hashlib.md5(input_str.encode())
    # Преобразуем первые 8 символов в целое число (достаточно для равномерности)
    hash_int = int(hash_obj.hexdigest()[:8], 16)
    # Остаток от деления на 100
    return 'test' if (hash_int % 100) < split_point else 'control'

# Применяем функцию ко всем user_id
groups = np.array([assign_group(uid, salt) for uid in user_ids])

# Выделяем группы
test_arpu = arpu_samples[groups == 'test']
control_arpu = arpu_samples[groups == 'control']

print(f"Размер теста: {len(test_arpu)}")
print(f"Размер контроля: {len(control_arpu)}")

# Проверка баланса
mean_test = np.mean(test_arpu)
mean_control = np.mean(control_arpu)
uplift = (mean_test - mean_control) / mean_control * 100
t_stat, p_value = ttest_ind(test_arpu, control_arpu, equal_var=False)

print(f"nСреднее в тесте: {mean_test:.4f}")
print(f"Среднее в контроле: {mean_control:.4f}")
print(f"Uplift: {uplift:.2f}%")
print(f"p-value (t-тест): {p_value:.4f}")

# Вычислим среднее контроля до добавления (оно же после добавления двух таких же изменится)
mean_control_before = np.mean(control_arpu)

# Добавляем в контроль два значения, равных его текущему среднему
control_modified = np.append(control_arpu, [mean_control_before, mean_control_before])

# Добавляем в тест два выброса по 400 000
test_modified = np.append(test_arpu, [400000, 400000])

print("n=== После добавления значений ===")
print(f"Размер теста: {len(test_modified)}")
print(f"Размер контроля: {len(control_modified)}")

# --- 4. Расчёт средних, стандартных ошибок и uplift ---
mean_t = np.mean(test_modified)
mean_c = np.mean(control_modified)
std_t = np.std(test_modified, ddof=1)
std_c = np.std(control_modified, ddof=1)
n_t = len(test_modified)
n_c = len(control_modified)

se_t = std_t / np.sqrt(n_t)
se_c = std_c / np.sqrt(n_c)

uplift = (mean_t / mean_c - 1) * 100  # в процентах

# --- 5. Доверительный интервал для uplift через дельта-метод ---
# Приближённая стандартная ошибка для R = mean_t/mean_c - 1 (в долях)
R = mean_t / mean_c
se_R = (1 / mean_c) * np.sqrt(se_t**2 + (R**2) * (se_c**2))

# Степени свободы (по аналогии с t-тестом Уэлча)
# Используем приближение Саттертуэйта для отношения? Для простоты возьмём min(n_t-1, n_c-1)
dof = min(n_t - 1, n_c - 1)

# Квантиль t-распределения
t_crit = t.ppf(0.975, df=dof)

# Доверительный интервал для R (в долях)
ci_R_lower = R - t_crit * se_R
ci_R_upper = R + t_crit * se_R

# Переводим в проценты для uplift (R*100% - 100%)
ci_uplift_lower = (ci_R_lower - 1) * 100
ci_uplift_upper = (ci_R_upper - 1) * 100

print(f"nСреднее теста: {mean_t:.2f}")
print(f"Среднее контроля: {mean_c:.2f}")
print(f"Effect = {uplift:.2f}%")
print(f"95% доверительный интервал: [{ci_uplift_lower:.2f}%, {ci_uplift_upper:.2f}%]")

# Для сравнения: t-тест для абсолютной разницы
t_stat, p_val = ttest_ind(test_modified, control_modified, equal_var=False)
print(f"nT-тест для разности средних: p-value = {p_val:.4f}")

# --- 6. График: точечная оценка uplift с аналитическим доверительным интервалом ---
plt.figure(figsize=(7, 5))
plt.errorbar(x=['Uplift'], y=[uplift], 
             yerr=[[uplift - ci_uplift_lower], [ci_uplift_upper - uplift]],
             fmt='o', color='red', capsize=10, markersize=10, label='Uplift (95% CI, delta method)')
plt.axhline(y=0, color='black', linestyle='--', linewidth=1, label='Ноль')
plt.ylabel('Относительное изменение, %')
plt.title('95% доверительный интервал')
plt.grid(axis='y', linestyle='--', alpha=0.5)
plt.legend()
plt.tight_layout()
plt.show()

n = len(arpu_samples)
user_ids = np.arange(n).astype(str)  # строковые идентификаторы

# Соль эксперимента
salt = "experiment_14"

# Функция для назначения группы на основе хэша
def assign_group(user_id, salt, split_point=50):
    """split_point - процент пользователей в тесте (0-100)"""
    # Создаём строку для хэширования
    input_str = str(user_id) + "_" + salt
    # Берём MD5-хэш (можно любой другой)
    hash_obj = hashlib.md5(input_str.encode())
    # Преобразуем первые 8 символов в целое число (достаточно для равномерности)
    hash_int = int(hash_obj.hexdigest()[:8], 16)
    # Остаток от деления на 100
    return 'test' if (hash_int % 100) < split_point else 'control'

# Применяем функцию ко всем user_id
groups = np.array([assign_group(uid, salt) for uid in user_ids])

# Выделяем группы
test_arpu = arpu_samples[groups == 'test']
control_arpu = arpu_samples[groups == 'control']

print(f"Размер теста: {len(test_arpu)}")
print(f"Размер контроля: {len(control_arpu)}")

# Проверка баланса
mean_test = np.mean(test_arpu)
mean_control = np.mean(control_arpu)
uplift = (mean_test - mean_control) / mean_control * 100
t_stat, p_value = ttest_ind(test_arpu, control_arpu, equal_var=False)

print(f"nСреднее в тесте: {mean_test:.4f}")
print(f"Среднее в контроле: {mean_control:.4f}")
print(f"Uplift: {uplift:.2f}%")
print(f"p-value (t-тест): {p_value:.4f}")

# Вычислим среднее контроля до добавления (оно же после добавления двух таких же изменится)
mean_control_before = np.mean(control_arpu)

# Добавляем в контроль два значения, равных его текущему среднему
control_modified = np.append(control_arpu, [mean_control_before, mean_control_before])

# Добавляем в тест два выброса по 400 000
test_modified = np.append(test_arpu, [400000, 400000])

print("n=== После добавления значений ===")
print(f"Размер теста: {len(test_modified)}")
print(f"Размер контроля: {len(control_modified)}")

# --- 4. Расчёт средних, стандартных ошибок и uplift ---
mean_t = np.mean(test_modified)
mean_c = np.mean(control_modified)
std_t = np.std(test_modified, ddof=1)
std_c = np.std(control_modified, ddof=1)
n_t = len(test_modified)
n_c = len(control_modified)

se_t = std_t / np.sqrt(n_t)
se_c = std_c / np.sqrt(n_c)

uplift = (mean_t / mean_c - 1) * 100  # в процентах

# --- 5. Доверительный интервал для uplift через дельта-метод ---
# Приближённая стандартная ошибка для R = mean_t/mean_c - 1 (в долях)
R = mean_t / mean_c
se_R = (1 / mean_c) * np.sqrt(se_t**2 + (R**2) * (se_c**2))

# Степени свободы (по аналогии с t-тестом Уэлча)
# Используем приближение Саттертуэйта для отношения? Для простоты возьмём min(n_t-1, n_c-1)
dof = min(n_t - 1, n_c - 1)

# Квантиль t-распределения
t_crit = t.ppf(0.975, df=dof)

# Доверительный интервал для R (в долях)
ci_R_lower = R - t_crit * se_R
ci_R_upper = R + t_crit * se_R

# Переводим в проценты для uplift (R*100% - 100%)
ci_uplift_lower = (ci_R_lower - 1) * 100
ci_uplift_upper = (ci_R_upper - 1) * 100

print(f"nСреднее теста: {mean_t:.2f}")
print(f"Среднее контроля: {mean_c:.2f}")
print(f"Uplift = {uplift:.2f}%")
print(f"95% доверительный интервал для uplift (дельта-метод): [{ci_uplift_lower:.2f}%, {ci_uplift_upper:.2f}%]")

# Для сравнения: t-тест для абсолютной разницы
t_stat, p_val = ttest_ind(test_modified, control_modified, equal_var=False)
print(f"nT-тест для разности средних: p-value = {p_val:.4f}")

# --- 6. График: точечная оценка uplift с аналитическим доверительным интервалом ---
plt.figure(figsize=(7, 5))
plt.errorbar(x=['Uplift'], y=[uplift], 
             yerr=[[uplift - ci_uplift_lower], [ci_uplift_upper - uplift]],
             fmt='o', color='red', capsize=10, markersize=10, label='Uplift (95% CI, delta method)')
plt.axhline(y=0, color='black', linestyle='--', linewidth=1, label='Ноль')
plt.ylabel('Относительное изменение, %')
plt.title('95% доверительный интервал')
plt.grid(axis='y', linestyle='--', alpha=0.5)
plt.legend()
plt.tight_layout()
plt.show()