{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cell-001",
   "metadata": {},
   "source": [
    "# Testowanie statystyczne hipotez\n",
    "\n",
    "Cel notatnika: przejść przez **logikę testowania hipotez** na kilku narracyjnych przykładach i zobaczyć, skąd biorą się p-wartości.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-002",
   "metadata": {},
   "source": [
    "## Cele\n",
    "\n",
    "- Umiesz postawić $H_0$ i $H_A$ (jedno- i dwustronnie).\n",
    "- Rozumiesz, czym jest p-wartość: $P(\\text{dane tak skrajne lub bardziej} \\mid H_0)$.\n",
    "- Rozróżniasz błąd I rodzaju ($\\alpha$) i II rodzaju ($\\beta$) oraz pojęcie mocy ($1-\\beta$)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-003",
   "metadata": {},
   "source": [
    "## Wprowadzenie\n",
    "\n",
    "Testowanie statystyczne hipotez jest jedną z najważniejszych technik statystyki inferencyjnej. Używa się jej niemal zawsze wtedy, gdy pracujemy z danymi pochodzącymi z badań eksperymentalnych lub obserwacyjnych.\n",
    "\n",
    "Istnieje bardzo wiele testów statystycznych, ale w tej części kursu skoncentrujemy się na samej logice testowania hipotez.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-004",
   "metadata": {},
   "source": [
    "Procedurę statystycznego testowania hipotez możemy opisać w kilku krokach:\n",
    "\n",
    "1. Wymyślamy hipotezę alternatywną (ta nas interesuje!).\n",
    "2. Przyjmujemy hipotezę zerową (i ewentualnie dodatkowe założenia).\n",
    "3. Decydujemy się na poziom istotności $\\alpha$.\n",
    "4. Konstruujemy rozkład statystyki z próby (przy założeniu $H_0$).\n",
    "5. Określamy region krytyczny (biorąc pod uwagę $\\alpha$ i kierunkowość hipotezy).\n",
    "6. Losujemy próbę, zbieramy dane, liczymy statystykę.\n",
    "7. Odrzucamy bądź nie $H_0$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-005",
   "metadata": {},
   "source": [
    "## Hipoteza zerowa i hipoteza alternatywna\n",
    "\n",
    "**Hipoteza zerowa** ($H_0$) — hipoteza poddana procedurze weryfikacyjnej. Zwykle zakłada, że różnica między analizowanymi parametrami lub rozkładami wynosi zero.\n",
    "\n",
    "Przykładowo, wnioskując o parametrach, hipotezę zerową możemy zapisać jako:\n",
    "\n",
    "$$\n",
    "H_0 : \\theta_{1} = \\theta_{2},\n",
    "$$\n",
    "\n",
    "gdzie $\\theta_1$ i $\\theta_2$ to parametry rozkładu populacji — na przykład średnie albo odchylenia standardowe.\n",
    "\n",
    "**Hipoteza alternatywna** ($H_1$ lub $H_A$) — hipoteza przeciwstawna do weryfikowanej. Możemy ją zapisać na trzy sposoby:\n",
    "\n",
    "- hipoteza dwustronna:\n",
    "\n",
    "$$\n",
    "H_{1}\\colon \\theta_{1}\\neq \\theta_{2}\n",
    "$$\n",
    "\n",
    "- hipoteza prawostronna:\n",
    "\n",
    "$$\n",
    "H_{1}\\colon \\theta_{1} > \\theta_{2}\n",
    "$$\n",
    "\n",
    "- hipoteza lewostronna:\n",
    "\n",
    "$$\n",
    "H_{1}\\colon \\theta_{1} < \\theta_{2}\n",
    "$$\n",
    "\n",
    "W zależności od sformułowania $H_A$ będziemy badać jeden ogon rozkładu (jednostronnie) albo dwa ogony (dwustronnie).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-006",
   "metadata": {},
   "source": [
    "## Błąd typu I i błąd typu II\n",
    "\n",
    "**Błąd pierwszego rodzaju** (alfa-błąd, ang. *false positive*) — błąd polegający na odrzuceniu hipotezy zerowej, która w rzeczywistości nie jest fałszywa. Prawdopodobieństwo popełnienia błędu pierwszego rodzaju oznaczamy symbolem $\\alpha$ i nazywamy poziomem istotności testu.\n",
    "\n",
    "**Błąd drugiego rodzaju** (beta-błąd, ang. *false negative*) — błąd polegający na nieodrzuceniu hipotezy zerowej, która jest w rzeczywistości fałszywa. Prawdopodobieństwo popełnienia błędu drugiego rodzaju oznaczamy symbolem $\\beta$.\n",
    "\n",
    "**Moc testu** wynosi $1-\\beta$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-007",
   "metadata": {},
   "source": [
    "**P-wartość** (ang. *p-value*) to prawdopodobieństwo otrzymania wyniku **co najmniej tak skrajnego jak zaobserwowany**, przy założeniu, że hipoteza zerowa $H_0$ jest prawdziwa.\n",
    "\n",
    "Ważne doprecyzowanie: p-wartość **nie** jest miarą wielkości efektu i **nie** jest „prawdopodobieństwem, że $H_0$ jest prawdziwa”.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-008",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.patches import Patch\n",
    "\n",
    "import seaborn as sns\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",
    "\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": "cell-009",
   "metadata": {},
   "source": [
    "## Testowanie hipotez: przykład I\n",
    "\n",
    "Wyobraźmy sobie, że interesuje nas częstość występowania **synestezji grafem–kolor** (np. litery lub cyfry wywołują stałe wrażenia barwne).\n",
    "\n",
    "Przyjmijmy, że w populacji ogólnej około $1.4\\%$ osób ma taką cechę.\n",
    "\n",
    "Badamy grupę $n=500$ osób, które od dzieciństwa intensywnie uczą się muzyki (np. w szkole muzycznej). Podejrzewamy, że w tej grupie synestezja może występować częściej.\n",
    "\n",
    "Niech $p$ oznacza prawdopodobieństwo synestezji w badanej grupie. Hipoteza alternatywna jest prawostronna:\n",
    "\n",
    "$$H_A: p > 0.014$$\n",
    "\n",
    "Hipoteza zerowa:\n",
    "\n",
    "$$H_0: p = 0.014$$\n",
    "\n",
    "W próbie okazało się, że $k=10$ osób ma synestezję. Czy uprawnia nas to do odrzucenia hipotezy zerowej?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-010",
   "metadata": {},
   "source": [
    "Najpierw zobrazujmy rozkład $\\mathrm{Bin}(500, 0.014)$ i zaznaczmy prawy ogon (bo test jest prawostronny).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-011",
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 500\n",
    "k_obs = 10\n",
    "p0 = 0.014\n",
    "\n",
    "k = np.arange(0, 26)\n",
    "pmf = stats.binom.pmf(k, n=n, p=p0)\n",
    "\n",
    "tail = k >= k_obs\n",
    "colors = np.where(tail, \"lightblue\", \"grey\")\n",
    "\n",
    "plt.figure(figsize=(9, 4))\n",
    "plt.bar(k, pmf, color=colors, edgecolor=\"white\")\n",
    "plt.title(\"Rozkład Bin(500, 0.014)\")\n",
    "plt.xlabel(\"Liczba sukcesów k\")\n",
    "plt.ylabel(\"P(X = k)\")\n",
    "\n",
    "legend_handles = [\n",
    "    Patch(facecolor=\"lightblue\", label=\"Prawdopodobieństwo 10 i więcej sukcesów\"),\n",
    "    Patch(facecolor=\"grey\", label=\"Pozostałe wyniki\"),\n",
    "]\n",
    "plt.legend(handles=legend_handles, loc=\"upper right\")\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-012",
   "metadata": {},
   "source": [
    "Teraz obliczamy dokładną p-wartość dla testu prawostronnego:\n",
    "\n",
    "$$p = P(X \\ge 10 \\mid H_0).$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-013",
   "metadata": {},
   "outputs": [],
   "source": [
    "p_value = stats.binom.sf(k_obs - 1, n=n, p=p0)\n",
    "\n",
    "print(\"Prawdopodobieństwo 10 i więcej sukcesów:\")\n",
    "print(p_value)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-014",
   "metadata": {},
   "source": [
    "W tym przykładzie p-wartość wynosi około $0.168$. Nie daje nam to powodu do odrzucenia $H_0$, jeżeli wybraliśmy poziom istotności $\\alpha = 0.05$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0572e7c0",
   "metadata": {},
   "source": [
    "## Interaktywne ilustracje (widżety)\n",
    "\n",
    "W tej części możesz zmieniać parametry testu i od razu obserwować, jak zmienia się rozkład dwumianowy oraz p-wartość.\n",
    "\n",
    "W wielu prostych eksperymentach percepcyjnych uczestnik ma do wyboru dwie odpowiedzi (np. „A” albo „B”). Taki schemat nazywa się **zadaniem z dwiema alternatywami** (ang. *two-alternative forced choice*, skrót 2AFC).\n",
    "\n",
    "Jeśli w każdej próbie uczestnik zgaduje, to prawdopodobieństwo odpowiedzi poprawnej wynosi $p_0=0.5$, a liczba poprawnych odpowiedzi $K$ w $n$ próbach ma rozkład $\\mathrm{Bin}(n, p_0)$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b977b14",
   "metadata": {},
   "source": [
    "### Widżet 1: p-wartość jako „pole ogona”\n",
    "\n",
    "Wybierz $n$, $p_0$ i obserwowany wynik $k$ (liczbę „sukcesów”, np. poprawnych odpowiedzi). Zobaczysz, które wyniki wchodzą do p-wartości w zależności od tego, czy test jest jedno- czy dwustronny.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a95b588b",
   "metadata": {},
   "outputs": [],
   "source": [
    "if not WIDGETY_DOSTEPNE:\n",
    "    print('Widżety nie są dostępne: zainstaluj pakiet ipywidgets, aby korzystać z części interaktywnych.')\n",
    "else:\n",
    "    n_slider = widgets.IntSlider(\n",
    "        value=16,\n",
    "        min=5,\n",
    "        max=200,\n",
    "        step=1,\n",
    "        description=\"n\",\n",
    "        continuous_update=False,\n",
    "    )\n",
    "    p0_slider = widgets.FloatSlider(\n",
    "        value=0.50,\n",
    "        min=0.05,\n",
    "        max=0.95,\n",
    "        step=0.01,\n",
    "        description=\"p0\",\n",
    "        readout_format=\".2f\",\n",
    "        continuous_update=False,\n",
    "    )\n",
    "    k_slider = widgets.IntSlider(\n",
    "        value=13,\n",
    "        min=0,\n",
    "        max=n_slider.value,\n",
    "        step=1,\n",
    "        description=\"k\",\n",
    "        continuous_update=False,\n",
    "    )\n",
    "    alpha_slider = widgets.FloatSlider(\n",
    "        value=0.05,\n",
    "        min=0.001,\n",
    "        max=0.20,\n",
    "        step=0.001,\n",
    "        description=\"α\",\n",
    "        readout_format=\".3f\",\n",
    "        continuous_update=False,\n",
    "    )\n",
    "    alternative_dropdown = widgets.Dropdown(\n",
    "        options=[\n",
    "            (\"p > p0 (prawostronnie)\", \"greater\"),\n",
    "            (\"p < p0 (lewostronnie)\", \"less\"),\n",
    "            (\"p ≠ p0 (dwustronnie)\", \"two-sided\"),\n",
    "        ],\n",
    "        value=\"greater\",\n",
    "        description=\"Hipoteza\",\n",
    "    )\n",
    "\n",
    "    def aktualizuj_k_max(change=None):\n",
    "        k_slider.max = int(n_slider.value)\n",
    "        if k_slider.value > k_slider.max:\n",
    "            k_slider.value = k_slider.max\n",
    "\n",
    "    n_slider.observe(aktualizuj_k_max, names=\"value\")\n",
    "    aktualizuj_k_max()\n",
    "\n",
    "    def pokaz_test_dwumianowy(n, p0, k_obs, alpha, alternative):\n",
    "        n = int(n)\n",
    "        k_obs = int(k_obs)\n",
    "        k_all = np.arange(0, n + 1)\n",
    "        pmf_all = stats.binom.pmf(k_all, n=n, p=p0)\n",
    "\n",
    "        if alternative == \"greater\":\n",
    "            opis = \"p > p0 (prawostronnie)\"\n",
    "            mask = k_all >= k_obs\n",
    "            p_value = stats.binom.sf(k_obs - 1, n=n, p=p0)\n",
    "        elif alternative == \"less\":\n",
    "            opis = \"p < p0 (lewostronnie)\"\n",
    "            mask = k_all <= k_obs\n",
    "            p_value = stats.binom.cdf(k_obs, n=n, p=p0)\n",
    "        else:\n",
    "            opis = \"p ≠ p0 (dwustronnie)\"\n",
    "            pmf_obs = float(pmf_all[k_obs])\n",
    "            mask = pmf_all <= pmf_obs + 1e-12\n",
    "            p_value = float(pmf_all[mask].sum())\n",
    "\n",
    "        p_hat = k_obs / n\n",
    "        decyzja = \"Odrzucamy H0\" if float(p_value) < float(alpha) else \"Nie odrzucamy H0\"\n",
    "\n",
    "        print(f\"n = {n}, k = {k_obs}, p0 = {p0:.2f}\")\n",
    "        print(f\"Estymata z próby: p̂ = {p_hat:.3f}\")\n",
    "        print(f\"p-wartość = {float(p_value):.6f}\")\n",
    "        print(f\"Decyzja przy α={float(alpha):.3f}: {decyzja}\")\n",
    "\n",
    "        colors = np.where(mask, \"skyblue\", \"lightgrey\")\n",
    "        colors[k_obs] = \"darkblue\"\n",
    "\n",
    "        plt.figure(figsize=(10, 4))\n",
    "        plt.bar(k_all, pmf_all, color=colors, edgecolor=\"white\")\n",
    "        plt.title(f\"Rozkład Bin({n}, {p0:.2f}) i p-wartość ({opis})\")\n",
    "        plt.xlabel(\"Liczba sukcesów k\")\n",
    "        plt.ylabel(\"P(X = k)\")\n",
    "        plt.tight_layout()\n",
    "        plt.show()\n",
    "\n",
    "    out = widgets.interactive_output(\n",
    "        pokaz_test_dwumianowy,\n",
    "        {\n",
    "            \"n\": n_slider,\n",
    "            \"p0\": p0_slider,\n",
    "            \"k_obs\": k_slider,\n",
    "            \"alpha\": alpha_slider,\n",
    "            \"alternative\": alternative_dropdown,\n",
    "        },\n",
    "    )\n",
    "\n",
    "    display(\n",
    "        widgets.VBox([n_slider, p0_slider, k_slider, alpha_slider, alternative_dropdown]),\n",
    "        out,\n",
    "    )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b81406b",
   "metadata": {},
   "source": [
    "### Widżet 2: błąd I rodzaju „w praktyce”\n",
    "\n",
    "W tej symulacji zakładamy, że $H_0$ jest prawdziwa (w zadaniu 2AFC: $p_0=0.5$), a następnie powtarzamy „całe badanie” wiele razy. Zobaczysz, jak często mimo to otrzymamy wynik „istotny” (p-wartość mniejsza niż $\\alpha$).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a217e3c",
   "metadata": {},
   "outputs": [],
   "source": [
    "if not WIDGETY_DOSTEPNE:\n",
    "    print('Widżety nie są dostępne: zainstaluj pakiet ipywidgets, aby korzystać z części interaktywnych.')\n",
    "else:\n",
    "    n_slider = widgets.IntSlider(\n",
    "        value=40,\n",
    "        min=5,\n",
    "        max=200,\n",
    "        step=1,\n",
    "        description=\"n\",\n",
    "        continuous_update=False,\n",
    "    )\n",
    "    alpha_slider = widgets.FloatSlider(\n",
    "        value=0.05,\n",
    "        min=0.001,\n",
    "        max=0.20,\n",
    "        step=0.001,\n",
    "        description=\"α\",\n",
    "        readout_format=\".3f\",\n",
    "        continuous_update=False,\n",
    "    )\n",
    "    repeats_slider = widgets.IntSlider(\n",
    "        value=2000,\n",
    "        min=100,\n",
    "        max=20000,\n",
    "        step=100,\n",
    "        description=\"Powtórzenia\",\n",
    "        continuous_update=False,\n",
    "    )\n",
    "    seed_input = widgets.IntText(value=SEED, description=\"Ziarno\")\n",
    "    alternative_dropdown = widgets.Dropdown(\n",
    "        options=[\n",
    "            (\"p > p0 (prawostronnie)\", \"greater\"),\n",
    "            (\"p < p0 (lewostronnie)\", \"less\"),\n",
    "            (\"p ≠ p0 (dwustronnie)\", \"two-sided\"),\n",
    "        ],\n",
    "        value=\"two-sided\",\n",
    "        description=\"Hipoteza\",\n",
    "    )\n",
    "\n",
    "    def tabela_p_wartosci(n: int, p0: float, alternative: str) -> np.ndarray:\n",
    "        k_values = np.arange(0, n + 1)\n",
    "\n",
    "        if alternative == \"greater\":\n",
    "            return np.asarray(stats.binom.sf(k_values - 1, n=n, p=p0), dtype=float)\n",
    "        if alternative == \"less\":\n",
    "            return np.asarray(stats.binom.cdf(k_values, n=n, p=p0), dtype=float)\n",
    "\n",
    "        pmf = np.asarray(stats.binom.pmf(k_values, n=n, p=p0), dtype=float)\n",
    "        eps = 1e-12\n",
    "        p_table = np.empty(n + 1, dtype=float)\n",
    "        for k_obs in k_values:\n",
    "            p_table[int(k_obs)] = float(pmf[pmf <= pmf[int(k_obs)] + eps].sum())\n",
    "        return p_table\n",
    "\n",
    "    def symuluj_blad_I_rodzaju(n, alpha, repeats, seed, alternative):\n",
    "        n = int(n)\n",
    "        repeats = int(repeats)\n",
    "        seed = int(seed)\n",
    "        p0 = 0.5\n",
    "\n",
    "        rng_local = np.random.default_rng(seed)\n",
    "        k_samples = rng_local.binomial(n=n, p=p0, size=repeats)\n",
    "        p_table = tabela_p_wartosci(n=n, p0=p0, alternative=alternative)\n",
    "        p_values = p_table[k_samples]\n",
    "\n",
    "        odsetek = float((p_values < float(alpha)).mean())\n",
    "        print(f\"Założenie symulacji: H0 jest prawdziwa (p0 = {p0:.1f}).\")\n",
    "        print(f\"Powtórzenia: {repeats}\")\n",
    "        print(f\"Odsetek przypadków z p < α: {odsetek:.3f}\")\n",
    "\n",
    "        plt.figure(figsize=(9, 4))\n",
    "        plt.hist(p_values, bins=20, range=(0, 1), color=\"lightgrey\", edgecolor=\"white\")\n",
    "        plt.axvline(float(alpha), color=\"crimson\", linestyle=\"--\", linewidth=2, label=\"α\")\n",
    "        plt.title(\"Rozkład p-wartości, gdy H0 jest prawdziwa\")\n",
    "        plt.xlabel(\"p-wartość\")\n",
    "        plt.ylabel(\"liczba symulacji\")\n",
    "        plt.legend()\n",
    "        plt.tight_layout()\n",
    "        plt.show()\n",
    "\n",
    "    out = widgets.interactive_output(\n",
    "        symuluj_blad_I_rodzaju,\n",
    "        {\n",
    "            \"n\": n_slider,\n",
    "            \"alpha\": alpha_slider,\n",
    "            \"repeats\": repeats_slider,\n",
    "            \"seed\": seed_input,\n",
    "            \"alternative\": alternative_dropdown,\n",
    "        },\n",
    "    )\n",
    "\n",
    "    display(\n",
    "        widgets.VBox([n_slider, alpha_slider, repeats_slider, seed_input, alternative_dropdown]),\n",
    "        out,\n",
    "    )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-015",
   "metadata": {},
   "source": [
    "## Testowanie hipotez: przykład II\n",
    "\n",
    "Interesuje nas, czy niedobór snu wiąże się z częstszym występowaniem **bardzo niskich wyników** w prostym teście uwagi u dzieci w wieku 8 lat.\n",
    "\n",
    "Za „bardzo niski wynik” uznajemy wynik poniżej ustalonego progu. W populacji dzieci w tym wieku około $5\\%$ osób znajduje się poniżej tego progu.\n",
    "\n",
    "Wybrano losową próbę $n=66$ dzieci, które przez dwa tygodnie (według dzienniczka snu) spały przeciętnie mniej niż 8 godzin na dobę. Okazało się, że $k=7$ z nich uzyskało wynik poniżej progu.\n",
    "\n",
    "$$H_A: p > 0.05$$\n",
    "$$H_0: p = 0.05$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-016",
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 66\n",
    "k_obs = 7\n",
    "p0 = 0.05\n",
    "\n",
    "k = np.arange(0, 16)\n",
    "pmf = stats.binom.pmf(k, n=n, p=p0)\n",
    "colors = np.where(k >= k_obs, \"lightblue\", \"grey\")\n",
    "\n",
    "plt.figure(figsize=(9, 4))\n",
    "plt.bar(k, pmf, color=colors, edgecolor=\"white\")\n",
    "plt.title(\"Rozkład Bin(66, 0.05)\")\n",
    "plt.xlabel(\"Liczba sukcesów k\")\n",
    "plt.ylabel(\"P(X = k)\")\n",
    "\n",
    "legend_handles = [\n",
    "    Patch(facecolor=\"lightblue\", label=\"Prawdopodobieństwo 7 i więcej sukcesów\"),\n",
    "    Patch(facecolor=\"grey\", label=\"Pozostałe wyniki\"),\n",
    "]\n",
    "plt.legend(handles=legend_handles, loc=\"upper right\")\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-017",
   "metadata": {},
   "outputs": [],
   "source": [
    "p_value = stats.binom.sf(k_obs - 1, n=n, p=p0)\n",
    "print(\"Prawdopodobieństwo 7 i więcej sukcesów:\")\n",
    "print(float(p_value))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-018",
   "metadata": {},
   "source": [
    "Tym razem p-wartość wynosi około $0.046$. Wartość ta uprawnia nas do odrzucenia hipotezy zerowej przy $\\alpha = 0.05$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-019",
   "metadata": {},
   "source": [
    "## Testowanie hipotez: przykład III — *Bond tasting Martini*\n",
    "\n",
    "James Bond słynny jest z tego, że zawsze zamawia Martini wstrząśnięte, nie zmieszane. Czy jednak rzeczywiście potrafi odróżnić te dwie metody przygotowania swojego ulubionego drinku?\n",
    "\n",
    "Q zaprojektował prosty eksperyment: podczas $n=16$ losowych prób Bond ma orzec, czy Martini było wstrząśnięte, czy zmieszane.\n",
    "Gdyby wybierał całkowicie losowo, miałby $50\\%$ szans na poprawne zgadnięcie.\n",
    "\n",
    "Eksperyment pokazał, że Bond zgadł w $k=13$ próbach na 16.\n",
    "\n",
    "1) Czy uprawnia nas to do twierdzenia, że odpowiedział lepiej, niż gdyby odpowiadał losowo?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-020",
   "metadata": {},
   "source": [
    "Hipotezy (prawostronnie):\n",
    "\n",
    "$$H_A: p > 0.50$$\n",
    "$$H_0: p = 0.50$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-021",
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 16\n",
    "k_obs = 13\n",
    "p0 = 0.5\n",
    "\n",
    "k = np.arange(0, n + 1)\n",
    "pmf = stats.binom.pmf(k, n=n, p=p0)\n",
    "colors = np.where(k >= k_obs, \"lightblue\", \"grey\")\n",
    "\n",
    "plt.figure(figsize=(9, 4))\n",
    "plt.bar(k, pmf, color=colors, edgecolor=\"white\")\n",
    "plt.title(\"Rozkład Bin(16, 0.5)\")\n",
    "plt.xlabel(\"Liczba sukcesów k\")\n",
    "plt.ylabel(\"P(X = k)\")\n",
    "\n",
    "legend_handles = [\n",
    "    Patch(facecolor=\"lightblue\", label=\"Prawdopodobieństwo 13 i więcej sukcesów\"),\n",
    "    Patch(facecolor=\"grey\", label=\"Pozostałe wyniki\"),\n",
    "]\n",
    "plt.legend(handles=legend_handles, loc=\"upper right\")\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-022",
   "metadata": {},
   "outputs": [],
   "source": [
    "p_value_hand = stats.binom.sf(k_obs - 1, n=n, p=p0)\n",
    "print(\"Prawdopodobieństwo 13 i więcej sukcesów (na piechotę):\")\n",
    "print(p_value_hand)\n",
    "\n",
    "result = stats.binomtest(k_obs, n=n, p=p0, alternative=\"greater\")\n",
    "print(\"Sprawdzenie funkcją binomtest:\")\n",
    "print(result.pvalue)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-023",
   "metadata": {},
   "source": [
    "2) Możemy jednak sformułować pytanie badawcze inaczej: czy uzyskany wynik uprawnia nas do twierdzenia, że Bond odpowiadał **lepiej lub gorzej** niż losowo?\n",
    "\n",
    "$$H_A: p \\neq 0.50$$\n",
    "$$H_0: p = 0.50$$\n",
    "\n",
    "W przypadku hipotezy dwustronnej dzielimy poziom istotności $\\alpha$ na dwa ogony (np. przy $\\alpha=0.05$ po $2.5\\%$ na ogon).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-024",
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = 0.05\n",
    "crit_low = int(stats.binom.ppf(alpha / 2, n=n, p=p0))\n",
    "crit_high = int(stats.binom.ppf(1 - alpha / 2, n=n, p=p0))\n",
    "\n",
    "print(\"Wartości krytyczne:\")\n",
    "print(crit_low, crit_high)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-025",
   "metadata": {},
   "outputs": [],
   "source": [
    "k = np.arange(0, n + 1)\n",
    "pmf = stats.binom.pmf(k, n=n, p=p0)\n",
    "colors = np.where((k <= crit_low) | (k >= crit_high), \"lightblue\", \"grey\")\n",
    "\n",
    "plt.figure(figsize=(9, 4))\n",
    "plt.bar(k, pmf, color=colors, edgecolor=\"white\")\n",
    "plt.title(\"Rozkład Bin(16, 0.5) — hipoteza dwustronna\")\n",
    "plt.xlabel(\"Liczba sukcesów k\")\n",
    "plt.ylabel(\"P(X = k)\")\n",
    "\n",
    "legend_handles = [\n",
    "    Patch(facecolor=\"lightblue\", label=\"Obszar krytyczny (dwustronnie)\"),\n",
    "    Patch(facecolor=\"grey\", label=\"Poza obszarem krytycznym\"),\n",
    "]\n",
    "plt.legend(handles=legend_handles, loc=\"upper right\")\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-026",
   "metadata": {},
   "outputs": [],
   "source": [
    "p_right = stats.binom.sf(k_obs - 1, n=n, p=p0)\n",
    "p_left = stats.binom.cdf(k_obs, n=n, p=p0)\n",
    "p_two_naive = 2 * min(p_left, p_right)\n",
    "\n",
    "p_two = stats.binomtest(k_obs, n=n, p=p0, alternative=\"two-sided\").pvalue\n",
    "\n",
    "print(\"p (2 * min ogonów) =\", p_two_naive)\n",
    "print(\"p (binomtest, two-sided) =\", p_two)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-027",
   "metadata": {},
   "source": [
    "### Dygresja: przedział ufności dla $p$\n",
    "\n",
    "Szacowane z próby prawdopodobieństwo sukcesu:\n",
    "\n",
    "$$\\hat{p} = \\frac{k}{n}.$$ \n",
    "\n",
    "Jak możemy skonstruować wokół $\\hat{p}$ przedział ufności? Poniżej kilka metod.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-028",
   "metadata": {},
   "outputs": [],
   "source": [
    "p_hat = k_obs / n\n",
    "print(\"p_hat =\", p_hat)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-029",
   "metadata": {},
   "source": [
    "#### (1) Przybliżenie normalne\n",
    "\n",
    "To podejście jest proste, ale bywa niedokładne.\n",
    "\n",
    "$$\\hat{p} \\pm z_{0.975}\\,\\sqrt{\\frac{\\hat{p}(1-\\hat{p})}{n}}$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-030",
   "metadata": {},
   "outputs": [],
   "source": [
    "z_0975 = stats.norm.ppf(0.975)\n",
    "se_hat = np.sqrt(p_hat * (1 - p_hat) / n)\n",
    "\n",
    "ci_low = p_hat - z_0975 * se_hat\n",
    "ci_high = p_hat + z_0975 * se_hat\n",
    "\n",
    "print(\"CI (normal approx):\", ci_low, ci_high)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-031",
   "metadata": {},
   "source": [
    "#### (2) Przedział dokładny (Clopper–Pearson) i Wilson\n",
    "\n",
    "Możemy skorzystać z gotowej funkcji `binomtest`, która zwraca też przedział ufności dla $p$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-032",
   "metadata": {},
   "outputs": [],
   "source": [
    "res = stats.binomtest(k_obs, n=n, p=p0, alternative=\"two-sided\")\n",
    "ci_exact = res.proportion_ci(confidence_level=0.95, method=\"exact\")\n",
    "ci_wilson = res.proportion_ci(confidence_level=0.95, method=\"wilson\")\n",
    "\n",
    "print(\"CI (exact):\", ci_exact.low, ci_exact.high)\n",
    "print(\"CI (Wilson):\", ci_wilson.low, ci_wilson.high)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-033",
   "metadata": {},
   "source": [
    "#### (3) Symulacja Monte Carlo\n",
    "\n",
    "Losujemy wiele razy dane z $\\mathrm{Bin}(n, \\hat{p})$, liczymy proporcję i bierzemy kwantyle.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-034",
   "metadata": {},
   "outputs": [],
   "source": [
    "repeats = 200_000\n",
    "rng_ci = np.random.default_rng(SEED + 1)\n",
    "\n",
    "sim_k = rng_ci.binomial(n=n, p=p_hat, size=repeats)\n",
    "sim_phat = sim_k / n\n",
    "\n",
    "ci_mc_low, ci_mc_high = np.quantile(sim_phat, [0.025, 0.975])\n",
    "print(\"CI (Monte Carlo):\", ci_mc_low, ci_mc_high)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-035",
   "metadata": {},
   "source": [
    "## Testowanie hipotez — p-wartość przy niesymetrycznych rozkładach\n",
    "\n",
    "Załóżmy, że przeciętny student ma $80\\%$ szans na zdanie egzaminu. Wybraliśmy próbę $n=50$ studentów kognitywistyki, przeegzaminowaliśmy i okazało się, że $k=35$ z nich zdało egzamin.\n",
    "\n",
    "Pytanie jest dwustronne:\n",
    "\n",
    "$$H_A: p \\neq 0.80$$\n",
    "$$H_0: p = 0.80$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-036",
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 50\n",
    "k_obs = 35\n",
    "p0 = 0.8\n",
    "alpha = 0.05\n",
    "\n",
    "crit_low = stats.binom.ppf(alpha / 2, n=n, p=p0)\n",
    "crit_high = stats.binom.ppf(1 - alpha / 2, n=n, p=p0)\n",
    "\n",
    "print(\"Wartości krytyczne:\")\n",
    "print(crit_low, crit_high)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-037",
   "metadata": {},
   "outputs": [],
   "source": [
    "k = np.arange(0, n + 1)\n",
    "pmf = stats.binom.pmf(k, n=n, p=p0)\n",
    "colors = np.where((k <= crit_low) | (k >= crit_high), \"lightblue\", \"grey\")\n",
    "\n",
    "plt.figure(figsize=(10, 4))\n",
    "plt.bar(k, pmf, color=colors, edgecolor=\"white\")\n",
    "plt.title(\"Rozkład Bin(50, 0.8) — hipoteza dwustronna\")\n",
    "plt.xlabel(\"Liczba sukcesów k\")\n",
    "plt.ylabel(\"P(X = k)\")\n",
    "\n",
    "legend_handles = [\n",
    "    Patch(facecolor=\"lightblue\", label=\"Obszar krytyczny\"),\n",
    "    Patch(facecolor=\"grey\", label=\"Poza obszarem krytycznym\"),\n",
    "]\n",
    "plt.legend(handles=legend_handles, loc=\"upper left\")\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-038",
   "metadata": {},
   "outputs": [],
   "source": [
    "p_left = stats.binom.cdf(k_obs, n=n, p=p0)\n",
    "p_right = stats.binom.sf(k_obs - 1, n=n, p=p0)\n",
    "\n",
    "print(\"P(X ≤ 35) =\", p_left)\n",
    "print(\"P(X ≥ 35) =\", p_right)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-039",
   "metadata": {},
   "outputs": [],
   "source": [
    "p_two_naive = 2 * min(p_left, p_right)\n",
    "p_two = stats.binomtest(k_obs, n=n, p=p0, alternative=\"two-sided\").pvalue\n",
    "\n",
    "k_all = np.arange(0, n + 1)\n",
    "pmf_all = stats.binom.pmf(k_all, n=n, p=p0)\n",
    "pmf_obs = pmf_all[k_obs]\n",
    "p_two_manual = pmf_all[pmf_all <= pmf_obs].sum()\n",
    "\n",
    "print(\"2 * min(ogonów) =\", p_two_naive)\n",
    "print(\"binomtest (two-sided) =\", p_two)\n",
    "print(\"manual (suma P(X=k) ≤ P(X=35)) =\", p_two_manual)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-040",
   "metadata": {},
   "source": [
    "Żeby lepiej zobrazować, o co chodzi, możemy zaznaczyć wszystkie wyniki o prawdopodobieństwie nie większym niż prawdopodobieństwo wyniku obserwowanego.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-041",
   "metadata": {},
   "outputs": [],
   "source": [
    "colors = np.where(pmf_all <= pmf_obs + 1e-12, \"skyblue\", \"lightgrey\")\n",
    "colors[k_obs] = \"darkblue\"\n",
    "\n",
    "plt.figure(figsize=(10, 4))\n",
    "plt.bar(k_all, pmf_all, color=colors, edgecolor=\"white\")\n",
    "plt.title(\"Rozkład Bin(50, 0.8) — wyniki ≥ tak skrajne jak k=35\")\n",
    "plt.xlabel(\"Liczba sukcesów k\")\n",
    "plt.ylabel(\"P(X = k)\")\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fee2f2e2-5741-452b-a1d9-44e36dce465a",
   "metadata": {},
   "outputs": [],
   "source": [
    "stats.binomtest(35, 50, 0.8, alternative=\"less\")\n",
    "# stats.binomtest(35, 50, 0.8, alternative=\"two-sided\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-042",
   "metadata": {},
   "source": [
    "## Błąd standardowy średniej i centralne twierdzenie graniczne (CTG)\n",
    "\n",
    "**Błąd standardowy** — odchylenie standardowe średnich z prób.\n",
    "\n",
    "Jeśli odchylenie standardowe w populacji wynosi $\\sigma$, to:\n",
    "\n",
    "$$SE = \\frac{\\sigma}{\\sqrt{n}}.$$\n",
    "\n",
    "**Centralne twierdzenie graniczne (jedno ze sformułowań)**: dla populacji o średniej $\\mu$ i wariancji $\\sigma^2$ rozkład średniej z próby ma średnią $\\mu$ i wariancję $\\sigma^2/n$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-043",
   "metadata": {},
   "outputs": [],
   "source": [
    "def blad_standardowy_normal(n: int, mean: float = 0.0, sd: float = 1.0, repeats: int = 1000, seed: int = 1234) -> None:\n",
    "    \"\"\"Symulacja CTG: populacja normalna, rozkład średnich i QQ-plot.\"\"\"\n",
    "    rng_local = np.random.default_rng(seed)\n",
    "\n",
    "    sample_pop = rng_local.normal(loc=mean, scale=sd, size=100)\n",
    "    samples = rng_local.normal(loc=mean, scale=sd, size=(repeats, n))\n",
    "    means = samples.mean(axis=1)\n",
    "\n",
    "    fig, axes = plt.subplots(1, 3, figsize=(13, 3))\n",
    "\n",
    "    axes[0].hist(sample_pop, density=True, color=\"lightgrey\", edgecolor=\"white\")\n",
    "    x = np.linspace(mean - 4 * sd, mean + 4 * sd, 400)\n",
    "    axes[0].plot(x, stats.norm.pdf(x, loc=mean, scale=sd), color=\"black\")\n",
    "    axes[0].set_title(\"(normalny) Rozkład populacji\")\n",
    "    axes[0].set_xlabel(\"x\")\n",
    "    axes[0].set_ylabel(\"gęstość\")\n",
    "\n",
    "    se = sd / np.sqrt(n)\n",
    "    axes[1].hist(means, density=True, color=\"lightgrey\", edgecolor=\"white\")\n",
    "    x2 = np.linspace(means.min(), means.max(), 400)\n",
    "    axes[1].plot(x2, stats.norm.pdf(x2, loc=mean, scale=se), color=\"red\", linewidth=2)\n",
    "    axes[1].set_title(\"Rozkład średniej z próby\")\n",
    "    axes[1].set_xlabel(\"$\\\\bar{X}$\")\n",
    "\n",
    "    stats.probplot(means, dist=\"norm\", plot=axes[2])\n",
    "    axes[2].set_title(\"Wykres kwantyl-kwantyl\")\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-044",
   "metadata": {},
   "outputs": [],
   "source": [
    "blad_standardowy_normal(35, repeats=1000, seed=1234)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-045",
   "metadata": {},
   "outputs": [],
   "source": [
    "def blad_standardowy_unif(n: int, a: float, b: float, repeats: int = 1000, seed: int = 1234) -> None:\n",
    "    \"\"\"Symulacja CTG: populacja jednostajna U(a,b), rozkład średnich i QQ-plot.\"\"\"\n",
    "    rng_local = np.random.default_rng(seed)\n",
    "\n",
    "    sample_pop = rng_local.uniform(a, b, size=100)\n",
    "    samples = rng_local.uniform(a, b, size=(repeats, n))\n",
    "    means = samples.mean(axis=1)\n",
    "\n",
    "    pop_mean = (a + b) / 2\n",
    "    pop_sd = (b - a) / np.sqrt(12)\n",
    "    se = pop_sd / np.sqrt(n)\n",
    "\n",
    "    fig, axes = plt.subplots(1, 3, figsize=(13, 3))\n",
    "\n",
    "    axes[0].hist(sample_pop, density=True, color=\"lightgrey\", edgecolor=\"white\")\n",
    "    x = np.linspace(a, b, 400)\n",
    "    axes[0].plot(x, stats.uniform.pdf(x, loc=a, scale=b - a), color=\"black\")\n",
    "    axes[0].set_title(\"(jednostajny) Rozkład populacji\")\n",
    "    axes[0].set_xlabel(\"x\")\n",
    "    axes[0].set_ylabel(\"gęstość\")\n",
    "\n",
    "    axes[1].hist(means, density=True, color=\"lightgrey\", edgecolor=\"white\")\n",
    "    x2 = np.linspace(means.min(), means.max(), 400)\n",
    "    axes[1].plot(x2, stats.norm.pdf(x2, loc=pop_mean, scale=se), color=\"blue\", linewidth=2)\n",
    "    axes[1].set_title(\"Rozkład średniej z próby\")\n",
    "    axes[1].set_xlabel(\"$\\\\bar{X}$\")\n",
    "\n",
    "    stats.probplot(means, dist=\"norm\", plot=axes[2])\n",
    "    axes[2].set_title(\"Wykres kwantyl-kwantyl\")\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-046",
   "metadata": {},
   "outputs": [],
   "source": [
    "blad_standardowy_unif(35, a=0, b=100, repeats=1000, seed=1234)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-047",
   "metadata": {},
   "source": [
    "## Testowanie hipotez: przykład IV (test $z$ dla średniej)\n",
    "\n",
    "Rozważmy **standaryzowany test uwagi**, którego wyniki zostały przeskalowane tak, aby w populacji miały w przybliżeniu rozkład normalny o średniej $\\mu_0 = 50$ i odchyleniu standardowym $\\sigma = 10$ (np. skala typu T).\n",
    "\n",
    "Zebrano próbę $n=166$ osób z pewnej badanej grupy. Średnia w próbie wyniosła $\\bar{X}=55.71$.\n",
    "\n",
    "W tym przykładzie traktujemy $\\sigma=10$ jako znane z normalizacji testu. Dzięki temu możemy użyć testu $z$ dla średniej i policzyć p-wartość na podstawie rozkładu normalnego.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-048",
   "metadata": {},
   "outputs": [],
   "source": [
    "mu0 = 50.0\n",
    "sigma = 10.0\n",
    "n = 166\n",
    "xbar_obs = 55.71\n",
    "\n",
    "se = sigma / np.sqrt(n)\n",
    "print(\"SE =\", se)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-049",
   "metadata": {},
   "source": [
    "Zobaczmy rozkład $\\bar{X}$ przy założeniu prawdziwości $H_0$. Dla ilustracji losujemy 1000 prób po 166 obserwacji z $\\mathcal{N}(50, 10)$ i liczymy średnie.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-050",
   "metadata": {},
   "outputs": [],
   "source": [
    "rng_sim = np.random.default_rng(1234)\n",
    "sim_samples = rng_sim.normal(loc=mu0, scale=sigma, size=(1000, n))\n",
    "sim_means = sim_samples.mean(axis=1)\n",
    "\n",
    "fig, axes = plt.subplots(1, 2, figsize=(10, 4))\n",
    "\n",
    "axes[0].hist(sim_means, bins=25, density=True, color=\"lightgrey\", edgecolor=\"white\")\n",
    "x = np.linspace(sim_means.min(), sim_means.max(), 400)\n",
    "axes[0].plot(x, stats.norm.pdf(x, loc=mu0, scale=se), color=\"black\")\n",
    "axes[0].set_title(\"Rozkład $\\\\bar{X}$ przy $H_0$\")\n",
    "axes[0].set_xlabel(\"$\\\\bar{X}$\")\n",
    "axes[0].set_ylabel(\"gęstość\")\n",
    "\n",
    "stats.probplot(sim_means, dist=\"norm\", plot=axes[1])\n",
    "axes[1].set_title(\"QQ-plot średnich\")\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-051",
   "metadata": {},
   "source": [
    "Wyznaczmy wartości krytyczne dla hipotezy dwustronnej i prawostronnej.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-052",
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = 0.05\n",
    "\n",
    "# Dwustronnie\n",
    "crit_l = stats.norm.ppf(alpha / 2, loc=mu0, scale=se)\n",
    "crit_u = stats.norm.ppf(1 - alpha / 2, loc=mu0, scale=se)\n",
    "\n",
    "# Prawostronnie\n",
    "crit_u_right = stats.norm.ppf(1 - alpha, loc=mu0, scale=se)\n",
    "\n",
    "print(\"Wartości krytyczne (dwustronnie):\", crit_l, crit_u)\n",
    "print(\"Wartość krytyczna (prawostronnie):\", crit_u_right)\n",
    "\n",
    "print(\"(Dla porównania) kwantyle z symulacji:\")\n",
    "print(np.quantile(sim_means, [0.025, 0.975]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-053",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_normal_critical(mu: float, se: float, alpha: float, observed: float, two_sided: bool) -> None:\n",
    "    x = np.linspace(mu - 5 * se, mu + 5 * se, 600)\n",
    "    y = stats.norm.pdf(x, loc=mu, scale=se)\n",
    "\n",
    "    plt.figure(figsize=(9, 4))\n",
    "    plt.plot(x, y, color=\"black\")\n",
    "\n",
    "    if two_sided:\n",
    "        l = stats.norm.ppf(alpha / 2, loc=mu, scale=se)\n",
    "        u = stats.norm.ppf(1 - alpha / 2, loc=mu, scale=se)\n",
    "        plt.fill_between(x[x <= l], y[x <= l], color=\"skyblue\", alpha=0.4)\n",
    "        plt.fill_between(x[x >= u], y[x >= u], color=\"skyblue\", alpha=0.4)\n",
    "        plt.axvline(l, color=\"grey\", linestyle=\":\")\n",
    "        plt.axvline(u, color=\"grey\", linestyle=\":\")\n",
    "        plt.title(\"Obszary krytyczne (dwustronnie)\")\n",
    "    else:\n",
    "        u = stats.norm.ppf(1 - alpha, loc=mu, scale=se)\n",
    "        plt.fill_between(x[x >= u], y[x >= u], color=\"skyblue\", alpha=0.4)\n",
    "        plt.axvline(u, color=\"grey\", linestyle=\":\")\n",
    "        plt.title(\"Obszar krytyczny (prawostronnie)\")\n",
    "\n",
    "    plt.axvline(observed, color=\"crimson\", linestyle=\"--\", linewidth=2)\n",
    "    plt.xlabel(\"$\\\\bar{X}$\")\n",
    "    plt.ylabel(\"gęstość\")\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-054",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_normal_critical(mu=mu0, se=se, alpha=0.05, observed=xbar_obs, two_sided=False)\n",
    "plot_normal_critical(mu=mu0, se=se, alpha=0.05, observed=xbar_obs, two_sided=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-055",
   "metadata": {},
   "outputs": [],
   "source": [
    "p_right = stats.norm.sf(xbar_obs, loc=mu0, scale=se)\n",
    "p_two = 2 * stats.norm.sf(abs((xbar_obs - mu0) / se))\n",
    "\n",
    "print(\"p (prawostronnie) =\", float(p_right))\n",
    "print(\"p (dwustronnie) =\", float(p_two))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-056",
   "metadata": {},
   "source": [
    "W ogólności takie postępowanie nazywamy **testem $z$**:\n",
    "\n",
    "$$\n",
    "z = \\frac{\\bar{X} - \\mu_0}{\\sigma/\\sqrt{n}}.\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-057",
   "metadata": {},
   "outputs": [],
   "source": [
    "z_stat = (xbar_obs - mu0) / se\n",
    "p_two_z = 2 * stats.norm.sf(abs(z_stat))\n",
    "\n",
    "print(\"z =\", z_stat)\n",
    "print(\"p (dwustronnie) =\", p_two_z)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-058",
   "metadata": {},
   "source": [
    "## Aproksymacja rozkładu dwumianowego rozkładem normalnym\n",
    "\n",
    "Rozkład normalny jest niezłym przybliżeniem rozkładu dwumianowego. Możemy się o tym przekonać wykorzystując:\n",
    "\n",
    "$$\\mu = np$$\n",
    "$$\\sigma = \\sqrt{np(1-p)}$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-059",
   "metadata": {},
   "outputs": [],
   "source": [
    "rng_norm = np.random.default_rng(SEED + 5)\n",
    "\n",
    "n = 100\n",
    "p = 0.5\n",
    "data = rng_norm.binomial(n=n, p=p, size=10_000)\n",
    "\n",
    "mu = n * p\n",
    "sd = np.sqrt(n * p * (1 - p))\n",
    "\n",
    "plt.figure(figsize=(9, 4))\n",
    "plt.hist(data, bins=20, density=True, color=\"lightgrey\", edgecolor=\"white\", label=\"Bin(n,p)\")\n",
    "\n",
    "x = np.linspace(data.min(), data.max(), 400)\n",
    "plt.plot(x, stats.norm.pdf(x, loc=mu, scale=sd), color=\"crimson\", linewidth=2, label=\"Normal\")\n",
    "\n",
    "plt.title(\"Aproksymacja: Bin(100, 0.5) przez normalny\")\n",
    "plt.xlabel(\"k\")\n",
    "plt.ylabel(\"gęstość\")\n",
    "plt.legend()\n",
    "plt.tight_layout()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-062",
   "metadata": {},
   "source": [
    "## Podsumowanie\n",
    "\n",
    "- W testowaniu hipotez porównujemy wynik z próby z tym, co byłoby typowe przy założeniu $H_0$.\n",
    "- Dla hipotez jednostronnych patrzymy na jeden ogon rozkładu, dla dwustronnych — na dwa.\n",
    "- W rozkładach dyskretnych dwustronna p-wartość może wymagać ostrożnej definicji „bardziej skrajnych” wyników.\n",
    "- CTG i błąd standardowy tłumaczą, skąd bierze się logika testu $z$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-063",
   "metadata": {},
   "source": [
    "## Ćwiczenia\n",
    "\n",
    "1) W przykładzie z synestezją zmień $k$ na 12 (zostaw $n=500$ i $p_0=0.014$). Czy decyzja przy $\\alpha=0.05$ się zmieni?\n",
    "\n",
    "2) W przykładzie Bonda zmień $k$ na 11 (zostaw $n=16$). Policz p-wartość jednostronną i dwustronną.\n",
    "\n",
    "3) W przykładzie z testem (gdzie $p_0=0.8$) spróbuj innego wyniku, np. $k=30$ przy $n=50$. Jak zmienia się p-wartość?\n",
    "\n",
    "4) W przykładzie testu $z$ zmień $\\bar{X}$ na 52.0 i sprawdź, jak wygląda p-wartość.\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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}