{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "03000000",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "# Rozkłady i symulacje\n",
    "\n",
    "Ten notatnik jest w 100% w Pythonie: używamy `numpy`, `scipy.stats` i `matplotlib`.\n",
    "\n",
    "W tej części kursu uczymy się, jak korzystać z rozkładów teoretycznych (jako modeli matematycznych) oraz jak wspierać się symulacją Monte Carlo, gdy rachunek „dokładny” jest niewygodny albo chcemy zbudować intuicję.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03000001",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "## Cele\n",
    "\n",
    "Po tym notatniku:\n",
    "\n",
    "- Rozpoznajesz, co znaczy „rozkład teoretyczny” i jakie ma parametry.\n",
    "- Umiesz w Pythonie policzyć: gęstość/masę (`pdf`/`pmf`), dystrybuantę (`cdf`), kwantyl (`ppf`) oraz wylosować próbę (`rvs`).\n",
    "- Umiesz policzyć prawdopodobieństwa ogonów rozkładu (np. $P(X \\ge x)$, $P(X \\le x)$) korzystając z `cdf`/`sf`.\n",
    "- Rozumiesz po co ustawiamy *seed* i jak to pomaga w reprodukowalności.\n",
    "- Potrafisz porównać wynik „z rozkładu” z estymacją z symulacji.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03000002",
   "metadata": {},
   "source": [
    "## Wymagania wstępne\n",
    "\n",
    "- Podstawy Pythona (zmienne, listy/tablice, pętle `for`).\n",
    "- Intuicja: czym jest prawdopodobieństwo i zmienna losowa.\n",
    "- Instalacja środowiska z tego repo (`uv sync`).\n",
    "- (Opcjonalnie) Obsługa widżetów w notatniku (pakiet `ipywidgets`).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03000003",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Rozkłady teoretyczne — przypomnienie\n",
    "\n",
    "**Rozkład teoretyczny**\n",
    "\n",
    "- Jest opisany przez funkcję matematyczną.\n",
    "- Jest określony przez parametry tej funkcji.\n",
    "\n",
    "W Pythonie do pracy z rozkładami teoretycznymi używamy głównie modułu `scipy.stats`.\n",
    "\n",
    "W `scipy.stats` większość rozkładów jest dostępna jako obiekt z metodami (dla przykładu: rozkład normalny `stats.norm`):\n",
    "\n",
    "- `pdf(...)` — funkcja gęstości (*density*) dla rozkładów ciągłych,\n",
    "- `pmf(...)` — funkcja masy prawdopodobieństwa dla rozkładów dyskretnych,\n",
    "- `cdf(...)` — dystrybuanta (*cumulative distribution function*),\n",
    "- `ppf(...)` — odwrotna dystrybuanta (kwantyle),\n",
    "- `rvs(...)` — generator liczb losowych.\n",
    "\n",
    "Dodatkowo bardzo przydatne są:\n",
    "\n",
    "- `sf(...)` (*survival function*) = $P(X > x)$ — „prawy ogon” bez odejmowania od 1,\n",
    "- `isf(...)` — odwrotność `sf`, przydatna np. do wyznaczania granicy obszaru krytycznego.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03000004",
   "metadata": {},
   "source": [
    "## Cztery podstawowe operacje na rozkładzie (pdf/cdf/ppf/rvs) — przykład: rozkład normalny\n",
    "\n",
    "Omówmy teraz zestaw czterech operacji dla rozkładu normalnego.\n",
    "\n",
    "W Pythonie (SciPy) będziemy pisać odpowiednio:\n",
    "\n",
    "- `stats.norm.pdf(...)`\n",
    "- `stats.norm.cdf(...)`\n",
    "- `stats.norm.ppf(...)`\n",
    "- `stats.norm.rvs(...)` (albo bezpośrednio `rng.normal(...)`)\n",
    "\n",
    "Zobaczmy definicje (intuicje) i spróbujmy powiedzieć sobie, co te funkcje robią.\n",
    "\n",
    "### `pdf`\n",
    "\n",
    "> Funkcja gęstości prawdopodobieństwa (gęstość zmiennej losowej) – nieujemna funkcja rzeczywista, określona dla rozkładu prawdopodobieństwa, taka że całka z tej funkcji, obliczona w odpowiednich granicach, jest równa prawdopodobieństwu wystąpienia danego zdarzenia losowego.\n",
    "\n",
    "Ta definicja ma ważną konsekwencję:\n",
    "\n",
    "- dla rozkładów **ciągłych** wartość `pdf(x)` **nie** jest prawdopodobieństwem „w punkcie” $x$;\n",
    "- prawdopodobieństwa dotyczą **przedziałów**, np. $P(80 \\le X \\le 90)$.\n",
    "\n",
    "`pdf` będzie nam służyć m.in. do rysowania krzywych rozkładu.\n",
    "\n",
    "### `cdf`\n",
    "\n",
    "> Dystrybuanta jest definiowana jako prawdopodobieństwo tego, że zmienna $X$ ma wartości mniejsze bądź równe $x$.\n",
    "\n",
    "Czyli `cdf(x)` to $P(X \\le x)$. Czasami interesują nas wartości **większe** od $x$ — wtedy wygodnie używać `sf(x)` zamiast pisać `1 - cdf(x)`.\n",
    "\n",
    "### `ppf`\n",
    "\n",
    "> Odwrotna dystrybuanta (wygodniej myśleć: *kwantyle*).\n",
    "\n",
    "Jeśli `cdf` odpowiada na pytanie „jakie jest $P(X \\le x)$?”, to `ppf` odpowiada na pytanie odwrotne:\n",
    "\n",
    "> „Jaka jest taka wartość $x$, że poniżej niej leży np. 80% rozkładu?”\n",
    "\n",
    "### `rvs`\n",
    "\n",
    "Funkcje losujące pozwalają nam generować próbę z danego rozkładu. Przy dużej liczbie losowań rozkład empiryczny będzie dążył do rozkładu teoretycznego.\n",
    "\n",
    "W dalszej części notatnika wiele zadań zrobimy **na dwa sposoby**:\n",
    "\n",
    "1) korzystając z rozkładu teoretycznego (dokładnie),\n",
    "2) estymując prawdopodobieństwo z losowań (Monte Carlo).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03000005",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from scipy import stats\n",
    "\n",
    "from IPython.display import display\n",
    "\n",
    "SEED = 20250116\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\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03000006",
   "metadata": {},
   "source": [
    "## `pdf`, `cdf` i `ppf` w praktyce (SciPy: `stats.norm`)\n",
    "\n",
    "Zacznijmy od narysowania krzywej ilustrującej rozkład normalny o parametrach $\\mu=100$ i $\\sigma=15$.\n",
    "\n",
    "Znamy regułę „trzech sigm”, a tu dla wygody weźmiemy zakres od $-4\\sigma$ do $4\\sigma$ (czyli w praktyce prawie cały rozkład).\n",
    "\n",
    "W kognitywistyce zmienna $X$ może oznaczać na przykład czas reakcji (ms).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03000007",
   "metadata": {},
   "outputs": [],
   "source": [
    "mu = 100\n",
    "sigma = 15\n",
    "\n",
    "x = np.linspace(-4, 4, num=1000) * sigma + mu\n",
    "pdf_values = stats.norm.pdf(x, loc=mu, scale=sigma)\n",
    "\n",
    "plt.figure(figsize=(8, 4))\n",
    "plt.plot(x, pdf_values, linewidth=2)\n",
    "plt.title(\"Rozkład normalny: gęstość (pdf)\")\n",
    "plt.xlabel(\"Wartość zmiennej X\")\n",
    "plt.ylabel(\"Gęstość\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03000008",
   "metadata": {},
   "source": [
    "Przyjrzyjmy się teraz, jak działa dystrybuanta.\n",
    "\n",
    "Chcemy dowiedzieć się, jakie jest prawdopodobieństwo, że $X \\le 80$.\n",
    "\n",
    "W rozkładzie normalnym odpowiada za to `stats.norm.cdf(80, loc=..., scale=...)`.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03000009",
   "metadata": {},
   "outputs": [],
   "source": [
    "x0 = 80\n",
    "\n",
    "p_x_le_80 = stats.norm.cdf(x0, loc=mu, scale=sigma)\n",
    "p_x_le_80\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0300000a",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8, 4))\n",
    "plt.plot(x, pdf_values, linewidth=2)\n",
    "plt.axvline(x=x0, linestyle=\"--\", linewidth=2)\n",
    "\n",
    "plt.title(\"Rozkład normalny: ile leży poniżej 80?\")\n",
    "plt.xlabel(\"Wartość zmiennej X\")\n",
    "plt.ylabel(\"Gęstość\")\n",
    "\n",
    "label = f\"{p_x_le_80 * 100:.2f}%\"\n",
    "plt.text(x0 - 25, 0, label, fontsize=12)\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0300000b",
   "metadata": {},
   "source": [
    "### Widżet: pole pod krzywą (dla rozkładu normalnego)\n",
    "\n",
    "Zmień parametry rozkładu oraz przedział $[a,b]$ i zobacz, jaki obszar pod krzywą odpowiada prawdopodobieństwu $P(a < X < b)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0300000c",
   "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",
    "    mu_slider = widgets.FloatSlider(value=100, min=0, max=200, step=1, description='μ')\n",
    "    sigma_slider = widgets.FloatSlider(value=15, min=1, max=60, step=1, description='σ')\n",
    "    przedzial_slider = widgets.FloatRangeSlider(\n",
    "        value=(80, 120),\n",
    "        min=0,\n",
    "        max=200,\n",
    "        step=1,\n",
    "        description='[a,b]',\n",
    "    )\n",
    "\n",
    "    def pokaz_pole(mu, sigma, przedzial):\n",
    "        a, b = przedzial\n",
    "        if a >= b:\n",
    "            print('Ustaw a < b.')\n",
    "            return\n",
    "\n",
    "        rozklad = stats.norm(loc=mu, scale=sigma)\n",
    "        p_between = float(rozklad.cdf(b) - rozklad.cdf(a))\n",
    "        print(f'P({a:.0f} < X < {b:.0f}) = {p_between:.4f}')\n",
    "\n",
    "        x_grid = np.linspace(mu - 4 * sigma, mu + 4 * sigma, 400)\n",
    "        pdf = rozklad.pdf(x_grid)\n",
    "\n",
    "        plt.figure(figsize=(8, 4))\n",
    "        plt.plot(x_grid, pdf, linewidth=2)\n",
    "        plt.fill_between(x_grid, 0, pdf, where=(x_grid >= a) & (x_grid <= b), alpha=0.3)\n",
    "        plt.title('Pole pod krzywą PDF na przedziale [a,b]')\n",
    "        plt.xlabel('Wartość zmiennej X')\n",
    "        plt.ylabel('Gęstość')\n",
    "        plt.show()\n",
    "\n",
    "    out = widgets.interactive_output(\n",
    "        pokaz_pole,\n",
    "        {'mu': mu_slider, 'sigma': sigma_slider, 'przedzial': przedzial_slider},\n",
    "    )\n",
    "    display(widgets.VBox([mu_slider, sigma_slider, przedzial_slider]), out)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0300000d",
   "metadata": {},
   "source": [
    "Jeżeli chcemy dowiedzieć się, **poniżej jakiej wartości** leży 80% rozkładu, używamy `ppf` (kwantyl).\n",
    "\n",
    "Czyli liczymy wartość $x$ spełniającą:\n",
    "\n",
    "$$\n",
    "P(X \\le x) = 0.8\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0300000e",
   "metadata": {},
   "outputs": [],
   "source": [
    "q80 = stats.norm.ppf(0.8, loc=mu, scale=sigma)\n",
    "q80\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0300000f",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8, 4))\n",
    "plt.plot(x, pdf_values, linewidth=2)\n",
    "plt.axvline(x=q80, linestyle=\":\", linewidth=3)\n",
    "\n",
    "plt.title(\"Rozkład normalny: 80. percentyl\")\n",
    "plt.xlabel(\"Wartość zmiennej X\")\n",
    "plt.ylabel(\"Gęstość\")\n",
    "\n",
    "plt.text(q80 - 18, max(pdf_values) * 0.6, \"80%\", fontsize=14)\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03000010",
   "metadata": {},
   "source": [
    "## Symulacje statystyczne (losowanie)\n",
    "\n",
    "### Po co w ogóle to robić?\n",
    "\n",
    "Jest kilka powodów, dla których „eksperymenty” polegające na losowaniu z rozkładów teoretycznych są interesujące:\n",
    "\n",
    "- to świetne narzędzie dydaktyczne — pozwala zobaczyć, jak „działają” zależności statystyczne,\n",
    "- czasami pozwala dowiedzieć się czegoś nowego, kiedy rachunek jest trudny lub kłopotliwy,\n",
    "- jest podstawą wielu praktycznych procedur (np. bootstrap, symulacyjne przedziały ufności).\n",
    "\n",
    "W tym notatniku będziemy rozwiązywać zadania na dwa sposoby:\n",
    "\n",
    "1) korzystając z rozkładów teoretycznych,\n",
    "2) estymując prawdopodobieństwa z losowań.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03000011",
   "metadata": {},
   "source": [
    "### Stały seed w Pythonie: reprodukowalność\n",
    "\n",
    "1. Komputer generuje liczby *pseudolosowe*. To znaczy: nie są one „prawdziwie losowe”, ale dla celów statystyki zwykle możemy je traktować jak losowe.\n",
    "\n",
    "2. Konsekwencja pseudolosowości: jeśli ustawimy ten sam stan początkowy algorytmu (ziarno, *seed*), to przy tym samym kodzie uzyskamy te same wyniki.\n",
    "\n",
    "W statystyce to zaleta — chcemy, żeby ktoś mógł łatwo odtworzyć naszą symulację.\n",
    "\n",
    "W NumPy najwygodniej użyć generatora:\n",
    "\n",
    "```python\n",
    "rng = np.random.default_rng(SEED)\n",
    "```\n",
    "\n",
    "Ważna dydaktyczna uwaga: żeby unikać „ukrytego stanu” między komórkami, często w kolejnych przykładach tworzymy generator od nowa, np. `np.random.default_rng(SEED + 1)`.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03000012",
   "metadata": {},
   "outputs": [],
   "source": [
    "rng1 = np.random.default_rng(SEED)\n",
    "rng2 = np.random.default_rng(SEED)\n",
    "\n",
    "losowania_1 = rng1.integers(low=0, high=10, size=5)\n",
    "losowania_2 = rng2.integers(low=0, high=10, size=5)\n",
    "\n",
    "losowania_1, losowania_2\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03000013",
   "metadata": {},
   "source": [
    "### Rzut kostką\n",
    "\n",
    "Pierwszy przykład będzie dotyczył rzucania kostką.\n",
    "\n",
    "Na początek wykonajmy symulację 6000 rzutów uczciwą kostką i policzmy, ile razy wypadł każdy wynik.\n",
    "\n",
    "W NumPy do losowania z wektora/tablicy użyjemy `rng.choice(...)`:\n",
    "\n",
    "```python\n",
    "rng.choice(wektor, size=..., replace=True)  # losowanie ze zwracaniem\n",
    "rng.choice(wektor, size=..., replace=False) # losowanie bez zwracania\n",
    "```\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03000014",
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.default_rng(SEED + 1)\n",
    "\n",
    "kostka = np.arange(1, 7)\n",
    "\n",
    "liczba_rzutow = 6000\n",
    "symulowane_rzuty = rng.choice(kostka, size=liczba_rzutow, replace=True)\n",
    "\n",
    "symulowane_rzuty[:10]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03000015",
   "metadata": {},
   "outputs": [],
   "source": [
    "rzuty_series = pd.Series(symulowane_rzuty, name=\"wynik\")\n",
    "\n",
    "licznosci = rzuty_series.value_counts().sort_index()\n",
    "licznosci\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03000016",
   "metadata": {},
   "source": [
    "Spodziewamy się około 1000 wystąpień każdej z sześciu możliwości (bo $6000 / 6 = 1000$).\n",
    "\n",
    "Teraz przejdźmy do pytania probabilistycznego.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03000017",
   "metadata": {},
   "source": [
    "#### Historia → pytanie\n",
    "\n",
    "Wyobraźmy sobie, że ktoś mówi:\n",
    "\n",
    "> „W 6000 rzutach kostką wypadło **więcej niż 1030 jedynek**”.\n",
    "\n",
    "Jak często takie zdarzenie powinno się zdarzać, jeśli zakładamy uczciwą kostkę?\n",
    "\n",
    "Żeby to policzyć, zdefiniujemy zmienną losową:\n",
    "\n",
    "$$\n",
    "X = \\#\\{\\text{jedynek w 6000 rzutach}\\}\n",
    "$$\n",
    "\n",
    "i policzymy prawy ogon rozkładu $P(X > 1030)$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03000018",
   "metadata": {},
   "source": [
    "#### Model probabilistyczny\n",
    "\n",
    "Zakładamy:\n",
    "\n",
    "- w każdym rzucie $P(\\text{jedynka}) = 1/6$,\n",
    "- rzuty są niezależne.\n",
    "\n",
    "Wtedy liczba jedynek w 6000 rzutach ma rozkład dwumianowy.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03000019",
   "metadata": {},
   "source": [
    "#### Rozkład zmiennej $X$\n",
    "\n",
    "$$\n",
    "X \\sim \\mathrm{Bin}(n=6000, p=1/6)\n",
    "$$\n",
    "\n",
    "To jest rozkład „liczby sukcesów” w $n$ niezależnych próbach.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0300001a",
   "metadata": {},
   "source": [
    "#### Kwantyle: gdzie leży 95% rozkładu?\n",
    "\n",
    "Zanim policzymy ogon, sprawdźmy, jakie wartości $X$ są typowe.\n",
    "\n",
    "Na przykład 95. percentyl $q_{0.95}$ to taka liczba, że około 95% rozkładu leży **na lewo** od niej:\n",
    "\n",
    "$$\n",
    "P(X \\le q_{0.95}) \\approx 0.95\n",
    "$$\n",
    "\n",
    "W SciPy liczymy to przez `ppf(0.95, ...)`.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0300001b",
   "metadata": {},
   "outputs": [],
   "source": [
    "n = liczba_rzutow\n",
    "p_jedynka = 1 / 6\n",
    "\n",
    "srednia = n * p_jedynka\n",
    "odch_std = (n * p_jedynka * (1 - p_jedynka)) ** 0.5\n",
    "\n",
    "q95 = stats.binom.ppf(0.95, n=n, p=p_jedynka)\n",
    "q99 = stats.binom.ppf(0.99, n=n, p=p_jedynka)\n",
    "\n",
    "srednia, odch_std, q95, q99\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0300001c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Dla rozkładów dyskretnych `ppf` zwraca liczbę całkowitą (ale jako float).\n",
    "q95_int = int(q95)\n",
    "\n",
    "# sf(k) = P(X > k)\n",
    "p_gt_q95_via_cdf = 1 - stats.binom.cdf(q95_int, n=n, p=p_jedynka)\n",
    "p_gt_q95_via_sf = stats.binom.sf(q95_int, n=n, p=p_jedynka)\n",
    "\n",
    "p_gt_q95_via_cdf, p_gt_q95_via_sf\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0300001d",
   "metadata": {},
   "source": [
    "#### Prawdopodobieństwo ogona „na piechotę”\n",
    "\n",
    "Interesuje nas:\n",
    "\n",
    "$$\n",
    "P(X > 1030) = P(X \\ge 1031)\n",
    "$$\n",
    "\n",
    "W `scipy.stats` dla rozkładu dwumianowego wygodnie użyć `sf(1030, ...)`, bo `sf(k)` liczy $P(X > k)$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0300001e",
   "metadata": {},
   "outputs": [],
   "source": [
    "p_x_gt_1030 = stats.binom.sf(1030, n=n, p=p_jedynka)  # P(X > 1030)\n",
    "p_x_gt_1030\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0300001f",
   "metadata": {},
   "source": [
    "#### Weryfikacja obliczenia: `sf` vs `1 - cdf`\n",
    "\n",
    "Dla ogonów często spotkasz dwie wersje tego samego obliczenia:\n",
    "\n",
    "- bezpośrednio: `sf(k)`,\n",
    "- przez dopełnienie: `1 - cdf(k)`.\n",
    "\n",
    "`sf` jest zwykle bezpieczniejszy numerycznie, ale oba podejścia powinny dać ten sam wynik (z dokładnością do zaokrągleń).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03000020",
   "metadata": {},
   "outputs": [],
   "source": [
    "p_x_gt_1030_from_cdf = 1 - stats.binom.cdf(1030, n=n, p=p_jedynka)\n",
    "\n",
    "p_x_gt_1030, p_x_gt_1030_from_cdf\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03000021",
   "metadata": {},
   "source": [
    "#### Interpretacja\n",
    "\n",
    "- `p_x_gt_1030` to prawdopodobieństwo zdarzenia „więcej niż 1030 jedynek” **w przyjętym modelu uczciwej kostki**.\n",
    "- Jeśli ta liczba jest mała, to takie zdarzenie jest rzadkie w tym modelu; jeśli duża — jest zupełnie typowe.\n",
    "\n",
    "Za chwilę zobaczymy to samo na symulacji Monte Carlo.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03000022",
   "metadata": {},
   "source": [
    "#### Symulacja: jak często zdarza się „więcej niż 1030”?\n",
    "\n",
    "Powtórzmy doświadczenie wiele razy i policzmy, w jakiej części replikacji liczba jedynek przekroczy 1030.\n",
    "\n",
    "Zrobimy to w Pythonie czytelną pętlą.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03000023",
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.default_rng(SEED + 2)\n",
    "\n",
    "liczba_replikacji = 1000\n",
    "jedynki = np.empty(liczba_replikacji, dtype=int)\n",
    "\n",
    "for i in range(liczba_replikacji):\n",
    "    rzuty = rng.choice(kostka, size=liczba_rzutow, replace=True)\n",
    "    liczba_jedynek = np.sum(rzuty == 1)\n",
    "    jedynki[i] = liczba_jedynek\n",
    "\n",
    "jedynki[:10]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03000024",
   "metadata": {},
   "outputs": [],
   "source": [
    "estymacja = np.mean(jedynki > 1030)\n",
    "\n",
    "p_teoretyczne = stats.binom.sf(1030, n=n, p=p_jedynka)  # P(X > 1030) = P(X >= 1031)\n",
    "\n",
    "estymacja, p_teoretyczne\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03000025",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8, 4))\n",
    "plt.hist(jedynki, bins=30, edgecolor=\"white\")\n",
    "plt.axvline(x=1030, linestyle=\"--\", linewidth=2, color=\"black\")\n",
    "\n",
    "plt.title(\"Rozkład liczby jedynek w 6000 rzutach (symulacja)\")\n",
    "plt.xlabel(\"Liczba jedynek\")\n",
    "plt.ylabel(\"Liczba replikacji\")\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03000026",
   "metadata": {},
   "source": [
    "### Rzut monetą\n",
    "\n",
    "Zróbmy kilka ćwiczeń z prawdopodobieństwa.\n",
    "\n",
    "Rozważmy dziesięciokrotny rzut symetryczną monetą. Użyjemy rozkładu dwumianowego, który ma dwa parametry:\n",
    "\n",
    "- liczba prób $n$,\n",
    "- prawdopodobieństwo sukcesu $p$ (tu: „orzeł”).\n",
    "\n",
    "W naszym przypadku: $n=10$, $p=0.5$.\n",
    "\n",
    "Ten sam rozkład dwumianowy opisuje też liczbę odpowiedzi poprawnych w serii prób zadania dwuwyborowego (2AFC).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03000027",
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 10\n",
    "p = 0.5\n",
    "\n",
    "# X = liczba orłów w 10 rzutach\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03000028",
   "metadata": {},
   "source": [
    "#### Prawdopodobieństwo: dokładnie 4 orły\n",
    "\n",
    "**Historia:** rzucamy 10 razy monetą.\n",
    "\n",
    "**Pytanie:** jakie jest $P(X = 4)$, jeśli moneta jest symetryczna?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03000029",
   "metadata": {},
   "outputs": [],
   "source": [
    "p_x_eq_4 = stats.binom.pmf(4, n=n, p=p)\n",
    "p_x_eq_4\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0300002a",
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.default_rng(SEED + 3)\n",
    "\n",
    "s = rng.binomial(n=n, p=p, size=10_000)\n",
    "\n",
    "estymacja_eq_4 = np.mean(s == 4)\n",
    "estymacja_eq_4\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0300002b",
   "metadata": {},
   "source": [
    "#### Prawdopodobieństwo: co najmniej 4 orły\n",
    "\n",
    "Teraz pytanie o ogon rozkładu:\n",
    "\n",
    "$$\n",
    "P(X \\ge 4)\n",
    "$$\n",
    "\n",
    "W rozkładach dyskretnych warto uważać na nierówności:\n",
    "\n",
    "- $P(X \\ge 4)$ to to samo co $P(X > 3)$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0300002c",
   "metadata": {},
   "outputs": [],
   "source": [
    "k_values = np.arange(4, n + 1)\n",
    "\n",
    "p_ge_4_sum = stats.binom.pmf(k_values, n=n, p=p).sum()\n",
    "\n",
    "# Wygodniej i stabilniej numerycznie:\n",
    "p_ge_4_sf = stats.binom.sf(3, n=n, p=p)  # P(X > 3) = P(X >= 4)\n",
    "\n",
    "p_ge_4_sum, p_ge_4_sf\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0300002d",
   "metadata": {},
   "outputs": [],
   "source": [
    "estymacja_ge_4 = np.mean(s >= 4)\n",
    "estymacja_ge_4\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0300002e",
   "metadata": {},
   "source": [
    "#### `pmf`, `cdf` i `sf` — trzy różne pytania\n",
    "\n",
    "Na tym samym przykładzie (10 rzutów, $X$ = liczba orłów) warto utrwalić różnicę:\n",
    "\n",
    "- `pmf(4)` to $P(X=4)$ (dokładnie 4),\n",
    "- `cdf(4)` to $P(X\\le 4)$ (co najwyżej 4),\n",
    "- `sf(3)` to $P(X>3)=P(X\\ge 4)$ (co najmniej 4).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0300002f",
   "metadata": {},
   "source": [
    "### Widżet: prawy ogon w rozkładzie dwumianowym\n",
    "\n",
    "Ustaw $n$, $p$ oraz próg $k$ i zobacz, jak wygląda prawdopodobieństwo ogona $P(X\\ge k)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03000030",
   "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(value=10, min=1, max=200, step=1, description='n')\n",
    "    p_slider = widgets.FloatSlider(value=0.5, min=0.01, max=0.99, step=0.01, description='p')\n",
    "    k_slider = widgets.IntSlider(value=4, min=0, max=n_slider.value, step=1, description='k')\n",
    "\n",
    "    def _update_k_max(change):\n",
    "        k_slider.max = change['new']\n",
    "        if k_slider.value > k_slider.max:\n",
    "            k_slider.value = k_slider.max\n",
    "\n",
    "    n_slider.observe(_update_k_max, names='value')\n",
    "\n",
    "    def pokaz_ogon(n, p, k):\n",
    "        k = int(min(max(k, 0), n))\n",
    "\n",
    "        k_values = np.arange(0, n + 1)\n",
    "        pmf_values = stats.binom.pmf(k_values, n=n, p=p)\n",
    "\n",
    "        p_x_ge_k = float(stats.binom.sf(k - 1, n=n, p=p))\n",
    "\n",
    "        print(f'P(X >= {k}) = {p_x_ge_k:.6f}')\n",
    "\n",
    "        kolory = np.array(['lightgrey'] * len(k_values), dtype=object)\n",
    "        kolory[k_values >= k] = '#4C78A8'\n",
    "\n",
    "        plt.figure(figsize=(8, 3))\n",
    "        plt.bar(k_values, pmf_values, color=kolory, edgecolor='white')\n",
    "        plt.title('Rozkład dwumianowy: ogon P(X >= k)')\n",
    "        plt.xlabel('k = liczba sukcesów')\n",
    "        plt.ylabel('P(X=k)')\n",
    "        plt.show()\n",
    "\n",
    "    out = widgets.interactive_output(\n",
    "        pokaz_ogon,\n",
    "        {'n': n_slider, 'p': p_slider, 'k': k_slider},\n",
    "    )\n",
    "    display(widgets.VBox([n_slider, p_slider, k_slider]), out)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03000031",
   "metadata": {},
   "outputs": [],
   "source": [
    "p_eq_4 = stats.binom.pmf(4, n=n, p=p)\n",
    "p_le_4 = stats.binom.cdf(4, n=n, p=p)\n",
    "p_ge_4 = stats.binom.sf(3, n=n, p=p)\n",
    "\n",
    "p_eq_4, p_le_4, p_ge_4\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03000032",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Weryfikacja na symulacji (te same losowania `s` z poprzednich komórek)\n",
    "estymacja_eq_4 = np.mean(s == 4)\n",
    "estymacja_le_4 = np.mean(s <= 4)\n",
    "estymacja_ge_4 = np.mean(s >= 4)\n",
    "\n",
    "estymacja_eq_4, estymacja_le_4, estymacja_ge_4\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03000033",
   "metadata": {},
   "source": [
    "Porównaj wyniki „dokładne” (z rozkładu) z wynikami z symulacji.\n",
    "\n",
    "Przy 10 000 losowań różnice powinny być niewielkie, ale nie będą idealnie równe — to naturalny błąd Monte Carlo.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03000034",
   "metadata": {},
   "source": [
    "### Rozkład hipergeometryczny\n",
    "\n",
    "Kolejne zadanie:\n",
    "\n",
    "Jako szef kampusu masz za zadanie skompletować delegację składającą się z 7 członków organizacji.\n",
    "\n",
    "Organizacja składa się z 18 kobiet oraz 15 mężczyzn.\n",
    "\n",
    "**Pytanie:** jakie jest prawdopodobieństwo, że w losowo wybranej delegacji będzie **więcej niż 4 mężczyzn**?\n",
    "\n",
    "Najpierw zróbmy symulację, a potem policzmy to samo „dokładnie” z rozkładu hipergeometrycznego.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03000035",
   "metadata": {},
   "outputs": [],
   "source": [
    "liczba_mezczyzn = 15\n",
    "liczba_kobiet = 18\n",
    "\n",
    "org = np.array([\"M\"] * liczba_mezczyzn + [\"K\"] * liczba_kobiet)\n",
    "\n",
    "org[:10], org.size\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03000036",
   "metadata": {},
   "source": [
    "Pojedyncze powtórzenie eksperymentu: losujemy 7 osób **bez zwracania** (`replace=False`).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03000037",
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.default_rng(SEED + 4)\n",
    "\n",
    "delegacja = rng.choice(org, size=7, replace=False)\n",
    "delegacja\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03000038",
   "metadata": {},
   "outputs": [],
   "source": [
    "liczba_M = int(np.sum(delegacja == \"M\"))\n",
    "liczba_M\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03000039",
   "metadata": {},
   "source": [
    "Teraz powtórzmy to 10 000 razy i sprawdźmy, jak często mamy więcej niż 4 mężczyzn.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0300003a",
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.default_rng(SEED + 5)\n",
    "\n",
    "liczba_replikacji = 10_000\n",
    "s = np.empty(liczba_replikacji, dtype=int)\n",
    "\n",
    "for i in range(liczba_replikacji):\n",
    "    delegacja = rng.choice(org, size=7, replace=False)\n",
    "    s[i] = np.sum(delegacja == \"M\")\n",
    "\n",
    "estymacja = np.mean(s > 4)\n",
    "estymacja\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0300003b",
   "metadata": {},
   "source": [
    "#### Rozwiązanie „dokładne”: rozkład hipergeometryczny\n",
    "\n",
    "Jeśli losujemy bez zwracania z populacji o znanej liczbie „sukcesów” (tu: mężczyźni), to liczba sukcesów w próbie ma rozkład hipergeometryczny.\n",
    "\n",
    "Oznaczenia w `scipy.stats.hypergeom`:\n",
    "\n",
    "- `M` — liczebność populacji,\n",
    "- `n` — liczba sukcesów w populacji,\n",
    "- `N` — wielkość próby (liczba losowań).\n",
    "\n",
    "U nas:\n",
    "\n",
    "- `M = 33`, `n = 15`, `N = 7`.\n",
    "\n",
    "Interesuje nas:\n",
    "\n",
    "$$\n",
    "P(X > 4) = P(X \\ge 5)\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0300003c",
   "metadata": {},
   "outputs": [],
   "source": [
    "M = liczba_mezczyzn + liczba_kobiet\n",
    "n_success = liczba_mezczyzn\n",
    "N_draws = 7\n",
    "\n",
    "p_gt_4 = stats.hypergeom.sf(4, M=M, n=n_success, N=N_draws)  # P(X > 4)\n",
    "\n",
    "# Możemy też jawnie zsumować pmf od 5 do 7:\n",
    "k_values = np.arange(5, N_draws + 1)\n",
    "p_gt_4_sum = stats.hypergeom.pmf(k_values, M=M, n=n_success, N=N_draws).sum()\n",
    "\n",
    "p_gt_4, p_gt_4_sum\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0300003d",
   "metadata": {},
   "source": [
    "Całkiem nieźle — wynik z symulacji powinien być blisko wyniku „dokładnego”. Różnica wynika z błędu Monte Carlo (im więcej replikacji, tym mniejszy).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0300003e",
   "metadata": {},
   "source": [
    "## Ćwiczenia\n",
    "\n",
    "1) Rozkład normalny:\n",
    "\n",
    "- Dla $X \\sim N(100, 15)$ policz $P(X > 130)$.\n",
    "- Znajdź 95. percentyl (taki $x$, że $P(X \\le x)=0.95$).\n",
    "\n",
    "2) Kostka:\n",
    "\n",
    "- Dla uczciwej kostki policz $P(X \\ge 1100)$, gdzie $X$ to liczba jedynek w 6000 rzutach.\n",
    "- Oszacuj to samo metodą Monte Carlo (np. 20 000 replikacji).\n",
    "\n",
    "3) Moneta:\n",
    "\n",
    "- Dla 20 rzutów symetryczną monetą policz $P(X \\ge 14)$, gdzie $X$ to liczba orłów.\n",
    "- Oszacuj to samo metodą Monte Carlo (np. 50 000 replikacji).\n",
    "\n",
    "4) Hipergeom:\n",
    "\n",
    "- Organizacja ma 12 kobiet i 8 mężczyzn, losujesz 5 osób. Jakie jest $P(X \\ge 4)$?\n",
    "\n",
    "5) Wariant (2AFC):\n",
    "\n",
    "- Załóż, że w 10 próbach zadania dwuwyborowego prawdopodobieństwo odpowiedzi poprawnej wynosi $p=0.7$. Policz $P(X \\ge 8)$.\n",
    "- Oszacuj to samo metodą Monte Carlo (np. 50 000 replikacji).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0300003f",
   "metadata": {},
   "source": [
    "## Recap\n",
    "\n",
    "- Rozkład teoretyczny to model opisany funkcją i parametrami.\n",
    "- W Pythonie: `pdf/pmf`, `cdf`, `ppf`, `rvs` (oraz `sf/isf`).\n",
    "- Umiesz liczyć prawdopodobieństwa ogonów (np. $P(X\\le x)$ i $P(X>x)$) i porównywać wynik „z rozkładu” z Monte Carlo.\n",
    "- *Seed* to narzędzie reprodukowalności.\n",
    "- Symulacje Monte Carlo pomagają budować intuicję i przybliżać prawdopodobieństwa.\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
}