Split System 2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 123,
   "id": "0308446a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 128,
   "id": "252be1e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.cluster.hierarchy import linkage, fcluster\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import stats\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from scipy.stats import mannwhitneyu, ttest_ind, levene\n",
    "from itertools import combinations\n",
    "import matplotlib\n",
    "import seaborn as sns\n",
    "from sklearn.cluster import KMeans"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0bf01c7c",
   "metadata": {},
   "source": [
    "# Условие"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a6f402e",
   "metadata": {},
   "source": [
    "**В этом домашнем задании вы должны написать свою сплит систему для проведения A/B тестов**\n",
    "\n",
    "В предоставленных данных есть 50000 юзеров, для каждого есть запись какой-то метрики за 2 года (730 дней).\n",
    "\n",
    "Данные: https://disk.yandex.ru/d/2mH8fuKCKLeozg\n",
    "\n",
    "\n",
    "### Ваша задача:\n",
    "\n",
    "1) Написать алгоритм, которые разбивает данные на 3 группы размера около 5000, которые будут считаться \"похожими\".\n",
    "\n",
    "- Для этого вам самостоятельно нужно решить, какие критерии схожести тут применимы и построить необходимые тесты для выбора этих групп.\n",
    "\n",
    "- Данные 3 группы должны считаться похожими на основе ваших метрик.\n",
    "\n",
    "- Разбиение затем будет протестировано на следующих 2х годах нами.\n",
    "\n",
    "2) Визуализировать распределения выбранных вами метрик для найденных групп и сравнить с аналогичными разультатами для случайных групп"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b289ee8c",
   "metadata": {},
   "source": [
    "### Нужно реализовать следующие функции (вы можете добавлять какие-то побочные, если считаете нужным):\n",
    "\n",
    "`generate_groups`  - функция, которая по данным и числу групп делает случайное разбиение данных на примерно одинаковые части по количеству групп\n",
    "\n",
    "`conduct_test_on_pair` - функция проводит статистические тесты между двумя данными группами и возвращает нужные статистики и результаты\n",
    "\n",
    "`conduct_tests` -- по списку групп и количеству необходимых похожих групп (зафиксировано в условии задачи как 3) перебирает возможные тройки групп пока не найдет похожие, критерий остановки за вами\n",
    "\n",
    "`find_groups` -- берет данные и число групп, на которое вы их разбиваете и ищет 3 группы удовлетворяющие вашим критериям схожести.\n",
    "\n",
    "Рекомендуем начинать с того, чтобы искать разбиение на 2 группы, а не 3, чтобы проверить, что сплит система работает адекватно. Задача с двумя группами зачастую сильно проще, чем с тремя."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f3a96384",
   "metadata": {},
   "source": [
    "Группы должны состоять из уникальных клиентов без дубликатов, один клиент может быть только в одной группе. То есть должны выбирать группы из датасета без повторений и у групп не должно быть пересечений."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1136d992",
   "metadata": {},
   "source": [
    " Про то, по каким параметрам оценивать группы: в этом задании вам требуется найти группы, в которых распределения **средних значений и дисперсии данной последовательности** по пользователям похожи. То есть вам нужно вычислить среднее значение и дисперсию данной последовательности для каждого пользователя, а затем проверить, что распределения результатов между группами схожи по вашим тестам.\n",
    "\n",
    "\n",
    "**Замечание**: Для разбиения на группы и сравнения похожести можете брать любую метрику и считать любой критерии, но мы будем сравнивать именно метрику «среднее» и дисперсию на следующих по времени данных"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c07b871a",
   "metadata": {},
   "source": [
    "#### Про систему оценивания: При оценивании мы будем проверять полученные вами группы как на этих данных, так и на следующих двух годах. При проверке будем учитывать и качество полученных групп при помощи метрик."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 129,
   "id": "63c1d214",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>2019-01-01</th>\n",
       "      <th>2019-01-02</th>\n",
       "      <th>2019-01-03</th>\n",
       "      <th>2019-01-04</th>\n",
       "      <th>2019-01-05</th>\n",
       "      <th>2019-01-06</th>\n",
       "      <th>2019-01-07</th>\n",
       "      <th>2019-01-08</th>\n",
       "      <th>2019-01-09</th>\n",
       "      <th>2019-01-10</th>\n",
       "      <th>...</th>\n",
       "      <th>2020-12-21</th>\n",
       "      <th>2020-12-22</th>\n",
       "      <th>2020-12-23</th>\n",
       "      <th>2020-12-24</th>\n",
       "      <th>2020-12-25</th>\n",
       "      <th>2020-12-26</th>\n",
       "      <th>2020-12-27</th>\n",
       "      <th>2020-12-28</th>\n",
       "      <th>2020-12-29</th>\n",
       "      <th>2020-12-30</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>-3.533155</td>\n",
       "      <td>-4.475755</td>\n",
       "      <td>-4.748074</td>\n",
       "      <td>-5.569953</td>\n",
       "      <td>-8.017896</td>\n",
       "      <td>-12.396771</td>\n",
       "      <td>-11.853961</td>\n",
       "      <td>-12.046498</td>\n",
       "      <td>-10.297042</td>\n",
       "      <td>-10.075553</td>\n",
       "      <td>...</td>\n",
       "      <td>-168.163849</td>\n",
       "      <td>-168.728405</td>\n",
       "      <td>-168.971555</td>\n",
       "      <td>-168.824042</td>\n",
       "      <td>-167.700081</td>\n",
       "      <td>-167.634364</td>\n",
       "      <td>-167.856116</td>\n",
       "      <td>-168.743465</td>\n",
       "      <td>-168.231590</td>\n",
       "      <td>-167.913551</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>-3.824475</td>\n",
       "      <td>-1.530769</td>\n",
       "      <td>0.038338</td>\n",
       "      <td>-1.475237</td>\n",
       "      <td>-1.602029</td>\n",
       "      <td>-0.992509</td>\n",
       "      <td>0.302906</td>\n",
       "      <td>1.789140</td>\n",
       "      <td>4.177589</td>\n",
       "      <td>9.273372</td>\n",
       "      <td>...</td>\n",
       "      <td>75.919063</td>\n",
       "      <td>72.909418</td>\n",
       "      <td>75.289867</td>\n",
       "      <td>77.770907</td>\n",
       "      <td>74.099941</td>\n",
       "      <td>72.411474</td>\n",
       "      <td>67.472791</td>\n",
       "      <td>67.084800</td>\n",
       "      <td>66.805609</td>\n",
       "      <td>72.078232</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>-1.616434</td>\n",
       "      <td>-4.658834</td>\n",
       "      <td>-6.746498</td>\n",
       "      <td>-9.979424</td>\n",
       "      <td>-11.683396</td>\n",
       "      <td>-9.507703</td>\n",
       "      <td>-7.614438</td>\n",
       "      <td>-9.735666</td>\n",
       "      <td>-10.754863</td>\n",
       "      <td>-11.058539</td>\n",
       "      <td>...</td>\n",
       "      <td>47.759469</td>\n",
       "      <td>46.596598</td>\n",
       "      <td>50.135737</td>\n",
       "      <td>48.424880</td>\n",
       "      <td>47.740658</td>\n",
       "      <td>46.676116</td>\n",
       "      <td>45.611939</td>\n",
       "      <td>43.331310</td>\n",
       "      <td>45.655325</td>\n",
       "      <td>39.770838</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>-0.157282</td>\n",
       "      <td>1.827617</td>\n",
       "      <td>-2.585284</td>\n",
       "      <td>-4.397189</td>\n",
       "      <td>0.113494</td>\n",
       "      <td>7.614615</td>\n",
       "      <td>5.456654</td>\n",
       "      <td>5.928451</td>\n",
       "      <td>6.953534</td>\n",
       "      <td>4.514409</td>\n",
       "      <td>...</td>\n",
       "      <td>318.364962</td>\n",
       "      <td>319.281629</td>\n",
       "      <td>322.614999</td>\n",
       "      <td>325.200268</td>\n",
       "      <td>324.745069</td>\n",
       "      <td>326.519830</td>\n",
       "      <td>326.467326</td>\n",
       "      <td>321.806523</td>\n",
       "      <td>321.286907</td>\n",
       "      <td>323.570847</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2.056969</td>\n",
       "      <td>10.187929</td>\n",
       "      <td>8.559520</td>\n",
       "      <td>8.579555</td>\n",
       "      <td>9.464262</td>\n",
       "      <td>9.038088</td>\n",
       "      <td>10.875990</td>\n",
       "      <td>9.883777</td>\n",
       "      <td>8.026472</td>\n",
       "      <td>13.040411</td>\n",
       "      <td>...</td>\n",
       "      <td>-38.079584</td>\n",
       "      <td>-38.515039</td>\n",
       "      <td>-38.031876</td>\n",
       "      <td>-36.610698</td>\n",
       "      <td>-37.753387</td>\n",
       "      <td>-35.470938</td>\n",
       "      <td>-35.454435</td>\n",
       "      <td>-35.026979</td>\n",
       "      <td>-35.633306</td>\n",
       "      <td>-35.460477</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 730 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   2019-01-01  2019-01-02  2019-01-03  2019-01-04  2019-01-05  2019-01-06  \\\n",
       "0   -3.533155   -4.475755   -4.748074   -5.569953   -8.017896  -12.396771   \n",
       "1   -3.824475   -1.530769    0.038338   -1.475237   -1.602029   -0.992509   \n",
       "2   -1.616434   -4.658834   -6.746498   -9.979424  -11.683396   -9.507703   \n",
       "3   -0.157282    1.827617   -2.585284   -4.397189    0.113494    7.614615   \n",
       "4    2.056969   10.187929    8.559520    8.579555    9.464262    9.038088   \n",
       "\n",
       "   2019-01-07  2019-01-08  2019-01-09  2019-01-10  ...  2020-12-21  \\\n",
       "0  -11.853961  -12.046498  -10.297042  -10.075553  ... -168.163849   \n",
       "1    0.302906    1.789140    4.177589    9.273372  ...   75.919063   \n",
       "2   -7.614438   -9.735666  -10.754863  -11.058539  ...   47.759469   \n",
       "3    5.456654    5.928451    6.953534    4.514409  ...  318.364962   \n",
       "4   10.875990    9.883777    8.026472   13.040411  ...  -38.079584   \n",
       "\n",
       "   2020-12-22  2020-12-23  2020-12-24  2020-12-25  2020-12-26  2020-12-27  \\\n",
       "0 -168.728405 -168.971555 -168.824042 -167.700081 -167.634364 -167.856116   \n",
       "1   72.909418   75.289867   77.770907   74.099941   72.411474   67.472791   \n",
       "2   46.596598   50.135737   48.424880   47.740658   46.676116   45.611939   \n",
       "3  319.281629  322.614999  325.200268  324.745069  326.519830  326.467326   \n",
       "4  -38.515039  -38.031876  -36.610698  -37.753387  -35.470938  -35.454435   \n",
       "\n",
       "   2020-12-28  2020-12-29  2020-12-30  \n",
       "0 -168.743465 -168.231590 -167.913551  \n",
       "1   67.084800   66.805609   72.078232  \n",
       "2   43.331310   45.655325   39.770838  \n",
       "3  321.806523  321.286907  323.570847  \n",
       "4  -35.026979  -35.633306  -35.460477  \n",
       "\n",
       "[5 rows x 730 columns]"
      ]
     },
     "execution_count": 129,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.read_csv('train_data.csv', delimiter=',')\n",
    "df.head() # Посомтрим, как выглядят наши данные"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b79446e5",
   "metadata": {},
   "source": [
    "## Напишите следующие функции"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "id": "f2c22ee2",
   "metadata": {},
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "Cannot take a larger sample than population when 'replace=False'",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "Input \u001b[1;32mIn [131]\u001b[0m, in \u001b[0;36m<cell line: 23>\u001b[1;34m()\u001b[0m\n\u001b[0;32m     23\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m i \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mrange\u001b[39m(group_num):\n\u001b[0;32m     24\u001b[0m     cluster_indices \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39mwhere(cluster_labels \u001b[38;5;241m==\u001b[39m i)[\u001b[38;5;241m0\u001b[39m] \u001b[38;5;66;03m# Пользователи, принадлежащие кластеру i\u001b[39;00m\n\u001b[1;32m---> 25\u001b[0m     selected_indices \u001b[38;5;241m=\u001b[39m \u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mrandom\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mchoice\u001b[49m\u001b[43m(\u001b[49m\u001b[43mcluster_indices\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mgroup_size\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mreplace\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mFalse\u001b[39;49;00m\u001b[43m)\u001b[49m \u001b[38;5;66;03m# Выборка размера group_size из кластера i\u001b[39;00m\n\u001b[0;32m     26\u001b[0m     group_data \u001b[38;5;241m=\u001b[39m df\u001b[38;5;241m.\u001b[39miloc[selected_indices]\n\u001b[0;32m     27\u001b[0m     groups\u001b[38;5;241m.\u001b[39mappend(group_data)\n",
      "File \u001b[1;32mmtrand.pyx:965\u001b[0m, in \u001b[0;36mnumpy.random.mtrand.RandomState.choice\u001b[1;34m()\u001b[0m\n",
      "\u001b[1;31mValueError\u001b[0m: Cannot take a larger sample than population when 'replace=False'"
     ]
    }
   ],
   "source": [
    "# Сначала мы сделаем кластеризацию, чтобы выявить похожие группы\n",
    "# M(x) и D(x) для каждого пользователя\n",
    "user_means = df.mean(axis=1)\n",
    "user_variances = df.var(axis=1)\n",
    "group_num=10\n",
    "\n",
    "# Массив признаков для каждого пользователя, состоящий из ключевых метрик: матожидания и дисперсии по времени\n",
    "features = np.column_stack((user_means, user_variances))\n",
    "\n",
    "# Кластеризация KMeans\n",
    "kmeans = KMeans(n_clusters=group_num, random_state=42)\n",
    "kmeans.fit(features)\n",
    "\n",
    "# Метки кластеров для каждого пользователя\n",
    "cluster_labels = kmeans.labels_\n",
    "\n",
    "# Число пользователей в каждой группе\n",
    "group_size = len(df) // group_num\n",
    "\n",
    "groups = []\n",
    "\n",
    "# Делим пользователей на группы\n",
    "for i in range(group_num):\n",
    "    cluster_indices = np.where(cluster_labels == i)[0] # Пользователи, принадлежащие кластеру i\n",
    "    selected_indices = np.random.choice(cluster_indices, group_size, replace=False) # Выборка размера group_size из кластера i\n",
    "    group_data = df.iloc[selected_indices]\n",
    "    groups.append(group_data)\n",
    "\n",
    "# Получим кластеризированные группы, далее увидим, какие кластеры выделились и осознаем, что наибольший кластер - это под индексом 2, остальные кластеры содержат мало наблюдений - значит они сильно отличаются от других кластеров, эти наблюдения стоит убрать из нашей выборки\n",
    "\n",
    "def conduct_tests_on_pair(df_1, df_2):\n",
    "    \"\"\"\n",
    "    Conducts statistical test or tests on two given groups to decide if they are similar\n",
    "\n",
    "    df_i : Dataframe of users of a single group\n",
    "\n",
    "    returns p-value or anything you need further\n",
    "    \"\"\"\n",
    "    means_df_1 = df_1.mean(axis=1)  # Сначала посчитаем среднее значение по времени для каждого пользователя в каждой группе, потом сравним средние значения пользователей двух групп между собой на равенство мат ожидания двух групп\n",
    "    means_df_2 = df_2.mean(axis=1)\n",
    "\n",
    "    var_df_1 = df_1.var(axis=1) # То же самое с дисперсиями\n",
    "    var_df_2 = df_2.var(axis=1)\n",
    "\n",
    "    stat_tt, p_tt = stats.ttest_ind(means_df_1, means_df_2)\n",
    "    stat_mw, p_mw = mannwhitneyu(means_df_1, means_df_2)\n",
    "    stat_lev, p_lev = levene(var_df_1, var_df_2) # На дисперсию\n",
    "\n",
    "    # print(f'Двусторонний t-тест, тестовая статистика: {stat_tt}, p-value: {p_tt}')\n",
    "    # print(f'Тест Mann-Witney, тестовая статистика: {stat_mw}, p-value: {p_mw}')\n",
    "    # print(f'Тест Левене, тестовая статистика: {stat_lev}, p-value: {p_lev}')\n",
    "    return stat_tt, p_tt, stat_mw, p_mw, stat_lev, p_lev\n",
    "\n",
    "\n",
    "# print(conduct_tests_on_pair(generated_groups[0], generated_groups[1]))\n",
    "\n",
    "def conduct_tests(dfs, group_num):\n",
    "    \"\"\"\n",
    "    Goal is to find out if there is group_num groups in the dfs list which are similar enough according to your criterion.\n",
    "    For every (or maybe not every) selection of group_num groups from the dfs list and conducts pairwise tests on selected subset using conduct_tests_on_pair.\n",
    "\n",
    "    dfs : list of dataframes of groups\n",
    "\n",
    "    returns list of selected groups or None if there are no similar groups\n",
    "    \"\"\"\n",
    "    group_combinations = combinations(dfs, group_num) # комбинации групп для проверки\n",
    "    for combination in group_combinations:\n",
    "        # Стат тесты для проверки похожести: левене для проверки схожести дисперсий\n",
    "        p_values_tt = []\n",
    "        p_values_mw = []\n",
    "        p_values_lev = []\n",
    "        for df1, df2 in combinations(combination, 2):\n",
    "            tt_stat, p_tt, mw_stat, p_mw, lev_stat, p_lev = conduct_tests_on_pair(df1, df2)\n",
    "            p_values_tt.append(p_tt)\n",
    "            p_values_mw.append(p_mw)\n",
    "            p_values_lev.append(p_lev)\n",
    "\n",
    "        # Будем требовать, чтобы все тесты не отвергались, уровень значимости возьмем побольше, чтобы группы были более похожими\n",
    "        if all(p > 0.2 for p in p_values_tt) and all(p > 0.2 for p in p_values_mw) and all(p > 0.2 for p in p_values_lev):\n",
    "            return list(combination)\n",
    "    return None\n",
    "\n",
    "def find_groups(df, group_num):\n",
    "    \"\"\"\n",
    "    Splits df into groups and conducts tests on subsets of the groups until necessary groups are found\n",
    "    \"\"\"\n",
    "    for i in range(100000):\n",
    "        dfs = generate_groups(df, group_num)\n",
    "        result = conduct_tests(dfs, 3)\n",
    "\n",
    "        if result is not None:\n",
    "            return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "id": "da17471c",
   "metadata": {},
   "outputs": [],
   "source": [
    "group1 = generated_groups[:2] + generated_groups[3:]  # Оставили только группы с адекватными значениями средних и дисперсий, которые похожи между собой, согласно кластеризации\n",
    "all_indices = [index for group in group1 for index in group.index]  # индексы наблюдений, которые надо удалить\n",
    "df1 = df.copy().drop(all_indices)  # удалим этих юзеров"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "id": "54ff8adb",
   "metadata": {},
   "outputs": [],
   "source": [
    "high_var_indices = df1[df1.var(axis=1) > 9**6].index  # Удалим данные с большой дисперсией, очень шумная дисперсия ломает тесты, поэтому удалим небольшое число шумных наблюдений\n",
    "high_var_indices.shape\n",
    "df2 = df1.copy().drop(high_var_indices)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "id": "abadb8dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_groups(df, group_num):\n",
    "    \"\"\"\n",
    "    Splits users into group_num groups of similar size\n",
    "\n",
    "    df : Dataframe with all user data\n",
    "    group_num : number of groups to split into\n",
    "\n",
    "    returns a list of dataframes of generated groups\n",
    "    \"\"\"\n",
    "    # permutations обеспечит случайность выбора наших групп, поэтому основные описательные статистики между группами будут равны при большом числе наблюдений, и группы будут похожими\n",
    "    n_users = len(df)\n",
    "    indexes_to_choose = np.random.permutation(n_users) # Создадим последовательность перестановок длины n\n",
    "    group_indexes = np.array_split(indexes_to_choose, group_num) # Поделим данные на равные группы, получим списки индексов для каждой группы (список списков)\n",
    "    groups = []\n",
    "    for indexes in group_indexes:\n",
    "        group = df.iloc[indexes]\n",
    "        groups.append(group)\n",
    "    return groups\n",
    "\n",
    "def conduct_tests_on_pair(df_1, df_2):\n",
    "    \"\"\"\n",
    "    Conducts statistical test or tests on two given groups to decide if they are similar\n",
    "\n",
    "    df_i : Dataframe of users of a single group\n",
    "\n",
    "    returns p-value or anything you need further\n",
    "    \"\"\"\n",
    "    means_df_1 = df_1.mean(axis=1)  # Сначала посчитаем среднее значение по времени для каждого пользователя в каждой группе, потом сравним средние значения пользователей двух групп между собой на равенство мат ожидания двух групп\n",
    "    means_df_2 = df_2.mean(axis=1)\n",
    "\n",
    "    var_df_1 = df_1.var(axis=1) # То же самое с дисперсиями\n",
    "    var_df_2 = df_2.var(axis=1)\n",
    "\n",
    "    stat_tt, p_tt = stats.ttest_ind(means_df_1, means_df_2)\n",
    "    stat_mw, p_mw = mannwhitneyu(means_df_1, means_df_2)\n",
    "    stat_lev, p_lev = levene(var_df_1, var_df_2)\n",
    "\n",
    "    # print(f'Двусторонний t-тест, тестовая статистика: {stat_tt}, p-value: {p_tt}')\n",
    "    # print(f'Тест Mann-Witney, тестовая статистика: {stat_mw}, p-value: {p_mw}')\n",
    "    # print(f'Тест Левене, тестовая статистика: {stat_lev}, p-value: {p_lev}')\n",
    "    return stat_tt, p_tt, stat_mw, p_mw, stat_lev, p_lev\n",
    "\n",
    "\n",
    "# print(conduct_tests_on_pair(generated_groups[0], generated_groups[1]))\n",
    "\n",
    "def conduct_tests(dfs, group_num):\n",
    "    \"\"\"\n",
    "    Goal is to find out if there is group_num groups in the dfs list which are similar enough according to your criterion.\n",
    "    For every (or maybe not every) selection of group_num groups from the dfs list and conducts pairwise tests on selected subset using conduct_tests_on_pair.\n",
    "\n",
    "    dfs : list of dataframes of groups\n",
    "\n",
    "    returns list of selected groups or None if there are no similar groups\n",
    "    \"\"\"\n",
    "    group_combinations = combinations(dfs, group_num) # комбинации групп для проверки\n",
    "    for combination in group_combinations:\n",
    "        # Стат тесты для проверки похожести: левене для проверки схожести дисперсий\n",
    "        p_values_tt = []\n",
    "        p_values_mw = []\n",
    "        p_values_lev = []\n",
    "        for df1, df2 in combinations(combination, 2):\n",
    "            tt_stat, p_tt, mw_stat, p_mw, lev_stat, p_lev = conduct_tests_on_pair(df1, df2)\n",
    "            p_values_tt.append(p_tt)\n",
    "            p_values_mw.append(p_mw)\n",
    "            p_values_lev.append(p_lev)\n",
    "\n",
    "        # Будем требовать, чтобы все тесты не отвергались, уровень значимости возьмем побольше, чтобы группы были более похожими\n",
    "        if all(p > 0.2 for p in p_values_tt) and all(p > 0.2 for p in p_values_mw) and all(p > 0.2 for p in p_values_lev):\n",
    "            return list(combination)\n",
    "    return None\n",
    "\n",
    "def find_groups(df, group_num):\n",
    "    \"\"\"\n",
    "    Splits df into groups and conducts tests on subsets of the groups until necessary groups are found\n",
    "    \"\"\"\n",
    "    for i in range(100000):\n",
    "        dfs = generate_groups(df, group_num)\n",
    "        result = conduct_tests(dfs, 3)\n",
    "\n",
    "        if result is not None:\n",
    "            return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "id": "e992a5ea",
   "metadata": {},
   "outputs": [],
   "source": [
    "# groups = generate_groups(df1, 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "id": "e8035f63",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAD4CAYAAAAQP7oXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAiRklEQVR4nO3dcXAU9d0/8PfeXRI8ELjsaWIkDk0gT2sfqjkP0VRj09xcZ9o+Nk+Hh2lGK5rJxKenZtTiIIjijAbDD2g60Th2mniD0BGlA+pMp5jnpJopsSVHLjNSq0lQKujFmLsAMQjhst/fH5hrLhduv+Am2YP36x9u976f2/cuHJ/bvdtdRQghQEREJMEy0wGIiCh9sGkQEZE0Ng0iIpLGpkFERNLYNIiISBqbBhERSbPNdIDz9dlnn810BF1OpxMDAwMzHUNXOuRMh4yAeXPetvk2PPU//fiv7/0Xjh//f6bNORFzGisvL8+w1+KeBhERSWPTICIiaWwaREQkjU2DiIiksWkQEZE0Ng0iIpLGpkFERNLYNIiISBqbBhERSWPTILqE/Lct7S4CQSbDpkFERNLYNIiISBqbBhERSWPTILqEjByPzXQESnNsGkREJE3qpxRdXV3w+/3QNA3l5eWoqKhIeF4IAb/fj1AohKysLPh8PhQUFKSsPXz4MH7/+99jZGQEVqsV1dXVWLRokaErR0RExtLd09A0DS0tLVi7di0aGhqwb98+HD16NGFMKBRCX18fGhsbUVNTg+bmZt3a7du3Y/ny5di0aRNWrFiB7du3T8HqERGRkXSbRm9vL3Jzc5GTkwObzYaSkhJ0dHQkjAkGgygtLYWiKCgqKsLw8DAGBwdT1iqKgq+++goAcPLkSTgcjilYPSIiMpLu4aloNApVVePTqqqip6cnaYzT6UwYE41GU9auXLkSdXV12LZtGzRNw9NPPz3p8gOBAAKBAACgvr4+YTlmZbPZmNMg6ZARMG9Oi9UCRVEwa1YWMjKcUE4qpsw5kVm350TpktNIuk1DCJE0T1EUqTGpaltbW7Fy5UrcdNNNaG9vxwsvvIDHH388abzH44HH44lPp8P9eNPlvsHpkDMdMgLmzamNahBC4NSp0zh+fADiMnPmnMis23OidMk5rfcIV1UVkUgkPh2JRJIOJamqmrDhxsakqn3nnXewbNkyAMDNN9+M3t7eb7YmREQ05XSbRmFhIcLhMPr7+xGLxdDe3g63250wxu12o62tDUIIdHd3w263w+FwpKzNzs7G+++/DwA4ePAgcnNzp2D1iIjISLqHp6xWK6qqqlBXVwdN01BWVob8/Hy0trYCALxeL4qLi9HZ2Yna2lpkZmbC5/OlrAWAe++9N/5T3IyMDNx7771TuJpERGQEqfM0XC4XXC5Xwjyv1xt/rCgKqqurpWsB4Nvf/jY2btx4PlmJiGiG8YxwIiKSxqZBRETS2DSIiEgamwYREUlj0yAiImlsGkREJI1Ng4iIpLFpEBGRNDYNIiKSxqZBRETS2DSIiEgamwYREUlj0yAiImlsGkREJI1Ng4iIpLFpEBGRNKmbMHV1dcXvsldeXo6KioqE54UQ8Pv9CIVCyMrKgs/nQ0FBQcrahoYGfPbZZwCAkydPwm63Y9OmTcatGRERGU63aWiahpaWFqxbtw6qqmLNmjVwu91YsGBBfEwoFEJfXx8aGxvR09OD5uZmbNiwIWXtQw89FK9/6aWXYLfbp2YNiYjIMLqHp3p7e5Gbm4ucnBzYbDaUlJSgo6MjYUwwGERpaSkURUFRURGGh4cxODgoVSuEwLvvvovvf//7xq4ZEREZTndPIxqNQlXV+LSqqujp6Uka43Q6E8ZEo1Gp2n/+85+YN28errrqqkmXHwgEEAgEAAD19fUJyzErm83GnAZJh4yAeXNarBYoioJZs7KQkeGEclIxZc6JzLo9J0qXnEbSbRpCiKR5iqJIjZGp3bdvX8q9DI/HA4/HE58eGBjQizzjnE4ncxokHTIC5s2pjWoQQuDUqdM4fnwA4jJz5pzIrNtzonTJmZeXZ9hr6R6eUlUVkUgkPh2JROBwOJLGjN9wY2P0akdHR7F//36UlJR8o5UgIqLpods0CgsLEQ6H0d/fj1gshvb2drjd7oQxbrcbbW1tEEKgu7sbdrsdDodDt/a9995DXl5ewiEsIiIyL93DU1arFVVVVairq4OmaSgrK0N+fj5aW1sBAF6vF8XFxejs7ERtbS0yMzPh8/lS1o7ROzRFRETmInWehsvlgsvlSpjn9XrjjxVFQXV1tXTtmPvuu082JxERmQDPCCciImlsGkREJI1Ng+gSMXRwZKYj0EWATYPoElFYuQw7rl040zEozbFpEBGRNDYNIiKSxqZBRETS2DSIiEgamwYREUlj0yAiImlsGkREJI1Ng4iIpLFpEBGRNDYNIiKSxqZBRETSpO6n0dXVBb/fD03TUF5ejoqKioTnhRDw+/0IhULIysqCz+dDQUGBbu2f//xn7NmzB1arFS6XC3feeadhK0ZERMbTbRqapqGlpQXr1q2DqqpYs2YN3G43FixYEB8TCoXQ19eHxsZG9PT0oLm5GRs2bEhZe/DgQQSDQWzevBkZGRk4fvz4lK4oERF9c7qHp3p7e5Gbm4ucnBzYbDaUlJSgo6MjYUwwGERpaSkURUFRURGGh4cxODiYsra1tRU/+9nPkJGRAQCYN2/eFKweEREZSXdPIxqNQlXV+LSqqujp6Uka43Q6E8ZEo9GUteFwGB988AF27NiBjIwM/PKXv8SiRYu+8QoREdHU0W0aQoikeYqiSI1JVatpGr788kvU1dXh0KFDaGhowHPPPZf02oFAAIFAAABQX1+f0JzMymazMadB0iEjYN6cFqsFiqJg1qwswHL2vWXGnBOZdXtOlC45jaTbNFRVRSQSiU9HIhE4HI6kMQMDA0ljYrHYOWuzs7OxbNkyKIqCRYsWwWKxYGhoCHPnzk14bY/HA4/HE58evxyzcjqdzGmQdMgImDenNqpBCIFTp04jSzv7Ic6MOScy6/acKF1y5uXlGfZaut9pFBYWIhwOo7+/H7FYDO3t7XC73Qlj3G432traIIRAd3c37HY7HA5HytqlS5fi4MGDAIDPPvsMsVgMl19+uWErRkRExtPd07BaraiqqkJdXR00TUNZWRny8/PR2toKAPB6vSguLkZnZydqa2uRmZkJn8+XshYAfvjDH+L555/Hr3/9a9hsNtx3331Jh6aIiMhcpM7TcLlccLlcCfO8Xm/8saIoqK6ulq4Fzh4LrK2tPZ+sREQ0w3hGOBERSWPTICIiaWwaREQkjU2DiIiksWkQEZE0Ng2iS8zK7OyZjkBpjE2DiIiksWkQEZE0Ng0iIpLGpkFERNLYNIiISBqbBhERSWPTICIiaWwaREQkjU2DiIiksWkQEZE0qZswdXV1we/3Q9M0lJeXo6KiIuF5IQT8fj9CoRCysrLg8/lQUFCQsvbVV1/FW2+9Fb8neGVl5aQ3ayIiIvPQbRqapqGlpQXr1q2DqqpYs2YN3G43FixYEB8TCoXQ19eHxsZG9PT0oLm5GRs2bNCt/clPfoLbb7996taOiIgMpXt4qre3F7m5ucjJyYHNZkNJSQk6OjoSxgSDQZSWlkJRFBQVFWF4eBiDg4NStURElD509zSi0ShUVY1Pq6qKnp6epDFOpzNhTDQa1a1988030dbWhoKCAtx1112YM2dO0vIDgQACgQAAoL6+PmE5ZmWz2ZjTIOmQETBvTovVAkVRMGtWFmBRAACZmZmmzDqeWbfnROmS00i6TUMIkTRPURSpMalqvV4vli9fDgB45ZVX8NJLL8Hn8yWN93g88Hg88emBgQG9yDPO6XQyp0HSISNg3pzaqAYhBE6dOo0s7ez7cWRkBAPR6AwnS82s23OidMmZl5dn2GvpHp5SVRWRSCQ+HYlE4HA4ksaM33BjY1LVzp8/HxaLBRaLBeXl5Th06NA3XhkiIppauk2jsLAQ4XAY/f39iMViaG9vh9vtThjjdrvR1tYGIQS6u7tht9vhcDhS1g4ODsbr9+/fj/z8fINXjYiIjKZ7eMpqtaKqqgp1dXXQNA1lZWXIz89Ha2srgLOHmYqLi9HZ2Yna2lpkZmbGDzOdqxYAtm/fjsOHD0NRFFxxxRWoqamZwtUkIiIjSJ2n4XK5ks6h8Hq98ceKoqC6ulq6FgAeeOCB88lJRAaJDWkzHYHSGM8IJ7oEvJ2VhWMZZ9/uO5cUznAaSmdsGkREJI1Ng4iIpLFpEBGRNDYNIiKSxqZBRETS2DSIiEgamwYREUlj0yAiImlsGkREJI1Ng4iIpLFpEBGRNDYNIiKSxqZBRETS2DSIiEgamwYREUmTuglTV1cX/H4/NE1DeXk5KioqEp4XQsDv9yMUCiErKws+nw8FBQVStW+88Qa2b9+O5uZmzJ0715CVIiKiqaG7p6FpGlpaWrB27Vo0NDRg3759OHr0aMKYUCiEvr4+NDY2oqamBs3NzVK1AwMDeO+99+B0Og1eLSIimgq6TaO3txe5ubnIycmBzWZDSUkJOjo6EsYEg0GUlpZCURQUFRVheHgYg4ODurVbt27FHXfcAUVRjF8zIiIynO7hqWg0ClVV49OqqqKnpydpzPi9BVVVEY1GU9YGg0FkZ2dj4cKFKZcfCAQQCAQAAPX19WmxV2Kz2ZjTIOmQETBvTovVAkVRYP36TwBQFMWUWccz6/acKF1yGkm3aQghkuZN3DM415hzzT99+jR27dqFdevW6Qb0eDzweDzx6YGBAd2ameZ0OpnTIOmQETBvTm1UgxACo1//Ofa+NGPW8cy6PSdKl5x5eXmGvZZu01BVFZFIJD4diUTgcDiSxozfcGNjYrHYpLWff/45+vv78cgjj8Tnr169Gs888wzmz5//TdeJiIimiO53GoWFhQiHw+jv70csFkN7ezvcbnfCGLfbjba2Nggh0N3dDbvdDofDcc7aa665Bs3NzWhqakJTUxNUVcXGjRvZMIiITE53T8NqtaKqqgp1dXXQNA1lZWXIz89Ha2srAMDr9aK4uBidnZ2ora1FZmYmfD5fyloiIkpPUudpuFwuuFyuhHlerzf+WFEUVFdXS9dO1NTUJBODiIhmGM8IJyIiaWwaREQkjU2D6BIgziT//J3oQrBpEBGRNDYNIiKSxqZBRETS2DSIiEgamwYREUlj0yAiImlsGkREJI1Ng4iIpLFpEBGRNDYNIiKSxqZBRETS2DSIiEgamwYREUmTuglTV1cX/H4/NE1DeXk5KioqEp4XQsDv9yMUCiErKws+nw8FBQUpa3fs2IFgMAhFUTBv3jz4fD5kZ2cbunJERGQs3T0NTdPQ0tKCtWvXoqGhAfv27cPRo0cTxoRCIfT19aGxsRE1NTVobm7Wrb399tuxefNmbNq0CS6XC3/84x+nYPWIiMhIuk2jt7cXubm5yMnJgc1mQ0lJCTo6OhLGBINBlJaWQlEUFBUVYXh4GIODgylr7XZ7vP706dNQFMXgVSMiIqPpHp6KRqNQVTU+raoqenp6ksY4nc6EMdFoVLf25ZdfRltbG+x2O9avXz/p8gOBAAKBAACgvr4+YTlmZbPZmNMg6ZARMG9Oi9UCRVGgWJT4BzNFUUyZdTyzbs+J0iWnkXSbhhDJd/yauFdwrjF6tZWVlaisrMTu3buxZ88erFixImm8x+OBx+OJTw8MDOhFnnFOp5M5DZIOGQHz5tRGNQghIDQBIUT8fWnGrOOZdXtOlC458/LyDHst3cNTqqoiEonEpyORCBwOR9KY8RtubIxMLQDccsst+Pvf/35BK0BERNNHt2kUFhYiHA6jv78fsVgM7e3tcLvdCWPcbjfa2toghEB3dzfsdjscDkfK2nA4HK8PBoOGdkIiIpoauoenrFYrqqqqUFdXB03TUFZWhvz8fLS2tgIAvF4viouL0dnZidraWmRmZsLn86WsBYA//OEPCIfD8eOrNTU1U7iaRERkBKnzNFwuF1wuV8I8r9cbf6woCqqrq6VrAWDVqlXnk5OIiEyAZ4QTEZE0Ng0iIpLGpkFERNLYNIiISBqbBhERSWPTICIiaWwaREQkjU2DiIiksWkQEZE0Ng0iIpLGpkFERNLYNIiISBqbBtElaGV29kxHoDTFpkFERNLYNIguQbEhbaYjUJqSup9GV1cX/H4/NE1DeXk5KioqEp4XQsDv9yMUCiErKws+nw8FBQUpa7dt24YDBw7AZrMhJycHPp8Ps2fPNnTliIjIWLp7GpqmoaWlBWvXrkVDQwP27duHo0ePJowJhULo6+tDY2Mjampq0NzcrFv7ve99D1u2bMHmzZtx1VVXYffu3VOwekSXlqGDIzMdgS5yuk2jt7cXubm5yMnJgc1mQ0lJCTo6OhLGBINBlJaWQlEUFBUVYXh4GIODgylrr7vuOlitVgBAUVERotHoFKweEREZSffwVDQahaqq8WlVVdHT05M0xul0JoyJRqNStQCwd+9elJSUTLr8QCCAQCAAAKivr09YjlnZbDbmNEg6ZARMlHPul3A658QnLVYLFEWBYlGgKAqAs7dnttlscM4zQd5zMM321JEuOY2k2zSEEEnzxv7x6Y2Rqd21axesVituvfXWSZfv8Xjg8Xji0wMDA3qRZ5zT6WROg6RDRsA8Oe8onI8/jMuhjWoQQkBoAkKI+PsyFothYODYzAXVYZbtqSddcubl5Rn2WrqHp1RVRSQSiU9HIhE4HI6kMeM33NgYvdq3334bBw4cQG1tbVIzISIi89FtGoWFhQiHw+jv70csFkN7ezvcbnfCGLfbjba2Nggh0N3dDbvdDofDkbK2q6sLr7/+OlavXo2srKypWTsiIjKU7uEpq9WKqqoq1NXVQdM0lJWVIT8/H62trQAAr9eL4uJidHZ2ora2FpmZmfD5fClrAaClpQWxWAxPPfUUAGDx4sWoqamZqvUkIiIDSJ2n4XK54HK5EuZ5vd74Y0VRUF1dLV0LAM8+++z55CQiIhPgGeFERCSNTYMojVzIyXsn+SMTMhCbBtFFzn7iBDx/3DHTMegiwaZBRETS2DSIiEgamwYREUlj0yAiImlsGkQXiRyXCzuXFM50DLrIsWkQEZE0Ng0iIpLGpkFERNLYNIiISBqbBhERSWPTICIiaWwaRGmksHIZcia51QDRdGHTIEojxzL4lqWZJXUTpq6uLvj9fmiahvLyclRUVCQ8L4SA3+9HKBRCVlYWfD4fCgoKUta+++672LlzJz799FNs2LABhYU8KYkIOHv588v/M3OmYxBNSvdji6ZpaGlpwdq1a9HQ0IB9+/bh6NGjCWNCoRD6+vrQ2NiImpoaNDc369bm5+dj1apV+M53vjMFq0V06bqQe24QydJtGr29vcjNzUVOTg5sNhtKSkrQ0dGRMCYYDKK0tBSKoqCoqAjDw8MYHBxMWbtgwQLk5eVNzVoREdGU0D08FY1GoapqfFpVVfT09CSNcTqdCWOi0ahUrZ5AIIBAIAAAqK+vT1iOWdlsNuY0SDpkBAzOOfdLOJ1zJn1qSFFgUZRJl2WxWKAI4PK5c+P1FuvXnwsVQPn6Dn6KouCP31uEbAU4c/iwMZkNdkn+vacJ3aYhhEiap0y4feS5xsjU6vF4PPB4PPHpgYGB86qfCU6nkzkNkg4ZAWNzDp0YAQZOTfpcphDQhJh0WTmaBgEFQydOxOu1Ue3sk+Ls+3T8+/Jcr2MGl+Lf+1Qy8qiO7uEpVVURiUTi05FIBA6HI2nM+A03NkamloiM9b/FV850BLqI6TaNwsJChMNh9Pf3IxaLob29HW63O2GM2+1GW1sbhBDo7u6G3W6Hw+GQqiWimSHOJB8JINKje3jKarWiqqoKdXV10DQNZWVlyM/PR2trKwDA6/WiuLgYnZ2dqK2tRWZmJnw+X8paANi/fz9efPFFnDhxAvX19Vi4cCEee+yxKVxVokvP6GlNeuwdGfPxhzPHpi4MXRSkztNwuVxwTTgL1ev1xh8rioLq6mrpWgC48cYbceONN55PViIyEE8UpAvBfzVEJrIyO1t3zLHz/DEJkZGk9jSIaHrs+Na3IAQw0nnAkNfrmyX3Fs9xubATimHLpYsX9zSIiEgamwbRRWbnksLzvhIuf0lFstg0iExkkvNhiUyFTYMozZyrsVzI3sLEixvyYoekh02D6CJ0vg2EP78lWfyXQjTNVmZnS/20diakujMg90IIYNMgmnaxIfmztL8JoxsTr2lFAM/TIJp2O5cUQlGA0wd4TgSlH+5pEE0joz79T3YIafz3EscyLHh54bdw9UeHzut1j2VYUp5xzkNUxKZBlCam4nsQ2UNOK7OzsXNJIQorlxmegdILD08RmcTQwRGo+sMMtXNJIeaf0fjrKZLGfylEM8BMJ/HJNIwd3/rWNCShdMA9DaJpJPvLqaGDI7j8PzPPWTv+PIwclwvijJjSvQUzNTmaWdzTIDKJVN8v5Lhc2LmkMD7Nw0k0U6T2NLq6uuD3+6FpGsrLy1FRUZHwvBACfr8foVAIWVlZ8Pl8KCgoSFn75ZdfoqGhAV988QWuuOIKPPTQQ5gzZ46hK0dkZuP3JsYuTT6msHIZlAwFn3d26r6O0RcbFOLsF99bo9Gk59isSPdfgKZpaGlpwdq1a9HQ0IB9+/bh6NGjCWNCoRD6+vrQ2NiImpoaNDc369a+9tprWLJkCRobG7FkyRK89tprxq8dkYlk3XBDwt5Cqj2LsZ++jv1iahDJP4PNdN0Q/+ntVPxnvjI7G0MHR/gzW0qgu6fR29uL3Nxc5OTkAABKSkrQ0dGBBQsWxMcEg0GUlpZCURQUFRVheHgYg4OD+OKLL85Z29HRgSeffBIAcNttt+HJJ5/EnXfeOQWrSDR1Jn4iHzo4gkV33QRxRkDJSPyPfuL5D2MNxAExaVMQAnh5YeovoAehABnG38kv1XKzbrgBQgD/897Zc0DGN0IA8V9jVR7+GLEhDbbLzza0rdFofHtl3XADgLMnOI41xq3RKHJcLvS+9Dc4f2D4KpFBdJtGNBqFqv77h4CqqqKnpydpjNPpTBgTjUZT1h4/fhwOhwMA4HA4cOLEiUmXHwgEEAgEAAD19fXIy8uTXbcZxZzGMXPG/wOAr/Pl5eUBeQD6+iYdm+rntM4Uz30T49+ps35p7Gv/ZexBXzjpOSe+3jazxs3My/v39gr/u+b/xj2Pvj78R3zSvH/v46VLTqPo7tOKSX42oUz4xHSuMTK1ejweD+rr61FfX49HH330vGpnCnMaJx0yAsxpNOY0lpE5dZuGqqqIRCLx6UgkEt9DGD9mYGAgaUyq2nnz5mFwcBAAMDg4iLlz536zNSEioimn2zQKCwsRDofR39+PWCyG9vZ2uN3uhDFutxttbW0QQqC7uxt2ux0OhyNlrdvtxjvvvAMAeOedd7B06dIpWD0iIjKS7ncaVqsVVVVVqKurg6ZpKCsrQ35+PlpbWwEAXq8XxcXF6OzsRG1tLTIzM+Hz+VLWAkBFRQUaGhqwd+9eOJ1OPPzww7phPR7PN1nXacOcxkmHjABzGo05jWVkTkVM9sUDERHRJHimDhERSWPTICIiaaa5YOGrr76Kt956K/4rqsrKSri+Ptt19+7d2Lt3LywWC+655x5cf/31AICPPvoITU1NGBkZQXFxMe655x4oioIzZ87gueeew0cffYTLL78cDz74IK680thbVb7xxhvYvn07mpub45nNlHPHjh0IBoNQFAXz5s2Dz+dD9tcnUZkp57Zt23DgwAHYbDbk5OTA5/Nh9uzZpsv57rvvYufOnfj000+xYcMGFBb++4Q2M+VMRe9yQFPp+eefR2dnJ+bNm4ctW7YASH0pofPdpkYZGBhAU1MTjh07BkVR4PF48OMf/9h0WUdGRrB+/XrEYjGMjo7ipptuwooVK6YnpzCJV155Rbz++utJ848cOSJWrVolRkZGxOeffy7uv/9+MTo6KoQQ4tFHHxUffvih0DRN1NXVic7OTiGEEHv27BG/+93vhBBC/PWvfxW/+c1vDM36xRdfiKefflr86le/EsePHzdlzuHh4fjjP/3pT/HlmC1nV1eXiMViQgghtm3bJrZt22bKnEeOHBGffvqpWL9+vejt7U2Yb6ac5zI6Oiruv/9+0dfXJ86cOSNWrVoljhw5MuXLHfOPf/xDHDp0SDz88MPxedu2bRO7d+8WQgixe/fub/R3b5RoNCoOHTokhBDi5MmTora2Vhw5csR0WTVNE1999ZUQQogzZ86INWvWiA8//HBacpr+8FRHRwdKSkqQkZGBK6+8Erm5uejt7cXg4CC++uorFBUVQVEUlJaWoqOjA8DZy5r84Ac/AADcdNNNOHjw4KQnGl6orVu34o477kjoxmbLabfb449Pnz4dz2q2nNdddx2sVisAoKioCNGvL8lhtpwLFiyY9Mxfs+U8l/GXA7LZbPFL+kyXa6+9NumCpB0dHbjtttsAnL2U0FieC9mmRnE4HPGLrV522WW4+uqrEY1GTZdVURTMmnX2dPvR0VGMjo5CUZRpyWmaw1MA8Oabb6KtrQ0FBQW46667MGfOHESjUSxevDg+Jjs7G9FoFFarNekSJWP/4Yy/fInVaoXdbsfQ0JAhJxAGg0FkZ2dj4cKFCfPNlhMAXn75ZbS1tcFut2P9+vWmzTlm7969KCkpMX3O8dIpp97lgKbbuS4ldCHbdCr09/fj448/xqJFi0yZVdM0rF69Gn19ffjRj36ExYsXT0vOaW0aTz31FI4dO5Y0/xe/+AW8Xi+WL18OAHjllVfw0ksvwefznfMTWKpPZpM9dz7HElPl3L17N9atWyedZ6ZyLl26FJWVlaisrMTu3buxZ88erFixwpQ5AWDXrl2wWq249dZbU+aZ6ZyyeaYy54WYiWVeqAvZpkY7deoUtmzZgrvvvjthr10203RktVgs2LRpE4aHh7F582Z88skn5xxrZM5pbRqPP/641Ljy8nJs3LgRQPJlTKLRKLKzsye9RMnYF71jz6mqitHRUZw8efK87tVxrpyffPIJ+vv78cgjj8SXuXr1ajzzzDOmyjnRLbfcgvr6eqxYscKUOd9++20cOHAATzzxRPw/MjPmnMxM5LwQMpcDmm5jlxJyOBwJlxK6kG1qpFgshi1btuDWW2/FsmXLTJ0VAGbPno1rr70WXV1d05LTNN9pjF2HCgD2798fP3Pc7Xajvb0dZ86cQX9/P8LhMBYtWgSHw4HLLrsM3d3dEEKgra0tfomSG264AW+//TYA4G9/+xu++93vGvKp6pprrkFzczOamprQ1NQEVVWxceNGzJ8/31Q5ASA87iqiwWAwfjzebDm7urrw+uuvY/Xq1cjKyorPN1vOc0mXnDKXA5pu57qU0IVsU6MIIfDCCy/g6quvxk9/+lPTZj1x4gSGh4cBnP0l1XvvvYerr756WnKa5ozwZ599FocPH4aiKLjiiitQU1MT/yS0a9cu/OUvf4HFYsHdd9+N4uJiAMChQ4fw/PPPY2RkBNdffz2qqqqgKApGRkbw3HPP4eOPP8acOXPw4IMPxu/pYaT77rsPzzzzTLybmynn5s2bEQ6HoSgKnE4nampq4p8gzJTzgQceQCwWi3/SXrx4MWpqakyXc//+/XjxxRdx4sQJzJ49GwsXLsRjjz1mupypdHZ2YuvWrfFL+vz85z+f8mWO+e1vf4v3338fQ0NDmDdvHlasWIGlS5eioaEBAwMD8UsJjf07ON9tapQPPvgATzzxBK655pr461ZWVmLx4sWmyvqvf/0LTU1N0DQNQgjcfPPNWL58OYaGhqY8p2maBhERmZ9pDk8REZH5sWkQEZE0Ng0iIpLGpkFERNLYNIiISBqbBhERSWPTICIiaf8fAJF9/X75lEAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "groups = find_groups(df2, 10)\n",
    "x1 = np.linspace(-1, 1, 100)\n",
    "y1 = groups[0].mean(axis=1)\n",
    "y2 = groups[1].mean(axis=1)\n",
    "y3 = groups[2].mean(axis=1)\n",
    "fig, ax = plt.subplots(1, 1)\n",
    "ax.hist(y1, density=True, histtype='bar', bins=300, alpha=0.3, color='mediumorchid')\n",
    "ax.hist(y2, density=True, histtype='bar', bins=300, alpha=0.7, color='cyan')\n",
    "ax.hist(y3, density=True, histtype='bar', bins=300, alpha=0.9, color='red')\n",
    "ax.axvline(np.mean(y1), ymax=1, color='black', alpha=0.8)\n",
    "ax.axvline(np.mean(y2), ymax=1, color='green', alpha=0.8)\n",
    "ax.axvline(np.mean(y3), ymax=1, color='yellow', alpha=0.8)\n",
    "# rect = matplotlib.patches.Rectangle((np.mean(y1)-3*np.std(y1),0), 6*np.std(y1), 0.3, color='black', alpha=0.2)\n",
    "# ax.add_patch(rect)\n",
    "plt.xlim(-5000, 3000)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "id": "4fed1ea4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD8CAYAAAB3u9PLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAVR0lEQVR4nO3df2zU9eHH8ee1By1Y7No7aYPgDBU2XeKg1mGaCXbcmsVshujmNCaLNoaZqgjObYj4Y3Nduk2Gc8rmAsON7A/dIjPfP+ZI5w8mjFCEuglxUGAGtNK1x2+B0t7n+8c+u3BSvAJ39Miej7/u8/m8e/fyfXe+7vO5+3yIBEEQIEn6n1c03AEkSYXBQpAkARaCJClkIUiSAAtBkhSyECRJAESzDViyZAkbN26kvLycRYsWnbQ9CAKWL1/Opk2bKCkpobm5mYkTJ+YlrCQpf7LuIVx33XUsWLDglNs3bdrEBx98wFNPPcXs2bNZunRpTgNKks6NrIVwxRVXUFZWdsrtGzZsYPr06UQiESZPnszhw4fZu3dvTkNKkvIv6yGjbJLJJPF4PL0ci8VIJpNUVFScNLatrY22tjYAWltbz/ahJUk5dNaFMNiVLyKRyKBjE4kEiUQivfz++++f7cPnXTwep6enZ7hjZGXO3Dofcp4PGcGcuTZu3Li83fdZ/8ooFotlTGJvb++geweSpMJ21oVQV1fH6tWrCYKArVu3Mnr0aAtBks5DWQ8ZPfnkk2zZsoWDBw9y1113cfPNN9Pf3w9AY2MjU6dOZePGjcyZM4eRI0fS3Nyc99CSpNzLWghz58792O2RSIQ777wzV3kkScPEM5UlSYCFIEkKWQiSJMBCkCSFLARJEmAhSJJCFoIkCbAQJEkhC0GSBFgIkqSQhSBJAiwESVLIQpAkARaCJClkIUiSAAtBkhSyECRJgIUgSQpZCJIkwEKQJIUsBEkSYCFIkkIWgiQJsBAkSSELQZIEWAiSpJCFIEkCLARJUshCkCQBFoIkKWQhSJIAC0GSFLIQJEmAhSBJClkIkiQAokMZ1NHRwfLly0mlUsycOZNZs2ZlbP/www956qmn6O3tZWBggK985Ss0NDTkI68kKU+yFkIqlWLZsmUsXLiQWCzGgw8+SF1dHePHj0+Pefnllxk/fjzz58/nwIED3HfffVx77bVEo0PqG0lSAch6yKizs5Pq6mqqqqqIRqPU19fT3t6eMSYSiXD06FGCIODo0aOUlZVRVOTRKEk6n2T9CJ9MJonFYunlWCzGtm3bMsZ86Utf4sc//jHf/OY3OXLkCPPmzRu0ENra2mhrawOgtbWVeDx+tvnzLhqNmjOHzJk750NGMOf5JGshBEFw0rpIJJKx/NZbb/HJT36SRx55hD179vD444/z6U9/mtGjR2eMSyQSJBKJ9HJPT8+Z5j5n4vG4OXPInLlzPmQEc+bauHHj8nbfWY/rxGIxent708u9vb1UVFRkjHn11VeZNm0akUiE6upqxo4dy/vvv5/7tJKkvMlaCDU1NXR1ddHd3U1/fz9r166lrq4uY0w8Hucf//gHAPv27eP9999n7Nix+UksScqLrIeMiouLaWpqoqWlhVQqRUNDAxMmTGDVqlUANDY2ctNNN7FkyRK+9a1vAXDbbbdx4YUX5je5JCmnhvS70NraWmprazPWNTY2pm9XVlaycOHC3CaTJJ1T/jZUkgRYCJKkkIUgSQIsBElSyEKQJAEWgiQpZCFIkgALQZIUshAkSYCFIEkKWQiSJMBCkCSFLARJEmAhSJJCFoIkCbAQJEmhIf0DOdL/gv9bt46eY8cAKC0poQz4yjXXDG8o6RyyEKRQz7FjHA9vFwUBPX19w5pHOtc8ZCRJAiwESVLIQpAkARaCJClkIUiSAAtBkhSyECRJgIUgSQpZCJIkwEKQJIUsBEkSYCFIkkIWgiQJsBAkSSELQZIEWAiSpNCQ/oGcjo4Oli9fTiqVYubMmcyaNeukMZs3b+a5555jYGCAMWPG8L3vfS/XWSVJeZS1EFKpFMuWLWPhwoXEYjEefPBB6urqGD9+fHrM4cOHWbp0KQ899BDxeJz9+/fnNbQkKfeyHjLq7OykurqaqqoqotEo9fX1tLe3Z4x54403mDZtGvF4HIDy8vL8pJUk5U3WPYRkMkksFksvx2Ixtm3bljGmq6uL/v5+HnvsMY4cOcL111/PjBkzTrqvtrY22traAGhtbU0XSCGLRqPmzKFCzllaUkJREAAQKSqitKSkYLNCYc/licx5/shaCEH4BjlRJBLJWB4YGGDnzp08/PDD9PX1sXDhQiZNmsS4ceMyxiUSCRKJRHq5p6fnTHOfM/F43Jw5VMg5jx47xvHwdsnIkRzt6yvYrFDYc3kic+bWR/+/mktZCyEWi9Hb25te7u3tpaKi4qQxY8aMobS0lNLSUi6//HLefffdvAaXJOVW1u8Qampq6Orqoru7m/7+ftauXUtdXV3GmLq6Ot555x0GBgY4duwYnZ2dXHzxxXkLLUnKvax7CMXFxTQ1NdHS0kIqlaKhoYEJEyawatUqABobGxk/fjxTpkzhgQceoKioiC984QtccskleQ8vScqdIZ2HUFtbS21tbca6xsbGjOUbbriBG264IXfJJEnnlGcqS5IAC0GSFLIQJEmAhSBJClkIkiTAQpAkhSwESRJgIUiSQhaCJAmwECRJIQtBkgRYCJKkkIUgSQIsBElSyEKQJAEWgiQpZCFIkgALQZIUshAkSYCFIEkKWQiSJMBCkCSFLARJEmAhSJJCFoIkCbAQJEkhC0GSBFgIkqSQhSBJAiwESVLIQpAkARaCJClkIUiSAAtBkhSyECRJwBALoaOjg/vuu497772XP/7xj6cc19nZyde//nXWrVuXq3ySpHMkayGkUimWLVvGggULWLx4MWvWrGH37t2Djvvd737HlClT8pFTkpRnWQuhs7OT6upqqqqqiEaj1NfX097eftK4P/3pT0ybNo0LL7wwL0ElSfkVzTYgmUwSi8XSy7FYjG3btp00Zv369Tz66KP84he/OOV9tbW10dbWBkBrayvxePxMc58z0WjUnDlUyDlLS0ooCgIAIkVFlJaUFGxWKOy5PJE5zx9ZCyEI3yAnikQiGcvPPfcct912G0VFH7/DkUgkSCQS6eWenp6h5hw28XjcnDlUyDmPHjvG8fB2yciRHO3rK9isUNhzeSJz5ta4cePydt9ZCyEWi9Hb25te7u3tpaKiImPM9u3b+dnPfgbAgQMH2LRpE0VFRXzuc5/LcVxJUr5kLYSamhq6urro7u6msrKStWvXMmfOnIwxzzzzTMbtq666yjKQpPNM1kIoLi6mqamJlpYWUqkUDQ0NTJgwgVWrVgHQ2NiY95CSpPzLWggAtbW11NbWZqw7VRHcfffdZ59KknTOeaayJAmwECRJIQtBkgRYCJKkkIUgSQIsBElSyEKQJAEWgiQpZCFIkgALQZIUshAkSYCFIEkKWQiSJMBCkCSFLARJEmAhSJJCFoIkCbAQJEkhC0GSBFgIkqSQhSBJAiwESVLIQpAkARaCJClkIUiSAAtBkhSyECRJgIUgSQpZCJIkwEKQJIUsBEkSYCFIkkIWgiQJsBAkSaHoUAZ1dHSwfPlyUqkUM2fOZNasWRnb//rXv/LSSy8BUFpayp133smll16a66ySpDzKuoeQSqVYtmwZCxYsYPHixaxZs4bdu3dnjBk7diyPPfYYTzzxBDfddBO/+tWv8hZYkpQfWQuhs7OT6upqqqqqiEaj1NfX097enjHmU5/6FGVlZQBMmjSJ3t7e/KSVJOVN1kNGyWSSWCyWXo7FYmzbtu2U41955RWmTp066La2tjba2toAaG1tJR6Pn27ecy4ajZozhwo5Z2lJCUVBAECkqIjSkpKCzQqFPZcnMuf5I2shBOEb5ESRSGTQsW+//Tavvvoq3//+9wfdnkgkSCQS6eWenp6h5hw28XjcnDlUyDmPHjvG8fB2yciRHO3rK9isUNhzeSJz5ta4cePydt9ZDxnFYrGMQ0C9vb1UVFScNO7dd9/l2Wef5dvf/jZjxozJbUpJUt5lLYSamhq6urro7u6mv7+ftWvXUldXlzGmp6eHJ554gnvuuSev7SVJyp+sh4yKi4tpamqipaWFVCpFQ0MDEyZMYNWqVQA0Njbyhz/8gUOHDrF06dL037S2tuY3uSQpp4Z0HkJtbS21tbUZ6xobG9O377rrLu66667cJpMknVOeqSxJAiwESVLIQpAkARaCJClkIUiSAAtBkhSyECRJgIUgSQpZCJIkwEKQJIUsBEkSYCFIkkIWgiQJsBAkSSELQZIEWAiSpJCFIEkCLARJUshCkCQBFoIkKWQhSJIAC0GSFLIQJEmAhSBJClkIkiTAQpAkhSwESRJgIUiSQhaCJAmwECRJIQtBkgRYCJKkkIUgSQIsBElSyEKQJAEQHcqgjo4Oli9fTiqVYubMmcyaNStjexAELF++nE2bNlFSUkJzczMTJ07MR15JUp5k3UNIpVIsW7aMBQsWsHjxYtasWcPu3bszxmzatIkPPviAp556itmzZ7N06dK8BZYk5UfWQujs7KS6upqqqiqi0Sj19fW0t7dnjNmwYQPTp08nEokwefJkDh8+zN69e/MWWpKUe1kPGSWTSWKxWHo5Fouxbdu2k8bE4/GMMclkkoqKioxxbW1ttLW1AdDa2sq4cePOKvy5Ys7cKtScD91663BHOG2FOpcfZc7zQ9Y9hCAITloXiUROewxAIpGgtbWV1tZW5s+ffzo5h405c8ucuXM+ZARz5lo+c2YthFgsRm9vb3q5t7f3pE/+sViMnp6ejx0jSSpsWQuhpqaGrq4uuru76e/vZ+3atdTV1WWMqaurY/Xq1QRBwNatWxk9erSFIEnnmazfIRQXF9PU1ERLSwupVIqGhgYmTJjAqlWrAGhsbGTq1Kls3LiROXPmMHLkSJqbm7M+cCKROPv054A5c8ucuXM+ZARz5lo+c0aCwb4AkCT9z/FMZUkSYCFIkkJDunTFUKxYsYI333yTaDRKVVUVzc3NXHDBBQCsXLmSV155haKiIu644w6mTJkCwI4dO3jmmWfo6+tj6tSp3HHHHUQiEY4fP87TTz/Njh07GDNmDHPnzmXs2LEAvPbaa7z44osA3HjjjVx33XWnlfNvf/sbv//973nvvff44Q9/SE1NDQDd3d3Mmzcv/TvkSZMmMXv27ILLWWjzeaIXXniBv/zlL1x44YUA3HrrrdTW1uY8cz5lu0zLuXD33XdTWlpKUVERxcXFtLa2cujQIRYvXsy///1vLrroIubNm0dZWRlw+nN7JpYsWcLGjRspLy9n0aJFADnNlKvne7Cchfi67Onp4ZlnnmHfvn1EIhESiQTXX3/98M9pkCMdHR1Bf39/EARBsGLFimDFihVBEATBrl27ggceeCDo6+sL9uzZE9xzzz3BwMBAEARBMH/+/OCf//xnkEqlgpaWlmDjxo1BEATByy+/HDz77LNBEATBG2+8Efz0pz8NgiAIDh48GNx9993BwYMHM26fjl27dgXvvfde8OijjwadnZ3p9Xv27Anuv//+Qf+mkHIW2nye6Pnnnw9eeumlQf9bcpU5nwYGBoJ77rkn+OCDD4Ljx48HDzzwQLBr1668P+5HNTc3B/v3789Yt2LFimDlypVBEATBypUrz+r9dSY2b94cbN++PeM9kstMuXq+B8tZiK/LZDIZbN++PQiCIPjwww+DOXPmBLt27Rr2Oc3ZIaPPfvazFBcXAzB58mSSySQA7e3t1NfXM2LECMaOHUt1dTWdnZ3s3buXI0eOMHnyZCKRCNOnT09fEmPDhg3pT6rXXHMNb7/9NkEQ0NHRwZVXXklZWRllZWVceeWVdHR0nFbO8ePHn9bZiIWWs9DmcyhymTmfhnKZluHS3t7OjBkzAJgxY0Y615nM7Zm44oor0p9U85EpV8/3YDlPZThzVlRUpC8AOmrUKC6++GKSyeSwz2nODhmd6JVXXqG+vh74z2UtJk2alN5WWVlJMpmkuLj4pEti/LdETrxcRnFxMaNHj+bgwYMnXUbjv/eVK93d3XznO99h1KhR3HLLLVx++eWDXrpjOHMW+nz++c9/ZvXq1UycOJFvfOMblJWV5TTzf3f782Eol2k5V1paWgD44he/SCKRYP/+/elzeyoqKjhw4ABwZq+HXMllpnw/34X8uuzu7mbnzp1cdtllwz6np1UIjz/+OPv27Ttp/S233MLVV18NwIsvvkhxcTHXXnstMPhlLT5u/am2nerY52Drh5LzoyoqKliyZAljxoxhx44d/OQnP2HRokUFl3M45vNEH5e5sbGRr371qwA8//zz/Pa3v6W5uTnvmXNlOB5zMI8//jiVlZXs37+fH/zgBx+7R3smc5tvhfZ8F/Lr8ujRoyxatIjbb7+d0aNHn9bjftz6U23LlvW0CuHhhx/+2O2vvfYab775Jo888kj6gT966YtkMkllZeWgl8SorKzM+JtYLMbAwAAffvghZWVlVFZWsmXLloz7uuKKK04752BGjBjBiBEjAJg4cSJVVVV0dXUVXM7hmM8zyTxz5kx+9KMf5TxzPg3lMi3nwn/noLy8nKuvvprOzk7Ky8vZu3cvFRUV7N27N/0p70zmNldymSmfz/cnPvGJ9O1Cel329/ezaNEirr32WqZNmwYM/5zm7DuEjo4OXnrpJb773e9SUlKSXl9XV8fatWs5fvw43d3ddHV1cdlll1FRUcGoUaPYunUrQRCwevXq9CUxrrrqKl577TUA1q1bx2c+8xkikQhTpkzhrbfe4tChQxw6dIi33nor/U372Tpw4ACpVAqAPXv20NXVRVVVVcHlLOT5PPGS5+vXr2fChAk5z5xPQ7lMS74dPXqUI0eOpG///e9/55JLLqGuro7XX38dgNdffz29B3kmc5srucyUz+e7EF+XQRDwy1/+kosvvpgvf/nL6fXDPac5O1P53nvvpb+/P91AJ/5s88UXX+TVV1+lqKiI22+/nalTpwKwfft2lixZQl9fH1OmTKGpqYlIJEJfXx9PP/00O3fupKysjLlz51JVVQX85/uJlStXAv/5mWRDQ8Np5Vy/fj2//vWvOXDgABdccAGXXnopDz30EOvWreOFF16guLiYoqIivva1r6UntpByFtp8nujnP/85//rXv4hEIlx00UXMnj07/Qk7l5nzaePGjfzmN79JX6blxhtvzPtjnmjPnj088cQTAAwMDPD5z3+eG2+8kYMHD7J48WJ6enqIx+Pcf//96ffa6c7tmXjyySfZsmULBw8epLy8nJtvvpmrr746Z5ly9XwPlnPz5s0F97p85513eOSRR7jkkkvSz8mtt97KpEmThnVOvXSFJAnwTGVJUshCkCQBFoIkKWQhSJIAC0GSFLIQJEmAhSBJCv0/xGm6+ngptAUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 1)\n",
    "rect1 = matplotlib.patches.Rectangle((np.mean(y1)-np.std(y1),0), 2*np.std(y1), 0.7, color='black', alpha=0.25)\n",
    "rect2 = matplotlib.patches.Rectangle((np.mean(y2)-np.std(y2),0), 2*np.std(y2), 0.7, color='red', alpha=0.3)\n",
    "rect3 = matplotlib.patches.Rectangle((np.mean(y3)-np.std(y3),0), 2*np.std(y3), 0.7, color='cyan', alpha=0.35)\n",
    "ax.add_patch(rect1)\n",
    "ax.add_patch(rect2)\n",
    "ax.add_patch(rect3)\n",
    "plt.xlim(-20000, 20000)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 121,
   "id": "d57d673f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0b0b83dcb2fa46408e15952936fc61f6",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/10000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEaCAYAAAAL7cBuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAo5UlEQVR4nO3de1RU5f4/8PcMd0QIGBARFBC8pWZHUDGDUrQsLaTS47GT19Kvtk55Ocu7mVoHS8VMLcNL1vFYmYLVz0shCil1RDnK8R4aikFyBw25DHx+f/h1vo5scIMMA/J+rcVazJ699/N5Zjb7zd7Pnj0aEREQERHdRWvuAoiIqGliQBARkSIGBBERKWJAEBGRIgYEEREpYkAQEZEiBkQT5+Pjg2XLlpm7jGZFr9djwoQJcHV1hUajwaFDh8xdUrOWnp4OjUaDw4cPN0p7Go0G//znP+97PYsXL4a/v38DVNRyMSDMYNy4cQgLC1N87u4/juTkZEyfPl3Veg8fPgyNRoP09PSGKLPZ2rlzJ/71r3/h22+/RVZWFvr3719tnqtXryqGx7Jly+Dj49PgNf3zn/+ERqOp0zIrVqyAhYUFZsyYUeM8JSUlsLe3R2pqqtH0YcOGwcLCAt9880296r2Tt7c3srKy0Ldv3/telynUtN3PmjULP//8s3mKekAwIJo4Nzc3tGrVytxlKCovLzd3CYp++eUXtGvXDv3794eHhwesra3NXVK9REdHY968efjss89QVlamOM++ffvQtm1b9OzZ0zAtIyMD8fHxmDVrFj755JP7rsPCwgIeHh6wsrK673U1JgcHB+h0OnOX0bwJNbqxY8fKoEGDFJ8DIJ9//rnhcYcOHWTp0qWGx7GxsdKrVy+xs7MTJycnCQoKkpSUFPn1118FgNFPaGioiIhUVVXJ+++/L76+vmJlZSV+fn4SFRVl1G5ubq68+OKLYm9vL+7u7rJgwQJ55ZVXjOoMDQ2VCRMmyIIFC8TDw0N0Op2IiGzbtk369Okjjo6O4urqKs8884ycP3/esNzt2rZt2yZDhgwROzs76dy5sxw6dEiuXr0qQ4cOFXt7e+nataskJibW+trdqy+hoaFGr0GHDh1qfJ3vnm/Lli3Vpr/11lsiIlJRUSFvvfWW+Pj4iI2NjXTr1k0+/vhjo3VGR0dLly5dxMbGRlxcXOTxxx+XjIwMOXjwYLX1jh07ttZ+xsfHi7u7u1RUVEi3bt1k27ZtivONGTNGZs6caTRt0aJFEhERIZmZmWJtbS1XrlyptS2Rmrcrkf97/3788Uejx3V9P2+/DhkZGUZtW1hYyJYtWwyP7/4bWL16tTzyyCPSqlUradOmjYwaNUoyMzONalHa7t966y3p2LGjUVuffvqpdO3aVaytraVdu3Yyf/58qaioMDwfGhoqEydOlCVLlkibNm3E2dlZxo4dKzdu3DDMc+rUKRkyZIg4OTmJvb29dOnSRT777LN7vsbNEQPCDOobEFlZWWJlZSXLly+XS5cuyZkzZ2Tbtm2Smpoqer1edu/eLQDk6NGjkpWVJXl5eSIisnbtWrG1tZUNGzbIhQsX5KOPPhIbGxvZuHGjoZ3hw4dLQECAxMfHy6lTp2TcuHHi6OhYLSAcHBxk8uTJcvr0aUlNTRURkc2bN8u3334raWlpkpKSIsOHDxd/f38pKysTkf/7I/bz85OYmBg5f/68hIeHS9u2bWXQoEGya9cuOX/+vERERIiXl5eUl5fX+Nrdqy95eXkyc+ZM8fHxkaysLMnOzlZcT0pKigCQnTt3GuYrKSmR2bNni5eXl2RlZUlWVpZcv37d8J716NFD9u/fL5cuXZIvvvhCnJycDO0eO3ZMLCwsZOvWrZKeni6pqakSHR0tGRkZUlZWJmvXrhUAhvUWFhbWsoWIjBo1SqZPny4iIsuXLzfs9O5UXl4uDz30kBw+fNgwTa/XS7t27WT37t0iIjJ06FBDyNWktu1KpOaAqOv7eT8B8cMPP8ilS5ckKSlJgoODJSQkxNDfmrb7uwPiu+++E61WK++++66cP39evvjiC3nooYdkwYIFhnlCQ0PFyclJ3nzzTTl79qzs3btXnJycZNGiRYZ5evToIaNHj5bTp0/LxYsXZc+ePfLtt9/W+ho3VwwIMxg7dqxYWFhIq1atqv3UFhC3d2q//vqr4np//PFHxee9vLzk73//u9G0N998U3x9fUVE5MKFCwJA4uLiDM+Xl5eLl5dXtYAICAiQysrKWvuXl5cnAAw7rts7lDv/0z969KgAkBUrVhim3e7ff//73xrXfa++iCj/53i3jIwMASAHDx40mr506dJqRx2XLl0SjUYjZ8+eNZr+9ttvyyOPPCIiIrt27RJHR0cpKipSbO/zzz8XtQfs2dnZYm1tLSdOnBARkczMTLG0tJRz584Zzbd//35p06aN0fsRGxsrbm5uhp3yl19+KV5eXqLX62ts717bVU0BUdf3s74BUVO9V69eFZGat/u7t4MBAwbISy+9ZDTP6tWrxdbW1vDPTGhoqPTo0cNonsmTJ0u/fv0Mjx0dHY3qfZBxDMJM+vbtixMnTlT7qU3Pnj3x1FNPoXv37hgxYgQ++OADZGRk1LpMcXExrl69ipCQEKPpoaGhSE9PR0lJCc6cOQMA6Nevn+F5KysrBAYGVltf7969odUabzYnTpzAiBEj4Ovri9atW6N9+/YAgMuXLxvN98gjjxh+9/DwMPTp7mnZ2dn17ospHDt2DCKCwMBAODg4GH7effdd/PLLLwCAwYMHw8/PD76+vvjzn/+MTz75BLm5ufVqb8uWLejatavh9Wrbti0GDx6M6Ohoo/l27dqF559/3uj92LBhA/7yl78Yxguef/55/PHHH9i7d2+N7dVnuwLu//1U69ChQ3jqqafg7e2N1q1bY8CAAQCqb1/3cvr0acVtp7S0FBcvXjRM69Wrl9E87dq1w7Vr1wyPZ82ahUmTJuGJJ57A4sWLkZKSUsceNR8MCDOxs7ODv79/tZ/aWFhYYO/evYiPj0dQUBB27tyJTp064bvvvrtne3dfQSMKN/FVc5XN3QPmJSUlGDJkCDQaDTZv3oyjR48iOTkZGo2m2iD2nYOct9tSmlZVVXXffWlIt+tJSkoyCvNTp04Zrh5ycHDAsWPHEBMTg06dOuHjjz+Gv78/jh8/Xqe2RAQbN25EamoqLC0tDT/79u3D1q1bDa+piGD37t0YMWKEYdkrV65g//79+PDDDw3LtWrVCgUFBbUOVtd3u6rr+3k7yO58vyorK2t9v69cuYJnnnkGPj4++OKLL3Ds2DHDlVn1uUiipm3nzul3X9Sg0WiMaly4cCEuXLiAkSNH4tSpU+jXrx8WLFhQ51qaAwZEM6PRaNCnTx/MmzcPiYmJCA0NxZYtWwD834ZdWVlpmN/R0RFeXl5ISEgwWk9iYiJ8fX1hb2+Pbt26AQB++uknw/N6vV7Vzu3s2bPIycnBO++8gyeffBJdu3ZFQUGBSXbaavqiltJrdXv63dN69+4N4NbO6u5A79ixo2E+CwsLhISEYMmSJTh+/Djatm2Lf/3rX7W2d7f4+HhcunQJR44cqXZ0WVFRgV27dgG4FVYlJSUYOHCgYdno6Gh07doVJ0+eNFpux44d2LNnD3777bca261tu2oo7u7uAIDMzEzDtBMnTtS6rSQnJ+PmzZtYvXo1HnvsMXTu3Nnov3lA/Wv78MMPK247dnZ28PPzq1Nf/Pz8MHXqVHz99ddYsmQJPvroozot31xYmrsAUi8pKQkHDhzAkCFD0LZtW/zyyy9ITU3FxIkTAQAdOnSAVqvFnj17MGrUKNjY2MDJyQlz587FzJkzERAQgCeeeALx8fH46KOPsG7dOgBAQEAAhg8fjmnTpmHDhg1wc3PDypUrUVxcfM+jig4dOsDGxgYffvghZs6cifT0dMyZM6fO1/yrda++qKXT6eDg4IDvv/8eDz/8MGxsbODs7AxfX1/8/vvv+OmnnxAQEAB7e3v4+/tjwoQJePXVV/Hee+8hODgYf/zxB44fP46cnBzMnj0bu3fvxqVLlxASEgI3NzccP34cGRkZhvD19fUFAHzzzTcYMGAA7Ozs4ODgUK2uDRs2IDQ0FMHBwdWee+655/DJJ5/gz3/+M2JiYvDss88ado56vR6bN2/GlClT0L17d6PlunfvDi8vL2zatAmLFi2qtt57bVcNxd/fHx06dMDixYsRFRWF3NxczJs3r9ZtJSAgABqNBitXrsSYMWNw8uRJLFmyxGiemrb7u82dOxfDhw9HZGQkIiIicOLECSxevBgzZ85UfSn0jRs3MHv2bLzwwgvw9fVFYWEh9u3bZ3ifHzhmG/1owep7FdOpU6dk6NCh0qZNG7G2tpb27dvLrFmzDANsIreuePH09BStVmt0met7770nPj4+YmlpKb6+voqXub7wwgtiZ2cnbm5usnDhQnnxxRdl2LBhhnluXwJ4tx07doi/v7/Y2NhIr1695NChQ0YDj3cPcoooDxJnZWUJAPnhhx9qfO3U9EXNILWIyNatWw3ruT0wXV5eLqNHjxZnZ2ejy1z1er0sX75cOnfuLFZWVuLq6iohISHy1VdfiYhIQkKCPPnkk6LT6cTGxkb8/f3lH//4h1RVVRnae+ONN8Td3V00Go3iZa7Xrl0TKyurapfP3vbdd9+JRqORCxcuiJ+fn6FtkVuD5ACqDWTfNmvWLGnfvr3iBQb32q5qGqSuz/v5888/y5/+9CextbWVnj17SmJi4j0HqdeuXSteXl5ia2srjz32mOzdu7daW0rbfU2XuXbp0kWsrKzE09NT5s2bp3iZ653uvHDh5s2bMnr0aMPlzm5ubjJy5EhVlxI3RxoRfqMcVVdZWYkuXbrgueeew8qVK81dDt3h5MmT6Nu3L3JzcxWPQogaCk8xEYBb52Kzs7Px6KOP4vr164iKikJ6ejrGjRtn7tLoLmVlZVi7di3DgUyOAUEAbh0xLFu2DGlpabCyskL37t1x8OBB9OjRw9yl0V369OmDPn36mLsMagF4iomIiBTxMlciIlLEgCAiIkUP1BjEnR/AqQudTlfv2yI0V+xzy8A+twz302dPT88an+MRBBERKWJAEBGRIgYEEREpYkAQEZEiBgQRESliQBARkSIGBBERKWJAEBGRIgYEEREpeqA+SU3U3FS++pzidIvobxq5EqLqeARBRESKGBBERKSIAUFERIoYEEREpIgBQUREihgQRESkiAFBRESKGBBERKSIAUFERIoYEEREpIgBQUREihgQRESkiDfrI5PhjeiImjcGxAOupp00YpIat5BG1BL7TGQKDAgza4r/ZTfFmoio8XEMgoiIFPEIgu5bjad0iKhZY0A0Eu5Eiai5YUBQk9HUxj5qC/WaampqfSC6HwwIavF4dEekjAFBVA8MFbpbfY44mzoGBNEDoDl99oOn4ZoPXuZKRESKeARBLca1Ef3NXQJRs8KAINV43p2oZWFAEDVBPE9PTQEDgpo87iyJzIOD1EREpKhRjiByc3Oxbt06FBYWQqPRICwsDM888wxu3LiBqKgo5OTkwM3NDdOnT4eDgwMAICYmBvHx8dBqtRg/fjx69erVGKUSUTNh6kt7G/LItbkeBTdKQFhYWOCvf/0r/Pz8cPPmTcyZMwc9e/bEoUOH0KNHD4SHhyM2NhaxsbF4+eWXcfXqVSQlJWHVqlUoKCjA0qVL8cEHH0Cr5QEPEVFjaZQ9rrOzM/z8/AAAdnZ2aNeuHfLz85GcnIzQ0FAAQGhoKJKTkwEAycnJ6N+/P6ysrODu7g4PDw+kpaU1RqlERPS/Gn2QOjs7G7/++iv8/f1RVFQEZ2dnALdCpLi4GACQn5+PgIAAwzIuLi7Iz8+vtq64uDjExcUBACIjI6HT6epVk6WlZb2XVetaHedvqHpqare2Pte11rpqqHbNVX99mKvPjbFt11VNtda1TlP3ua511me7a6j3xlTvc6MGRGlpKVauXIlx48bB3t6+xvlERNX6wsLCEBYWZnicm5tbr7p0Ol29lzUVU9ej1+vN1ueGarepvWe1MVeffxveR3F6Uzz33VCvkam37YZcd0Ot6372YZ6enjU+12gBodfrsXLlSjz++OPo27cvAMDJyQkFBQVwdnZGQUEBHB0dAQCurq7Iy8szLJufnw8XF5fGKrVJaK6DWtS8cbujOzVKQIgIPv74Y7Rr1w7Dhg0zTA8MDERCQgLCw8ORkJCAoKAgw/Q1a9Zg2LBhKCgoQFZWFvz9/Ruj1BaDt52ghsJQeXA1SkCcP38eiYmJaN++Pf7+978DAEaPHo3w8HBERUUhPj4eOp0OM2bMAAB4e3sjODgYM2bMgFarxcSJE3kFE1ED4m1TSI1GCYguXbrgq6++Unxu0aJFitMjIiIQERFhyrKIqBloqDB7EI50Gvu27rzVRgPjf2ZE9KBgQBARmUlTP6rhiX0iIlLEI4h64qkkMgdud9SYGBAPiJa442iJfSbza0nbHQOimWlJGydRfTwIn/FpKn/nHIMgIiJFDAgiIlLEgCAiIkUcgyAik2gq59Gp/hgQRI2AO0tqjniKiYiIFPEI4h74n1/D42tKSrhdND0MCCK6J+68WyaeYiIiIkUMCCIiUsSAICIiRRyDIKIWjeMrNeMRBBERKWJAEBGRIgYEEREpYkAQEZEiBgQRESliQBARkSIGBBERKWJAEBGRIn5Q7n/xwzJERMZ4BEFERIoYEEREpIgBQUREihgQRESkiAFBRESKGBBERKRIVUCsWLECR48ehV6vN3U9RETURKgKiM6dO2Pnzp147bXXEB0djfPnz5u6LiIiMjONiIjamTMyMvDjjz/i8OHDsLCwQGhoKAYMGAAPD49al1u/fj1SUlLg5OSElStXAgC++uorHDhwAI6OjgCA0aNH409/+hMAICYmBvHx8dBqtRg/fjx69eqlqr7MzEy1XTGi0+lwbUT/ei1LRGRubWKSkJubW69lPT09a3yuTp+k9vb2xl/+8hc8+uij2Lx5M3bs2IFvv/0W/v7++Otf/wofHx/F5Z544gk8/fTTWLdundH0Z599Fs89Z/wJ5qtXryIpKQmrVq1CQUEBli5dig8++ABaLYdLiIgak+qAyMzMRGJiIo4cOQJLS0s8/vjjmD17NhwdHfH999/j/fffrxYAt3Xr1g3Z2dmq2klOTkb//v1hZWUFd3d3eHh4IC0tDZ06dVJbKhERNQBVATFnzhzk5OQgODgYf/vb3xAQEGD0/LBhw7B37946N75//34kJibCz88Pr7zyChwcHJCfn2+0fhcXF+Tn5ysuHxcXh7i4OABAZGQkdDpdnWsAAEtL3pKKiJovS0vLeu//al2vmpnCw8MRGBhY6460pqOHmgwZMgQvvvgiAODLL7/EZ599hqlTp6IOQyIICwtDWFiY4XF9z8GZ4oUlImoser3eJGMQqk7s29nZVTtFlJmZidTU1HoVBAAPPfQQtFottFotBg0ahIsXLwIAXF1dkZeXZ5gvPz8fLi4u9W6HiIjqR1VAbNq0CXZ2dkbTbG1tsWnTpno3XFBQYPj96NGj8Pb2BgAEBgYiKSkJFRUVyM7ORlZWFvz9/evdDhER1Y+qU0xFRUVwdnY2mubs7IzCwkJVjaxevRpnzpzB9evXMWXKFIwcORKnT59Geno6NBoN3Nzc8NprrwG4daVUcHAwZsyYAa1Wi4kTJ/IKJiIiM1AVEG3atMGpU6fQvXt3w7TTp0/D3d1dVSNvvvlmtWkDBw6scf6IiAhERESoWjcREZmGqoB46aWXsGLFCgwcOBBt2rTBtWvXcPDgQUydOtXU9RERkZmoOncTFBSEBQsWoLS0FCkpKSgtLcX8+fMRFBRk6vqIiMhMVH8AwN/fn4PFREQtiKqA0Ov1OHToENLT01FaWmr03Ouvv26SwoiIyLxUBcTatWtx+fJl9O7dG05OTqauiYiImgBVAXHy5EmsXbsWrVq1MnU9RETURKgapNbpdKioqDB1LURE1ISoOoIICQnB+++/j6FDh+Khhx4yeu7Oz0YQEdGDQ1VA7Nu3DwCwfft2o+kajQZr165t+KqIiMjsVAVEXe/USkREzZ/qmxzp9XqcPXsWSUlJAIDS0tJql7wSEdGDQ9URxJUrV7B8+XJYWVkhLy8P/fv3x5kzZ5CQkIDp06ebukYiIjIDVUcQ0dHRGDVqFFavXm340qBu3brh3LlzJi2OiIjMR1VAXL16FY8//rjRNFtbW5SXl5ukKCIiMj9VAeHm5oZLly4ZTUtLS4OHh4dJiiIiIvNTNQYxatQoREZGYvDgwdDr9YiJicEPP/yAyZMnm7o+IiIyE1VHEL1798bcuXNRXFyMbt26IScnB7NmzcIjjzxi6vqIiMhMVN/u28/PD35+fqashYiImhBVAfHll1/W+NyoUaMarBgiImo6VAVEXl6e0ePCwkKcOXMGffr0MUlRRERkfqoCQum7p0+cOIHDhw83eEFERNQ0qL7Vxt169uyJ5OTkhqyFiIiaEFVHENeuXTN6XFZWhsOHD0On05mkKCIiMj9VAfG3v/3N6LG1tTV8fX0xbdo0kxRFRETmd99XMRER0YOp3mMQRET0YFN1BPE///M/qlb20Ucf3VcxRETUdKgKiKFDhyIhIQFDhw6Fm5sbcnJysH//foSEhKBjx46mrpGIiMxAVUAkJCRg/vz5cHFxMUx79NFH8e6772L48OEmK46IiMxH1RhEfn4+bG1tjabZ2toiPz/fJEUREZH5qTqCCAwMxPLly/HCCy/AxcUFeXl5iI2NRe/evU1dHxERmYmqgHj11VexY8cOREdHIz8/Hy4uLujXrx9eeuklU9dHRERmoiogrK2tMWbMGIwZM8bU9RARUROh+vsgUlNTceTIERQVFWHOnDm4ePEibt68ie7du5uyPiIiMhNVAbF3717s2bMHgwYNws8//wzg1lHFli1bsGzZsnsuv379eqSkpMDJyQkrV64EANy4cQNRUVHIycmBm5sbpk+fDgcHBwBATEwM4uPjodVqMX78ePTq1aue3SMiovpSdRXTnj17sHDhQoSHh0OrvbVIu3btkJmZqaqRJ554AvPmzTOaFhsbix49emDNmjXo0aMHYmNjAQBXr15FUlISVq1ahfnz52PTpk2oqqqqQ5eIiKghqAqImzdvVrtzq16vh6WlujNU3bp1Mxwd3JacnIzQ0FAAQGhoqOHW4cnJyejfvz+srKzg7u4ODw8PpKWlqWqHiIgajqo9fNeuXREbG4uIiAjDtL179+Lhhx+ud8NFRUVwdnYGADg7O6O4uBjArc9cBAQEGOZzcXGp8fMWcXFxiIuLAwBERkbW+/bjaoOOiKgpsrS0NMnXL6jaM06YMAHLly/HgQMHUFpaijfeeAP29vaYPXt2gxckIqrnDQsLQ1hYmOFxbm5uvdrk91oQUXOm1+vrvf/z9PSs8bl7BkRVVRV+++03LFmyBFeuXEFOTg5cXV3h7+9vGI+oDycnJxQUFMDZ2RkFBQVwdHQEALi6uhp9B/btz10QEVHjuuceXqvV4r333oO1tTX8/f0RHByMTp063Vc4ALc+nZ2QkADg1r2egoKCDNOTkpJQUVGB7OxsZGVlwd/f/77aIiKiulM9BnHhwgV06tSpXo2sXr0aZ86cwfXr1zFlyhSMHDkS4eHhiIqKQnx8PHQ6HWbMmAEA8Pb2RnBwMGbMmAGtVouJEyfedxgREVHdaUTFSf+NGzfiyJEjCAwMhKurKzQajeG5UaNGmbTAulB72e3ddDodro3o38DVEBE1jjYxSSYZg6jxX/N9+/YZfi8pKUFQUBA0Gg3y8/ORl5dn+CEiogdTjaeYtm/fjqeffhoAcPz4cWzdurXRiiIiIvOrMSA8PDzw2WefwcvLC3q9HgcPHlS8BHXgwIEmLZCIiMyjxoB444038M033+DIkSOorKxEYmKi4nwMCCKiB1ONAeHp6YkpU6YAAJYsWYJFixY1WlFERGR+qq4fZTgQEbU8/IABEREpYkAQEZEiBgQRESliQBARkSIGBBERKWJAEBGRIgYEEREpYkAQEZEiBgQRESliQBARkSIGBBERKWJAEBGRIgYEEREpYkAQEZEiBgQRESliQBARkSIGBBERKWJAEBGRIgYEEREpYkAQEZEiBgQRESliQBARkSIGBBERKWJAEBGRIgYEEREpYkAQEZEiBgQRESliQBARkSIGBBERKbI0dwHTpk2Dra0ttFotLCwsEBkZiRs3biAqKgo5OTlwc3PD9OnT4eDgYO5SiYhaFLMHBAC89dZbcHR0NDyOjY1Fjx49EB4ejtjYWMTGxuLll182Y4VERC1PkzzFlJycjNDQUABAaGgokpOTzVwREVHL0ySOIN555x0AwODBgxEWFoaioiI4OzsDAJydnVFcXKy4XFxcHOLi4gAAkZGR0Ol09Wrf0rJJvAxERPViaWlZ7/1frett8DXW0dKlS+Hi4oKioiIsW7YMnp6eqpcNCwtDWFiY4XFubm69ajDFC0tE1Fj0en2993+17XPNforJxcUFAODk5ISgoCCkpaXByckJBQUFAICCggKj8QkiImocZg2I0tJS3Lx50/B7amoq2rdvj8DAQCQkJAAAEhISEBQUZM4yiYhaJLOeYioqKsKKFSsAAJWVlRgwYAB69eqFjh07IioqCvHx8dDpdJgxY4Y5yyQiapE0IiLmLqKhZGZm1ms5nU6HayP6N3A1RESNo01M0oM5BkFERE0TA4KIiBQxIIiISBEDgoiIFDEgiIhIEQOCiIgUMSCIiEgRA4KIiBQxIIiISBEDgoiIFDEgiIhIEQOCiIgUMSCIiEgRA4KIiBQxIIiISBEDgoiIFDEgiIhIEQOCiIgUMSCIiEgRA4KIiBQxIIiISBEDgoiIFDEgiIhIEQOCiIgUMSCIiEgRA4KIiBQxIIiISBEDgoiIFDEgiIhIEQOCiIgUMSCIiEgRA4KIiBQxIIiISBEDgoiIFDEgiIhIkaW5C6jNiRMnsGXLFlRVVWHQoEEIDw83d0lERC1Gkz2CqKqqwqZNmzBv3jxERUXhyJEjuHr1qrnLIiJqMZpsQKSlpcHDwwNt2rSBpaUl+vfvj+TkZHOXRUTUYjTZU0z5+flwdXU1PHZ1dcUvv/xiNE9cXBzi4uIAAJGRkfD09Kx3e97/71i9lyUiMrf72f/VpMkeQYhItWkajcbocVhYGCIjIxEZGXlfbc2ZM+e+lm+O2OeWgX1uGUzV5yYbEK6ursjLyzM8zsvLg7OzsxkrIiJqWZpsQHTs2BFZWVnIzs6GXq9HUlISAgMDzV0WEVGL0WTHICwsLDBhwgS88847qKqqwpNPPglvb2+TtBUWFmaS9TZl7HPLwD63DKbqs0aUTvYTEVGL12RPMRERkXkxIIiISFGTHYMwhXvdukNEsGXLFvznP/+BjY0Npk6dCj8/P/MU20Du1ecff/wRu3fvBgDY2tpi0qRJ8PHxafxCG5DaW7SkpaVh/vz5mD59Ovr169e4RTYwNX0+ffo0Pv30U1RWVqJ169Z4++23G7/QBnSvPpeUlGDNmjXIy8tDZWUlhg8fjieffNI8xTaA9evXIyUlBU5OTli5cmW1502y/5IWorKyUl5//XX5/fffpaKiQmbNmiUZGRlG8xw/flzeeecdqaqqkvPnz8vcuXPNVG3DUNPnc+fOyfXr10VEJCUlpUX0+fZ8ixcvlnfffVd++uknM1TacNT0+caNG/Lmm29KTk6OiIgUFhaao9QGo6bPO3fulM8//1xERIqKimTcuHFSUVFhjnIbxOnTp+XixYsyY8YMxedNsf9qMaeY1Ny649ixYwgJCYFGo0GnTp3wxx9/oKCgwEwV3z81fe7cuTMcHBwAAAEBAUafPWmO1N6iZe/evejbty8cHR3NUGXDUtPnw4cPo2/fvtDpdAAAJycnc5TaYNT0WaPRoLS0FCKC0tJSODg4QKttvru8bt26Gf5WlZhi/9V8X606Urp1R35+frV5bv8B1TRPc6Kmz3eKj4/Ho48+2hilmYza9/no0aMYMmRIY5dnEmr6nJWVhRs3bmDx4sWYPXs2EhISGrvMBqWmz08//TR+++03TJ48GTNnzsT48eObdUDciyn2Xy1mDEJU3LpDzTzNSV36c+rUKRw8eBBLliwxdVkmpabPn376KcaMGfPA7CzU9LmyshK//vorFi5ciPLycixYsAABAQEmuX9PY1DT55MnT6JDhw5YtGgRrl27hqVLl6JLly6wt7dvrDIblSn2Xy0mINTcusPV1RW5ubm1ztOcqL1dyeXLl7FhwwbMnTsXrVu3bswSG5yaPl+8eBEffPABAKC4uBj/+c9/oNVq0adPn0attaGo3bZbt24NW1tb2NraomvXrrh8+XKzDQg1fT548CDCw8Oh0Wjg4eEBd3d3ZGZmwt/fv7HLbRSm2H89GP9CqaDm1h2BgYFITEyEiODChQuwt7dv1gGhps+5ublYsWIFXn/99Wa7s7iTmj6vW7fO8NOvXz9MmjSp2YYDoH7bPnfuHCorK1FWVoa0tDS0a9fOTBXfPzV91ul0+O9//wsAKCwsRGZmJtzd3c1RbqMwxf6rRX2SOiUlBVu3bjXcuiMiIgLff/89AGDIkCEQEWzatAknT56EtbU1pk6dio4dO5q56vtzrz5//PHH+Pe//204d2lhYXHfd8c1t3v1+U7r1q1D7969m/1lrmr6/M033+DgwYPQarUYOHAgnn32WXOWfN/u1ef8/HysX7/eMFD7/PPPIyQkxJwl35fVq1fjzJkzuH79OpycnDBy5Ejo9XoAptt/taiAICIi9VrMKSYiIqobBgQRESliQBARkSIGBBERKWJAEBGRIgYEkQmtW7cOX3zxhbnLIKoXBgQRESliQBARkaIWcy8monuZNm0awsLCkJiYiMLCQgQFBWHSpEmYPXs2Xn75ZfTu3RvArRvfvfbaa5g/fz78/PywatUqnD17FuXl5fDx8cGkSZPg7e1dbf2HDh3CgQMHsHTpUsO0kSNHYs2aNfDw8EBFRQW2b9+On376CXq9HkFBQRg3bhysra1RXFyM9evX49y5c9BoNPD29sbixYsfmBsOUtPErYvoDocPH8b8+fPx4YcfIisrC7t27cJjjz2GI0eOGOY5efIkWrdubfi2rl69emHNmjXYuHEjfH19sWbNmnq1vW3bNmRlZeH999/HmjVrkJ+fj6+//hoA8N1338HFxQUbN25EdHQ0Ro8e3azvNEzNAwOC6A5PPfUUdDodHBwcMGLECBw5cgQDBgzAsWPHUFZWBuBWiAwYMMCwzMCBA2FnZwcrKyu89NJLuHz5MkpKSurUrojgwIEDGDt2LBwcHGBnZ4eIiAhDMFlYWKCwsBC5ubmwtLRE165dGRBkcjzFRHSHO79wxc3NDfn5+fDw8EC7du1w/Phx9O7dG8eOHcN7770HAKiqqsL27dvx888/o7i42LDTLi4urtP3DhQXF6OsrAxz5swxTBMRVFVVAQCee+457NixA8uWLQMAhIWF1fhd20QNhQFBdIc776efm5sLFxcXAMBjjz2Gw4cPo6qqCl5eXvDw8ABw62ji2LFjWLhwIdzc3FBSUoLx48crrtvGxgbl5eWGx4WFhYbfW7duDWtra6xatcrQ5p3s7Ozwyiuv4JVXXkFGRgbefvttdOzYET169GiIbhMp4ikmojvs378feXl5uHHjBmJiYhAcHAzgVkCkpqbihx9+MDq9dPPmTVhaWsLBwQFlZWXYvn17jevu0KEDMjIykJ6ejvLycnz11VeG57RaLQYNGoRPP/0URUVFAG59heSJEycAAMePH8fvv/8OEYGdnR20Wi0HqMnkeARBdIcBAwZg2bJlKCgoQGBgIF544QUAgLOzMzp16oQzZ85g+vTphvlDQ0Nx8uRJTJkyBQ4ODhg1apThOwnu5unpiRdffBFLly6FtbU1Ro8ejbi4OMPzY8aMwddff4358+fj+vXrcHFxweDBg9GrVy9kZWVh8+bNKC4uRqtWrTBkyBA8/PDDpn0xqMXj90EQ/a9p06Zh8uTJ6Nmzp7lLIWoSeIxKRESKGBBERKSIp5iIiEgRjyCIiEgRA4KIiBQxIIiISBEDgoiIFDEgiIhI0f8HtxTSYpHaCGkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "0.045"
      ]
     },
     "execution_count": 121,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "simulations = 10000\n",
    "n_s = 100\n",
    "res = []\n",
    "from tqdm.auto import tqdm\n",
    "\n",
    "# Запуск симуляций A/A теста\n",
    "for i in tqdm(range(simulations)):\n",
    "  s1 = y1.sample(n_s, replace = False).values\n",
    "  s2 = y3.sample(n_s, replace = False).values\n",
    "  res.append(stats.ttest_ind(s1, s2, equal_var = True)[1]) # сохраняем pvalue\n",
    "\n",
    "plt.hist(res, bins = 50)\n",
    "plt.style.use('ggplot')\n",
    "plt.xlabel('pvalues')\n",
    "plt.ylabel('frequency')\n",
    "plt.title(\"Histogram of ttest A/A simulations \")\n",
    "plt.show()\n",
    "\n",
    "# Проверяем, что количество ложноположительных случаев не превышает альфа\n",
    "sum(np.array(res) < 0.05) / simulations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "id": "03b0306d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Uniform\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.054322475244723556"
      ]
     },
     "execution_count": 122,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.stats import kstest, uniform\n",
    "\n",
    "def is_uniform(sample, significance_level):\n",
    "    # Perform the Kolmogorov-Smirnov test\n",
    "    _, p_value = kstest(sample, uniform.cdf)\n",
    "\n",
    "    return p_value\n",
    "\n",
    "significance_level = 0.05\n",
    "# Example usage\n",
    "p_value = is_uniform(res, significance_level)\n",
    "\n",
    "if p_value > significance_level:\n",
    "    print(\"Uniform\")  # Sample distribution is likely uniform\n",
    "else:\n",
    "    print(\"Not Uniform\")  # Sample distribution is not likely uniform\n",
    "\n",
    "p_value"
   ]
  }
 ],
 "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}