{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "01137ceb",
   "metadata": {},
   "source": [
    "# Wielkość efektu i moc testów\n",
    "\n",
    "Na tym wykładzie przyjrzymy się dwóm pojęciom:\n",
    "\n",
    "- **wielkość efektu**: jak duża jest różnica / zależność (w jednostkach oryginalnych lub w postaci standaryzowanej),\n",
    "- **moc testu**: jak duża jest szansa, że przy danym planie badania wykryjemy efekt o interesującej nas wielkości.\n",
    "\n",
    "Dla intuicji będziemy używać prostego przykładu z pogranicza sztucznej inteligencji i kognitywistyki: **etykieta „wygenerowane przez AI” a zaufanie do tekstu**. Wynikiem będzie ocena zaufania na skali 0-10.\n",
    "\n",
    "W kilku miejscach pojawią się też **widżety** (interaktywne suwaki), które pozwalają zobaczyć „na żywo”, jak na moc wpływają: liczebność, poziom istotności, zmienność i wielkość efektu.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "37293bd3",
   "metadata": {},
   "source": [
    "## Cele\n",
    "\n",
    "Po zajęciach:\n",
    "\n",
    "- Rozumiesz pojęcia: poziom istotności $\\alpha$, błąd II rodzaju $\\beta$, moc $1-\\beta$.\n",
    "- Potrafisz policzyć moc w prostym modelu (znana wariancja → rozkład normalny).\n",
    "- Potrafisz policzyć moc testu t (niecentralny rozkład t).\n",
    "- Widzisz związek między wielkością efektu, wariancją, liczebnością i mocą.\n",
    "- Umiesz planować liczebność próby dla zadanego celu (np. moc 0.80)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "66dcb94d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from scipy import integrate\n",
    "from scipy import stats\n",
    "\n",
    "from IPython.display import display\n",
    "\n",
    "SEED = 20250116\n",
    "rng = np.random.default_rng(SEED)\n",
    "\n",
    "sns.set_theme(style=\"whitegrid\")\n",
    "plt.rcParams[\"figure.figsize\"] = (9, 4)\n",
    "plt.rcParams[\"figure.dpi\"] = 120\n",
    "\n",
    "try:\n",
    "    import ipywidgets as widgets\n",
    "\n",
    "    WIDGETY_DOSTEPNE = True\n",
    "except ImportError:\n",
    "    widgets = None\n",
    "    WIDGETY_DOSTEPNE = False"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "94161952",
   "metadata": {},
   "source": [
    "## Znana wariancja\n",
    "\n",
    "Moc testu to prawdopodobieństwo odrzucenia hipotezy zerowej, kiedy hipoteza zerowa jest fałszywa. Jest to więc miara czułości naszego testu.\n",
    "Moc to $1 - \\beta$, gdzie $\\beta$ jest prawdopodobieństwem popełnienia błędu II rodzaju (czyli nieodrzucenia hipotezy zerowej, kiedy w populacji istnieje efekt).\n",
    "\n",
    "Tak jak do kontrolowania $\\alpha$ (czyli prawdopodobieństwa odrzucenia hipotezy zerowej, gdy jest ona prawdziwa) musimy znać rozkład statystyki z próby przy założeniu, że $H_0$ jest prawdziwa, podobnie, żeby oszacować moc i $\\beta$, musimy opisać rozkład statystyki z próby przy założeniu, że $H_0$ jest fałszywa.\n",
    "W praktyce oznacza to, że musimy przyjąć jakąś **punktową hipotezę alternatywną**.\n",
    "\n",
    "Skoro zazwyczaj hipoteza zerowa mówi o braku różnic (przy dwóch grupach będzie to $\\mu_1 - \\mu_2 = 0$), to hipoteza alternatywna może mieć postać $\\mu_1 - \\mu_2 = d$, gdzie $d \\neq 0$.\n",
    "\n",
    "**Uwaga (wielkość efektu).** Różnicę między średnimi $d$ często wyraża się również w formie standaryzowanej (jako stosunek do odchylenia standardowego w populacji):\n",
    "\n",
    "$$\n",
    " d_{\\text{Cohena}} = \\frac{\\mu_1 - \\mu_2}{\\sigma}.\n",
    "$$\n",
    "\n",
    "### Przykład: etykieta AI a zaufanie do tekstu\n",
    "\n",
    "Wyobraźmy sobie proste badanie. Dwie grupy studentów czytają ten sam krótki tekst poradnikowy. Różni się tylko etykieta nad tekstem:\n",
    "\n",
    "- **Grupa „człowiek”** widzi informację: „tekst napisany przez człowieka”.\n",
    "- **Grupa „AI”** widzi informację: „tekst wygenerowany przez AI”.\n",
    "\n",
    "Po lekturze każda osoba ocenia, jak bardzo ufa tekstowi, na skali od 0 do 10 punktów.\n",
    "\n",
    "Interesuje nas, czy etykieta „AI” obniża zaufanie do tekstu. Dla wygody ustawiamy różnicę jako:\n",
    "\n",
    "$$\n",
    "\\Delta = \\mu_{\\text{człowiek}} - \\mu_{\\text{AI}}.\n",
    "$$\n",
    "\n",
    "Dla prostoty zaczniemy od wersji „międzyosobniczej”: porównujemy dwie niezależne grupy osób, z których każda widzi tylko jedną etykietę.\n",
    "\n",
    "- Pytanie: czy średnie zaufanie w grupie „człowiek” jest **większe** niż w grupie „AI”?\n",
    "- Hipotezy (jednostronnie): $H_0: \\mu_1-\\mu_2=0$ vs. $H_A: \\mu_1-\\mu_2>0$.\n",
    "- Ustalamy poziom istotności: $\\alpha=0.025$ (to dla dla uproszczenia - robimy test jednostronny ale z alfą tak jakbyśmy robili dwustronny)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ba05c1cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Liczebność w każdej grupie (zakładamy dwie próby równoliczne)\n",
    "n = 20\n",
    "\n",
    "# Poziom istotności. Zostajemy przy teście jednostronnym (górny ogon).\n",
    "alpha = 0.025\n",
    "\n",
    "# Dla intuicji: ocena zaufania do tekstu na skali 0-10.\n",
    "mean_human_label = 6.0\n",
    "mean_ai_label = 5.0\n",
    "\n",
    "# Odchylenie standardowe w obu populacjach/grupach (pkt) – w tej części traktujemy je jako znane.\n",
    "sigma = 2.0\n",
    "\n",
    "# Różnica średnich, którą chcemy wykryć: człowiek minus AI (pkt)\n",
    "diff = mean_human_label - mean_ai_label\n",
    "\n",
    "cohen_d = diff / sigma\n",
    "print(\"Różnica średnich (pkt):\", float(diff))\n",
    "print(\"Standaryzowana wielkość efektu (d Cohena):\", float(cohen_d))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6637766c",
   "metadata": {},
   "source": [
    "### Najmniejszy interesujący efekt (SESOI)\n",
    "\n",
    "W planowaniu mocy potrzebujemy przyjąć konkretną alternatywę: „jak duży efekt chcemy wykryć?”\n",
    "\n",
    "W praktyce często warto zacząć od **SESOI** (ang. *smallest effect size of interest*), czyli najmniejszego efektu, który byłby dla nas merytorycznie istotny.\n",
    "\n",
    "W naszym przykładzie bardzo mała różnica, np. `0.1` punktu na skali zaufania 0-10, może być statystycznie wykrywalna w ogromnej próbie, ale trudno uznać ją za praktycznie ważną. Dlatego dobrze jest myśleć o efekcie w jednostkach naturalnych (punkty skali) i dopiero pomocniczo w jednostkach standaryzowanych (np. $d$ Cohena).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b5fc2b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "sesoi_points = 0.5\n",
    "sesoi_d = sesoi_points / sigma\n",
    "\n",
    "print(f\"SESOI: {sesoi_points:.1f} pkt\")\n",
    "print(f\"To odpowiada d Cohena ≈ {sesoi_d:.3f} przy σ = {sigma:.1f} pkt\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95eb02d1",
   "metadata": {},
   "source": [
    "#### Uwaga: test jednostronny i wybór $\\alpha$\n",
    "\n",
    "W tym przykładzie rozważamy hipotezę kierunkową („etykieta AI obniża zaufanie”), więc test jest jednostronny.\n",
    "W literaturze spotkasz oba podejścia:\n",
    "\n",
    "- test jednostronny z $\\alpha = 0.05$ (mniej konserwatywny),\n",
    "- test dwustronny z $\\alpha = 0.05$ (bez zakładania kierunku).\n",
    "\n",
    "My ustawiamy $\\alpha = 0.025$ w teście jednostronnym.\n",
    "To wybór **konserwatywny**: trudniej wtedy o odrzucenie $H_0$, a moc jest mniejsza niż przy $\\alpha = 0.05$ jednostronnie.\n",
    "\n",
    "Dodatkowa korzyść dydaktyczna: przy tej wartości $\\alpha$ dostajemy tę samą wartość krytyczną\n",
    "co w 95% przedziale ufności (czyli w przybliżeniu $z \\approx 1.96$).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e384698",
   "metadata": {},
   "source": [
    "### Błąd standardowy różnicy średnich\n",
    "\n",
    "Na początek obliczmy błąd standardowy różnicy dwóch średnich.\n",
    "\n",
    "Na wariancję różnicy dwóch średnich składają się wariancje generowane przez obie średnie:\n",
    "\n",
    "$$\n",
    "\\frac{\\sigma_1^2}{n} + \\frac{\\sigma_2^2}{n}\n",
    "$$\n",
    "\n",
    "Stąd błąd standardowy (czyli pierwiastek z wariancji) w naszym przypadku to $\\sigma\\sqrt{2/n}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22bf8f2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "standard_error = sigma * np.sqrt(2 / n)\n",
    "print(\"Błąd standardowy różnicy średnich:\", standard_error)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f88908d",
   "metadata": {},
   "source": [
    "### Rozkład statystyki pod $H_0$ oraz pod (punktową) $H_A$\n",
    "\n",
    "Znając błąd standardowy, możemy nakreślić rozkład różnicy średnich z próby przy założeniu hipotezy zerowej oraz – dla porównania – rozkład tej samej statystyki przy założeniu, że prawdziwa jest punktowa hipoteza alternatywna $\\mu_1-\\mu_2=\\text{diff}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "417b42d8",
   "metadata": {},
   "outputs": [],
   "source": [
    "x_min = min(-4 * standard_error, diff - 4 * standard_error)\n",
    "x_max = max(4 * standard_error, diff + 4 * standard_error)\n",
    "x = np.linspace(x_min, x_max, 500)\n",
    "\n",
    "pdf_h0 = stats.norm.pdf(x, loc=0.0, scale=standard_error)\n",
    "pdf_ha = stats.norm.pdf(x, loc=diff, scale=standard_error)\n",
    "\n",
    "plt.plot(x, pdf_h0, linewidth=3, label=r\"$H_0: N(0,\\,se)$\")\n",
    "plt.plot(x, pdf_ha, linewidth=2, label=rf\"$H_A: N({diff:.0f},\\,se)$\")\n",
    "\n",
    "plt.title(\"Rozkład różnicy średnich: $H_0$ vs. punktowa $H_A$\")\n",
    "plt.xlabel(\"różnica średnich z próby (pkt)\")\n",
    "plt.ylabel(\"gęstość\")\n",
    "plt.legend()\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c7b2301",
   "metadata": {},
   "source": [
    "### Region krytyczny\n",
    "\n",
    "Jeśli statystyka policzona z naszych danych (różnica między średnimi) przekroczy wartość krytyczną, odrzucimy $H_0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70b3bb8f",
   "metadata": {},
   "outputs": [],
   "source": [
    "critical_value = stats.norm.isf(alpha, loc=0.0, scale=standard_error)\n",
    "print(\"Wartość krytyczna:\", float(critical_value))\n",
    "\n",
    "plt.plot(x, pdf_h0, linewidth=3, label=r\"$H_0: N(0,\\,se)$\")\n",
    "plt.plot(x, pdf_ha, linewidth=2, label=rf\"$H_A: N({diff},\\,se)$\")\n",
    "plt.axvline(critical_value, linestyle=\"--\", color=\"black\", label=\"wartość krytyczna\")\n",
    "\n",
    "plt.title(\"Wartość krytyczna wyznacza region odrzucenia $H_0$\")\n",
    "plt.xlabel(\"różnica średnich z próby\")\n",
    "plt.ylabel(\"gęstość\")\n",
    "plt.legend()\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a15364ce",
   "metadata": {},
   "source": [
    "### p-wartość obliczana „ręcznie” (dla pojedynczej obserwacji statystyki)\n",
    "\n",
    "Załóżmy, że z danych otrzymaliśmy różnicę średnich równą `obs_diff`.\n",
    "\n",
    "- Pod $H_0$ nasza statystyka ma rozkład $N(0, se)$.\n",
    "- p-wartość dla testu jednostronnego (górny ogon) to:\n",
    "\n",
    "$$\n",
    " p = P_{H_0}(D \\ge d_{\\text{obs}}).\n",
    "$$\n",
    "\n",
    "Policzymy to najpierw numerycznie, całkując gęstość (bez użycia gotowej funkcji),\n",
    "a potem zweryfikujemy wyniki funkcją `stats.norm.sf`.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd555d21",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Przykładowy wynik z danych: obserwowana różnica średnich (pkt)\n",
    "obs_diff = 1.0\n",
    "p_value = stats.norm.sf(obs_diff, loc=0.0, scale=standard_error)\n",
    "print(\"p (stats.norm.sf):\", float(p_value))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6119c903",
   "metadata": {},
   "source": [
    "### Decyzja i interpretacja\n",
    "\n",
    "- Odrzucamy $H_0$, gdy `obs_diff` leży w regionie krytycznym, czyli `obs_diff >= critical_value`.\n",
    "- Równoważnie: odrzucamy $H_0$, gdy p-wartość $\\le \\alpha$.\n",
    "\n",
    "Pamiętaj: nawet jeśli odrzucamy $H_0$, to wciąż warto od razu raportować **wielkość efektu**\n",
    "(np. różnicę średnich lub $d$ Cohena) oraz **niepewność** (np. przedział ufności).\n",
    "Moc odpowiada na pytanie: „czy przy takim planie badania mamy szansę wykryć różnicę, która nas interesuje?”"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5970a05f",
   "metadata": {},
   "outputs": [],
   "source": [
    "reject_h0 = obs_diff >= critical_value\n",
    "print(\"Czy odrzucamy H0 (region krytyczny)?\", reject_h0)\n",
    "print(\"Czy p <= alpha?\", bool(p_value <= alpha))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3884dcea",
   "metadata": {},
   "outputs": [],
   "source": [
    "z_crit_975 = stats.norm.isf(0.025)\n",
    "ci_low = obs_diff - z_crit_975 * standard_error\n",
    "ci_high = obs_diff + z_crit_975 * standard_error\n",
    "\n",
    "print(\"95% przedział ufności dla różnicy średnich (pkt):\", (float(ci_low), float(ci_high)))\n",
    "print(\"d Cohena (na podstawie obs_diff i znanego sigma):\", float(obs_diff / sigma))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5cb0d637",
   "metadata": {},
   "source": [
    "### Moc (analitycznie)\n",
    "\n",
    "Moc testu, czyli prawdopodobieństwo odrzucenia hipotezy zerowej, kiedy prawdziwa jest hipoteza alternatywna, to powierzchnia pod rozkładem dla $H_A$ na prawo od wartości krytycznej.\n",
    "\n",
    "$$\n",
    "\\text{power} = P_{H_A}(D \\ge \\text{critical}).\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ac9372e",
   "metadata": {},
   "outputs": [],
   "source": [
    "power_known_sigma = stats.norm.sf(critical_value, loc=diff, scale=standard_error)\n",
    "beta_known_sigma = 1 - power_known_sigma\n",
    "\n",
    "print(\"Moc:\", float(power_known_sigma))\n",
    "print(\"Beta (błąd II rodzaju):\", float(beta_known_sigma))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "735768a9",
   "metadata": {},
   "source": [
    "### Moc (symulacyjnie)\n",
    "\n",
    "Zobaczmy, czy symulacja daje to samo (w granicach błędu losowego).\n",
    "\n",
    "Symulujemy wiele „mini-badań”: losujemy dwie próby z populacji różniących się średnią o `diff`, liczymy różnicę średnich i sprawdzamy, czy przekroczyliśmy wartość krytyczną."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1706f7e8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def simulate_diff_of_means(\n",
    "    n_simulations: int,\n",
    "    n_per_group: int,\n",
    "    mu1: float,\n",
    "    mu2: float,\n",
    "    sigma: float,\n",
    "    rng: np.random.Generator,\n",
    ") -> np.ndarray:\n",
    "    \"\"\"Symuluje rozkład różnicy średnich (średnia1 - średnia2).\"\"\"\n",
    "\n",
    "    diffs = np.empty(n_simulations, dtype=float)\n",
    "\n",
    "    for i in range(n_simulations):\n",
    "        sample1 = rng.normal(loc=mu1, scale=sigma, size=n_per_group)\n",
    "        sample2 = rng.normal(loc=mu2, scale=sigma, size=n_per_group)\n",
    "        diffs[i] = sample1.mean() - sample2.mean()\n",
    "\n",
    "    return diffs\n",
    "\n",
    "\n",
    "n_sim = 50_000\n",
    "simulated_diffs_ha = simulate_diff_of_means(\n",
    "    n_sim, n, mu1=mean_human_label, mu2=mean_ai_label, sigma=sigma, rng=rng\n",
    ")\n",
    "simulated_power = float(np.mean(simulated_diffs_ha >= critical_value))\n",
    "\n",
    "print(\"Moc (symulacja):\", simulated_power)\n",
    "print(\"Moc (analitycznie):\", float(power_known_sigma))\n",
    "print(\"Różnica:\", float(simulated_power - power_known_sigma))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "917df592",
   "metadata": {},
   "source": [
    "### Widżet: moc jako pole pod krzywą (znana wariancja)\n",
    "\n",
    "Poniżej możesz zmieniać parametry planu badania i od razu zobaczyć na wykresie:\n",
    "\n",
    "- obszar błędu I rodzaju ($\\alpha$) pod rozkładem dla $H_0$,\n",
    "- obszar mocy ($1-\\beta$) pod rozkładem dla $H_A$.\n",
    "\n",
    "Uwaga: to jest model uproszczony (znana $\\sigma$, normalność). W kolejnej części wrócimy do testu t.\n",
    "\n",
    "Skala osi x jest stała, żeby przy zmianie `n`, `σ` i `Δ` było widać rzeczywistą zmianę szerokości i położenia rozkładów, a nie tylko efekt autoskalowania wykresu.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea86ad6a",
   "metadata": {},
   "outputs": [],
   "source": [
    "if not WIDGETY_DOSTEPNE:\n",
    "    print(\n",
    "        \"Widżety nie są dostępne. Jeśli pracujesz w Jupyterze, doinstaluj pakiet ipywidgets, \"\n",
    "        \"aby korzystać z części interaktywnych.\"\n",
    "    )\n",
    "else:\n",
    "    n_slider = widgets.IntSlider(\n",
    "        value=int(n),\n",
    "        min=5,\n",
    "        max=200,\n",
    "        step=1,\n",
    "        description=\"n\",\n",
    "        continuous_update=False,\n",
    "    )\n",
    "    alpha_slider = widgets.FloatSlider(\n",
    "        value=float(alpha),\n",
    "        min=0.001,\n",
    "        max=0.20,\n",
    "        step=0.001,\n",
    "        description=\"α\",\n",
    "        readout_format=\".3f\",\n",
    "        continuous_update=False,\n",
    "    )\n",
    "    sigma_slider = widgets.FloatSlider(\n",
    "        value=max(3.0, float(sigma)),\n",
    "        min=3.0,\n",
    "        max=8.0,\n",
    "        step=0.1,\n",
    "        description=\"σ (pkt)\",\n",
    "        readout_format=\".1f\",\n",
    "        continuous_update=False,\n",
    "    )\n",
    "    diff_slider = widgets.FloatSlider(\n",
    "        value=max(2.0, float(diff)),\n",
    "        min=2.0,\n",
    "        max=4.0,\n",
    "        step=0.1,\n",
    "        description=\"Δ (pkt)\",\n",
    "        readout_format=\".1f\",\n",
    "        continuous_update=False,\n",
    "    )\n",
    "    test_dropdown = widgets.Dropdown(\n",
    "        options=[\n",
    "            (\"Jednostronny (μ1 > μ2)\", \"greater\"),\n",
    "            (\"Dwustronny (μ1 ≠ μ2)\", \"two-sided\"),\n",
    "        ],\n",
    "        value=\"greater\",\n",
    "        description=\"Test\",\n",
    "    )\n",
    "\n",
    "    # Stała skala osi x: zostawiamy celowo węższą, żeby zmiany σ i Δ były\n",
    "    # wyraźniej widoczne na wykładzie, zamiast zachowywać ten sam wygląd po przeskalowaniu.\n",
    "    fixed_x_min = -7.0\n",
    "    fixed_x_max = 12.0\n",
    "    fixed_x = np.linspace(fixed_x_min, fixed_x_max, 1000)\n",
    "\n",
    "    def pokaz_moc_normal(n, alpha, sigma, diff, test):\n",
    "        n = int(n)\n",
    "        alpha = float(alpha)\n",
    "        sigma = float(sigma)\n",
    "        diff = float(diff)\n",
    "\n",
    "        se = sigma * np.sqrt(2 / n)\n",
    "        d_cohena = diff / sigma if sigma > 0 else float(\"nan\")\n",
    "\n",
    "        if test == \"greater\":\n",
    "            crit_high = float(stats.norm.isf(alpha, loc=0.0, scale=se))\n",
    "            crit_low = None\n",
    "            beta = float(stats.norm.cdf(crit_high, loc=diff, scale=se))\n",
    "            power = 1 - beta\n",
    "        else:\n",
    "            crit_high = float(stats.norm.isf(alpha / 2, loc=0.0, scale=se))\n",
    "            crit_low = -crit_high\n",
    "            beta = float(\n",
    "                stats.norm.cdf(crit_high, loc=diff, scale=se)\n",
    "                - stats.norm.cdf(crit_low, loc=diff, scale=se)\n",
    "            )\n",
    "            power = 1 - beta\n",
    "\n",
    "        print(f\"n w grupie = {n} (razem {2 * n})\")\n",
    "        print(f\"σ = {sigma:.1f} pkt, Δ = {diff:.1f} pkt, d Cohena ≈ {d_cohena:.3f}\")\n",
    "        print(f\"se = {se:.2f} pkt\")\n",
    "\n",
    "        if crit_low is None:\n",
    "            print(f\"wartość krytyczna = {crit_high:.2f} pkt (górny ogon)\")\n",
    "        else:\n",
    "            print(f\"wartości krytyczne = {crit_low:.2f} pkt i {crit_high:.2f} pkt\")\n",
    "\n",
    "        print(f\"moc ≈ {power:.3f} (β ≈ {beta:.3f})\")\n",
    "\n",
    "        x = fixed_x\n",
    "\n",
    "        pdf_h0 = stats.norm.pdf(x, loc=0.0, scale=se)\n",
    "        pdf_ha = stats.norm.pdf(x, loc=diff, scale=se)\n",
    "\n",
    "        plt.figure(figsize=(9, 4))\n",
    "        plt.plot(x, pdf_h0, color=\"red\", linewidth=2, label=r\"$H_0$\")\n",
    "        plt.plot(x, pdf_ha, color=\"blue\", linewidth=2, label=r\"$H_A$\")\n",
    "\n",
    "        if crit_low is None:\n",
    "            mask_alpha = x >= crit_high\n",
    "            plt.fill_between(\n",
    "                x[mask_alpha],\n",
    "                pdf_h0[mask_alpha],\n",
    "                color=\"red\",\n",
    "                alpha=0.20,\n",
    "                label=\"α (błąd I rodzaju)\",\n",
    "            )\n",
    "            mask_power = x >= crit_high\n",
    "            plt.fill_between(\n",
    "                x[mask_power],\n",
    "                pdf_ha[mask_power],\n",
    "                color=\"blue\",\n",
    "                alpha=0.20,\n",
    "                label=\"moc (1−β)\",\n",
    "            )\n",
    "            plt.axvline(crit_high, color=\"black\", linestyle=\"--\", linewidth=2)\n",
    "        else:\n",
    "            mask_alpha_left = x <= crit_low\n",
    "            mask_alpha_right = x >= crit_high\n",
    "            plt.fill_between(\n",
    "                x[mask_alpha_left],\n",
    "                pdf_h0[mask_alpha_left],\n",
    "                color=\"red\",\n",
    "                alpha=0.20,\n",
    "                label=\"α (błąd I rodzaju)\",\n",
    "            )\n",
    "            plt.fill_between(\n",
    "                x[mask_alpha_right],\n",
    "                pdf_h0[mask_alpha_right],\n",
    "                color=\"red\",\n",
    "                alpha=0.20,\n",
    "            )\n",
    "\n",
    "            mask_power_left = x <= crit_low\n",
    "            mask_power_right = x >= crit_high\n",
    "            plt.fill_between(\n",
    "                x[mask_power_left],\n",
    "                pdf_ha[mask_power_left],\n",
    "                color=\"blue\",\n",
    "                alpha=0.20,\n",
    "                label=\"moc (1−β)\",\n",
    "            )\n",
    "            plt.fill_between(\n",
    "                x[mask_power_right],\n",
    "                pdf_ha[mask_power_right],\n",
    "                color=\"blue\",\n",
    "                alpha=0.20,\n",
    "            )\n",
    "\n",
    "            plt.axvline(crit_low, color=\"black\", linestyle=\"--\", linewidth=2)\n",
    "            plt.axvline(crit_high, color=\"black\", linestyle=\"--\", linewidth=2)\n",
    "\n",
    "        plt.axvline(0, color=\"red\", linestyle=\":\", linewidth=1.5, alpha=0.8)\n",
    "        plt.axvline(diff, color=\"blue\", linestyle=\":\", linewidth=1.5, alpha=0.8)\n",
    "        plt.xlim(fixed_x_min, fixed_x_max)\n",
    "        plt.title(\"Moc testu jako pole pod krzywą\")\n",
    "        plt.xlabel(\"różnica średnich z próby (pkt; stała skala)\")\n",
    "        plt.ylabel(\"gęstość\")\n",
    "        plt.legend()\n",
    "        plt.tight_layout()\n",
    "        plt.show()\n",
    "\n",
    "    out = widgets.interactive_output(\n",
    "        pokaz_moc_normal,\n",
    "        {\n",
    "            \"n\": n_slider,\n",
    "            \"alpha\": alpha_slider,\n",
    "            \"sigma\": sigma_slider,\n",
    "            \"diff\": diff_slider,\n",
    "            \"test\": test_dropdown,\n",
    "        },\n",
    "    )\n",
    "\n",
    "    controls = widgets.VBox(\n",
    "        [\n",
    "            widgets.HBox([n_slider, alpha_slider]),\n",
    "            widgets.HBox([sigma_slider, diff_slider]),\n",
    "            test_dropdown,\n",
    "        ]\n",
    "    )\n",
    "\n",
    "    display(controls, out)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e70c5628",
   "metadata": {},
   "source": [
    "### Co wpływa na moc?\n",
    "\n",
    "Żeby mówić o mocy, musimy określić:\n",
    "\n",
    "- hipotezę i statystykę testową,\n",
    "- poziom istotności ($\\alpha$),\n",
    "- wariancję/odchylenie standardowe w populacji ($\\sigma^2$ / $\\sigma$),\n",
    "- spodziewany efekt (np. różnicę średnich przy $H_A$).\n",
    "\n",
    "**Czynniki wpływające na moc testu**\n",
    "\n",
    "- Im większa różnica, na wykryciu której nam zależy, tym łatwiej ją wykryć. Zatem, im większa różnica między średnimi przewidywana przez hipotezę alternatywną (im większy spodziewany efekt), tym większa moc testu.\n",
    "- Im mniejszy „szum” w naszych danych, tym łatwiej wykryć interesującą nas różnicę. Zatem, im mniejsza wariancja, tym większa moc testu (bo od wariancji zależy błąd standardowy, a im mniejszy błąd standardowy, tym większa moc). Na wariancję wyników mamy niewielki wpływ, ale warto pamiętać, że precyzja pomiaru może ją zmniejszać.\n",
    "- Im bardziej jesteśmy skłonni akceptować „fałszywe alarmy”, tym łatwiej będzie nam wykryć prawdziwe różnice. Zatem, zwiększając poziom $\\alpha$, zwiększamy równocześnie moc testu: większa $\\alpha$ oznacza niższą wartość krytyczną. Coś za coś: rośnie wtedy prawdopodobieństwo błędu I rodzaju.\n",
    "- Wreszcie, im więcej pomiarów, tym łatwiej wykryć interesującą nas różnicę. Im większa liczebność, tym większa moc testu (bo od liczebności zależy błąd standardowy, a im mniejszy błąd standardowy, tym większa moc)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10829875",
   "metadata": {},
   "source": [
    "### (Dygresja) Widżet: dlaczego plan wewnątrzosobniczy bywa mocniejszy?\n",
    "\n",
    "W głównym przykładzie porównujemy dwie niezależne grupy. W części badań nad oceną technologii można jednak zastosować plan wewnątrzosobniczy: ta sama osoba ocenia wiele krótkich tekstów, z których część ma etykietę „człowiek”, a część etykietę „AI”.\n",
    "\n",
    "Wtedy porównujemy **dwie średnie oceny u tej samej osoby**, a nie dwie niezależne grupy.\n",
    "\n",
    "Jeśli oceny tej samej osoby w dwóch warunkach są dodatnio skorelowane, to zmienność różnicy maleje, a moc rośnie.\n",
    "\n",
    "Poniżej porównujemy (w uproszczonym modelu) moc planu:\n",
    "\n",
    "- międzyosobniczego (dwie niezależne grupy),\n",
    "- wewnątrzosobniczego (ten sam uczestnik w obu warunkach),\n",
    "\n",
    "w zależności od korelacji $\\rho$ między warunkami."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e9539bf2",
   "metadata": {},
   "outputs": [],
   "source": [
    "if not WIDGETY_DOSTEPNE:\n",
    "    print(\n",
    "        \"Widżety nie są dostępne. Jeśli pracujesz w Jupyterze, doinstaluj pakiet ipywidgets, \"\n",
    "        \"aby korzystać z części interaktywnych.\"\n",
    "    )\n",
    "else:\n",
    "    n_slider = widgets.IntSlider(\n",
    "        value=int(n),\n",
    "        min=5,\n",
    "        max=200,\n",
    "        step=1,\n",
    "        description=\"n\",\n",
    "        continuous_update=False,\n",
    "    )\n",
    "    alpha_slider = widgets.FloatSlider(\n",
    "        value=float(alpha),\n",
    "        min=0.001,\n",
    "        max=0.20,\n",
    "        step=0.001,\n",
    "        description=\"α\",\n",
    "        readout_format=\".3f\",\n",
    "        continuous_update=False,\n",
    "    )\n",
    "    sigma_slider = widgets.FloatSlider(\n",
    "        value=float(sigma),\n",
    "        min=0.2,\n",
    "        max=5.0,\n",
    "        step=0.1,\n",
    "        description=\"σ (pkt)\",\n",
    "        readout_format=\".1f\",\n",
    "        continuous_update=False,\n",
    "    )\n",
    "    diff_slider = widgets.FloatSlider(\n",
    "        value=float(diff),\n",
    "        min=0.0,\n",
    "        max=3.0,\n",
    "        step=0.1,\n",
    "        description=\"Δ (pkt)\",\n",
    "        readout_format=\".1f\",\n",
    "        continuous_update=False,\n",
    "    )\n",
    "    rho_slider = widgets.FloatSlider(\n",
    "        value=0.50,\n",
    "        min=0.00,\n",
    "        max=0.90,\n",
    "        step=0.05,\n",
    "        description=\"ρ\",\n",
    "        readout_format=\".2f\",\n",
    "        continuous_update=False,\n",
    "    )\n",
    "\n",
    "    def pokaz_plan_within_between(n, alpha, sigma, diff, rho):\n",
    "        n = int(n)\n",
    "        alpha = float(alpha)\n",
    "        sigma = float(sigma)\n",
    "        diff = float(diff)\n",
    "        rho = float(rho)\n",
    "\n",
    "        # Międzyosobniczy: dwie niezależne próby\n",
    "        se_between = sigma * np.sqrt(2 / n)\n",
    "        crit_between = float(stats.norm.isf(alpha, loc=0.0, scale=se_between))\n",
    "        power_between = float(stats.norm.sf(crit_between, loc=diff, scale=se_between))\n",
    "\n",
    "        # Wewnątrzosobniczy (sparowany): przybliżenie przez korelację między warunkami\n",
    "        se_within = sigma * np.sqrt(2 * (1 - rho) / n)\n",
    "        crit_within = float(stats.norm.isf(alpha, loc=0.0, scale=se_within))\n",
    "        power_within = float(stats.norm.sf(crit_within, loc=diff, scale=se_within))\n",
    "\n",
    "        print(f\"n = {n}, α = {alpha:.3f}, ρ = {rho:.2f}\")\n",
    "        print(f\"se (międzyosobniczy) = {se_between:.2f} pkt → moc ≈ {power_between:.3f}\")\n",
    "        print(f\"se (wewnątrzosobniczy) = {se_within:.2f} pkt → moc ≈ {power_within:.3f}\")\n",
    "\n",
    "        rho_grid = np.linspace(0.0, 0.9, 100)\n",
    "        se_grid = sigma * np.sqrt(2 * (1 - rho_grid) / n)\n",
    "        crit_grid = stats.norm.isf(alpha, loc=0.0, scale=se_grid)\n",
    "        power_grid = stats.norm.sf(crit_grid, loc=diff, scale=se_grid)\n",
    "\n",
    "        plt.figure(figsize=(9, 4))\n",
    "        plt.plot(rho_grid, power_grid, color=\"blue\", linewidth=2, label=\"wewnątrzosobniczy\")\n",
    "        plt.axhline(\n",
    "            power_between,\n",
    "            color=\"red\",\n",
    "            linestyle=\"--\",\n",
    "            linewidth=2,\n",
    "            label=\"międzyosobniczy\",\n",
    "        )\n",
    "        plt.axvline(rho, color=\"black\", linestyle=\":\", linewidth=2)\n",
    "        plt.scatter([rho], [power_within], color=\"black\", zorder=5)\n",
    "\n",
    "        plt.title(\"Moc a korelacja między warunkami (plan wewnątrzosobniczy)\")\n",
    "        plt.xlabel(\"korelacja między warunkami (ρ)\")\n",
    "        plt.ylabel(\"moc\")\n",
    "        plt.ylim(0.0, 1.0)\n",
    "        plt.legend()\n",
    "        plt.tight_layout()\n",
    "        plt.show()\n",
    "\n",
    "    out = widgets.interactive_output(\n",
    "        pokaz_plan_within_between,\n",
    "        {\n",
    "            \"n\": n_slider,\n",
    "            \"alpha\": alpha_slider,\n",
    "            \"sigma\": sigma_slider,\n",
    "            \"diff\": diff_slider,\n",
    "            \"rho\": rho_slider,\n",
    "        },\n",
    "    )\n",
    "\n",
    "    controls = widgets.VBox(\n",
    "        [\n",
    "            widgets.HBox([n_slider, alpha_slider, rho_slider]),\n",
    "            widgets.HBox([sigma_slider, diff_slider]),\n",
    "        ]\n",
    "    )\n",
    "\n",
    "    display(controls, out)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "153f5064",
   "metadata": {},
   "source": [
    "## Nieznana wariancja\n",
    "\n",
    "Obliczenie mocy w przykładzie z poprzedniej części było prostsze niż zazwyczaj. Założyliśmy bowiem, że znamy wariancję w populacji. Wtedy korzystamy z faktu, że rozkład różnicy średnich to rozkład normalny. Zmiana hipotezy zerowej na alternatywną powoduje, że jeden parametr tego rozkładu — jego średnia — zmienia się z 0 na d, ale nie wpływa to na drugi parametr — błąd standardowy.\n",
    "\n",
    "Nie znając wariancji w populacji, a zazwyczaj jej nie znamy, musimy skorzystać z rozkładu t: jeśli błąd standardowy liczymy, biorąc wariancję policzoną z danych, różnica między statystyką a parametrem podzielona przez ten błąd ma taki właśnie rozkład.\n",
    "\n",
    "Jak w takim wypadku wygląda rozkład statystyki przy założeniu hipotezy alternatywnej? Z pomocą prostej symulacji łatwo to sprawdzić."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3fa70d63",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_sim = 10_000\n",
    "\n",
    "diffs_h0 = simulate_diff_of_means(\n",
    "    n_sim, n, mu1=mean_ai_label, mu2=mean_ai_label, sigma=sigma, rng=rng\n",
    ")\n",
    "diffs_ha = simulate_diff_of_means(\n",
    "    n_sim, n, mu1=mean_human_label, mu2=mean_ai_label, sigma=sigma, rng=rng\n",
    ")\n",
    "\n",
    "plt.hist(diffs_h0, bins=60, density=True, color=\"red\", alpha=0.5, label=r\"$H_0$\")\n",
    "plt.hist(diffs_ha, bins=60, density=True, color=\"blue\", alpha=0.5, label=r\"$H_A$\")\n",
    "\n",
    "plt.title(\"Rozkład różnicy średnich (symulacja)\")\n",
    "plt.xlabel(\"różnica średnich z próby (pkt)\")\n",
    "plt.ylabel(\"gęstość\")\n",
    "plt.legend()\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4efe727e",
   "metadata": {},
   "source": [
    "Łatwo możemy sprawdzić, że te histogramy odpowiadają krzywym normalnym o odpowiednich parametrach (średnia równa 0 albo `diff`, a odchylenie standardowe równe `se`). Możemy też dorysować wartość krytyczną."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c1feb00",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(diffs_h0, bins=60, density=True, color=\"red\", alpha=0.5, label=r\"$H_0$\")\n",
    "plt.hist(diffs_ha, bins=60, density=True, color=\"blue\", alpha=0.5, label=r\"$H_A$\")\n",
    "\n",
    "plt.plot(x, pdf_h0, color=\"red\", linewidth=2)\n",
    "plt.plot(x, pdf_ha, color=\"blue\", linewidth=2)\n",
    "\n",
    "plt.axvline(critical_value, linestyle=\"--\", color=\"black\", label=\"wartość krytyczna\")\n",
    "\n",
    "plt.title(\"Symulacja + krzywe normalne\")\n",
    "plt.xlabel(\"różnica średnich z próby\")\n",
    "plt.ylabel(\"gęstość\")\n",
    "plt.xlim(-5, 6)\n",
    "plt.legend()\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24fa21c6",
   "metadata": {},
   "source": [
    "### Standaryzacja: dzielimy przez błąd standardowy\n",
    "\n",
    "W kolejnym kroku nakreślmy te same histogramy, ale standaryzując różnice między średnimi, czyli dzieląc je przez błąd standardowy.\n",
    "Dla $H_0$ możemy dorysować standardowy $N(0,1)$, ponieważ standaryzujemy przez `se`. \n",
    "\n",
    "Zmianie ulegają oba rozkłady: zarówno ich średnie, jak i odchylenia standardowe należy podzielić przez tę samą wartość.Dlatego w przypadku hipotezy zerowej mamy teraz do czynienia ze standardowym $N(0, 1)$, a w przypadku hipotezy alternatywnej — z $N(\\text{diff}/se, 1)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "439f92c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "standardized_h0 = diffs_h0 / standard_error\n",
    "standardized_ha = diffs_ha / standard_error\n",
    "x_std = np.linspace(-5, 11, 500)\n",
    "mean_under_ha_standardized = diff / standard_error\n",
    "\n",
    "plt.hist(standardized_h0, bins=60, density=True, color=\"red\", alpha=0.5, label=r\"$H_0$\")\n",
    "plt.hist(standardized_ha, bins=60, density=True, color=\"blue\", alpha=0.5, label=r\"$H_A$\")\n",
    "\n",
    "plt.plot(x_std, stats.norm.pdf(x_std, loc=0.0, scale=1.0), color=\"red\", linewidth=2, label=r\"$H_0: N(0,1)$\")\n",
    "plt.plot(\n",
    "    x_std,\n",
    "    stats.norm.pdf(x_std, loc=mean_under_ha_standardized, scale=1.0),\n",
    "    color=\"blue\",\n",
    "    linewidth=2,\n",
    "    label=rf\"$H_A: N({mean_under_ha_standardized:.2f},1)$\",\n",
    ")\n",
    "plt.axvline(stats.norm.isf(alpha), linestyle=\"--\", color=\"black\", label=\"wartość krytyczna\")\n",
    "plt.title(\"Standaryzacja: $H_0$ vs $H_A$\")\n",
    "plt.xlabel(r\"$(\\bar{x}_1-\\bar{x}_2)/se$\")\n",
    "plt.ylabel(\"gęstość\")\n",
    "plt.xlim(-5, 11)\n",
    "plt.legend()\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "82868203",
   "metadata": {},
   "source": [
    "### Wartość krytyczna po standaryzacji\n",
    "\n",
    "Zmianie ulega również wartość krytyczna: dzielimy ją przez `se`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "908b8574",
   "metadata": {},
   "outputs": [],
   "source": [
    "critical_value_standardized = critical_value / standard_error\n",
    "print(\"Wartość krytyczna (standaryzowana):\", critical_value_standardized)\n",
    "print(\"Dla porównania, norm.isf(alpha):\", stats.norm.isf(alpha))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4aa6eb33",
   "metadata": {},
   "source": [
    "## Moc testu $t$\n",
    "\n",
    "Dla przypomnienia: różnica między średnimi podzielona przez błąd standardowy to statystyka $z$, jeśli znamy wariancję w populacji, albo statystyka $t$, jeśli szacujemy wariancję na podstawie danych z próby.\n",
    "\n",
    "W tym drugim przypadku, żeby zasymulować rozkład odpowiadający hipotezie zerowej, przy każdym losowaniu musimy na nowo szacować błąd standardowy, korzystając z wariancji policzonych w wylosowanych próbach."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "39f468d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "def simulate_t_statistic_two_sample(\n",
    "    n_simulations: int,\n",
    "    n_per_group: int,\n",
    "    mu1: float,\n",
    "    mu2: float,\n",
    "    sigma: float,\n",
    "    rng: np.random.Generator,\n",
    ") -> np.ndarray:\n",
    "    \"\"\"Symuluje statystykę t dla różnicy średnich (wariancje z próby).\"\"\"\n",
    "\n",
    "    t_values = np.empty(n_simulations, dtype=float)\n",
    "\n",
    "    for i in range(n_simulations):\n",
    "        sample1 = rng.normal(loc=mu1, scale=sigma, size=n_per_group)\n",
    "        sample2 = rng.normal(loc=mu2, scale=sigma, size=n_per_group)\n",
    "\n",
    "        mean_diff = sample1.mean() - sample2.mean()\n",
    "        var1 = sample1.var(ddof=1)\n",
    "        var2 = sample2.var(ddof=1)\n",
    "\n",
    "        se_hat = np.sqrt((var1 + var2) / n_per_group)\n",
    "        t_values[i] = mean_diff / se_hat\n",
    "\n",
    "    return t_values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df13f6f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_sim = 10_000\n",
    "df = 2 * n - 2\n",
    "\n",
    "t_h0 = simulate_t_statistic_two_sample(\n",
    "    n_sim, n, mu1=mean_ai_label, mu2=mean_ai_label, sigma=sigma, rng=rng\n",
    ")\n",
    "\n",
    "plt.hist(t_h0, bins=60, density=True, color=\"red\", alpha=0.5, label=r\"$H_0$ (symulacja)\")\n",
    "\n",
    "x_t = np.linspace(-5, 11, 400)\n",
    "plt.plot(x_t, stats.t.pdf(x_t, df=df), color=\"black\", linewidth=2, label=rf\"$t$ (df={df})\")\n",
    "\n",
    "plt.title(\"Statystyka t pod $H_0$\")\n",
    "plt.xlabel(\"t\")\n",
    "plt.ylabel(\"gęstość\")\n",
    "plt.xlim(-5, 11)\n",
    "plt.legend()\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1fb6bb12",
   "metadata": {},
   "source": [
    "Teraz dorysujemy histogram, który przedstawia przybliżenie rozkładu statystyki t przy założeniu hipotezy alternatywnej."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e0a35ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "t_ha = simulate_t_statistic_two_sample(\n",
    "    n_sim, n, mu1=mean_human_label, mu2=mean_ai_label, sigma=sigma, rng=rng\n",
    ")\n",
    "\n",
    "plt.hist(t_h0, bins=60, density=True, color=\"red\", alpha=0.5, label=r\"$H_0$ (symulacja)\")\n",
    "plt.hist(t_ha, bins=60, density=True, color=\"blue\", alpha=0.5, label=r\"$H_A$ (symulacja)\")\n",
    "\n",
    "plt.plot(x_t, stats.t.pdf(x_t, df=df), color=\"black\", linewidth=2, label=rf\"$t$ (df={df})\")\n",
    "\n",
    "plt.title(\"Statystyka t: $H_0$ vs $H_A$ (symulacja)\")\n",
    "plt.xlabel(\"t\")\n",
    "plt.ylabel(\"gęstość\")\n",
    "plt.xlim(-5, 11)\n",
    "plt.legend()\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "65b3b4fe",
   "metadata": {},
   "source": [
    "Jeśli się bliżej przyjrzeć nowemu rozkładowi, nie jest to po prostu rozkład $t$ przesunięty w prawo. Jest on trochę niższy i nieznacznie skośny.\n",
    "\n",
    "Okazuje się, że rozkład t, z jakim mieliśmy do tej pory do czynienia, jest tylko szczególnym przypadkiem ogólniejszego rozkładu: **niecentralnego rozkładu t**.\n",
    "\n",
    "Parametr niecentralności dla rozkładu t przy założeniu hipotezy alternatywnej wynosi `diff/se`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "264df3c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "noncentrality = diff / standard_error\n",
    "\n",
    "plt.hist(t_h0, bins=60, density=True, color=\"red\", alpha=0.5, label=r\"$H_0$ (symulacja)\")\n",
    "plt.hist(t_ha, bins=60, density=True, color=\"blue\", alpha=0.5, label=r\"$H_A$ (symulacja)\")\n",
    "\n",
    "plt.plot(x_t, stats.t.pdf(x_t, df=df), color=\"red\", linewidth=2, label=rf\"centralny t (df={df})\")\n",
    "plt.plot(\n",
    "    x_t,\n",
    "    stats.nct.pdf(x_t, df=df, nc=noncentrality), # nct to non-central t\n",
    "    color=\"blue\",\n",
    "    linewidth=2,\n",
    "    label=rf\"niecentralny t (nc={noncentrality:.2f})\",\n",
    ")\n",
    "\n",
    "plt.title(\"Niecentralny rozkład t jako model pod $H_A$\")\n",
    "plt.xlabel(\"t\")\n",
    "plt.ylabel(\"gęstość\")\n",
    "plt.xlim(-5, 11)\n",
    "plt.legend()\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f1b1775",
   "metadata": {},
   "outputs": [],
   "source": [
    "sample1 = rng.normal(loc=mean_human_label, scale=sigma, size=n)\n",
    "sample2 = rng.normal(loc=mean_ai_label, scale=sigma, size=n)\n",
    "\n",
    "mean_diff_obs = sample1.mean() - sample2.mean()\n",
    "var1 = sample1.var(ddof=1)\n",
    "var2 = sample2.var(ddof=1)\n",
    "se_hat = np.sqrt((var1 + var2) / n)\n",
    "t_obs = mean_diff_obs / se_hat\n",
    "\n",
    "print(\"Różnica średnich (obs, pkt):\", float(mean_diff_obs))\n",
    "print(\"t (obs):\", float(t_obs))\n",
    "\n",
    "p_t_manual, p_t_manual_error = integrate.quad(\n",
    "    lambda t: stats.t.pdf(t, df=df),\n",
    "    t_obs,\n",
    "    np.inf,\n",
    ")\n",
    "p_t_function = stats.t.sf(t_obs, df=df)\n",
    "\n",
    "print(\"p (całkowanie, górny ogon):\", float(p_t_manual))\n",
    "print(\"p (stats.t.sf):\", float(p_t_function))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e83faace",
   "metadata": {},
   "outputs": [],
   "source": [
    "t_critical = stats.t.isf(alpha, df=df)\n",
    "reject_h0_t = t_obs >= t_critical\n",
    "\n",
    "print(\"Wartość krytyczna t:\", float(t_critical))\n",
    "print(\"Czy odrzucamy H0?\", bool(reject_h0_t))\n",
    "print(\"Czy p <= alpha?\", bool(p_t_function <= alpha))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "24655cb7",
   "metadata": {},
   "outputs": [],
   "source": [
    "ttest = stats.ttest_ind(sample1, sample2, equal_var=True, alternative=\"greater\")\n",
    "print(\"Test t (SciPy): stat=\", float(ttest.statistic), \"p=\", float(ttest.pvalue))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c220c98e",
   "metadata": {},
   "source": [
    "### Wartość krytyczna i moc testu t\n",
    "\n",
    "Żeby obliczyć moc testu t, wyznaczamy wartość krytyczną z centralnego rozkładu t, a następnie liczymy, jaka część niecentralnego rozkładu t znajduje się na prawo od tej wartości."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f3f7b33",
   "metadata": {},
   "outputs": [],
   "source": [
    "t_critical = stats.t.isf(alpha, df=df)\n",
    "power_t = stats.nct.sf(t_critical, df=df, nc=noncentrality)\n",
    "\n",
    "print(\"Wartość krytyczna t:\", float(t_critical))\n",
    "print(\"Moc testu t:\", float(power_t))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "93be3c8b",
   "metadata": {},
   "source": [
    "Może wydawać się dziwne, że do obliczenia mocy testu t korzystamy z `se`, skoro test t stosujemy wtedy,\n",
    "kiedy tego błędu nie znamy.\n",
    "\n",
    "Należy jednak pamiętać, że przy planowaniu mocy `se` opisuje **parametry populacji** (model),\n",
    "a nie konkretną zrealizowaną próbę.\n",
    "\n",
    "Parametr niecentralności możemy też przepisać w sposób, który od razu pokazuje rolę wielkości efektu\n",
    "(d Cohena) i liczebności:\n",
    "\n",
    "$$\n",
    "\\frac{\\text{diff}}{se} = \\frac{\\text{diff}}{\\sigma}\\sqrt{\\frac{n}{2}}.\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f33a20c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "noncentrality_alt = (diff / sigma) * np.sqrt(n / 2)\n",
    "print(\"Niecentralność (przepisana):\", float(noncentrality_alt))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47d283cb",
   "metadata": {},
   "source": [
    "### Widżet: moc testu t (nieznana wariancja)\n",
    "\n",
    "W teście t nie znamy odchylenia standardowego w populacji, więc szacujemy je na podstawie danych.\n",
    "Wtedy statystyka testowa ma (pod $H_0$) rozkład t, a pod hipotezą alternatywną — **niecentralny rozkład t**.\n",
    "\n",
    "Poniższy widżet pokazuje, jak zmienia się moc w zależności od:\n",
    "\n",
    "- liczebności w grupie ($n$),\n",
    "- poziomu istotności ($\\alpha$),\n",
    "- wielkości efektu (tu: różnicy w punktach skali oraz d Cohena).\n",
    "\n",
    "Uwaga: to model dla dwóch **niezależnych** prób o równej liczebności i założeniu równości wariancji.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "01843185",
   "metadata": {},
   "outputs": [],
   "source": [
    "if not WIDGETY_DOSTEPNE:\n",
    "    print(\n",
    "        \"Widżety nie są dostępne. Jeśli pracujesz w Jupyterze, doinstaluj pakiet ipywidgets, \"\n",
    "        \"aby korzystać z części interaktywnych.\"\n",
    "    )\n",
    "else:\n",
    "    n_slider = widgets.IntSlider(\n",
    "        value=int(n),\n",
    "        min=5,\n",
    "        max=200,\n",
    "        step=1,\n",
    "        description=\"n\",\n",
    "        continuous_update=False,\n",
    "    )\n",
    "    alpha_slider = widgets.FloatSlider(\n",
    "        value=float(alpha),\n",
    "        min=0.001,\n",
    "        max=0.20,\n",
    "        step=0.001,\n",
    "        description=\"α\",\n",
    "        readout_format=\".3f\",\n",
    "        continuous_update=False,\n",
    "    )\n",
    "    sigma_slider = widgets.FloatSlider(\n",
    "        value=float(sigma),\n",
    "        min=0.2,\n",
    "        max=5.0,\n",
    "        step=0.1,\n",
    "        description=\"σ (pkt)\",\n",
    "        readout_format=\".1f\",\n",
    "        continuous_update=False,\n",
    "    )\n",
    "    diff_slider = widgets.FloatSlider(\n",
    "        value=float(diff),\n",
    "        min=0.0,\n",
    "        max=3.0,\n",
    "        step=0.1,\n",
    "        description=\"Δ (pkt)\",\n",
    "        readout_format=\".1f\",\n",
    "        continuous_update=False,\n",
    "    )\n",
    "    test_dropdown = widgets.Dropdown(\n",
    "        options=[\n",
    "            (\"Jednostronny (μ1 > μ2)\", \"greater\"),\n",
    "            (\"Dwustronny (μ1 ≠ μ2)\", \"two-sided\"),\n",
    "        ],\n",
    "        value=\"greater\",\n",
    "        description=\"Test\",\n",
    "    )\n",
    "\n",
    "    def pokaz_moc_t(n, alpha, sigma, diff, test):\n",
    "        n = int(n)\n",
    "        alpha = float(alpha)\n",
    "        sigma = float(sigma)\n",
    "        diff = float(diff)\n",
    "\n",
    "        df = 2 * n - 2\n",
    "        d_cohena = diff / sigma if sigma > 0 else float(\"nan\")\n",
    "        noncentrality = d_cohena * np.sqrt(n / 2)\n",
    "\n",
    "        if test == \"greater\":\n",
    "            t_crit = float(stats.t.isf(alpha, df=df))\n",
    "            power = float(stats.nct.sf(t_crit, df=df, nc=noncentrality))\n",
    "            crits = [t_crit]\n",
    "        else:\n",
    "            t_crit = float(stats.t.isf(alpha / 2, df=df))\n",
    "            power = float(\n",
    "                stats.nct.sf(t_crit, df=df, nc=noncentrality)\n",
    "                + stats.nct.cdf(-t_crit, df=df, nc=noncentrality)\n",
    "            )\n",
    "            crits = [-t_crit, t_crit]\n",
    "\n",
    "        print(f\"n w grupie = {n} (razem {2 * n}), df = {df}\")\n",
    "        print(f\"σ = {sigma:.1f} pkt, Δ = {diff:.1f} pkt, d Cohena ≈ {d_cohena:.3f}\")\n",
    "        print(f\"parametr niecentralności ≈ {noncentrality:.3f}\")\n",
    "        print(f\"moc ≈ {power:.3f}\")\n",
    "\n",
    "        x_min = -6.0\n",
    "        x_max = max(6.0, float(noncentrality) + 8.0)\n",
    "        x = np.linspace(x_min, x_max, 1000)\n",
    "\n",
    "        pdf_h0 = stats.t.pdf(x, df=df)\n",
    "        pdf_ha = stats.nct.pdf(x, df=df, nc=noncentrality)\n",
    "\n",
    "        plt.figure(figsize=(9, 4))\n",
    "        plt.plot(x, pdf_h0, color=\"red\", linewidth=2, label=r\"$H_0$ (t)\")\n",
    "        plt.plot(x, pdf_ha, color=\"blue\", linewidth=2, label=r\"$H_A$ (niecentralny t)\")\n",
    "\n",
    "        if test == \"greater\":\n",
    "            mask_alpha = x >= crits[0]\n",
    "            plt.fill_between(\n",
    "                x[mask_alpha],\n",
    "                pdf_h0[mask_alpha],\n",
    "                color=\"red\",\n",
    "                alpha=0.20,\n",
    "                label=\"α (błąd I rodzaju)\",\n",
    "            )\n",
    "            plt.fill_between(\n",
    "                x[mask_alpha],\n",
    "                pdf_ha[mask_alpha],\n",
    "                color=\"blue\",\n",
    "                alpha=0.20,\n",
    "                label=\"moc\",\n",
    "            )\n",
    "            plt.axvline(crits[0], color=\"black\", linestyle=\"--\", linewidth=2)\n",
    "        else:\n",
    "            mask_left = x <= crits[0]\n",
    "            mask_right = x >= crits[1]\n",
    "\n",
    "            plt.fill_between(\n",
    "                x[mask_left],\n",
    "                pdf_h0[mask_left],\n",
    "                color=\"red\",\n",
    "                alpha=0.20,\n",
    "                label=\"α (błąd I rodzaju)\",\n",
    "            )\n",
    "            plt.fill_between(\n",
    "                x[mask_right],\n",
    "                pdf_h0[mask_right],\n",
    "                color=\"red\",\n",
    "                alpha=0.20,\n",
    "            )\n",
    "\n",
    "            plt.fill_between(\n",
    "                x[mask_left],\n",
    "                pdf_ha[mask_left],\n",
    "                color=\"blue\",\n",
    "                alpha=0.20,\n",
    "                label=\"moc\",\n",
    "            )\n",
    "            plt.fill_between(\n",
    "                x[mask_right],\n",
    "                pdf_ha[mask_right],\n",
    "                color=\"blue\",\n",
    "                alpha=0.20,\n",
    "            )\n",
    "\n",
    "            plt.axvline(crits[0], color=\"black\", linestyle=\"--\", linewidth=2)\n",
    "            plt.axvline(crits[1], color=\"black\", linestyle=\"--\", linewidth=2)\n",
    "\n",
    "        plt.title(\"Moc testu t\")\n",
    "        plt.xlabel(\"t\")\n",
    "        plt.ylabel(\"gęstość\")\n",
    "        plt.legend()\n",
    "        plt.tight_layout()\n",
    "        plt.show()\n",
    "\n",
    "    out = widgets.interactive_output(\n",
    "        pokaz_moc_t,\n",
    "        {\n",
    "            \"n\": n_slider,\n",
    "            \"alpha\": alpha_slider,\n",
    "            \"sigma\": sigma_slider,\n",
    "            \"diff\": diff_slider,\n",
    "            \"test\": test_dropdown,\n",
    "        },\n",
    "    )\n",
    "\n",
    "    controls = widgets.VBox(\n",
    "        [\n",
    "            widgets.HBox([n_slider, alpha_slider]),\n",
    "            widgets.HBox([sigma_slider, diff_slider]),\n",
    "            test_dropdown,\n",
    "        ]\n",
    "    )\n",
    "\n",
    "    display(controls, out)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75811722",
   "metadata": {},
   "source": [
    "### Zadanie (planowanie liczebności)\n",
    "\n",
    "W badaniach nad postrzeganiem technologii często interesuje nas, czy efekt jest na tyle duży, żeby dało się go wykryć w rozsądnym badaniu.\n",
    "Zrobimy więc proste planowanie liczebności.\n",
    "\n",
    "Załóżmy, że chcemy wykryć efekt etykiety „AI” rzędu **d Cohena = 0.5** (np. różnica średnich około 1 punktu przy odchyleniu standardowym 2 punkty).\n",
    "\n",
    "**Pytanie.** Co najmniej ilu osób potrzebujemy w każdej grupie, aby mieć nie mniej niż 80% szans na wykrycie takiego efektu w teście **dwustronnym** przy $\\alpha=0.05$?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2a1c88a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.stats.power import TTestIndPower\n",
    "\n",
    "analysis = TTestIndPower()\n",
    "\n",
    "required_n = analysis.solve_power(\n",
    "    effect_size=0.5,\n",
    "    alpha=0.05,\n",
    "    power=0.80,\n",
    "    alternative=\"two-sided\",\n",
    ")\n",
    "\n",
    "print(\"Wynik (ciągły):\", float(required_n))\n",
    "print(\"Minimalnie (zaokrąglamy w górę):\", int(np.ceil(required_n)))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c0bd43e",
   "metadata": {},
   "source": [
    "## Dodatkowe symulacje: selekcja wyników i przedziały ufności\n",
    "\n",
    "Poniższe widżety pokazują dwie rzeczy, które często są dla studentów zaskakujące:\n",
    "\n",
    "1) \"istotne\" wyniki nie są reprezentatywne dla wszystkich możliwych wyników badania (selekcja wyników),\n",
    "2) 95% przedział ufności ma sens długookresowy: w wielu powtórzeniach około 95% takich przedziałów zawiera prawdziwą wartość."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a67595c",
   "metadata": {},
   "source": [
    "### Widżet: selekcja wyników – zawyżanie estymat\n",
    "\n",
    "Symulujemy wiele powtórzeń tego samego badania. Następnie porównujemy rozkład estymowanej różnicy dla wszystkich symulacji oraz tylko dla tych, które dały wynik istotny statystycznie.\n",
    "\n",
    "To ilustruje, dlaczego w literaturze (zwłaszcza przy niskiej mocy) \"istotne\" oszacowania efektu często bywają zawyżone."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f5a656f",
   "metadata": {},
   "source": [
    "## Krótka powtórka\n",
    "\n",
    "- **Wielkość efektu** mówi: jak duża jest różnica/zależność (w jednostkach oryginalnych lub standaryzowanych).\n",
    "- **Moc** mówi: jak duża jest szansa wykrycia efektu o pewnej wielkości przy danym $\\alpha$, zmienności i liczebności.\n",
    "- Moc zawsze liczymy dla konkretnej alternatywy (np. d Cohena = 0.5).\n",
    "- Dla testu t naturalnym narzędziem teoretycznym jest **niecentralny** rozkład t.\n",
    "\n",
    "## Ćwiczenia (do samodzielnego wykonania)\n",
    "\n",
    "1. Zmień `n` (np. 10, 20, 50, 100) i sprawdź, jak rośnie moc.\n",
    "2. Zmień `sigma` (większy „szum”) i sprawdź, jak spada moc.\n",
    "3. Zmień `alpha` (np. 0.10 vs 0.01) i zobacz kompromis: moc rośnie, ale rośnie też ryzyko błędu I rodzaju.\n",
    "4. Skorzystaj z widżetów i sprawdź, jak zmienia się moc, gdy przełączasz test jednostronny/dwustronny.\n",
    "5. Narysuj krzywą mocy: moc w funkcji `n` dla kilku wartości efektu (np. d = 0.2, 0.5, 0.8).\n",
    "\n",
    "Wskazówka: możesz wykorzystać `analysis.power(effect_size=..., nobs1=..., alpha=..., ratio=1.0, alternative=...)`.\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.14.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}