Дополнительный датасет: рестораны TripAdvisor (HoReCa)

#@title Импорты, загрузка данных, проверка sha256 и чистка
%matplotlib inline
import hashlib, os, io, urllib.request, warnings
import numpy as np
import pandas as pd
from scipy import stats
import matplotlib.pyplot as plt
import lab_style

warnings.filterwarnings("ignore", category=FutureWarning)
warnings.filterwarnings("ignore", category=UserWarning)

SHA256 = "73fb6c804f801485fd84be8c2814181f306ae925e09030309a2a5744b2aeb1fe"
URL = "https://raw.githubusercontent.com/Ananyadevaraju/TripAdvisor---Sentiment-Analysis/master/TA_restaurants_curated.csv"
LOCAL = "TA_restaurants_curated.csv"

def load_raw_bytes():
    # Сначала локальная копия рядом с блокнотом, иначе — скачивание по URL.
    if os.path.exists(LOCAL):
        data = open(LOCAL, "rb").read()
        src = "локальная копия"
    else:
        with urllib.request.urlopen(URL) as resp:
            data = resp.read()
        src = "скачано по URL"
    digest = hashlib.sha256(data).hexdigest()
    assert digest == SHA256, f"sha256 не совпал: {digest}"
    return data, src

raw, data_src = load_raw_bytes()
df = pd.read_csv(io.BytesIO(raw))
# Проверка колонок (структура файла):
need = {"Name", "City", "Rating", "Number of Reviews"}
assert need <= set(df.columns), f"нет колонок: {need - set(df.columns)}"

# Чистка: рейтинг > 0 (−1 и пусто = нет данных) И число отзывов >= 2, обе не пусты.
raw_rows = len(df)
d = (df[(df["Rating"] > 0) & (df["Number of Reviews"] >= 2)]
        .dropna(subset=["Rating", "Number of Reviews"])
        .copy())
d["ln_reviews"] = np.log(d["Number of Reviews"])

# Единый тёмный стиль «Тёмная обжарка», акцент horeca (жёлтый).
C = lab_style.apply_dark_style("horeca")
L = lab_style.LABELS

def rus(x, nd=None):
    # Русская десятичная запятая для чисел в выводе.
    if nd is not None:
        x = round(float(x), nd)
    return str(x).replace(".", ",")

print(f"Источник данных: {data_src}, sha256 проверен.")
print(f"Строк до чистки: {raw_rows}")
print(f"Строк после чистки (рабочий набор): {len(d)}")
print(f"Городов в наборе: {d['City'].nunique()}")

#@title Таблица: размер набора до/после чистки и охват городов
summary = pd.DataFrame({
    "Показатель": ["Строк в исходном файле",
                   "Строк после чистки (рабочий набор)",
                   "Отброшено при чистке",
                   "Городов",
                   "Диапазон Rating (балл)",
                   "Диапазон Number of Reviews (шт.)"],
    "Значение": [len(df),
                 len(d),
                 len(df) - len(d),
                 d["City"].nunique(),
                 f"{d['Rating'].min():.1f} … {d['Rating'].max():.1f}",
                 f"{int(d['Number of Reviews'].min())} … {int(d['Number of Reviews'].max())}"]
})
summary

#@title Сводные статистики по обеим величинам
desc = d[["Rating", "Number of Reviews", "ln_reviews"]].describe(
    percentiles=[0.1, 0.25, 0.5, 0.75, 0.9]).round(3)
desc

#@title Rating: распределение (столбцы по 9 дискретным значениям) и ECDF
vc = d["Rating"].value_counts().sort_index()
fig, (ax1, ax2) = lab_style.styled_fig(1, 2, figsize=(12, 4.5))

ax1.bar(vc.index, vc.values, width=0.35, color=C["accent"], edgecolor=C["axis"])
ax1.set_xlabel(L["rating"]); ax1.set_ylabel("Число ресторанов (шт.)")
ax1.set_title("Распределение рейтинга (дискретное, шаг 0,5)")
ax1.set_xticks(np.arange(1.0, 5.5, 0.5))

xs = np.sort(d["Rating"].values)
ecdf = np.arange(1, len(xs) + 1) / len(xs)
ax2.step(xs, ecdf, where="post", color=C["second"])
ax2.set_xlabel(L["rating"]); ax2.set_ylabel("ECDF $F_n(x)$")
ax2.set_title("Эмпирическая функция распределения рейтинга")
plt.show()

#@title Number of Reviews: гистограмма в обычном и логарифмическом масштабе
fig, (ax1, ax2) = lab_style.styled_fig(1, 2, figsize=(12, 4.5))

ax1.hist(d["Number of Reviews"], bins=60, color=C["accent"], edgecolor=C["axis"])
ax1.set_xlabel(L["reviews"]); ax1.set_ylabel("Число ресторанов (шт.)")
ax1.set_title("Число отзывов — обычный масштаб (тяжёлый хвост)")

ax2.hist(d["Number of Reviews"], bins=np.logspace(np.log10(2), np.log10(d["Number of Reviews"].max()), 50),
         color=C["second"], edgecolor=C["axis"])
ax2.set_xscale("log")
ax2.set_xlabel(L["reviews"] + " — лог-ось"); ax2.set_ylabel("Число ресторанов (шт.)")
ax2.set_title("Число отзывов — логарифмическая ось X")
plt.show()

#@title Боксплоты и лог-преобразование: ln(reviews) ≈ нормально (гистограмма + Q–Q)
fig, axes = lab_style.styled_fig(1, 3, figsize=(14, 4.3))
ax1, ax2, ax3 = axes

bp = ax1.boxplot([d["Rating"]], patch_artist=True, tick_labels=["Rating"])
for b in bp["boxes"]:
    b.set_facecolor(C["panel"]); b.set_edgecolor(C["accent"])
ax1.set_ylabel(L["rating"]); ax1.set_title("Боксплот рейтинга")

ax2.hist(d["ln_reviews"], bins=50, density=True, color=C["second"], edgecolor=C["axis"])
mu, sd = d["ln_reviews"].mean(), d["ln_reviews"].std()
xx = np.linspace(d["ln_reviews"].min(), d["ln_reviews"].max(), 200)
ax2.plot(xx, stats.norm.pdf(xx, mu, sd), color=C["accent"], lw=2, label="норм. модель")
ax2.set_xlabel(L["log_reviews"]); ax2.set_ylabel("Плотность")
ax2.set_title("ln(число отзывов) ≈ нормально"); ax2.legend()

stats.probplot(d["ln_reviews"], dist="norm", plot=ax3)
ax3.get_lines()[0].set_color(C["second"]); ax3.get_lines()[0].set_markersize(2)
ax3.get_lines()[1].set_color(C["accent"])
ax3.set_xlabel("Теоретические квантили"); ax3.set_ylabel("Выборочные квантили")
ax3.set_title("Q–Q график ln(reviews)")
plt.show()

rmin, rmax = d['Rating'].min(), d['Rating'].max(); rmean = d['Rating'].mean()
print("Нормализация Rating (справочно):")
print(f"  min-max в [0;1]: напр. среднее {rus(rmean,2)} балла -> {rus((rmean-rmin)/(rmax-rmin),3)}")
print(f"  z-score: mean=0, sd=1 по построению (исходные mean={rus(rmean,3)}, sd={rus(d['Rating'].std(),3)})")
print(f"ln(reviews): mean={rus(mu,2)}  sd={rus(sd,2)}")

#@title Rating: среднее и 95% доверительный интервал
n = len(d); mean_r = d["Rating"].mean(); sd_r = d["Rating"].std(ddof=1)
se = sd_r / np.sqrt(n); z = 1.96
lo, hi = mean_r - z*se, mean_r + z*se
print(f"n = {n}")
print(f"Среднее Rating        = {rus(mean_r,4)} балла")
print(f"Исправл. СКО (s)      = {rus(sd_r,4)}")
print(f"Стандартная ошибка SE = {rus(se,5)}")
print(f"95% ДИ для среднего   = [{rus(lo,4)}; {rus(hi,4)}] балла")
print(f"Асимметрия (skew)     = {rus(stats.skew(d['Rating']),2)}  (левый перекос)")
print("(ДИ трактуем как неопределённость среднего при взгляде на набор как на ВЫБОРКУ из более")
print(" широкой совокупности ресторанов; для него самого как полной популяции среднее известно точно.)")

#@title Number of Reviews: сводка и доверительный интервал для среднего ln(reviews)
rv = d["Number of Reviews"]
print(f"min={int(rv.min())}  max={int(rv.max())}  медиана={int(rv.median())}  "
      f"p90={int(rv.quantile(0.9))}  среднее={rus(rv.mean(),2)} шт.")
print(f"Асимметрия Number of Reviews = {rus(stats.skew(rv),2)}  (тяжёлый правый хвост)")
print(f"Асимметрия ln(reviews)       = {rus(stats.skew(d['ln_reviews']),2)}  (почти нормально)")

n = len(d); m = d["ln_reviews"].mean(); s = d["ln_reviews"].std(ddof=1)
se = s/np.sqrt(n)
print(f"\nln(reviews): среднее = {rus(m,4)}  s = {rus(s,4)}")
print(f"95% ДИ для среднего ln(reviews) = [{rus(m-1.96*se,4)}; {rus(m+1.96*se,4)}]")
print(f"Обратный перевод exp(среднего ln) = {rus(np.exp(m),2)} шт. (геом. среднее ≈ медиана {int(rv.median())})")
print("(ДИ — неопределённость среднего при взгляде на набор как на выборку, см. пояснение выше.)")

#@title Rating: эмпирическая PMF и CDF
vc = d["Rating"].value_counts().sort_index()
pmf = vc / vc.sum(); cdf = pmf.cumsum()
fig, (ax1, ax2) = lab_style.styled_fig(1, 2, figsize=(12, 4.3))
ax1.stem(pmf.index, pmf.values, basefmt=" ")
ax1.set_xlabel(L["rating"]); ax1.set_ylabel("Вероятность $p(x)$")
ax1.set_title("Эмпирическая PMF рейтинга")
ax2.step(cdf.index, cdf.values, where="post", color=C["accent"])
ax2.set_xlabel(L["rating"]); ax2.set_ylabel("$F(x)$")
ax2.set_title("Эмпирическая CDF рейтинга")
plt.show()
for k, v in pmf.items():
    print(f"  P(Rating={rus(k,1)}) = {rus(v,4)}")

#@title Number of Reviews: проверка логнормальной модели (Q–Q по ln + K–S)
mu, sd = d["ln_reviews"].mean(), d["ln_reviews"].std(ddof=1)
# K-S согласие ln(reviews) с нормальным: стандартизуем и сравним со стандартным N(0,1).
z_ln = (d["ln_reviews"] - mu) / sd
ks = stats.kstest(z_ln, stats.norm.cdf)
fig, (ax1, ax2) = lab_style.styled_fig(1, 2, figsize=(12, 4.3))

ax1.hist(d["Number of Reviews"],
         bins=np.logspace(np.log10(2), np.log10(d["Number of Reviews"].max()), 50),
         density=True, color=C["second"], edgecolor=C["axis"], label="эмпирика")
xx = np.logspace(np.log10(2), np.log10(d["Number of Reviews"].max()), 300)
ax1.plot(xx, stats.lognorm.pdf(xx, s=sd, scale=np.exp(mu)), color=C["accent"], lw=2,
         label="логнормальная модель")
ax1.set_xscale("log")
ax1.set_xlabel(L["reviews"] + " — лог-ось"); ax1.set_ylabel("Плотность")
ax1.set_title("Логнормальная модель числа отзывов"); ax1.legend()

stats.probplot(d["ln_reviews"], dist="norm", sparams=(mu, sd), plot=ax2)
ax2.get_lines()[0].set_color(C["second"]); ax2.get_lines()[0].set_markersize(2)
ax2.get_lines()[1].set_color(C["accent"])
ax2.set_xlabel("Теоретические квантили"); ax2.set_ylabel("Выборочные квантили")
ax2.set_title("Q–Q: ln(reviews) против нормального")
plt.show()
print(f"Оценки: mu = {rus(mu,3)}, sigma = {rus(sd,3)}")
print(f"K–S статистика = {rus(ks.statistic,4)}  (D мал → форма близка к логнормальной)")

#@title Корреляции, эффект-сайз и p-value на полном наборе
r_pr = stats.pearsonr(d["Rating"], d["Number of Reviews"])
r_prln = stats.pearsonr(d["Rating"], d["ln_reviews"])
rho = stats.spearmanr(d["Rating"], d["Number of Reviews"])
R2 = rho.statistic**2
print(f"Пирсон r(Rating, Number)      = {rus(r_pr.statistic,4)}")
print(f"Пирсон r(Rating, ln reviews)  = {rus(r_prln.statistic,4)}")
print(f"Спирмен rho_s(Rating, Number) = {rus(rho.statistic,4)}")
print(f"p-value на полном наборе (n={len(d)}) = {rho.pvalue:.3e}")
print(f"Эффект-сайз: |rho_s| = {rus(abs(rho.statistic),4)} — ранговая связь ничтожна")
print(f"  (rho_s^2 = {rus(R2,4)}, т.е. {rus(R2*100,2)}%; это НЕ доля дисперсии как r^2 Пирсона, а мера крошечности ранговой связи)")

#@title Диаграмма рассеяния Rating vs ln(reviews) с усреднением по бинам
fig, ax = lab_style.styled_fig(figsize=(9, 5.2))
idx = np.random.RandomState(0).choice(len(d), size=8000, replace=False)
ax.scatter(d["ln_reviews"].values[idx], d["Rating"].values[idx],
           s=6, alpha=0.08, color=C["second"], label="рестораны (подвыборка)")
# среднее Rating по бинам ln(reviews):
bins = np.linspace(d["ln_reviews"].min(), d["ln_reviews"].max(), 14)
cat = pd.cut(d["ln_reviews"], bins)
gm = d.groupby(cat)["Rating"].mean()
centers = [(b.left + b.right) / 2 for b in gm.index]
ax.plot(centers, gm.values, "-o", color=C["accent"], lw=2, label="среднее Rating по бинам")
ax.set_xlabel(L["log_reviews"]); ax.set_ylabel(L["rating"])
ax.set_title("Rating против ln(числа отзывов): связь почти горизонтальна")
ax.legend()
plt.show()

#@title Перестановочный тест по |rho_s| (подвыборка n=2000)
rng = np.random.RandomState(42)
sub = d.sample(2000, random_state=1)
obs = abs(stats.spearmanr(sub["Rating"], sub["Number of Reviews"]).statistic)
a = sub["Rating"].values; b = sub["Number of Reviews"].values
B = 2000; cnt = 0
for _ in range(B):
    perm = rng.permutation(b)
    if abs(stats.spearmanr(a, perm).statistic) >= obs:
        cnt += 1
p_perm = (cnt + 1) / (B + 1)
print(f"Подвыборка n=2000: |rho_s| набл. = {rus(obs,4)}")
print(f"Перестановочный p (B={B}) = {rus(p_perm,4)}")
print("При n=2000 тест ещё отвергает H0, но на меньших n (300/1000) результат")
print("неустойчив — то значимо, то нет (см. таблицу ниже). Доминирует шум выборки.")

#@title Устойчивость значимости: p-value на 5 случайных подвыборках при разных n
rng = np.random.RandomState(7)
rows = []
for n in [300, 1000, 5000]:
    ps = []
    for _ in range(5):
        s = d.sample(n, random_state=rng.randint(1_000_000))
        ps.append(stats.spearmanr(s["Rating"], s["Number of Reviews"]).pvalue)
    rows.append([n, min(ps), max(ps), sum(p < 0.05 for p in ps)])
tab = pd.DataFrame(rows, columns=["n", "p_min", "p_max", "значимо из 5 (p<0,05)"])
tab["p_min"] = tab["p_min"].map(lambda x: rus(x, 4))
tab["p_max"] = tab["p_max"].map(lambda x: rus(x, 4))
print(tab.to_string(index=False))

#@title График: типичное p-value падает с ростом n (эффект фиксирован)
rng = np.random.RandomState(11)
ns = [100, 300, 1000, 3000, 10000, 30000, 100000]
med_p = []
for n in ns:
    n = min(n, len(d))
    ps = [stats.spearmanr(s["Rating"], s["Number of Reviews"]).pvalue
          for s in (d.sample(n, random_state=rng.randint(1_000_000)) for _ in range(9))]
    med_p.append(np.median(ps))
fig, ax = lab_style.styled_fig(figsize=(9, 5))
ax.plot(ns, med_p, "-o", color=C["accent"], label="медианное p (9 подвыборок)")
ax.axhline(0.05, color=C["second"], ls="--", lw=2, label="порог значимости 0,05")
ax.set_xscale("log"); ax.set_yscale("log")
ax.set_xlabel("Объём выборки n (шт. ресторанов)")
ax.set_ylabel("p-value (лог-ось)")
ax.set_title("Значимость — функция n: эффект тот же, p падает с ростом n")
ax.legend()
plt.show()

#@title chi-квадрат тест независимости на биннированных величинах
rbin = pd.cut(d["Rating"], bins=[0, 3.0, 3.5, 4.0, 4.5, 5.0])
nbin = pd.qcut(d["Number of Reviews"], q=5, duplicates="drop")
ct = pd.crosstab(rbin, nbin)
chi2, p_chi, dof, _ = stats.chi2_contingency(ct)
print(f"chi^2 = {rus(chi2,1)}  dof = {dof}  p = {p_chi:.3e}")
print("chi^2 тоже отвергает независимость — но лишь из-за огромного n;")
print("сила связи (Крамер V) мала:")
V = np.sqrt(chi2 / (len(d) * (min(ct.shape) - 1)))
print(f"  Крамер V = {rus(V,4)}  (намного меньше 1 — связь слабая, практической зависимости нет)")

#@title RMSE оценки среднего падает как 1/sqrt(n) — сверка с теорией
rng = np.random.RandomState(3)
fig, (ax1, ax2) = lab_style.styled_fig(1, 2, figsize=(12, 4.6))
for col, ax, color, name in [("Rating", ax1, C["accent"], "Rating"),
                             ("ln_reviews", ax2, C["second"], "ln(reviews)")]:
    vals = d[col].values
    true_mean = vals.mean(); sigma = vals.std(ddof=0)
    ns = [10, 30, 100, 300, 1000, 3000]
    rmse = []
    for n in ns:
        means = [rng.choice(vals, size=n, replace=True).mean() for _ in range(400)]
        rmse.append(np.sqrt(np.mean((np.array(means) - true_mean) ** 2)))
    ax.plot(ns, rmse, "o-", color=color, label="RMSE (пересэмплинг)")
    ax.plot(ns, sigma / np.sqrt(ns), "--", color=C["muted"], label="теория σ/√n")
    ax.set_xscale("log"); ax.set_yscale("log")
    ax.set_xlabel("Объём выборки n (шт.)"); ax.set_ylabel("RMSE среднего")
    ax.set_title(f"Состоятельность среднего: {name}"); ax.legend()
plt.show()
print("Эмпирический RMSE ложится на теоретическую 1/√n — оценки состоятельны.")

#@title Выборки с возвращением и без; FPC; расчёт объёма по точности
rng = np.random.RandomState(5)
vals = d["Rating"].values
N = len(vals); n = 1000
with_rep = [rng.choice(vals, size=n, replace=True).mean() for _ in range(2000)]
without_rep = [rng.choice(vals, size=n, replace=False).mean() for _ in range(2000)]
fpc = (N - n) / (N - 1)
print(f"N = {N}, n = {n}")
print(f"С возвращением:  среднее оценок = {rus(np.mean(with_rep),4)}  дисперсия = {rus(np.var(with_rep),6)}")
print(f"Без возвращения: среднее оценок = {rus(np.mean(without_rep),4)}  дисперсия = {rus(np.var(without_rep),6)}")
print(f"Отношение дисперсий (без/с) = {rus(np.var(without_rep)/np.var(with_rep),4)}  ~ FPC = {rus(fpc,4)}")

sigma = vals.std(ddof=0); z = 1.96; E = 0.01
n_need = (z * sigma / E) ** 2
print(f"\nДля точности E={rus(E,2)} балла (95%) нужно n = {int(np.ceil(n_need))} ресторанов")

fig, ax = lab_style.styled_fig(figsize=(9, 5))
ax.hist(with_rep, bins=40, alpha=0.7, color=C["accent"], label="с возвращением")
ax.hist(without_rep, bins=40, alpha=0.6, color=C["second"], label="без возвращения")
ax.axvline(vals.mean(), color=C["muted"], ls="--", lw=2, label="среднее генеральной")
ax.set_xlabel("Оценка среднего " + L["rating"]); ax.set_ylabel("Частота (шт.)")
ax.set_title("Распределение оценок среднего: два способа отбора (n=1000)")
ax.legend()
plt.show()

#@title Статичный пример: ECDF выборки против генеральной совокупности
rng = np.random.RandomState(0)
vals = np.sort(d["Rating"].values)
pop_ecdf = np.arange(1, len(vals) + 1) / len(vals)
fig, ax = lab_style.styled_fig(figsize=(9, 5))
ax.step(vals, pop_ecdf, where="post", color=C["muted"], lw=2.5,
        label=f"генеральная (N={len(vals)})")
for n, col in [(50, C["accent"]), (500, C["second"])]:
    s = np.sort(rng.choice(d["Rating"].values, size=n, replace=False))
    e = np.arange(1, n + 1) / n
    ax.step(s, e, where="post", color=col, alpha=0.9, label=f"выборка n={n}")
ax.set_xlabel(L["rating"]); ax.set_ylabel("ECDF $F(x)$")
ax.set_title("Выборочная ECDF приближает генеральную (интерактив — в marimo)")
ax.legend()
plt.show()