Bootstrap-t-student/bootstrap-t.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Projekt - Test t studenta\n",
    "\n",
    "- Marcin Kostrzewski\n",
    "- Krystian Wasilewski\n",
    "- Mateusz Tylka\n",
    "\n",
    "## Test t studenta\n",
    "\n",
    "Metoda statystyczna służącą do porównania dwóch średnich między sobą gdy znamy liczbę badanych próbek, średnią arytmetyczną oraz wartość odchylenia standardowego lub wariancji.\n",
    "Jest to jeden z mniej skomplikowanych i bardzo często wykorzystywanych testów statystycznych używanych do weryfikacji hipotez. Dzięki niemu możemy dowiedzieć się czy dwie różne średnie są różne niechcący (w wyniku przypadku) czy są różne istotnie statystycznie (np. z uwagi na naszą manipulację eksperymentalna).\n",
    "Wyróżniamy 3 wersję testu t:\n",
    "\n",
    "1. test t Studenta dla jednej próby\n",
    "2. test t Studenta dla prób niezależnych\n",
    "3. test t Studenta dla prób zależnych\n",
    "\n",
    "Wszystkie rodzaje testów są testami parametrycznymi, a co za tym idzie nasze mierzone zmienne ilościowe powinny mieć rozkład normalny."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from math import sqrt\n",
    "from scipy import stats\n",
    "from scipy.stats import sem\n",
    "from scipy.stats import t\n",
    "import matplotlib.pyplot as plt\n",
    "from statistics import mean, stdev\n",
    "from scipy.stats import ttest_ind, ttest_1samp, ttest_rel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset = pd.read_csv('experiment_data.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_p(t_stat, df):\n",
    "    \"\"\"Funkcja oblicza wartość *p* na podstawie statystyki testowej i stopni swobody\"\"\"\n",
    "    return (1.0 - t.cdf(abs(t_stat), df)) * 2.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_cv(df, alpha=0.05):\n",
    "    \"\"\"Funkcja oblicza wartość krytyczną (critical value)\"\"\"\n",
    "    return t.ppf(1.0 - alpha, df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "def t_test(sample_1, sample_2=None, df_fn=df_single, t_stat_fn=t_stat_single, population_mean=None, alpha=0.05):\n",
    "    \"\"\"\n",
    "    Funkcja przeprowadza test T-studenta dla dwóch zmiennych.\n",
    "    liczba kolumn wynosi 1, test jest przeprowadzany dla jednej zmiennej.\n",
    "    @param df_fn - funkcja obliczająca stopnie swobody\n",
    "    @param t_stat_fn - funkcja obliczająca statystykę T\n",
    "    \"\"\"\n",
    "    t_stat_list = get_t_stats(sample_1, sample_2, t_stat_fn, population_mean=population_mean)\n",
    "    t_stat_sum = sum(t_stat_list)\n",
    "\n",
    "    data_size = sample_1.shape[0]\n",
    "\n",
    "    t_stat = t_stat_sum / data_size\n",
    "    # TODO: dolna i górna opcja dają inne wyniki z jakiegoś powodu (???)\n",
    "    t_stat = mean(t_stat_list)\n",
    "\n",
    "    if sample_2 is None:\n",
    "        df = df_fn(sample_1)\n",
    "    else:\n",
    "        df = df_fn(sample_1, sample_2)\n",
    "    cv = calculate_cv(df, alpha)\n",
    "    p = calculate_p(t_stat, df)\n",
    "    return t_stat, df, cv, p, t_stat_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_t_stats(sample_1, sample_2=None, t_stat_fn=t_stat_single, population_mean=None):\n",
    "    \"\"\"Funkcja oblicza listę statystyk testowych dla każdej próbki bootstrapowej wybranej na podstawie danych sample_1 i sample_2\"\"\"\n",
    "    t_stat_list = []\n",
    "\n",
    "    # One sample test\n",
    "    if t_stat_fn==t_stat_single:\n",
    "        if not population_mean:\n",
    "            raise Exception(\"population_mean not provided\")\n",
    "        for bootstrap in generate_bootstraps(sample_1):\n",
    "            stat = t_stat_fn(bootstrap, population_mean)\n",
    "            t_stat_list.append(stat)\n",
    "        return t_stat_list\n",
    "\n",
    "    # Two sample test\n",
    "    for bootstrap_1, bootstrap_2 in zip(generate_bootstraps(sample_1), generate_bootstraps(sample_2)):\n",
    "        stat = t_stat_fn(bootstrap_1, bootstrap_2)\n",
    "        t_stat_list.append(stat)\n",
    "    return t_stat_list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Test Shapiro Wilka\n",
    "\n",
    "Wszystkie rodzaje testów są testami parametrycznymi, a co za tym idzie nasze mierzone zmienne ilościowe powinny mieć rozkład normalny."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9606528878211975\n",
      "2.666284970587185e-10\n",
      "Dane nie mają rozkładu normalnego.\n"
     ]
    }
   ],
   "source": [
    "# TODO: Test Shapiro Wilka sprawdzający czy nasze dane mają rozkład normalny\n",
    "x = dataset['Height'].to_numpy()\n",
    "shapiro_test = stats.shapiro(x)\n",
    "print(shapiro_test.statistic)\n",
    "print(shapiro_test.pvalue)\n",
    "\n",
    "if shapiro_test.pvalue > shapiro_test.statistic:\n",
    "    print(\"Dane mają rozkład normalny.\")\n",
    "else:\n",
    "    print(\"Dane nie mają rozkładu normalnego.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Testowanie hipotez metodą bootstrap\n",
    "\n",
    "**Bootstrap** – metoda szacowania (estymacji) wyników poprzez wielokrotne losowanie ze zwracaniem z próby. Polega ona na utworzeniu nowego rozkładu wyników, na podstawie posiadanych danych, poprzez wielokrotne losowanie wartości z posiadanej próby. Metoda ze zwracaniem polega na tym, że po wylosowaniu danej wartości, “wraca” ona z powrotem do zbioru.\n",
    "\n",
    "Metoda bootstrapowa znajduje zastosowanie w sytuacji, w której nie znamy rozkładu z populacji z której pochodzi próbka lub w przypadku rozkładów małych lub asymetrycznych. W takim wypadku, dzięki tej metodzie, wyniki testów parametrycznych i analiz opartych o modele liniowe są bardziej precyzyjne. Zazwyczaj losuje się wiele próbek, np. 2000 czy 5000."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "def generate_bootstraps(data, n_bootstraps=100):\n",
    "    data_size = data.shape[0]\n",
    "    for _ in range(n_bootstraps):\n",
    "        indices =  np.random.choice(len(data), size=data_size)\n",
    "        yield data.iloc[indices, :]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Test t studenta dla jednej próby\n",
    "\n",
    "**Test t Studenta dla jednej próby** wykorzystujemy gdy chcemy porównać średnią “teoretyczną” ze średnią, którą faktycznie możemy zaobserwować w naszej bazie danych. Średnia teoretyczna to średnia pochodząca z innych badań lub po prostu bez większych uzasadnień pochodząca z naszej głowy.\n",
    "\n",
    "Wyobraźmy sobie, że mamy dane z takimi zmiennymi jak wzrost pewnej grupy ludzi. Dzięki testowi t Studenta dla jednej próby możemy dowiedzieć się np. czy wzrost naszego młodszego brata wynoszący 155cm odbiega znacząco od średniej wzrostu tej grupy. Hipoteza zerowa w takim badaniu wyglądałaby następująco H0: Badana próba została wylosowana z populacji, w której wzrost osób wynosi średnio 155cm. Z kolei hipoteza alternatywna będzie brzmiała H1: Badana próba nie została wylosowana z populacji gdzie średni wzrost wynosi 155cm\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "def t_stat_single(sample, population_mean):\n",
    "    \"\"\"Funkcja oblicza wartość statystyki testowej dla jednej próbki\"\"\"\n",
    "    if sample.empty:\n",
    "        raise Exception(\"Empty sample\")\n",
    "    sample = sample['Height'].values.tolist()\n",
    "    sample_size = len(sample)\n",
    "    return (mean(sample) - population_mean) / (stdev(sample) / sqrt(sample_size))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "def df_single(sample_1):\n",
    "    \"\"\"Funkcja oblicza stopnie swobody dla jednej próbki\"\"\"\n",
    "    # TODO: I have no clue what to return from here\n",
    "    return len(sample_1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "def bootstrap_one_sample(sample, population_mean):\n",
    "    return t_test(\n",
    "        sample_1=sample,\n",
    "        df_fn=df_single,\n",
    "        t_stat_fn=t_stat_single,\n",
    "        population_mean=population_mean\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sprawdzenie czy osoba o wzroście 165cm pasuje do populacji (nie jest odmieńcem)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "t: 6.854929920812628, df: 500, cv: 1.6479068539295045, p: 2.1091128843409024e-11\n",
      "\n"
     ]
    }
   ],
   "source": [
    "#TODO: poprawić kod aby można było podawać kolumny\n",
    "\n",
    "t_stat, df, cv, p, _ = bootstrap_one_sample(dataset, 165)\n",
    "pretty_print_full_stats(t_stat, df, cv, p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "TODO: Wniosek"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Test t studenta dla prób niezależnych\n",
    "\n",
    "**Test t Studenta dla prób niezależnych** jest najczęściej stosowaną metodą statystyczną w celu porównania średnich z dwóch niezależnych od siebie grup. Wykorzystujemy go gdy chcemy porównać dwie grupy pod względem jakiejś zmiennej ilościowej. Na przykład gdy chcemy porównać średni wzrost kobiet i mężczyzn w danej grupie.\n",
    "\n",
    "Zazwyczaj dwie średnie z różnych od siebie grup będą się różnić. Test t Studenta powie nam jednak czy owe różnice są istotne statystycznie – czy nie są przypadkowe. Hipoteza zerowa takiego testu będzie brzmiała H0: Średni wzrost w grupie mężczyzn jest taki sam jak średni w grupie kobiet. Hipoteza alternatywna z kolei H1: Kobiety będą różnić się od mężczyzn pod wzrostu.\n",
    "Jeśli wynik testu t Studenta będzie istotny na poziomie p < 0,05 możemy odrzucić hipotezę zerową na rzecz hipotezy alternatywnej.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 159,
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "def t_stat_ind(sample_1, sample_2):\n",
    "    \"\"\"Funkcja oblicza wartość statystyki testowej dla dwóch próbek niezależnych\"\"\"\n",
    "    if sample_1.empty or sample_2.empty:\n",
    "        raise Exception(\"Empty sample\")\n",
    "    sample_1 = sample_1[0].values.tolist()\n",
    "    sample_2 = sample_2[0].values.tolist()\n",
    "    sed = sqrt(sem(sample_1)**2 + sem(sample_2)**2)\n",
    "    return (mean(sample_1) - mean(sample_2)) / sed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 162,
   "metadata": {},
   "outputs": [],
   "source": [
    "def df_ind(sample_1, sample_2):\n",
    "    \"\"\"Funkcja oblicza stopnie swobody dla dwóch próbek niezależnych\"\"\"\n",
    "    return len(sample_1) + len(sample_2) - 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 167,
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "def bootstrap_independent(sample_1, sample_2):\n",
    "    return t_test(\n",
    "        sample_1=sample_1,\n",
    "        sample_2=sample_2,\n",
    "        df_fn=df_ind,\n",
    "        t_stat_fn=t_stat_ind\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TODO: Wyciągnąć wysokości kobiet i mężczyzn oraz poprawić kod aby można było podawać kolumny\n",
    "t_stat, df, cv, p, _ = bootstrap_independent(dataset, dataset)\n",
    "pretty_print_full_stats(t_stat, df, cv, p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "TODO: Wniosek"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Test t studenta dla prób zależnych\n",
    "\n",
    "W odróżnieniu od testu t – Studenta dla prób niezależnych, gdzie porównujemy dwie grupy, ten rodzaj testu stosujemy gdy poddajemy analizie tą samą pojedynczą grupę, ale dwukrotnie w czasie. Na przykład gdy chcemy porównać średnie wagi grupy osób przed dietą oraz po diecie, aby sprawdzić czy dieta spowodowała istotne zmiany statystyczne.\n",
    "\n",
    "Hipoteza zerowa takiego testu będzie brzmiała H0: Średnia waga osób po diecie jest taka sama jak przed dietą. Hipoteza alternatywna z kolei H1: Dieta znacząco wpłynęła na średnią wagę danej grupy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 160,
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "def t_stat_dep(sample_1, sample_2, mu=0):\n",
    "    \"\"\"Funkcja oblicza wartość statystyki testowej dla dwóch próbek zależnych\"\"\"\n",
    "    if sample_1.empty or sample_2.empty:\n",
    "        raise Exception(\"Empty sample\")\n",
    "    sample_1 = sample_1[0].values.tolist()\n",
    "    sample_2 = sample_2[0].values.tolist()\n",
    "    differences = [x_1 - x_2 for x_1, x_2 in zip(sample_1, sample_2)]\n",
    "    sample_size = len(sample_1)\n",
    "    return (mean(differences) - mu) / (stdev(differences) / sqrt(sample_size))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 161,
   "metadata": {},
   "outputs": [],
   "source": [
    "def df_dep(sample_1, sample_2):\n",
    "    \"\"\"Funkcja oblicza stopnie swobody dla dwóch próbek zależnych\"\"\"\n",
    "    l1, l2 = len(sample_1), len(sample_2)\n",
    "    if l1 != l2:\n",
    "        raise Exception(\"Samples aren't of equal length\")\n",
    "    return l1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 168,
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "def bootstrap_dependent(sample_1, sample_2):\n",
    "    return t_test(\n",
    "        sample_1=sample_1,\n",
    "        sample_2=sample_2,\n",
    "        df_fn=df_dep,\n",
    "        t_stat_fn=t_stat_dep\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TODO: Wyciągnąć wagi przed dietą i po oraz poprawić kod aby można było podawać kolumny\n",
    "t_stat, df, cv, p, _ = bootstrap_dependent(dataset, dataset)\n",
    "pretty_print_full_stats(t_stat, df, cv, p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "TODO: Wniosek"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Wykresy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 171,
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "def draw_distribution(stats):\n",
    "    \"\"\"\n",
    "    Funkcja rysuje rozkład statystyki testowej\n",
    "    @param stats: lista statystyk testowych\n",
    "    \"\"\"\n",
    "    plt.hist(stats)\n",
    "    plt.xlabel('Test statistic value')\n",
    "    plt.ylabel('Frequency')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Testy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Statystyka testowa dla jednej próby:\n",
      "1.414213562373095 - z naszej funkcji\n",
      "[1.41421356] - z gotowej biblioteki\n",
      "\n",
      "Statystyka testowa dla dwóch prób niezależnych:\n",
      "-3.0 - z naszej funkcji\n",
      "[-3.] - z gotowej biblioteki\n",
      "\n",
      "Statystyka testowa dla dwóch prób zależnych:\n",
      "-1.6329931618554525 - z naszej funkcji\n",
      "[-1.63299316] - z gotowej biblioteki\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Testy dla samych statystyk testowych\n",
    "def pretty_print_stats(t_stat_selfmade, t_stat_lib, suffix):\n",
    "    print(f'Statystyka testowa dla {suffix}:')\n",
    "    print(t_stat_selfmade, '- z naszej funkcji')\n",
    "    print(t_stat_lib, '- z gotowej biblioteki')\n",
    "    print()\n",
    "    \n",
    "dummy = pd.DataFrame([1, 2, 3, 4, 5])\n",
    "dummy2 = pd.DataFrame([4, 5, 6, 7, 8])\n",
    "dummy3 = pd.DataFrame([1, 3 , 3, 4, 6])\n",
    "\n",
    "t_stat_selfmade = t_stat_single(dummy, 2)\n",
    "t_stat_lib, _ = ttest_1samp(dummy, 2)\n",
    "pretty_print_stats(t_stat_selfmade, t_stat_lib, 'jednej próby')\n",
    "\n",
    "t_stat_selfmade = t_stat_ind(dummy, dummy2)\n",
    "t_stat_lib, _ = ttest_ind(dummy, dummy2)\n",
    "pretty_print_stats(t_stat_selfmade, t_stat_lib, 'dwóch prób niezależnych')\n",
    "\n",
    "t_stat_selfmade = t_stat_dep(dummy, dummy3)\n",
    "t_stat_lib, _ = ttest_rel(dummy, dummy3)\n",
    "pretty_print_stats(t_stat_selfmade, t_stat_lib, 'dwóch prób zależnych')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'pandas.core.frame.DataFrame'>\n",
      "Statystyki dla jednej próby:\n",
      "t: 1.8073147056683616, df: 5, cv: 2.015048372669157, p: 0.13052275003443325\n",
      "\n",
      "Statystyki dla dwóch prób zależnych:\n",
      "t: 3.0790273716290404, df: 5, cv: 2.015048372669157, p: 0.027500015466573435\n",
      "\n",
      "Statystyki dla dwóch prób niezależnych:\n",
      "t: 2.8109511013364576, df: 8, cv: 1.8595480375228421, p: 0.02280961069987497\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Testy z bootstrappowaniem\n",
    "\n",
    "def pretty_print_full_stats(t_stat, df, cv, p):\n",
    "    print(f't: {t_stat}, df: {df}, cv: {cv}, p: {p}\\n')\n",
    "\n",
    "print(type(dummy))\n",
    "\n",
    "print('Statystyki dla jednej próby:')\n",
    "t_stat, df, cv, p, _ = bootstrap_one_sample(dummy, 2)\n",
    "pretty_print_full_stats(t_stat, df, cv, p)\n",
    "\n",
    "print('Statystyki dla dwóch prób zależnych:')\n",
    "t_stat, df, cv, p, _ = bootstrap_dependent(dummy2, dummy3)\n",
    "pretty_print_full_stats(t_stat, df, cv, p)\n",
    "\n",
    "print('Statystyki dla dwóch prób niezależnych:')\n",
    "t_stat, df, cv, p, _ = bootstrap_independent(dummy2, dummy3)\n",
    "pretty_print_full_stats(t_stat, df, cv, p)"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "11938c6bc6919ae2720b4d5011047913343b08a43b18698fd82dedb0d4417594"
  },
  "kernelspec": {
   "display_name": "Python 3.9.1 64-bit",
   "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.9.1"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}