Лабораторная работа 1: Fifty Shades of Coffee

#@title 1. Данные и единицы измерения
from pathlib import Path
import hashlib
import platform
import numpy as np
import pandas as pd
import scipy
from scipy import stats, optimize
import matplotlib
import matplotlib.pyplot as plt
import seaborn as sns

SEED = 42
ALPHA = 0.05
REPOSITORY = "Mironiks007/fifty-shades-of-coffee-stat"
NOTEBOOK_PATH = "stat-lab1-s1.ipynb"  # целевой путь после загрузки в корень репозитория
NOTEBOOK_URL = f"https://github.com/{REPOSITORY}/blob/main/{NOTEBOOK_PATH}"
TEAM = [
    {"name": "Петров Мирон", "github": "Mironiks007", "cat_id": 502585, "coeff": 0.5},
    {"name": "Каплунова Анна", "github": "kaplunovaania", "cat_id": 501572, "coeff": 0.5},
]

columns = ["moisture", "acidity", "body"]
labels = {"moisture": "Влажность (доля; 0.11 = 11%)", "acidity": "Кислотность (балл)", "body": "Полнота (балл)"}
sns.set_theme(style="whitegrid")
# единый стиль графиков
PALETTE = {"main": "#5cb5a6", "accent": "#c8884a", "ink": "#eceae6"}  # dark restyle: fills=amber(accent), lines=teal(main), ink=light
NAMES = {"moisture": "Влажность", "acidity": "Кислотность", "body": "Полнота"}
ESTIMATOR_RU = {"mean": "среднее", "variance": "дисперсия", "var": "дисперсия", "std": "СКО"}
SAMPLING_RU = {"with replacement": "с возвращением", "without replacement": "без возвращения"}
plt.rcParams.update({
    "figure.dpi": 110, "figure.titlesize": 14, "figure.titleweight": "bold",
    "axes.titlesize": 12, "axes.labelsize": 11, "legend.fontsize": 9,
    "axes.prop_cycle": matplotlib.cycler(color=[PALETTE["main"], PALETTE["accent"], PALETTE["ink"]]),
})
pd.set_option("display.max_columns", 20)

# --- наш тёмный векторный стиль «Тёмная обжарка» (перебивает светлую тему seaborn) ---
plt.rcParams.update({
    "font.family": "serif",
    "font.serif": ["DejaVu Serif"],
    "figure.facecolor": "#141418",
    "figure.dpi": 120,
    "axes.facecolor": "#141418",
    "savefig.facecolor": "#141418",
    "text.color": "#eceae6",
    "axes.labelcolor": "#9a97a0",
    "axes.titlecolor": "#eceae6",
    "figure.titlesize": 14, "figure.titleweight": "bold",
    "xtick.color": "#9a97a0",
    "ytick.color": "#9a97a0",
    "axes.edgecolor": "#6d6b75",
    "grid.color": "#30303a",
    "grid.alpha": 1,
    "axes.axisbelow": True,
    "legend.facecolor": "#1e1e24",
    "legend.edgecolor": "#6d6b75",
    "legend.labelcolor": "#eceae6",
    # boxplot: без явных props линии по умолчанию тёмные — поднимаем на тёмном фоне
    "boxplot.boxprops.color": "#9a97a0",
    "boxplot.whiskerprops.color": "#9a97a0",
    "boxplot.capprops.color": "#9a97a0",
    "boxplot.medianprops.color": "#c8884a",
    "boxplot.flierprops.markeredgecolor": "#9a97a0",
})

# тёмный экспорт + вектор (как в s2): инлайн-картинки с нашим фоном
try:
    ip = get_ipython()
    ip.run_line_magic("config", "InlineBackend.figure_format = 'retina'")
    ip.run_line_magic("config", "InlineBackend.print_figure_kwargs = {'facecolor': '#141418', 'bbox_inches': 'tight'}")
except Exception:
    pass


def show_md(text):
    try:
        from IPython.display import display, Markdown
        display(Markdown(text))
    except ImportError:
        print(text)

def show_table(frame, digits=5):
    try:
        from IPython.display import display
        display(frame.round(digits))
    except ImportError:
        print(frame.round(digits).to_string())

def fmt(value, digits=4):
    return f"{float(value):.{digits}f}".replace(".", ",")

def holm_adjust(p_values):
    p = np.asarray(p_values, dtype=float)
    order = np.argsort(p)
    adjusted = np.empty(len(p))
    adjusted[order] = np.minimum(1, np.maximum.accumulate(p[order] * (len(p)-np.arange(len(p)))))
    return adjusted

expected_sha256 = "ae6fdd58f4a7f77397c11c2657e22be8aa8f746bdc5a74cf479e981bd3c25d5e"
candidates = [Path("coffee_ratings.csv"), Path("data/coffee_ratings.csv"),
              Path("../coffee_ratings.csv"), Path("../data/coffee_ratings.csv")]
local_path = next((p for p in candidates if p.is_file()), None)
if local_path is None:
    from urllib.request import urlopen
    urls = [f"https://raw.githubusercontent.com/{REPOSITORY}/main/coffee_ratings.csv",
            "https://raw.githubusercontent.com/rfordatascience/tidytuesday/main/data/2020/2020-07-07/coffee_ratings.csv"]
    download_errors = []
    for url in urls:
        try:
            with urlopen(url, timeout=30) as response:
                payload = response.read()
            if hashlib.sha256(payload).hexdigest() != expected_sha256:
                raise ValueError("Контрольная сумма CSV не совпала")
            local_path = Path("coffee_ratings.csv")
            local_path.write_bytes(payload)
            break
        except Exception as error:
            download_errors.append(f"{url}: {type(error).__name__}: {error}")
    if local_path is None:
        raise FileNotFoundError("Положите исходный coffee_ratings.csv рядом с ноутбуком. " + " | ".join(download_errors))
assert hashlib.sha256(local_path.read_bytes()).hexdigest() == expected_sha256, "Это другой CSV: требуется повторная проверка происхождения"
raw = pd.read_csv(local_path)
assert set(columns + ["species"]).issubset(raw.columns)
show_table(raw["species"].value_counts().to_frame("records"))

#@title 3. Подготовка данных и EDA (критерий 4)
arabica = raw.loc[raw["species"].eq("Arabica")].copy()
complete = arabica.dropna(subset=columns)
sensory_base = complete.loc[complete[["acidity", "body"]].gt(0).all(axis=1)].copy()
data = sensory_base.loc[sensory_base["moisture"].gt(0)].copy()
assert len(data) >= 1000
assert data[columns].gt(0).all().all()
flow_table = pd.DataFrame([
    ["Исходный файл", len(raw)], ["Arabica", len(arabica)],
    ["Полные три показателя", len(complete)], ["Положительные сенсорные баллы", len(sensory_base)],
    ["Основной анализ: также положительная влажность", len(data)],
], columns=["stage", "n"])
show_table(flow_table.set_index("stage"), 0)
show_table(pd.DataFrame({"missing": raw[columns].isna().sum(), "zero": raw[columns].eq(0).sum()}), 0)
show_md(f"После отбора исключили **{len(arabica)-len(data)} из {len(arabica)}** записей Arabica ({fmt(100*(len(arabica)-len(data))/len(arabica),1)}%). Далее проверяем, как это повлияло на результаты.")

#@title 3.1 Как исключение нулей влияет на результаты
zero_moisture = sensory_base.loc[sensory_base["moisture"].eq(0)].copy()
sensitivity_rows = []
for label, subset in [("positive moisture", data), ("zero moisture only", zero_moisture), ("zeros retained", sensory_base)]:
    sensitivity_rows.append({"sample": label, "n": len(subset), "moisture_recorded_mean": subset["moisture"].mean(),
        "acidity_mean": subset["acidity"].mean(), "body_mean": subset["body"].mean(),
        "acidity_body_rho": stats.spearmanr(subset["acidity"], subset["body"]).statistic})
sensitivity_table = pd.DataFrame(sensitivity_rows).set_index("sample")
show_table(sensitivity_table)
country_selection = sensory_base.assign(country=sensory_base["country_of_origin"].fillna("unknown"),
                                         zero=sensory_base["moisture"].eq(0)).groupby("country").agg(
                                             n=("zero","size"), zeros=("zero","sum"))
country_selection["zero_fraction"] = country_selection["zeros"] / country_selection["n"]
show_table(country_selection.sort_values("n", ascending=False).head(10))
selection_effects = pd.DataFrame({
    "positive_mean": data[["acidity","body"]].mean(), "zero_group_mean": zero_moisture[["acidity","body"]].mean(),
    "all_sensory_mean": sensory_base[["acidity","body"]].mean(),
})
selection_effects["difference_zero_minus_positive"] = selection_effects["zero_group_mean"] - selection_effects["positive_mean"]
selection_effects["difference_in_all_sensory_SD"] = selection_effects["difference_zero_minus_positive"] / sensory_base[["acidity","body"]].std(ddof=1)
show_table(selection_effects)
missing_moisture_n = len(zero_moisture)
mean_bound_low = data["moisture"].sum() / len(sensory_base)
mean_bound_high = (data["moisture"].sum() + missing_moisture_n) / len(sensory_base)
show_md(f"Если оставляем нулевую влажность, в анализ входит **{len(sensory_base)}** строк. Корреляция кислотности и полноты составляет **{fmt(sensitivity_table.loc['zeros retained','acidity_body_rho'],3)}**, после отбора - **{fmt(sensitivity_table.loc['positive moisture','acidity_body_rho'],3)}**. Направление и величина связи остаются близкими. Если считать {missing_moisture_n} нулей неизвестными измерениями, получаем границы среднего [{fmt(mean_bound_low)}; {fmt(mean_bound_high)}]. Среднее по записям с положительной влажностью равно **{fmt(data['moisture'].mean())}**, или **{fmt(100*data['moisture'].mean(),2)}%**.")

#@title 3.2 Методы анализа
normality_rows = []
interval_rows = []
for column in columns:
    values = data[column].to_numpy()
    normality = stats.normaltest(values)
    normality_rows.append({"variable": column, "statistic": normality.statistic, "p_value": normality.pvalue})
    n = len(values)
    mean = values.mean()
    standard_error = values.std(ddof=1) / np.sqrt(n)
    t_low, t_high = stats.t.interval(0.95, df=n-1, loc=mean, scale=standard_error)
    bootstrap_rng = np.random.default_rng(SEED)
    bootstrap_means = bootstrap_rng.choice(values, size=(2000, n), replace=True).mean(axis=1)
    b_low, b_high = np.quantile(bootstrap_means, [0.025, 0.975])
    interval_rows.append({"variable": column, "n": n, "mean": mean, "t_low": t_low, "t_high": t_high, "bootstrap_low": b_low, "bootstrap_high": b_high})
normality_table = pd.DataFrame(normality_rows).set_index("variable")
normality_table["p_holm"] = holm_adjust(normality_table["p_value"])
interval_table = pd.DataFrame(interval_rows).set_index("variable")
print("Ориентировочный normaltest (округление не учтено в его калибровке):")
print(normality_table.to_string(float_format=lambda value: f"{value:.5g}"))
print("Средние и 95% доверительные интервалы:")
print(interval_table.round(5).to_string())

#@title Разведочный анализ и распределение
x = data["moisture"].to_numpy()
q1, q3 = np.quantile(x, [0.25, 0.75])
iqr = q3 - q1
flagged = (x < q1 - 1.5 * iqr) | (x > q3 + 1.5 * iqr)
summary = pd.Series({"n": len(x), "mean": x.mean(), "median": np.median(x), "variance_ddof1": x.var(ddof=1), "std_ddof1": x.std(ddof=1), "min": x.min(), "q1": q1, "q3": q3, "max": x.max(), "skewness": stats.skew(x, bias=False), "IQR_flags": flagged.sum()})
print(summary.round(5).to_string())

fig, axes = plt.subplots(2, 2, figsize=(14, 8))
axes[0, 0].hist(x, bins=np.arange(x.min()-.005, x.max()+.015, .01), density=True, color=PALETTE["accent"], alpha=0.75)
grid = np.linspace(x.min(), x.max(), 300)
axes[0, 0].plot(grid, stats.norm.pdf(grid, x.mean(), x.std(ddof=1)), color=PALETTE["main"], label="Нормальное приближение")
axes[0, 0].set(xlabel=labels["moisture"], ylabel="Плотность", title="Гистограмма")
axes[0, 0].legend()
axes[0, 1].step(np.sort(x), np.arange(1, len(x)+1) / len(x), where="post", color=PALETTE["main"])
axes[0, 1].set(xlabel=labels["moisture"], ylabel="F_n(x)", title="Эмпирическая функция распределения", ylim=(0, 1.02))
axes[1, 0].boxplot(x, vert=False, medianprops={"color": PALETTE["accent"], "linewidth": 2}, flierprops={"markeredgecolor": PALETTE["ink"], "markersize": 4})
axes[1, 0].set(xlabel=labels["moisture"], title="Ящик с усами (выбросы оставлены)", yticks=[])
stats.probplot(x, dist="norm", plot=axes[1, 1])
axes[1, 1].get_lines()[0].set(color=PALETTE["accent"], markersize=3)
axes[1, 1].get_lines()[1].set(color=PALETTE["main"])
axes[1, 1].set(xlabel="Теоретические квантили", ylabel="Упорядоченные значения")
axes[1, 1].set_title("Q-Q график против нормального")
fig.suptitle(f"Влажность зёрен: Arabica, n={len(x)}")
fig.tight_layout()
plt.show()

show_md(
    f"**Вывод по графикам влажности**\n\n"
    f"На гистограмме видим узкий пик около 11%. Центральная половина наблюдений лежит между "
    f"**{fmt(100*q1,1)} и {fmt(100*q3,1)}%**; на CDF этому участку соответствует резкий подъём. "
    f"Крупные ступени показывают, что одни и те же записанные значения встречаются много раз\n\n"
    f"Асимметрия равна **{fmt(stats.skew(x,bias=False),2)}**: распределение сильнее смещено в сторону малой влажности. "
    f"При этом редкие большие значения тоже есть: весь диапазон составляет "
    f"{fmt(100*x.min(),1)}–{fmt(100*x.max(),1)}%. Boxplot отмечает {int(flagged.sum())} наблюдения по правилу 1,5 IQR, "
    f"но само это правило не устанавливает ошибку измерения, поэтому эти строки сохраняем\n\n"
    "На Q–Q-графике видим выраженные отклонения от прямой на обоих концах и ступени в середине. "
    "Нормальная кривая плохо передаёт одновременно узкий центральный пик и крайние значения"
)

#@title Независимость и оценки
estimate = interval_table.loc["moisture"]
show_table(estimate.to_frame("estimate"))
show_md(f"Получили среднее **{fmt(estimate['mean'],5)}**. Приближённый 95% t-интервал: **[{fmt(estimate['t_low'],5)}; {fmt(estimate['t_high'],5)}]**. Bootstrap-интервал: **[{fmt(estimate['bootstrap_low'],5)}; {fmt(estimate['bootstrap_high'],5)}]**. Границы двух интервалов близки.")

#@title Разведочный анализ и распределение
x = data["acidity"].to_numpy()
q1, q3 = np.quantile(x, [0.25, 0.75])
iqr = q3 - q1
flagged = (x < q1 - 1.5 * iqr) | (x > q3 + 1.5 * iqr)
summary = pd.Series({"n": len(x), "mean": x.mean(), "median": np.median(x), "variance_ddof1": x.var(ddof=1), "std_ddof1": x.std(ddof=1), "min": x.min(), "q1": q1, "q3": q3, "max": x.max(), "skewness": stats.skew(x, bias=False), "IQR_flags": flagged.sum()})
print(summary.round(5).to_string())

fig, axes = plt.subplots(2, 2, figsize=(14, 8))
axes[0, 0].hist(x, bins="fd", density=True, color=PALETTE["accent"], alpha=0.75)
grid = np.linspace(x.min(), x.max(), 300)
axes[0, 0].plot(grid, stats.norm.pdf(grid, x.mean(), x.std(ddof=1)), color=PALETTE["main"], label="Нормальное приближение")
axes[0, 0].set(xlabel=labels["acidity"], ylabel="Плотность", title="Гистограмма")
axes[0, 0].legend()
axes[0, 1].step(np.sort(x), np.arange(1, len(x)+1) / len(x), where="post", color=PALETTE["main"])
axes[0, 1].set(xlabel=labels["acidity"], ylabel="F_n(x)", title="Эмпирическая функция распределения", ylim=(0, 1.02))
axes[1, 0].boxplot(x, vert=False, medianprops={"color": PALETTE["accent"], "linewidth": 2}, flierprops={"markeredgecolor": PALETTE["ink"], "markersize": 4})
axes[1, 0].set(xlabel=labels["acidity"], title="Ящик с усами (выбросы оставлены)", yticks=[])
stats.probplot(x, dist="norm", plot=axes[1, 1])
axes[1, 1].get_lines()[0].set(color=PALETTE["accent"], markersize=3)
axes[1, 1].get_lines()[1].set(color=PALETTE["main"])
axes[1, 1].set(xlabel="Теоретические квантили", ylabel="Упорядоченные значения")
axes[1, 1].set_title("Q-Q график против нормального")
fig.suptitle(f"Кислотность: Arabica, n={len(x)}")
fig.tight_layout()
plt.show()

show_md(
    f"**Вывод по графикам кислотности**\n\n"
    f"Оценки сосредоточены около 7,5 балла: среднее равно **{fmt(x.mean(),2)}**, "
    f"медиана — **{fmt(np.median(x),2)}**. Центральная половина наблюдений лежит между "
    f"{fmt(q1,2)} и {fmt(q3,2)} балла. На CDF наиболее быстрый рост видим в этой области, "
    f"а ступени соответствуют повторяющимся оценкам\n\n"
    f"Асимметрия небольшая: **{fmt(stats.skew(x,bias=False),2)}**. Гистограмма в целом похожа на колокол, "
    f"но boxplot отмечает {int(flagged.sum())} крайних наблюдений с обеих сторон. "
    "На Q–Q-графике центральные точки близки к прямой, а на концах отклоняются от неё. "
    "Поэтому внешнее сходство с нормальной кривой ещё не означает согласия всего распределения с этой моделью"
)

#@title Независимость и оценки
estimate = interval_table.loc["acidity"]
show_table(estimate.to_frame("estimate"))
show_md(f"Получили среднее **{fmt(estimate['mean'],5)}**. Приближённый 95% t-интервал: **[{fmt(estimate['t_low'],5)}; {fmt(estimate['t_high'],5)}]**. Bootstrap-интервал: **[{fmt(estimate['bootstrap_low'],5)}; {fmt(estimate['bootstrap_high'],5)}]**. Границы двух интервалов близки.")

#@title Разведочный анализ и распределение
x = data["body"].to_numpy()
q1, q3 = np.quantile(x, [0.25, 0.75])
iqr = q3 - q1
flagged = (x < q1 - 1.5 * iqr) | (x > q3 + 1.5 * iqr)
summary = pd.Series({"n": len(x), "mean": x.mean(), "median": np.median(x), "variance_ddof1": x.var(ddof=1), "std_ddof1": x.std(ddof=1), "min": x.min(), "q1": q1, "q3": q3, "max": x.max(), "skewness": stats.skew(x, bias=False), "IQR_flags": flagged.sum()})
print(summary.round(5).to_string())

fig, axes = plt.subplots(2, 2, figsize=(14, 8))
axes[0, 0].hist(x, bins="fd", density=True, color=PALETTE["accent"], alpha=0.75)
grid = np.linspace(x.min(), x.max(), 300)
axes[0, 0].plot(grid, stats.norm.pdf(grid, x.mean(), x.std(ddof=1)), color=PALETTE["main"], label="Нормальное приближение")
axes[0, 0].set(xlabel=labels["body"], ylabel="Плотность", title="Гистограмма")
axes[0, 0].legend()
axes[0, 1].step(np.sort(x), np.arange(1, len(x)+1) / len(x), where="post", color=PALETTE["main"])
axes[0, 1].set(xlabel=labels["body"], ylabel="F_n(x)", title="Эмпирическая функция распределения", ylim=(0, 1.02))
axes[1, 0].boxplot(x, vert=False, medianprops={"color": PALETTE["accent"], "linewidth": 2}, flierprops={"markeredgecolor": PALETTE["ink"], "markersize": 4})
axes[1, 0].set(xlabel=labels["body"], title="Ящик с усами (выбросы оставлены)", yticks=[])
stats.probplot(x, dist="norm", plot=axes[1, 1])
axes[1, 1].get_lines()[0].set(color=PALETTE["accent"], markersize=3)
axes[1, 1].get_lines()[1].set(color=PALETTE["main"])
axes[1, 1].set(xlabel="Теоретические квантили", ylabel="Упорядоченные значения")
axes[1, 1].set_title("Q-Q график против нормального")
fig.suptitle(f"Полнота напитка: Arabica, n={len(x)}")
fig.tight_layout()
plt.show()

show_md(
    f"**Вывод по графикам полноты**\n\n"
    f"Среднее и медиана близки: **{fmt(x.mean(),2)}** и **{fmt(np.median(x),2)}** балла. "
    f"Центральная половина наблюдений лежит между {fmt(q1,2)} и {fmt(q3,2)} балла. "
    f"Разброс меньше, чем у кислотности: стандартные отклонения равны "
    f"{fmt(x.std(ddof=1),3)} и {fmt(data['acidity'].std(ddof=1),3)} балла соответственно\n\n"
    f"Распределение почти симметрично: асимметрия равна **{fmt(stats.skew(x,bias=False),2)}**. "
    "На гистограмме и CDF видим заметные скопления на отдельных оценках. "
    f"Boxplot отмечает {int(flagged.sum())} крайних наблюдения, а концы Q–Q-графика расходятся с прямой. "
    "Нормальная кривая передаёт общий вид середины, но хуже описывает крайние значения и частоты отдельных баллов"
)

#@title Независимость и оценки
estimate = interval_table.loc["body"]
show_table(estimate.to_frame("estimate"))
show_md(f"Получили среднее **{fmt(estimate['mean'],5)}**. Приближённый 95% t-интервал: **[{fmt(estimate['t_low'],5)}; {fmt(estimate['t_high'],5)}]**. Bootstrap-интервал: **[{fmt(estimate['bootstrap_low'],5)}; {fmt(estimate['bootstrap_high'],5)}]**. Границы двух интервалов близки.")

#@title 4.4 Нормализация и стандартизация
source_features = data[columns].copy()
scaling_parameters = pd.DataFrame({
    "min": source_features.min(),
    "max": source_features.max(),
    "mean": source_features.mean(),
    "std_ddof1": source_features.std(ddof=1),
})
scaling_parameters["range"] = scaling_parameters["max"] - scaling_parameters["min"]
assert scaling_parameters["range"].gt(0).all()
assert scaling_parameters["std_ddof1"].gt(0).all()
data_minmax = (source_features - scaling_parameters["min"]) / scaling_parameters["range"]
data_zscore = (source_features - scaling_parameters["mean"]) / scaling_parameters["std_ddof1"]

scaling_checks = pd.DataFrame({
    "minmax_min": data_minmax.min(),
    "minmax_max": data_minmax.max(),
    "z_mean": data_zscore.mean(),
    "z_std_ddof1": data_zscore.std(ddof=1),
    "skew_original": source_features.skew(),
    "skew_minmax": data_minmax.skew(),
    "skew_zscore": data_zscore.skew(),
})
print("Параметры преобразований (исходные шкалы):")
print(scaling_parameters.round(6).to_string())
print("\nПроверка диапазона, стандартизации и сохранения асимметрии:")
print(scaling_checks.round(6).to_string())
assert np.allclose(data_minmax.min(), 0)
assert np.allclose(data_minmax.max(), 1)
assert np.allclose(data_zscore.mean(), 0, atol=1e-12)
assert np.allclose(data_zscore.std(ddof=1), 1)
assert np.allclose(source_features.skew(), data_zscore.skew())
assert np.allclose(source_features.corr(), data_minmax.corr())
assert np.allclose(source_features.corr(), data_zscore.corr())
assert data[columns].equals(source_features)

fig, axes = plt.subplots(1, 2, figsize=(14, 4.5))
for axis, transformed, title, ylabel in [
    (axes[0], data_minmax, "Min-max нормализация", "Значение в [0, 1]"),
    (axes[1], data_zscore, "Z-стандартизация (ddof=1)", "Стандартизованное значение"),
]:
    axis.boxplot([transformed[column].to_numpy() for column in columns])
    axis.set_xticks(np.arange(1, len(columns) + 1))
    axis.set_xticklabels(columns)
    axis.set(title=title, ylabel=ylabel)
axes[1].axhline(0, color=PALETTE["main"], linestyle="--", linewidth=1)
fig.suptitle(f"Сравнение масштабирования: Arabica, n={len(data)}, выбросы по IQR оставлены")
fig.tight_layout()
plt.show()

show_md(
    "**Вывод по графикам преобразований**\n\n"
    "После min–max все показатели лежат в диапазоне от 0 до 1. При этом коробка влажности выглядит узкой: "
    "её границы растягивают редкие крайние значения, а основная масса остаётся около 10–12%\n\n"
    "После z-score средние равны 0, а стандартные отклонения - 1, поэтому сравниваем положение наблюдений "
    "относительно среднего каждого показателя. Крайние значения и различия формы сохранились: "
    "преобразование шкалы не сделало распределения нормальными и не устранило зависимость показателей. "
    "Для объяснения свойств кофе дальше используем исходные единицы"
)

#@title 5.3 Нормальный закон как модель для сравнения
distribution_tables = {}
distribution_math_rows = []
fig, axes = plt.subplots(1, 3, figsize=(14, 4.5))
for axis, column in zip(axes, columns):
    counts = data[column].value_counts().sort_index()
    support = counts.index.to_numpy(dtype=float)
    probabilities = counts.to_numpy() / len(data)
    distribution_table = pd.DataFrame({
        "value": support, "count": counts.to_numpy(),
        "probability": probabilities, "CDF": probabilities.cumsum(),
    })
    distribution_tables[column] = distribution_table
    mean = np.sum(support * probabilities)
    variance = np.sum((support - mean)**2 * probabilities)
    skewness = np.sum((support - mean)**3 * probabilities) / variance**1.5
    excess = np.sum((support - mean)**4 * probabilities) / variance**2 - 3
    modes = support[counts.to_numpy() == counts.max()]
    distribution_math_rows.append({
        "variable": column, "N": len(data), "unique_values": len(support),
        "support_min": support.min(), "support_max": support.max(),
        "modes": ", ".join(f"{value:g}" for value in modes),
        "mode_probability": probabilities.max(),
        "mean_emp": mean, "variance_emp_ddof0": variance,
        "skewness_emp": skewness, "excess_emp": excess,
    })
    assert counts.sum() == len(data)
    assert np.isclose(probabilities.sum(), 1)
    assert np.all(np.diff(distribution_table["CDF"]) >= 0)
    assert np.isclose(distribution_table["CDF"].iloc[-1], 1)
    assert np.isclose(mean, data[column].mean())
    assert np.isclose(variance, data[column].var(ddof=0))
    axis.vlines(support, 0, probabilities, color=PALETTE["main"], alpha=.8)
    axis.scatter(support, probabilities, color=PALETTE["accent"], s=18)
    axis.set(title=f"{NAMES[column]}: функция вероятностей", xlabel="Значение (исходная шкала)", ylabel="Вероятность")
    print(f"\n{column}: пять наиболее частых значений")
    print(distribution_table.sort_values(["count", "value"], ascending=[False, True]).head(5).round(6).to_string(index=False))
fig.suptitle(f"Равновероятный выбор строки из {len(data)} записей Arabica")
fig.tight_layout()
plt.show()
distribution_math_table = pd.DataFrame(distribution_math_rows).set_index("variable")
print("\nМатематические характеристики эмпирических распределений:")
print(distribution_math_table.round(6).to_string())

show_md(
    f"**Вывод по графикам вероятностей**\n\n"
    f"Самый заметный пик у влажности: значение 0,11 встречается в "
    f"**{fmt(100*distribution_math_table.loc['moisture','mode_probability'],1)}%** строк. "
    f"Для кислотности и полноты наиболее частая оценка - 7,5 балла, её доли равны "
    f"{fmt(100*distribution_math_table.loc['acidity','mode_probability'],1)}% и "
    f"{fmt(100*distribution_math_table.loc['body','mode_probability'],1)}%. "
    "Высота каждой линии показывает вероятность записанного значения при равновероятном выборе строки\n\n"
    f"У кислотности и полноты асимметрия близка к нулю, но эмпирический эксцесс положителен: "
    f"{fmt(distribution_math_table.loc['acidity','excess_emp'],2)} и "
    f"{fmt(distribution_math_table.loc['body','excess_emp'],2)}. "
    f"У влажности он заметно больше: {fmt(distribution_math_table.loc['moisture','excess_emp'],2)}. "
    "Это согласуется с заметной ролью крайних значений на предыдущих графиках\n\n"
)

#@title 5.4 Подбор теоретических моделей и проверка согласия
"""Проверяем приближение: непрерывное распределение и простое округление до 0,01."""
import numpy as np
from scipy import stats, optimize

def rounded_fit(x, family, step=0.01):
    """Оцениваем параметры распределения по значениям, записанным на сетке."""
    x = np.asarray(x, dtype=float)
    support, counts = np.unique(x, return_counts=True)
    center, unit = x.mean(), x.std(ddof=0)
    dist = getattr(stats, family)
    loc0, scale0 = dist.fit(x)
    def negative_log_likelihood(theta):
        loc = center + unit * theta[0]
        scale = unit * np.exp(theta[1])
        lower = (support - step / 2 - loc) / scale
        upper = (support + step / 2 - loc) / scale
        probability = np.where(lower >= 0,
                               dist.sf(lower) - dist.sf(upper),
                               dist.cdf(upper) - dist.cdf(lower))
        return -np.dot(counts, np.log(np.maximum(probability, np.finfo(float).tiny)))
    start = [(loc0-center)/unit, np.log(scale0/unit)]
    result = optimize.minimize(negative_log_likelihood, start, method='L-BFGS-B',
                               bounds=[(-10,10),(-6,4)], options={'ftol':1e-10})
    if not result.success:
        result = optimize.minimize(negative_log_likelihood, start, method='Powell',
                                   bounds=[(-10,10),(-6,4)])
    if not result.success or not np.isfinite(result.fun):
        raise RuntimeError('Distribution fit did not converge: ' + result.message)
    loc = center + unit * result.x[0]
    scale = unit * np.exp(result.x[1])
    return float(loc), float(scale), float(-result.fun)

def rounded_cdf_distance(x, family, loc, scale, step=0.01):
    support, counts = np.unique(x, return_counts=True)
    right = np.cumsum(counts) / len(x)
    left = right - counts / len(x)
    dist = getattr(stats, family)
    model_right = dist.cdf((support + step/2 - loc)/scale)
    model_left = dist.cdf((support - step/2 - loc)/scale)
    return max(np.max(np.abs(right-model_right)), np.max(np.abs(left-model_left)))


fitted_model_rows = []
fitted_parameters = {}
model_rng = np.random.default_rng(SEED + 610)
model_repetitions = 499
for column in columns:
    values = data[column].to_numpy()
    for family in ["norm", "logistic"]:
        loc, scale, log_likelihood = rounded_fit(values, family)
        observed_D = rounded_cdf_distance(values, family, loc, scale)
        extreme = 0
        for _ in range(model_repetitions):
            latent = (model_rng.normal(loc, scale, len(values)) if family == "norm" else
                      model_rng.logistic(loc, scale, len(values)))
            simulated = np.round(latent, 2)
            sim_loc, sim_scale, _ = rounded_fit(simulated, family)
            extreme += rounded_cdf_distance(simulated, family, sim_loc, sim_scale) >= observed_D - 1e-12
        p_mc = (extreme+1)/(model_repetitions+1)
        lower, upper = (0,1) if column == "moisture" else (0,10)
        dist = getattr(stats, family)
        fitted_parameters[(column,family)] = (loc,scale)
        fitted_model_rows.append({"variable":column, "family":family, "loc":loc, "scale":scale,
            "log_likelihood":log_likelihood, "AIC":4-2*log_likelihood,
            "CDF_distance":observed_D, "p_MC":p_mc,
            "latent_mass_outside_physical_scale":dist.cdf((lower-loc)/scale)+dist.sf((upper-loc)/scale)})
fitted_model_table = pd.DataFrame(fitted_model_rows)
fitted_model_table["p_Holm_6"] = holm_adjust(fitted_model_table["p_MC"])
fitted_model_table["reject_model"] = fitted_model_table["p_Holm_6"] < ALPHA
show_table(fitted_model_table.set_index(["variable","family"]), 6)
fig, axes = plt.subplots(1,3,figsize=(14, 4.5))
for axis, column in zip(axes, columns):
    values = data[column].to_numpy()
    grid = np.round(np.arange(values.min()-.01, values.max()+.015, .01),2)
    axis.step(grid, np.searchsorted(np.sort(values),grid,side="right")/len(values), where="post", color=PALETTE["ink"],label="Эмпирическая ФР")
    for family in ["norm","logistic"]:
        loc,scale = fitted_parameters[(column,family)]
        axis.step(grid,getattr(stats,family).cdf((grid+.005-loc)/scale),where="post",label={"norm": "нормальное", "logistic": "логистическое"}.get(family, family))
    axis.set(title=NAMES[column],xlabel=labels[column],ylabel="F(x)")
    axis.legend()
fig.tight_layout()
plt.show()
rejected_count = int(fitted_model_table["reject_model"].sum())
show_md(
    "**Вывод по сравнению моделей**\n\n"
    "Для влажности модельные CDF заметно расходятся с эмпирической около центрального скачка. "
    "Для кислотности и полноты общий ход кривых похож, но отдельные скачки наблюдаемых частот модели передают хуже. "
    "Внешней близости кривых недостаточно, чтобы принять модель\n\n"
    "В таблице логистическая модель имеет меньший AIC для каждого показателя, однако это только выбор "
    "между двумя предложенными приближениями. "
    f"После поправки Холма отвергаем **{rejected_count} из {len(fitted_model_table)}** сочетаний "
    "семейства и простого округления до 0,01\n\n"
    "Не делаем из этого вывод, что скрытые величины заведомо не могут иметь нормальное или логистическое распределение "
    "при другом механизме получения записанных оценок. Для экспериментов с повторными выборками "
    "оставляем эмпирическое распределение"
)

#@title 6. Связи между показателями
pairs = [("moisture", "acidity"), ("moisture", "body"), ("acidity", "body")]
correlation_rows = []
fig, axes = plt.subplots(1, 3, figsize=(14, 4.5))
for axis, (first, second) in zip(axes, pairs):
    pearson = stats.pearsonr(data[first], data[second])
    spearman = stats.spearmanr(data[first], data[second])
    correlation_rows.append({"pair": first + " / " + second, "pearson_r": pearson.statistic, "pearson_p": pearson.pvalue, "spearman_rho": spearman.statistic, "spearman_p": spearman.pvalue})
    axis.scatter(data[first], data[second], s=16, alpha=0.22, color=PALETTE["accent"])
    axis.set(xlabel=labels[first], ylabel=labels[second], title=f"ρ Спирмена = {spearman.statistic:.3f}")
fig.tight_layout()
plt.show()
correlation_table = pd.DataFrame(correlation_rows).set_index("pair")
correlation_table["pearson_p_holm"] = holm_adjust(correlation_table["pearson_p"])
correlation_table["spearman_p_holm"] = holm_adjust(correlation_table["spearman_p"])
print(correlation_table.to_string(float_format=lambda value: f"{value:.5g}"))

show_md(
    "**Вывод по диаграммам рассеяния**\n\n"
    f"Самая заметная связь - между кислотностью и полнотой: облако точек вытянуто вверх, "
    f"коэффициент Спирмена равен **{fmt(correlation_table.loc['acidity / body','spearman_rho'],3)}**. "
    "Образцы с большей оценкой кислотности обычно получают и большую оценку полноты, "
    "но разброс точек показывает, что соответствие не однозначное\n\n"
    f"Связи влажности с кислотностью и полнотой слабее и направлены вниз: "
    f"{fmt(correlation_table.loc['moisture / acidity','spearman_rho'],3)} и "
    f"{fmt(correlation_table.loc['moisture / body','spearman_rho'],3)}. "
    "При большом числе наблюдений даже слабая связь может быть статистически значимой, "
    "поэтому учитываем величину коэффициента. Скопления точек на сетке связаны с повторяющимися записанными значениями; "
)

#@title 6.1 Зависимость и проверка гипотезы (критерии 2 и 7)
# пары оценок одного образца
criterion2_a = data["acidity"].to_numpy()
criterion2_b = data["body"].to_numpy()
criterion2_n = len(criterion2_a)

criterion2_event_a = criterion2_a >= 7.5
criterion2_event_b = criterion2_b >= 7.5
criterion2_pa = criterion2_event_a.mean()
criterion2_pb = criterion2_event_b.mean()
criterion2_pab = (criterion2_event_a & criterion2_event_b).mean()
print("Иллюстрация для эмпирического распределения:")
print(f"P(acidity >= 7.5) = {criterion2_pa:.6f}")
print(f"P(body >= 7.5) = {criterion2_pb:.6f}")
print(f"P(оба события) = {criterion2_pab:.6f}")
print(f"Произведение вероятностей = {criterion2_pa * criterion2_pb:.6f}")

# коэффициент Спирмена
criterion2_rank_a = stats.rankdata(criterion2_a, method="average")
criterion2_rank_b = stats.rankdata(criterion2_b, method="average")
criterion2_rank_a = criterion2_rank_a - criterion2_rank_a.mean()
criterion2_rank_b = criterion2_rank_b - criterion2_rank_b.mean()
criterion2_scale = np.linalg.norm(criterion2_rank_a) * np.linalg.norm(criterion2_rank_b)
criterion2_rho = np.dot(criterion2_rank_a, criterion2_rank_b) / criterion2_scale

criterion2_rng = np.random.default_rng(42)
criterion2_repetitions = 9999
criterion2_extreme = 0
for _ in range(criterion2_repetitions):
    criterion2_permuted = criterion2_rng.permutation(criterion2_rank_b)
    criterion2_rho_null = np.dot(criterion2_rank_a, criterion2_permuted) / criterion2_scale
    if abs(criterion2_rho_null) >= abs(criterion2_rho) - 1e-12:
        criterion2_extreme += 1
criterion2_p_mc = (criterion2_extreme + 1) / (criterion2_repetitions + 1)
criterion2_p_bonf = min(1.0, 3 * criterion2_p_mc)

print("\nПерестановочная проверка независимости:")
print(f"n = {criterion2_n}; Spearman rho = {criterion2_rho:.6f}")
print(f"Перестановок = {criterion2_repetitions}; не менее экстремальных = {criterion2_extreme}")
print(f"Monte Carlo p-value = {criterion2_p_mc:.6f}")
print(f"С поправкой Бонферрони на три пары = {criterion2_p_bonf:.6f}")
print("H0 независимости", "отвергается" if criterion2_p_bonf < 0.05 else "не отвергается", "при alpha = 0.05")

show_md(f"При α={ALPHA} **{'отвергаем' if criterion2_p_bonf<ALPHA else 'не отвергаем'} H0**. Получили ρs={fmt(criterion2_rho,4)}, p_MC={fmt(criterion2_p_mc,4)}; с поправкой на три пары - {fmt(criterion2_p_bonf,4)}. Вывод делаем при предположении обменности строк. Разрешение p_MC при 9999 перестановках составляет 0,0001")

#@title 6.2 Проверка по странам и фермам
def normalize_label(series):
    return series.astype("string").str.strip().str.lower().str.replace(r"\s+", " ", regex=True)

group_audit_rows = []
for field in ["country_of_origin","farm_name","owner","certification_body"]:
    labels_clean = normalize_label(data[field]).dropna()
    counts = labels_clean.value_counts()
    group_audit_rows.append({"field":field, "missing":int(data[field].isna().sum()),
        "distinct_labels":len(counts), "labels_with_multiple_rows":int(counts.gt(1).sum()),
        "largest_group":int(counts.max())})
group_audit = pd.DataFrame(group_audit_rows).set_index("field")
show_table(group_audit,0)
identity_fields = ["country_of_origin","farm_name","lot_number","grading_date"]
identity_complete = data.dropna(subset=identity_fields)
duplicate_identity_n = int(identity_complete.duplicated(identity_fields, keep=False).sum())
print("Строк с полным потенциальным ключом партии:",len(identity_complete),
      "; строк с совпадающим ключом:",duplicate_identity_n,"(кандидаты для проверки, не автоматические дубликаты)")

country_data = data.dropna(subset=["country_of_origin"]).reset_index(drop=True)
country_groups = list(country_data.groupby("country_of_origin",sort=True).indices.values())
within_country_rows = []
for country, subset in country_data.groupby("country_of_origin",sort=True):
    if len(subset)>=20:
        within_country_rows.append({"country":country,"n":len(subset),
            "rho":stats.spearmanr(subset["acidity"],subset["body"]).statistic})
within_country_table = pd.DataFrame(within_country_rows).set_index("country")
show_table(within_country_table.sort_values("n",ascending=False),4)
rank_a = stats.rankdata(country_data["acidity"].to_numpy()).astype(float)
rank_b = stats.rankdata(country_data["body"].to_numpy()).astype(float)
for positions in country_groups:
    rank_a[positions] -= rank_a[positions].mean()
    rank_b[positions] -= rank_b[positions].mean()
within_scale = np.linalg.norm(rank_a)*np.linalg.norm(rank_b)
within_rho = np.dot(rank_a,rank_b)/within_scale
within_rng = np.random.default_rng(SEED+270)
within_extreme = 0
within_repetitions = 1999
for _ in range(within_repetitions):
    permuted_b = rank_b.copy()
    for positions in country_groups:
        permuted_b[positions] = within_rng.permutation(rank_b[positions])
    within_extreme += abs(np.dot(rank_a,permuted_b)/within_scale) >= abs(within_rho)-1e-12
within_p = (within_extreme+1)/(within_repetitions+1)

farm_label = normalize_label(data["farm_name"])
country_label = normalize_label(data["country_of_origin"])
generic_pattern = r"^(?:various|varios|varias|several|multiple|mixed|unknown|none|n/?a|not specified|\d+)(?:$|\s)"
farm_usable = farm_label.notna() & country_label.notna() & farm_label.ne("") & ~farm_label.str.contains(generic_pattern,regex=True,na=False)
farm_data = data.loc[farm_usable].copy().reset_index(drop=True)
farm_data["farm_cluster"] = (country_label[farm_usable]+" | "+farm_label[farm_usable]).to_numpy()
farm_groups = list(farm_data.groupby("farm_cluster",sort=True).indices.values())
farm_values = farm_data[columns].to_numpy()
farm_rho = stats.spearmanr(farm_values[:,1],farm_values[:,2]).statistic
bootstrap_rng = np.random.default_rng(SEED+271)
farm_boot = np.empty((1000,4))
iid_farm_boot = np.empty((1000,4))
for repetition in range(1000):
    chosen = bootstrap_rng.integers(len(farm_groups),size=len(farm_groups))
    cluster_positions = np.concatenate([farm_groups[j] for j in chosen])
    iid_positions = bootstrap_rng.integers(len(farm_values),size=len(farm_values))
    for destination, positions in [(farm_boot,cluster_positions),(iid_farm_boot,iid_positions)]:
        values = farm_values[positions]
        destination[repetition,:3] = values.mean(axis=0)
        destination[repetition,3] = stats.spearmanr(values[:,1],values[:,2]).statistic
farm_cluster_ci = np.quantile(farm_boot,[.025,.975],axis=0)
farm_iid_ci = np.quantile(iid_farm_boot,[.025,.975],axis=0)
farm_bootstrap_table = pd.DataFrame({
    "point_estimate":np.r_[farm_values.mean(axis=0),farm_rho],
    "iid_low":farm_iid_ci[0],"iid_high":farm_iid_ci[1],
    "farm_cluster_low":farm_cluster_ci[0],"farm_cluster_high":farm_cluster_ci[1],
},index=columns+["acidity_body_rho"])
show_table(farm_bootstrap_table,5)
show_md(f"После центрирования рангов внутри стран получили статистику {fmt(within_rho,3)} и p_MC={fmt(within_p,4)} при n={len(country_data)}. В проверку по фермам вошли {len(farm_data)} строк и {len(farm_groups)} кластеров. Для этого поднабора получили ρs={fmt(farm_rho,3)} и 95% кластерный bootstrap-интервал [{fmt(farm_cluster_ci[0,3],3)}; {fmt(farm_cluster_ci[1,3],3)}]. Положительная связь сохраняется при обоих способах учёта группирования")

#@title 7.2 Теоретическая опора
sizes = np.array([20, 50, 100, 250, 500, 1000, 5000])
repetitions = 1000
simulation_rng = np.random.default_rng(SEED)
simulation_rows = []
fig, axes = plt.subplots(2, 3, figsize=(14, 8))
for column_index, column in enumerate(columns):
    values = data[column].to_numpy()
    mean_target = values.mean()
    variance_target = values.var(ddof=0)
    sigma = np.sqrt(variance_target)
    fourth_moment = np.mean((values - mean_target)**4)
    for size in sizes:
        samples = simulation_rng.choice(values, size=(repetitions, size), replace=True)
        estimates_by_kind = {
            "mean": samples.mean(axis=1),
            "variance": samples.var(axis=1, ddof=1),
        }
        for estimator, estimates in estimates_by_kind.items():
            target = mean_target if estimator == "mean" else variance_target
            scale = sigma if estimator == "mean" else variance_target
            epsilon = 0.25 * scale
            errors = estimates - target
            theory_mse = (variance_target / size if estimator == "mean" else
                          (fourth_moment - (size - 3) / (size - 1) * variance_target**2) / size)
            simulation_rows.append({
                "variable": column, "estimator": estimator, "n": size,
                "target": target, "bias": errors.mean(),
                "RMSE": np.sqrt(np.mean(errors**2)),
                "normalized_RMSE": np.sqrt(np.mean(errors**2)) / scale,
                "theory_normalized_RMSE": np.sqrt(theory_mse) / scale,
                "epsilon": epsilon,
                "P_error_gt_epsilon": np.mean(np.abs(errors) > epsilon),
            })
    variable_results = pd.DataFrame(simulation_rows)
    for row_index, estimator in enumerate(["mean", "variance"]):
        subset = variable_results.loc[
            variable_results["variable"].eq(column) & variable_results["estimator"].eq(estimator)
        ]
        axis = axes[row_index, column_index]
        axis.loglog(subset["n"], subset["normalized_RMSE"], "o-", label="Моделирование")
        axis.loglog(subset["n"], subset["theory_normalized_RMSE"], "--", label="Теория (н.о.р.)")
        axis.set(xlabel="Объём выборки n", ylabel="Нормированная RMSE", title=f"{NAMES[column]}: {ESTIMATOR_RU.get(estimator, estimator)}")
        axis.legend()
fig.tight_layout()
plt.show()
simulation_table = pd.DataFrame(simulation_rows)
print(f"Повторений для каждого объёма: {repetitions}")
print("Средние: нормировка на sigma_emp; дисперсии: на variance_emp.")
print(simulation_table.drop(columns=["theory_normalized_RMSE"]).round(6).to_string(index=False))

summary_rows = []
for column in columns:
    for estimator in ["mean","variance"]:
        subset = simulation_table.loc[simulation_table["variable"].eq(column)&simulation_table["estimator"].eq(estimator)]
        first,last = subset.iloc[0],subset.iloc[-1]
        summary_rows.append({"variable":column,"estimator":estimator,
            "RMSE_small_n":first["RMSE"],"RMSE_large_n":last["RMSE"],
            "P_large_error_small_n":first["P_error_gt_epsilon"],"P_large_error_large_n":last["P_error_gt_epsilon"]})
consistency_summary = pd.DataFrame(summary_rows).set_index(["variable","estimator"])
show_table(consistency_summary,6)
show_md(
    "**Вывод по графикам состоятельности**\n\n"
    "Во всех шести панелях ошибка уменьшается с ростом n, а результаты повторных выборок близки "
    "к теоретическим пунктирным кривым. Для среднего скорость снижения RMSE соответствует зависимости "
    "от обратного квадратного корня из n\n\n"
    "Для дисперсии влажности нормированная ошибка больше, чем для дисперсий дегустационных баллов. "
    "Это согласуется с большим четвёртым стандартизованным моментом влажности: её оценка дисперсии "
    "сильнее зависит от попадания крайних наблюдений в выборку\n\n"
    "Полученные кривые иллюстрируют состоятельность среднего и исправленной дисперсии в выбранной iid-модели. "
    "Выборки до 5000 наблюдений получаем с возвращением из исходного набора. "
    "Нулевая наблюдаемая частота больших ошибок означает, что в 1000 повторениях такие ошибки не встретились"
)

#@title 8.2 Почему выбран такой размер
population = data[columns].copy()
population_n = len(population)
precision_targets = pd.Series({"moisture":.003,"acidity":.04,"body":.04})
planning_variance = population.var(ddof=0)
planning_z = stats.norm.ppf(.975)
planning_a = planning_z**2 * planning_variance / precision_targets**2
sample_size_plan = pd.DataFrame({"target_half_width":precision_targets,
    "minimum_n_WR":np.ceil(planning_a).astype(int),
    "minimum_n_WOR":np.ceil(population_n*planning_a/(population_n-1+planning_a)).astype(int)})
sample_n = int(np.ceil(sample_size_plan["minimum_n_WR"].max()/100)*100)
show_table(sample_size_plan,4)
assert 1 < sample_n < population_n
sampling_designs = [("with replacement", True), ("without replacement", False)]
sampling_rng = np.random.default_rng(SEED + 500)
sample_positions = {
    label: sampling_rng.choice(population_n, size=sample_n, replace=replace)
    for label, replace in sampling_designs
}
samples_by_design = {
    label: population.iloc[positions].copy()
    for label, positions in sample_positions.items()
}
sampling_design_table = pd.DataFrame([
    {"design": label, "draws": sample_n, "unique_rows": len(np.unique(positions)),
     "repeated_draws": sample_n - len(np.unique(positions))}
    for label, positions in sample_positions.items()
]).set_index("design")
assert sampling_design_table.loc["without replacement", "unique_rows"] == sample_n
population_stats = pd.DataFrame({
    "N": population_n, "mean": population.mean(), "median": population.median(),
    "variance_N": population.var(ddof=0), "variance_N_minus_1": population.var(ddof=1),
    "std_N": population.std(ddof=0), "q1": population.quantile(.25),
    "q3": population.quantile(.75), "min": population.min(), "max": population.max(),
})
sampling_comparison_rows = []
for label, replace in sampling_designs:
    sample = samples_by_design[label]
    for column in columns:
        values = sample[column].to_numpy()
        target_mean = population_stats.loc[column, "mean"]
        target_variance = population_stats.loc[column, "variance_N"]
        s2 = values.var(ddof=1)
        variance_estimate = s2 if replace else (population_n - 1) / population_n * s2
        correction = 1 if replace else np.sqrt(1 - sample_n / population_n)
        se_estimate = np.sqrt(s2 / sample_n) * correction
        low, high = values.mean() + np.array([-1, 1]) * stats.norm.ppf(.975) * se_estimate
        sampling_comparison_rows.append({
            "design": label, "variable": column, "n": len(values),
            "mean": values.mean(), "median": np.median(values), "s2_ddof1": s2,
            "variance_estimate_for_N": variance_estimate, "std_ddof1": np.sqrt(s2),
            "q1": np.quantile(values, .25), "q3": np.quantile(values, .75),
            "min": values.min(), "max": values.max(),
            "mean_error": values.mean() - target_mean,
            "variance_error": variance_estimate - target_variance,
            "SE_estimate": se_estimate, "CI95_low": low, "CI95_high": high,
            "covers_population_mean": bool(low <= target_mean <= high),
        })
sampling_comparison = pd.DataFrame(sampling_comparison_rows)
print(f"Конечная генеральная совокупность: N={population_n}; n={sample_n} ({sample_n/population_n:.1%} N)")
print("\nСхемы отбора:")
print(sampling_design_table.to_string())
print("\nИзвестные статистики генеральной совокупности:")
print(population_stats.round(6).to_string())
print("\nСтатистики отдельных выборок и ошибки относительно генеральной совокупности:")
show_table(sampling_comparison.set_index(["design","variable"])[["mean","median","s2_ddof1","variance_estimate_for_N","std_ddof1","q1","q3","min","max"]],6)
show_table(sampling_comparison.set_index(["design","variable"])[["mean_error","variance_error","SE_estimate","CI95_low","CI95_high","covers_population_mean"]],6)

# Сравниваем два способа отбора при одинаковом объёме выборки
sampling_repetitions = 1000
sampling_simulation_rng = np.random.default_rng(SEED + 501)
population_values = population.to_numpy()
population_mean = population_values.mean(axis=0)
population_variance = population_values.var(axis=0, ddof=0)
sampling_simulation_rows = []
sampling_mean_draws = {}
for label, replace in sampling_designs:
    means = np.empty((sampling_repetitions, len(columns)))
    variance_estimates = np.empty_like(means)
    covers = np.empty_like(means, dtype=bool)
    for repetition in range(sampling_repetitions):
        positions = sampling_simulation_rng.choice(population_n, size=sample_n, replace=replace)
        selected = population_values[positions]
        means[repetition] = selected.mean(axis=0)
        s2 = selected.var(axis=0, ddof=1)
        variance_estimates[repetition] = s2 if replace else (population_n - 1)/population_n*s2
        correction = 1 if replace else np.sqrt(1 - sample_n/population_n)
        half_width = stats.norm.ppf(.975) * np.sqrt(s2/sample_n) * correction
        covers[repetition] = np.abs(means[repetition] - population_mean) <= half_width
    sampling_mean_draws[label] = means
    exact_correction = 1 if replace else np.sqrt((population_n-sample_n)/(population_n-1))
    exact_se = np.sqrt(population_variance/sample_n) * exact_correction
    for j, column in enumerate(columns):
        errors = means[:, j] - population_mean[j]
        sampling_simulation_rows.append({
            "design": label, "variable": column,
            "mean_bias": errors.mean(), "mean_RMSE": np.sqrt(np.mean(errors**2)),
            "mean_SE_theory": exact_se[j],
            "variance_bias": (variance_estimates[:, j] - population_variance[j]).mean(),
            "variance_RMSE": np.sqrt(np.mean((variance_estimates[:, j] - population_variance[j])**2)),
            "CI95_coverage": covers[:, j].mean(),
        })
sampling_simulation = pd.DataFrame(sampling_simulation_rows)
print(f"\nСравнение схем по {sampling_repetitions} повторениям (n={sample_n}):")
print(sampling_simulation.round(6).to_string(index=False))
sampling_se_ratio = np.sqrt((population_n - sample_n)/(population_n - 1))
print(f"\nТочное отношение SE без возвращения / SE с возвращением: {sampling_se_ratio:.4f}")

fig, axes = plt.subplots(2, 3, figsize=(14, 8))
colors = {"with replacement": PALETTE["accent"], "without replacement": PALETTE["main"]}
for j, column in enumerate(columns):
    values = population[column].to_numpy()
    sorted_values = np.sort(values)
    axes[0, j].step(sorted_values, np.arange(1, population_n+1)/population_n,
                    where="post", color=PALETTE["ink"], linewidth=2, label="Генеральная совокупность")
    for label, _ in sampling_designs:
        sorted_sample = np.sort(samples_by_design[label][column].to_numpy())
        axes[0, j].step(sorted_sample, np.arange(1, sample_n+1)/sample_n,
                        where="post", color=colors[label], alpha=.8, label=SAMPLING_RU[label])
    axes[0, j].set(title=NAMES[column] + ": эмпирическая ФР", xlabel=labels[column], ylabel="F(x)")
    axes[0, j].legend()
    standardized_errors = {
        label: (sampling_mean_draws[label][:, j] - values.mean())/values.std(ddof=0)
        for label, _ in sampling_designs
    }
    bins = np.histogram_bin_edges(np.concatenate(list(standardized_errors.values())), bins=25)
    for label, errors in standardized_errors.items():
        axes[1, j].hist(errors, bins=bins, density=True, alpha=.45, color=colors[label], label=SAMPLING_RU[label])
    axes[1, j].axvline(0, color=PALETTE["ink"], linestyle="--")
    axes[1, j].set(title=NAMES[column] + ": ошибка выборочного среднего", xlabel="(X̄ выборки − μ) / σ генеральной", ylabel="Плотность")
    axes[1, j].legend()
fig.suptitle(f"Генеральная N={population_n}, выборка n={sample_n}, повторений {sampling_repetitions}")
fig.tight_layout()
plt.show()

show_md(f"По расчёту точности выбрали **n={sample_n} при N={population_n}**. Отношение стандартных ошибок без возвращения и с возвращением равно {fmt(sampling_se_ratio,3)}. Таким образом, отбор без возвращения уменьшает теоретическую стандартную ошибку на {fmt(100*(1-sampling_se_ratio),1)}%. В {sampling_repetitions} повторениях интервалы покрыли средние генеральной совокупности в {fmt(100*sampling_simulation['CI95_coverage'].min(),1)}–{fmt(100*sampling_simulation['CI95_coverage'].max(),1)}% случаев. Для оценки параметров этой конечной таблицы предпочитаем отбор без возвращения")

show_md(
    "**Что показывают графики отбора**\n\n"
    "В верхнем ряду CDF обеих отдельных выборок в целом повторяют форму распределения полной таблицы, "
    "но на отдельных участках отличаются от неё. Это ожидаемый разброс при отборе части строк\n\n"
    "В нижнем ряду ошибки средних по повторным выборкам сосредоточены около нуля. "
    "У отбора без возвращения распределения ошибок уже, что согласуется с меньшими RMSE в таблице "
    "и поправкой на конечную совокупность. Это преимущество по повторным отборам: "
    "одна конкретная выборка без возвращения всё равно может дать большую ошибку, чем одна выборка с возвращением"
)

#@title 9. Интерактивный генератор выборок (доп. задание)
def generate_sample(column="acidity", n=100, replace=False, seed=SEED):
    """Выборка объёма n из генеральной совокупности data[column]."""
    population = data[column].to_numpy()
    if not replace and n > len(population):
        raise ValueError(f"Без возвращения нельзя взять больше {len(population)} записей")
    rng = np.random.default_rng(seed)
    return rng.choice(population, size=n, replace=replace)

def show_generator(column="acidity", n=100, replace=False, seed=SEED):
    population = data[column].to_numpy()
    sample = generate_sample(column, n, replace, seed)
    scheme = "с возвращением" if replace else "без возвращения"

    fig, axes = plt.subplots(1, 2, figsize=(14, 4.5))
    grid = np.sort(population)
    axes[0].step(grid, np.arange(1, len(grid)+1) / len(grid), where="post",
                 color=PALETTE["ink"], linewidth=2, label="Генеральная совокупность")
    axes[0].step(np.sort(sample), np.arange(1, n+1) / n, where="post",
                 color=PALETTE["main"], label=f"Выборка, n={n}")
    axes[0].set(title="Эмпирические функции распределения", xlabel=labels[column], ylabel="F(x)")
    axes[0].legend()
    bins = np.histogram_bin_edges(population, bins="auto")
    axes[1].hist(population, bins=bins, density=True, alpha=.35, color=PALETTE["ink"], label="Генеральная совокупность")
    axes[1].hist(sample, bins=bins, density=True, alpha=.7, color=PALETTE["accent"], label="Выборка")
    axes[1].set(title="Гистограммы", xlabel=labels[column], ylabel="Плотность")
    axes[1].legend()
    fig.suptitle(f"{NAMES[column]}: выборка {scheme}, seed={seed}")
    fig.tight_layout()
    plt.show()

    table = pd.DataFrame({
        "генеральная": [population.mean(), population.var(ddof=0)],
        "выборка": [sample.mean(), sample.var(ddof=1)],
    }, index=["среднее", "дисперсия"])
    table["ошибка"] = table["выборка"] - table["генеральная"]
    show_table(table, 6)

try:
    import ipywidgets as widgets
    widgets.interact(
        show_generator,
        column=widgets.Dropdown(options=[(NAMES[c], c) for c in columns], value="acidity", description="Величина"),
        n=widgets.IntSlider(value=100, min=10, max=1000, step=10, description="n"),
        replace=widgets.Checkbox(value=False, description="с возвращением"),
        seed=widgets.IntText(value=SEED, description="seed"),
    )
except ImportError:
    show_md("ipywidgets не установлен: показываем статичный пример")

#@title 9. Интерактивный генератор выборок (доп. задание)
# пример для статичного просмотра (GitHub)
show_generator("acidity", n=100, replace=False, seed=SEED)

#@title 10. Итоги и выводы
show_md(f"Для основного анализа оставили {len(data)} записей Arabica из {len(raw)} строк исходного файла. Средняя влажность составила **{fmt(100*data['moisture'].mean(),2)}%**, кислотность - **{fmt(data['acidity'].mean(),4)}** балла, полнота - **{fmt(data['body'].mean(),4)}** балла")
show_md(f"Выявили положительную связь кислотности и полноты: **ρs={fmt(criterion2_rho,3)}**, p_MC={fmt(criterion2_p_mc,4)}. Она также сохранилась при проверках по странам и фермам")
show_md(f"Проверили нормальную и логистическую модели с допущением о простом округлении до 0,01 для трёх показателей. После поправки Холма отвергнуто **{int(fitted_model_table['reject_model'].sum())} из {len(fitted_model_table)}** моделей. Это результат для выбранного способа описания записанных значений, фактический механизм усреднения оценок по CSV не установлен. В экспериментах с выборками использовали эмпирические распределения. Показали уменьшение ошибки оценок с ростом n и сравнили два способа отбора выборки объёмом {sample_n} из {population_n} записей")
show_md("Свои выводы отношу к подготовленному набору Arabica. На результаты могут влиять исключение нулевой влажности, состав оценённых образцов и общие условия дегустации. Для интервалов использовали указанные в работе предположения о независимости и группировании наблюдений")

#@title 10.3 Участники и cat_mark_table
team_table = pd.DataFrame(TEAM)
team_table["notebook_link"] = NOTEBOOK_URL
team_table["comment"] = team_table["name"] + " / " + team_table["github"]
cat_mark_table = team_table[["cat_id","notebook_link","coeff","comment"]]
show_table(cat_mark_table)
if team_table[["cat_id","coeff"]].isna().any().any():
    show_md("Для завершения таблицы участников нам нужно заполнить `cat_id` и `coeff` в `TEAM`.")