code_languages_recognition.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "5b9ee5f7-3780-4f79-a79f-8c303f868875",
   "metadata": {},
   "source": [
    "# \n",
    "\n",
    "# \n",
    "   \n",
    " \n",
    "\n",
    "\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f786091a-a77c-44cd-8512-e0a1f2b45e3c",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "La modélisation des langues est un domaine où l'apprentissage non supervisé est très appliqué, l'analyse des textes et la detection des langues écrites sont parmi les applications très connues. Dans ce projet, on souhaite pouvoir identifier la langue d’un texte donné (français ou anglais). \n",
    "\n",
    "Pour cela, on va utiliser un tableau de données de textes où on a déjà labélisé chaque texte par sa langue. Dans un premier temps, l’objectif est de construire plusieurs modèles caractérisants les différentes langues, basé sur la fréquence d’apparition des symboles (lettres) dans chaque langue, pour ensuite procéder à un comparatif des diffférents modèles.\n",
    "\n",
    "Chaque exercice du projet consistera à :\n",
    "- Choisir un modèle.\n",
    "- Estimer ses paramètres.\n",
    "- Le programmer et commenter les résultats."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18d3d71b",
   "metadata": {},
   "source": [
    "## EXO 1 : Construire les données textes :\n",
    "\n",
    "### Question 1 :\n",
    " Construire le dataset des textes en français et en anglais, avec un label -1 et 1 pour chacune des langues."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 218,
   "id": "97cad0d7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import string\n",
    "import re\n",
    "from unicodedata import normalize\n",
    "from os import listdir\n",
    "from sklearn.model_selection import train_test_split, KFold\n",
    "\n",
    "from random import sample \n",
    "from sklearn.utils import shuffle\n",
    "os.chdir('/Users/mohamedennatiqi/Desktop/main dir/Bureaux/Unsupervised learning/projet/V2/data/')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 219,
   "id": "a7699f4e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def strip_accents(text):\n",
    "    import unicodedata\n",
    "    return ''.join(c for c in unicodedata.normalize('NFD', text)\n",
    "                  if unicodedata.category(c) != 'Mn')\n",
    "\n",
    "def cleaning_text(text):\n",
    "    special_char = {'~', ':', \"'\", '+', '[', '\\\\', '@', '^',\n",
    "                    '{', '%', '(', '-', '\"', '*', '|', ',', \n",
    "                    '&', '<', '`', '}', '.',\"’\", '_', '=', ']', \n",
    "                    '!', '>', ';', '?', '#', '$', ')', '/'}\n",
    "    for char in special_char:\n",
    "        text = text.replace(char, \"\")\n",
    "    return(strip_accents(text.upper()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 220,
   "id": "23d62caa",
   "metadata": {},
   "outputs": [],
   "source": [
    "import glob\n",
    "import os\n",
    "\n",
    "os.chdir(r'/Users/mohamedennatiqi/Desktop/main dir/Bureaux/Unsupervised learning/projet/V2/data')\n",
    "myFiles = glob.glob('*.txt')\n",
    "\n",
    "text_data = pd.DataFrame(columns= [\"text\", \"language\"])\n",
    "texts = []\n",
    "\n",
    "def read_text_file(file_path):\n",
    "    with open(file_path, 'r') as f:\n",
    "        return(f.read())\n",
    "        \n",
    "        \n",
    "for file in myFiles:\n",
    "    if file.endswith(\".txt\"):\n",
    "        file_path = f\"{'/Users/mohamedennatiqi/Desktop/main dir/Bureaux/Unsupervised learning/projet/V2/data'}/{file}\"\n",
    "        texts.append(read_text_file(file_path))\n",
    "        \n",
    "text_data.text, text_data.language = texts, [1 if (\"english\" in file_name) else -1 for file_name in myFiles]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 221,
   "id": "77a0dcb0",
   "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>text</th>\n",
       "      <th>language</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>I saw this DVD in my friends house and thought...</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Inhabituel, n'est-ce rien ? Seules les fréquen...</td>\n",
       "      <td>-1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Je cherche maintenant dans le passé avec une t...</td>\n",
       "      <td>-1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>What was with all the Turkish actors? No offen...</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>I am a big fan of Arnold Vosloo. Finally seein...</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>65</th>\n",
       "      <td>This film had a lot of promise, and the plot w...</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>66</th>\n",
       "      <td>All I could think of while watching this movie...</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>67</th>\n",
       "      <td>Florence venait d’avouer quelque chose qu’elle...</td>\n",
       "      <td>-1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>68</th>\n",
       "      <td>Restent toujours de minutieuses précautions à ...</td>\n",
       "      <td>-1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>69</th>\n",
       "      <td>If you look at Corey Large's information here ...</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>70 rows × 2 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                                                 text  language\n",
       "0   I saw this DVD in my friends house and thought...         1\n",
       "1   Inhabituel, n'est-ce rien ? Seules les fréquen...        -1\n",
       "2   Je cherche maintenant dans le passé avec une t...        -1\n",
       "3   What was with all the Turkish actors? No offen...         1\n",
       "4   I am a big fan of Arnold Vosloo. Finally seein...         1\n",
       "..                                                ...       ...\n",
       "65  This film had a lot of promise, and the plot w...         1\n",
       "66  All I could think of while watching this movie...         1\n",
       "67  Florence venait d’avouer quelque chose qu’elle...        -1\n",
       "68  Restent toujours de minutieuses précautions à ...        -1\n",
       "69  If you look at Corey Large's information here ...         1\n",
       "\n",
       "[70 rows x 2 columns]"
      ]
     },
     "execution_count": 221,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "text_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 222,
   "id": "b324731e-3533-4fb4-8560-9d670638ea61",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-----------------------\n",
      "Number of words\t\t:\t 51288\n",
      "Number of unique words\t:\t 103\n",
      "Nmber of english texts\t:\t35\n",
      "Nmber of french texts\t:\t35\n",
      "-----------------------\n"
     ]
    }
   ],
   "source": [
    "def show_statistics(data):\n",
    "        \n",
    "        texts = \"\"\n",
    "        for text in data.text.values:\n",
    "            texts += text\n",
    "        word_list = ' '.join(texts).split()\n",
    "        number_of_words = len(word_list)\n",
    "        number_of_unique_words = len(set(word_list))\n",
    "        number_en = data[data.language == 1].language.sum()\n",
    "        number_fr = np.abs(data[data.language == -1].language.sum())\n",
    "\n",
    "        print('-----------------------')\n",
    "        print(f'Number of words\\t\\t:\\t {number_of_words}')\n",
    "        print(f'Number of unique words\\t:\\t {number_of_unique_words}')\n",
    "        print(f'Nmber of english texts\\t:\\t{number_en}')\n",
    "        print(f'Nmber of french texts\\t:\\t{number_fr}')\n",
    "        print('-----------------------')\n",
    "\n",
    "        \n",
    "show_statistics(text_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b67e824-97c1-43e9-9ad9-58d51a3bb5cd",
   "metadata": {},
   "source": [
    "\n",
    "### Question 2 :\n",
    " Construire un tableau $X = (x_{ij})_{i= 1...n, j=1...p}$ telque $x_{ij}=log(1+f_{ij})$ où les $f_{ij}$ sont les fréquences du symbole $j$ dans le texte $i$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 223,
   "id": "3a8ab8b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "alphabet = list('abcdefghijklmnopqrstuvwxyz ')\n",
    "N = len(alphabet)\n",
    "\n",
    "def to_sentences(doc):\n",
    "    return doc.strip().split('$')\n",
    " \n",
    "def nb_lettres(text):\n",
    "    freq = [text.count(x) for x in alphabet]\n",
    "    return freq\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 224,
   "id": "5ee95a23",
   "metadata": {},
   "outputs": [],
   "source": [
    "nb_texts = text_data.shape[0]\n",
    "X = np.zeros((nb_texts, N))\n",
    "y = np.zeros(nb_texts)\n",
    "k = 0\n",
    "for text in text_data.text:\n",
    "    text_cleaned = cleaning_text(text)\n",
    "    X[k,:] = nb_lettres(text)\n",
    "    k=k+1\n",
    "    y = np.array(text_data.language.values)\n",
    "    \n",
    "X = np.log(X+1)\n",
    "y = np.array(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 225,
   "id": "2b4a90b7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(70,)"
      ]
     },
     "execution_count": 225,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 226,
   "id": "aa75103f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(70, 27)"
      ]
     },
     "execution_count": 226,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "70e7182b-8b2b-4a45-9871-5d048a578711",
   "metadata": {},
   "source": [
    "\n",
    "### Question 3 :\n",
    "\n",
    "Pour chacune des deux classes (anglais-français)on représente les log-fréquences sous forme d'un histogramme des symboles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 227,
   "id": "8a867000",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2AAAAE/CAYAAAAg1aCvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAmZ0lEQVR4nO3de7yUZbn/8e83QBaJ4VawtpEuT1tBE8wFm1J3poYa5rYiK/HYLrKdPy11F2X1ItO2lhn5U3+mZmWah7A8hJpnEcJ0gYtSwfKwFJIQIRE8FIfr98fzLBqW6zCD67lnzazP+/XixayZe57reg4z81xz3889jggBAAAAAIr3lmonAAAAAAB9BQUYAAAAACRCAQYAAAAAiVCAAQAAAEAiFGAAAAAAkAgFGAAAAAAkQgEGoCbZfrvtmbZX2f5+N23fYvtm259NlV9qthtth+3+1c6lUrZbbR+0ic8N2zv3dE69ne3jbc8qs+3XbF9edE5Fsr2P7T/bXm37iGrn08b2Jba/Ue08ANSWmvugBtD72W6V9JmIuKvAMJMlvSjpbdH9DxqeLenuiLiswHyAXikivtPTy7Q9VdLOEXF0Ty+7E2dKujAifpgoXlki4sRq5wCg9lCAAahV20t6vLPiy3b/iFgrSRHx1aSZAehp20t6rKMHbFuSI2J92pQAYNMwBBFAMrYH2p5m+/n83zTbA0se/7LtJfljn+lseJntn0o6TtKX8yFJB9meanu67atsvyzpeNtDbP84X+ZfbJ9lu1++jH62z7P9ou2nbX+hdAhf+2Fx+fKvKvl7nO3f2X7J9nzb+5c8dp/tb9uenQ+RvMP20JLH9y157iLbx5dsn/NsP2d7aT68aVAn23Kj/CVNaPd4V+u+s+37ba/Mn39dJzEa8u25PM/14Xzo58dtz23X9jTbN7btH9sX274t3z+zbb8j399/s73Q9l7two2x/Xj++E9sN5Qs+7O2n7S9Ih9Kum0n+Xa6/WwPtf2bfD1W2H7AdoefgbZ3s31n3u4J20eWPPZT2xfZnpHv29/b3qmj5eTtf2n7r/m2nml793KXZXt8Hn9lvj3vt/2ZTuL8MD+WXrY91/Z+JY9tOHY726edLHNb2zfYXmb7Gdsn5/cfIulrkj6R79/5lTy/JKfrbV+Zr/tjtps6Wc5TknaUdEseb6Cz19jZtmdLelXSjrZPsL0gX97Ttj9Xsoz9bS/Oj9MXnL0uTih5fJDt79t+Nt/es0qOne724Vn57bKPMQB9G28MAFI6Q9I4SaMljZI0VtLXpQ0ndadKOkjSzpLe39lCIuJ4SVdL+m5EDC4Z6vifkqZL2jJ//GeS1ubL20vSeEltJ7CflXRYfn+TpInlroTtd0qaIeksSVtJOl3SDbaHlTQ7StIJkraRtFneRra3k3SbpP8raVi+LVry55wr6d/y+3aW9E5J3+wkje7y72rdvy3pDkn/Iml4nktHjpM0RNK7JG0t6URJr0m6WdIOtkeUtD1a0s9L/j5S2b4dKunvkuZImpf/PV3S+e1iTZJ0sKSd8m3QdlwcIOl/8+X9q6RnJV3bSb5dbb/TJC1Wts3frqyAeEPvqe3NJd0p6RfK9t2nJF1cetKd3/ctZdvvSWVDXDtzm6Rd8mXNU3ZclupwWc4K9umSvqps2z8h6X1dxHk4X++t8tx/6ZIitkRn+3QjeeFwi6T5yrbjgZK+aPvgiLhd0nckXZe//kZV8vySZocr25dbKjumLuxoxSJiJ0nPSfpwHu/v+UPHKBuKvIWy4+IFZa+Jtyl77f3A9ntKFvWOfN3fKem/JF1k+1/yx86TtLeybbyVpC9LautR624ftinrGAMACjAAKU2SdGZEvBARy5SdeB6TP3akpJ9ExGMR8Wr+WKXmRMSN+VCkt0k6VNIXI+KViHhB0g8kfbIk3rSIWBQRK5Sd5JfraEm3RsStEbE+Iu6U1CzpQyVtfhIRf4qI1yRdr+zkWMq2wV0RcU1ErImI5RHRYtvKiqovRcSKiFil7CT3k+pYp/nnPRpdrfsaZUO6to2I1yOis8kc1ig7Sd85ItZFxNyIeDk/Ab4u3w7Ki5NGSb8pee6v8/avS/q1pNcj4sqIWJc/t30P2IUl63K2ssKkbXtdERHz8rhflfRe242lTy5j+61RVsBtn2/3BzoZvnqYpNaI+ElErI2IeZJu0MYF7q8i4qF8iOvV+ue+fYOIuCIiVuW5T5U0yvaQMpb1IUmPRcSv8scukPTXLuJclR9LayPi+5IGStq1g6Yd7tMO2o2RNCwizoyIf0TE05IuU+fH46Y8f1b+GlqnrHh/QyHXjZ/m7xdr8306IyKeisz9yr5k2K+k/Rpl7z9rIuJWSasl7ZoXi5+WdEpE/CXfLr9rK/TK2Ielyy/nGAPQx1GAAUhpW2XfVLd5Nr+v7bFFJY+V3i5X6XO2lzRA0pJ8SNBLkn6k7FvsjuKV5tWd7SV9vG25+bL3VXby1ab0ZPlVSYPz2++S9FQHyxwm6a2S5pYs8/b8/o50lX936/5lSZb0UD7069OdxPi5pN9KutbZsNDv2h6QP/YzSUflhc8xkq4v6ZmQpKUlt1/r4O/B2lj7dSk9LjasW0SslrRcWS9Gqe623/eU9TDdkQ9Pm9LJOm8v6d/b7dtJynpP2nS2bzfibJjoObafcjYstjV/aGhJs86WtdH+zU/kF3eSc9sQ0AX5MLmXlPX0DO2gaVf7tNT2krZttx2+pqxnpxzlPL/9uje4slk8N3qPsH2o7Qfz4X8vKStiS7fB8rbrQktiDs7bNKiD12WZ+7BNuccYgD6OSTgApPS8Nr6Yfrv8Pklaomw4XJt3bcLyS79tXqRs6NvQdiddbZa0i7Fdu8dfUXZC36b0BHyRpJ9HxKZMa79I2dDL9l5UVpjsHhF/KWM5XeXf5bpHxF+V9RbJ9r6S7rI9MyKebNdujbKeyG/lPU63KhsK9+OIeND2P5T1MByV/3sz2q9L23HRdswoz3dzZT047bdRl9sv7xE7TdJpeY/dvbYfjoi72zVdJOn+iPjgm1mZ3FHKhsUepOzEfYikvykrfruz0eshL3SHd9TQ2fVeX1E2zO+xiFhvu8M4Xe3Tdk0XSXomInbpJL/uena6e35P2JCDs2tJb5B0rKSbImKNs2sSy9nWL0p6Xdnw1/bXs5W9Dys4xgD0cfSAASjKAGcX/Lf96y/pGklftz0sv8blm5LaJra4XtIJtkfYfqs6v/apLBGxRNkQpO/bfpuz3wLbyXbbtWXXSzrZ9vD8OpD231a3SPqk7QHOJgcoHYJ2laQP2z44/4a8wdlF/h2eILdztaSDbB9pu7/trW2PzodNXqbsupVtpOxas3bXzJTqNP/u1t3ZJBptuf5N2YnsuvYBbH/A9rudTd7xsrIhVqXtrlR23c7aLoYxlusL+bpspaynpG1ikF8oOy5G5yfZ35H0+4hoLX1yd9vP9mHOJh9xvi7rOlpnZcMo/832Mfm+H2B7jDe+3q1cWygrhJcrK+YrmQ5+hqR32z4if+18QRt/CdA+zlpJyyT1t/1NZUNw36CMfdrmIUkv2/6Kswkq+tnew/aY/PGlkhrd+SQT3T2/p22mbNjlMklrbR+q7LrHbuXHzhWSznc2cUg/2+/Nj7ey92EFxxiAPo4CDEBRblXWI9H2b6qySSuaJf1B0h+VXdB+liRFxG3KrnO5V9kwnjn5ckqHtVXqWGUnZo8rKzSm65/DBC9TNhRrfp7Hr9o99xvKvhH/m7Ieg1+0PRARi5R9K/41ZSd8iyT9j8p4T42I55QNjTpN0gplhV7btS9fUbbuD+bDne5Sx9fxlJN/V+s+RtLvba9WNvnBKRHxTAcx3pE/72VJCyTdr38WzFI2nG0PbTz5xqb6hbKi8en8X9txcbeyfXGDsl6hndT5dUhdbb9d8r9XKzu2Lo6I+9ovIO/FGJ/HeF7ZMLlzlZ3cV+pKZcMn/6JsPzxY7hMj4kVJH5f0XWUn/yOVvXY6ej38VtlEEX/K472uzofwdrdP2+Kvk/RhZdekPaOsl+hyZT1AkvTL/P/ltudtwvN7VL7fTlb2xcTflPVc3VzBIk5X9p7UoqyYPVfZ67mSfVjWMQYA5vpQAL1R3uPwqKSBnQwh7Ol4jcpOFAekiFcPnE3T/YKk90TEn6udTz3Le5oWS5oUEfdWO596lfde3SHpkLyIBIAeRw8YgF7D9kdsb5YPqTtX0i0UQ73a5yU9TPFVjHyI65b5ULivKbvuqOxeNFQm/0KhX/5vhyqnA6COUYAB6E0+p2xI31PKrp34fHXTQWdst0o6RdlQShTjvcpeCy8qG853RGQ/a4BijJC0Utl1X5syCysAlIUhiAAAAACQCD1gAAAAAJAIBRgAAAAAJFLIDzEPHTo0Ghsbi1g0AAAAAPR6c+fOfTEihrW/v5ACrLGxUc3NzUUsGgAAAAB6PdvPdnQ/QxABAAAAIBEKMAAAAABIhAIMAAAAABIp5BqwjqxZs0aLFy/W66+/nipkzWhoaNDw4cM1YMCAaqcCAAAAoEDJCrDFixdriy22UGNjo2ynCtvrRYSWL1+uxYsXa4cddqh2OgAAAAAKlGwI4uuvv66tt96a4qsd29p6663pGQQAAAD6gKTXgFF8dYztAgAAAPQNTMJRodbWVu2xxx6SpObmZp188smdtr3vvvt02GGHpUoNAAAAQC+X7Bqw9hqnzOjR5bWeM6FHl1eOpqYmNTU1JY8LAAAAoDb1qR6wq666SmPHjtXo0aP1uc99TuvWrdPgwYN1xhlnaNSoURo3bpyWLl0qSXrqqac0btw4jRkzRt/85jc1ePDgNyyvtIfr/vvv1+jRozV69GjttddeWrVqlSRp9erVmjhxonbbbTdNmjRJEZFuhQEAAAD0Kn2mAFuwYIGuu+46zZ49Wy0tLerXr5+uvvpqvfLKKxo3bpzmz5+v//iP/9Bll10mSTrllFN0yimn6OGHH9a2227b7fLPO+88XXTRRWppadEDDzygQYMGSZIeeeQRTZs2TY8//riefvppzZ49u9D1BAAAANB7VW0IYmp333235s6dqzFjxkiSXnvtNW2zzTbabLPNNvRi7b333rrzzjslSXPmzNGNN94oSTrqqKN0+umnd7n8ffbZR6eeeqomTZqkj370oxo+fLgkaezYsRtujx49Wq2trdp3332LWEUAAACgV6rk8qNqXFqUUp/pAYsIHXfccWppaVFLS4ueeOIJTZ06VQMGDNgwC2G/fv20du3aTVr+lClTdPnll+u1117TuHHjtHDhQknSwIEDN7R5M8sHAAAAUPv6TAF24IEHavr06XrhhRckSStWrNCzzz7baftx48bphhtukCRde+213S7/qaee0rvf/W595StfUVNT04YCDAAAAADa9JkCbOTIkTrrrLM0fvx47bnnnvrgBz+oJUuWdNp+2rRpOv/88zV27FgtWbJEQ4YM6XL506ZN0x577KFRo0Zp0KBBOvTQQ3t6FQAAAADUOJczK5/tVkmrJK2TtDYiupx7vampKZqbmze6b8GCBRoxYsSmZ5rYq6++qkGDBsm2rr32Wl1zzTW66aabCotXa9sHAAAAKFdfvAbM9tyO6qZKJuH4QES82IM59Wpz587VSSedpIjQlltuqSuuuKLaKQEAAACocX1mFsRK7bfffpo/f3610wAAAABQR8q9Biwk3WF7ru3JRSYEAAAAAPWq3B6wfSLiedvbSLrT9sKImFnaIC/MJkvSdttt18NpAgAAAEDtK6sHLCKez/9/QdKvJY3toM2lEdEUEU3Dhg3r2SwBAAAAoA50W4DZ3tz2Fm23JY2X9GjRiQEAAABAvSmnB+ztkmbZni/pIUkzIuL2YtMqxgUXXKARI0Zo0qRJhSz/+OOP1/Tp0wtZNgAAAIDa1+01YBHxtKRRPR55atc/bFz58lZ22+Tiiy/Wbbfdph122GHDfWvXrlX//kwGCQAAAKB45c6CWPNOPPFEPf300zr88MM1ZMgQTZ48WePHj9exxx6rZcuW6WMf+5jGjBmjMWPGaPbs2ZKkqVOn6tOf/rT2339/7bjjjrrgggs2LO/KK6/UnnvuqVGjRumYY47ZcP/MmTP1vve9TzvuuCO9YQAAAAA20me6fi655BLdfvvtuvfee3XhhRfqlltu0axZszRo0CAdddRR+tKXvqR9991Xzz33nA4++GAtWLBAkrRw4ULde++9WrVqlXbddVd9/vOf15/+9CedffbZmj17toYOHaoVK1ZsiLNkyRLNmjVLCxcu1OGHH66JEydWa5UBAAAA9DJ9pgBr7/DDD9egQYMkSXfddZcef/zxDY+9/PLLWrVqlSRpwoQJGjhwoAYOHKhtttlGS5cu1T333KOJEydq6NChkqStttpqw3OPOOIIveUtb9HIkSO1dOnShGsEAAAAoLfrswXY5ptvvuH2+vXrNWfOnA0FWamBAwduuN2vXz+tXbtWESHbHS63tH1E9GDGAAAAAGpdn7kGrCvjx4/XhRdeuOHvlpaWLtsfeOCBuv7667V8+XJJ2mgIIgAAAAB0hgJM2fT0zc3N2nPPPTVy5EhdcsklXbbffffddcYZZ+j973+/Ro0apVNPPTVRpgAAAABqmYsYJtfU1BTNzc0b3bdgwQKNGDGix2PVC7YPAAAA6lXjlBllt209Z0KBmaRje25ENLW/nx4wAAAAAEiEAgwAAAAAEqEAAwAAAIBEkhZgTMveMbYLAAAA0DckK8AaGhq0fPlyio12IkLLly9XQ0NDtVMBAAAAULBkP8Q8fPhwLV68WMuWLUsVsmY0NDRo+PDh1U4DAAAAQMGSFWADBgzQDjvskCocAAAAAPQ6TMIBAAAAAIlQgAEAAABAIhRgAAAAAJBIsmvAAAAAAKBbU4dU2H5lMXkUhB4wAAAAAEiEAgwAAAAAEqEAAwAAAIBEuAYMNaVxyoyK2reeM6GgTAAAAIDK0QMGAAAAAIlQgAEAAABAIhRgAAAAAJAIBRgAAAAAJEIBBgAAAACJUIABAAAAQCIUYAAAAACQCAUYAAAAACRCAQYAAAAAifSvdgIAAAD1onHKjLLbtp4zocBMAPRWFGCob1OHVNB2ZXF5AAAAAGIIIgAAAAAkQwEGAAAAAIlQgAEAAABAIhRgAAAAAJAIBRgAAAAAJEIBBgAAAACJUIABAAAAQCJ96nfA+HFEAAAAANVUdg+Y7X62H7H9myITAgAAAIB6VUkP2CmSFkh6W0G5AADQazBqAgBQhLJ6wGwPlzRB0uXFpgMAAAAA9avcIYjTJH1Z0vriUgEAAACA+tbtEETbh0l6ISLm2t6/i3aTJU2WpO22266n8gNQhxjaBQAA+qpyesD2kXS47VZJ10o6wPZV7RtFxKUR0RQRTcOGDevhNAEAAACg9nVbgEXEVyNieEQ0SvqkpHsi4ujCMwMAAACAOsMPMQMAAABAIhX9EHNE3CfpvkIyAQAAAIA6V1EBBgAAUEuY9AdAb8MQRAAAAABIhB4wAEBNoCcDAFAP6AEDAAAAgEQowAAAAAAgEYYgAgCAJCoZRioxlBRAfaIHDAAAAAASoQADAAAAgEQowAAAAAAgEQowAAAAAEiEAgwAAAAAEqEAAwAAAIBEKMAAAAAAIBEKMAAAAABIhB9iBvo4fhgVAAAgHXrAAAAAACAResAAAACqYeqQCtuvLCYPAElRgAEAAEiVFUQUQwA2EQUYgN6NEyKg7+L1D6AOUYABQB2oZDIVJlIBAKB6mIQDAAAAABKhAAMAAACARBiCCACoP1w7BADopegBAwAAAIBE6AHrDN+eAgAAAOhh9IABAAAAQCL0gAEAAPQVjPABqo4eMAAAAABIhAIMAAAAABKhAAMAAACARCjAAAAAACARCjAAAAAASIQCDAAAAAASYRp6vCmNU2aU3bb1nAkFZgKgbJVMQy0xFXU52KYAgDLRAwYAAAAAiVCAAQAAAEAiFGAAAAAAkAgFGAAAAAAkQgEGAAAAAIlQgAEAAABAIhRgAAAAAJAIBRgAAAAAJNJtAWa7wfZDtufbfsz2t1IkBgAAAAD1pn8Zbf4u6YCIWG17gKRZtm+LiAcLzg0AAABdaJwyo6L2rQ0FJQKgbN0WYBERklbnfw7I/0WRSQEAereKT/rOmVBQJgAA1JayrgGz3c92i6QXJN0ZEb/voM1k2822m5ctW9bDaQIAAABA7SurAIuIdRExWtJwSWNt79FBm0sjoikimoYNG9bDaQIAAABA7SvnGrANIuIl2/dJOkTSo4VkhE3GkCAAAACgdytnFsRhtrfMbw+SdJCkhQXnBQAAAAB1p5wesH+V9DPb/ZQVbNdHxG+KTQsAAAAA6k85syD+QdJeCXIBAAAAgLpW1iQcAAAAAIA3jwIMAAAAABKpaBZE4E2ZOqTC9iuLyQMAAACoEgowoBuVTO/P1P4AAADoCkMQAQAAACARCjAAAAAASIQhiAAAAEAVcblD30IBViBeTNgUHDcAAAD1iyGIAAAAAJAIPWB9WSXTwjMlPAAAAPCm0QMGAAAAAInQAwYAAIBeq5JroyWuj0bvRw8YAAAAACRCDxgAALWI63gBoCbRAwYAAAAAidADBgAoHr01AABIogcMAAAAAJKhAAMAAACARCjAAAAAACARCjAAAAAASIRJOAAA6CUq+cHZ1oYCEwEAFIYeMAAAAABIhAIMAAAAABKhAAMAAACARCjAAAAAACARJuEAAABA2SqaLOacCQVmAtQmesAAAAAAIBEKMAAAAABIhAIMAAAAABKhAAMAAACARCjAAAAAACARCjAAAAAASIQCDAAAAAASoQADAAAAgEQowAAAAAAgEQowAAAAAEiEAgwAAAAAEqEAAwAAAIBEKMAAAAAAIBEKMAAAAABIhAIMAAAAABLp310D2++SdKWkd0haL+nSiPhh0YkBAACgxk0dUmH7lcXkAfQi3RZgktZKOi0i5tneQtJc23dGxOMF5wYAAAAAdaXbIYgRsSQi5uW3V0laIOmdRScGAAAAAPWmomvAbDdK2kvS7wvJBgAAAADqWNkFmO3Bkm6Q9MWIeLmDxyfbbrbdvGzZsp7MEQAAAADqQlkFmO0ByoqvqyPiVx21iYhLI6IpIpqGDRvWkzkCAAAAQF3otgCzbUk/lrQgIs4vPiUAAAAAqE/l9IDtI+kYSQfYbsn/fajgvAAAAACg7nQ7DX1EzJLkBLkAQN1onDKj7Lat50woMBMAANCbVDQLIgAAAABg01GAAQAAAEAiFGAAAAAAkAgFGAAAAAAkQgEGAAAAAIlQgAEAAABAIhRgAAAAAJAIBRgAAAAAJNLtDzEDAABo6pAK268sJg8AqHH0gAEAAABAIhRgAAAAAJAIBRgAAAAAJMI1YEBPquQaCa6PAAAA6HPoAQMAAACARCjAAAAAACARCjAAAAAASIRrwIBaxu/yAAAA1BR6wAAAAAAgEQowAAAAAEiEAgwAAAAAEqEAAwAAAIBEKMAAAAAAIBEKMAAAAABIhAIMAAAAABKhAAMAAACARPgh5t6CH9QFAAAA6h49YAAAAACQCD1gACpTSW8tPbUAAAAboQcMAAAAABKhAAMAAACARCjAAAAAACARCjAAAAAASIQCDAAAAAASoQADAAAAgEQowAAAAAAgEQowAAAAAEiEH2IGAAAANtXUIRW0XVlcHqgZ9IABAAAAQCIUYAAAAACQCAUYAAAAACRCAQYAAAAAiXRbgNm+wvYLth9NkRAAAAAA1KtyesB+KumQgvMAAAAAgLrXbQEWETMlrUiQCwAAAADUNa4BAwAAAIBEeqwAsz3ZdrPt5mXLlvXUYgEAAACgbvRYARYRl0ZEU0Q0DRs2rKcWCwAAAAB1gyGIAAAAAJBIOdPQXyNpjqRdbS+2/V/FpwUAAAAA9ad/dw0i4lMpEgEAAACAescQRAAAAABIhAIMAAAAABKhAAMAAACARCjAAAAAACARCjAAAAAASKTbWRABoNY1TplRUfvWcyYUlAkAAOjr6AEDAAAAgEQowAAAAAAgEQowAAAAAEiEAgwAAAAAEqEAAwAAAIBEmAURAAAAKFHJ7LmtDQUmgrpEDxgAAAAAJEIPGAC0N3VIBW1XFpcHAACoOxRgAAAAQK2o5EtCiS8KeyEKMAAAAACdY2RIj+IaMAAAAABIhAIMAAAAABKhAAMAAACARCjAAAAAACARCjAAAAAASIQCDAAAAAASYRp6AAD6qMYpM8pu29pQYCIA0IfQAwYAAAAAiVCAAQAAAEAiFGAAAAAAkAgFGAAAAAAkQgEGAAAAAIlQgAEAAABAIhRgAAAAAJAIBRgAAAAAJEIBBgAAAACJ9K92AgAAAADSaZwyo6L2rQ0FJdJH0QMGAAAAAIlQgAEAAABAIhRgAAAAAJAIBRgAAAAAJEIBBgAAAACJUIABAAAAQCJMQw8A1TZ1SAVtVxaXBwAAKBw9YAAAAACQSFkFmO1DbD9h+0nbU4pOCgAAAADqUbcFmO1+ki6SdKikkZI+ZXtk0YkBAAAAQL0p5xqwsZKejIinJcn2tZL+U9LjRSYGAAAAVIzratHLlTME8Z2SFpX8vTi/DwAAAABQAUdE1w3sj0s6OCI+k/99jKSxEfF/2rWbLGly/ueukp7o+XSTGirpxTqOV42YrGPtx6tGTNax9uNVI2a9x6tGTNax9uNVIybrWPvxqhGzGutYhO0jYlj7O8sZgrhY0rtK/h4u6fn2jSLiUkmXbnJ6vYzt5ohoqtd41YjJOtZ+vGrEZB1rP141YtZ7vGrEZB1rP141YrKOtR+vGjGrsY4plTME8WFJu9jewfZmkj4p6eZi0wIAAACA+tNtD1hErLV9kqTfSuon6YqIeKzwzAAAAACgzpQzBFERcaukWwvOpbdJPZyyGsM3WUfi1UJM1rH241UjZr3Hq0ZM1rH241UjJutY+/GqEbNuLmvqSLeTcAAAAAAAekY514ABAAAAAHoABViV2W60/Wi180jF9lTbp1c7jyLYPtn2AttXVzuXnlbN49T27+o1ZrW2q+3VqWMC2JjtLW3/d7XzAJAeBRjQc/5b0ociYlK1E6knEfG+vhATQMecqcfzlS2VfW4A6GPq8Q3tTbN9o+25th/Lf2C6aP1t/8z2H2xPt/3WogPaPjaPN9/2zwuOdYbtJ2zfpexHugtn+2jbD9lusf0j2/0KjneJpB0l3Wz7S0XGKon5DdsLbd9p+5oEPYv9bF+Wvy7usD2o4HiSqtNbU6WYO9p+xPaY1LGLkPfuLbR9ue1HbV9t+yDbs23/2fbYAuMuSHms2j41X8dHbX+xyFh5vLZtm+xzo/R9PMX7Tcl+vFjSPG38e6RFxNvc9oz8M/FR258oMl7uHEk75Z9T3ys6WPsed9un255aYLxzS3v48hEwpxUU68u2T85v/8D2PfntA21fVVDMMfnrryE/fh6zvUcRsUpiftv2KSV/n9223gXGPDE/RltsP2P73iLj9RUUYB37dETsLalJ0sm2ty443q6SLo2IPSW9rIK/EbO9u6QzJB0QEaMkndLNU95MrL2V/XbcXpI+Kqnwk0vbIyR9QtI+ETFa0jpJhfZKRcSJyn6g/AMR8YMiY0mS7SZJH9M/t2uKHyvcRdJFEbG7pJfy+OgBtneVdIOkEyLi4Wrn04N2lvRDSXtK2k3SUZL2lXS6pK8VGDfZsZq/x50g6d8ljZP0Wdt7FRWvRLLPjWq8j+d2lXRlROwVEc8WHOsQSc9HxKiI2EPS7QXHk6Qpkp6KiNER8T8J4qV2rbLP4jZHSvplQbFmStovv90kabDtAcrebx4oImD+Xn2zpLMkfVfSVRFR9JDyH0s6TpLyXuFPSir0soeIuCQ/lxojabGk84uM11dQgHXsZNvzJT2o7Fu3XQqOtygiZue3r1L2hlGkAyRNj4gXJSkiVhQYaz9Jv46IVyPiZaX5Ee8DJe0t6WHbLfnfOyaIm9K+km6KiNciYpWkWxLEfCYiWvLbcyU1JojZFwyTdJOko0u2b714JiL+GBHrJT0m6e7Ipt79o4o9flIeq/sqe497JSJWS/qV/nkiWKSUnxvVeB+XpGcj4sFEsf4o6aC812a/iFiZKG7diohHJG1je1vboyT9LSKeKyjcXEl7295C0t8lzVFWiO2nggqw3JmSPpjH+m6BcSRJEdEqaXn+Jc94SY9ExPKi4+Z+KOmeiEhxvlH3yvodsL7E9v6SDpL03oh41fZ9khoKDtv+twCK/m0AJ4hRKvVvHVjSzyLiq4njpuQqxPx7ye11kpIMQewDVkpaJGkfZUVKPSk9ZtaX/L1exX7+pDxWq/FalNJ/blTjN2teSRUoIv6U9/R9SNL/2r4jIs5MFT+Rtdr4i/eiz20kabqkiZLeoaxHrBARscZ2q7Le6N9J+oOkD0jaSdKCouJK2krSYEkDlG3PFMfs5ZKOV7ZNr0gQT7aPl7S9pJNSxOsL6AF7oyHKvqV51fZuyoaUFG072+/Nb39K0qyC490t6ci2oZW2tyow1kxJH7E9KP9m6sMFxmpzt6SJtreRsvWzvX2CuCnNkvThfOz5YEkTqp0QNtk/JB0h6VjbR1U5F1RupqQjbL/V9uaSPqJiv3Fvk/Jzoxrv40nZ3lbSqxFxlaTzJL0nQdhVkrZIEKfNUmU9UlvbHijpsAQxr1U2TG6ismKsSDOVDW+eqew1eKKklij2B28vlfQNZcMAzy0wTqlfKxsyO0bSb4sOln8xcbqyURrri47XV9AD9ka3SzrR9h8kPaFsGGLRFkg6zvaPJP1Z0v8rMlhEPGb7bEn3214n6RFl36YUEWue7esktUh6VglOTCLicdtfl3RHPkZ6jaQv5PHrQkQ8bPtmSfOVrVezsp4U9Iyk3/ZHxCu2D5N0p+1XIuKmlPGx6fL3uJ9Keii/6/J86FXRkn1uVON9vAreLel7ttcr+8z4fNEBI2J5PinNo5JuK/o6sLyX6ExJv5f0jKSFRcbLYz6WF+1/iYglBYd7QNn17XPy99TXVeCxavtYSWsj4hfOJvr6ne0DIuKeomJKUkT8I58I46WIWFdkrNxJynr67rUtSc0R8ZkEceuai/1iAEBRbA+OiNX57GczJU2OiHnVzqvW5T3D8yKi3npNUSdsN0r6TT5ZRDXiT5W0OiLOq0Z8oC/Lv1ieJ+njEfHnaueDTcMQRKB2XZpPMjJP0g0UX29ePgxpjrIhSAAA9Bq2R0p6UtmERhRfNYweMAAAAABIhB4wAAAAAEiEAgwAAAAAEqEAAwAAAIBEKMAAAAAAIBEKMAAAAABIhAIMAAAAABL5/1FAYdz0BTqcAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1080x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fr_freq  = np.mean(X[y==-1,:],axis=0)\n",
    "en_freq = np.mean(X[y== 1,:],axis=0) \n",
    "\n",
    "fr_f = pd.DataFrame(columns=[\"alphabet\",\"frequency\"])\n",
    "fr_f['alphabet'] = alphabet\n",
    "fr_f[\"frequency\"] = fr_freq\n",
    "\n",
    "en_f = pd.DataFrame(columns=[\"alphabet\",\"frequency\"])\n",
    "en_f['alphabet'] = alphabet\n",
    "en_f[\"frequency\"] = en_freq\n",
    "\n",
    "\n",
    "rg = np.arange(len(alphabet))\n",
    "y1 = list(en_f[\"frequency\"].values)\n",
    "y2 = list(fr_f[\"frequency\"].values)\n",
    "width = 0.40\n",
    "  \n",
    "plt.figure(figsize=(15,5))    \n",
    "plt.bar(rg-0.2, y1, width) \n",
    "plt.bar(rg+0.2, y2, width) \n",
    "plt.xticks(rg, alphabet) \n",
    "plt.legend([\"english\",\"french\"])\n",
    "plt.title(\"Log fréquence des symboles en anglais et en français\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e897ee7-6745-4bf6-8475-1c94f5ec914b",
   "metadata": {},
   "source": [
    "### *Commentaire :*\n",
    "\n",
    "On remarque très clairement que chaque langue est caractérisée par une destribution des fréquences bien distinguée. Pour l'anglais par exemple on remarque un grand nombre d'occurences des lettre (w, y, k) tandis que le français est caractérisé par une faible fréquence de ces lettres en particulier."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2e197ba-2c42-4d71-8e37-7e46acb6e25b",
   "metadata": {
    "tags": []
   },
   "source": [
    "## EXO 2 : Classifieur de Bayes naïf"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d72d14b-cd39-4077-8d9f-eb554c57ad7f",
   "metadata": {},
   "source": [
    "-----------------------------------------\n",
    " Supposons que : $f_k(x|\\theta_k) = \\mathcal{N}_{p}(x,\\mu_k,\\Sigma_k = \\sigma_k^2 I_p)$\n",
    " \n",
    "-----------------------------------------\n",
    "### Question 1 ###\n",
    "\n",
    " \n",
    " C'est à dire que la distribution des fréquences des symboles dans chaque classe $k$ est donnée par une distribution normale (avec $p$ le nombre des différents symboles). \n",
    " \n",
    " On calcule d'abord le maximum de vraisemblance pour les $n$ individus $x = (x_1,x_2, . . ,x_n)$\n",
    " \n",
    " $$ f_k(x|\\theta_k) = \\prod_{i=1}^n f_k(x_i|\\theta_k) $$ \n",
    " $$ f_k(x_i|\\theta_k) = \\frac{1}{(2\\pi)^{p/2}|\\Sigma_k|^{1/2}} e^{-\\frac{1}{2}(x_i-\\mu_k)^T\\Sigma_k^{-1}(x_i-\\mu_k)}$$\n",
    "\n",
    "On sait que : $|\\Sigma_k| = |\\sigma_k^2 I_p| = \\sigma_k^{2p}$ et $\\Sigma_k^{-1} = \\frac{1}{\\sigma_k^2}I_p$ \n",
    " \n",
    " \n",
    "Donc : $$ f_k(x_i|\\theta_k) = \\frac{1}{(2\\pi)^{p/2}\\sigma_k^{p}} e^{-\\frac{1}{2\\sigma_k^2}||x_i-\\mu_k||_{_2}^2}$$\n",
    "\n",
    "$$\\implies log f_k(x|\\theta_k)  = \\sum_{i=1}^n log(f_k(x_i|\\theta_k))$$\n",
    "$$ = \\sum_{i=1}^n - \\frac{p}{2}* log(2\\pi) - p log(\\sigma_k) - \\frac{1}{2\\sigma_k^2}||x_i-\\mu_k||_{_2}^2$$\n",
    "$$ = - \\frac{np}{2} log(2\\pi) - np  log(\\sigma_k)  - \\frac{1}{2\\sigma_k^2} \\sum_{i=1}^n||x_i-\\mu_k||_{_2}^2 $$\n",
    "\n",
    " On calcule l'estimateur de   $\\mu_k$ : $\\hat \\mu_k$ : \n",
    "\n",
    "\n",
    "$$\\frac{\\partial}{\\partial \\mu_k }log f_k(x|\\theta_k)  = 0 $$\n",
    "$$- \\frac{1}{2\\sigma_k^2} \\sum_{i=1}^n (x_i-\\mu_k) = 0 $$\n",
    "$$ \\implies \\mu_k = \\frac{1}{n}\\sum_{i=1}^n x_i$$\n",
    "\n",
    "On calcule l'estimateur de   $\\sigma_k$: $\\hat\\sigma_k$ :\n",
    "$$\\frac{\\partial}{\\partial \\sigma_k }log f_k(x|\\theta_k)  = 0 $$\n",
    "$$ \\frac{1}{\\sigma_k^3} \\sum_{i=1}^n ||x_i-\\mu_k||_{_2}^2 - \\frac{np}{\\sigma_k}= 0$$\n",
    "$$ \\implies \\sigma_k = \\sqrt{\\frac{\\sum_{i=1}^n ||x_i-\\mu_k||_{_2}^2}{np}}$$\n",
    "\n",
    "On considère une variable Z qui prends des valeurs dans $\\{1,..,K\\}$ qui représente la classe d'un individu, pour qu'on puisse déterminer une classe pour un nouveau individu, on utilise les paramètres $(\\sigma_k,\\mu_k)$ ,  si on suppose que la densité de chacune des classes est une loi normale.\n",
    "\n",
    "Pour un nouveau individu $X_i$, on a :\n",
    "\n",
    "$$\\hat{Z} = \\underset{k \\in \\{1,..,K\\}}{Argmax}  \\mathbb{P}(Z_i = k| X_i) = \\underset{k}{Argmax}\\ \\ \\frac{\\mathbb{P}(X_i|Z_i = k) \\mathbb{P}(Z_i = k)}{\\mathbb{P}(X_i) } = \\underset{k}{Argmax} \\ \\ \\mathbb{P}(X_i|Z_i = k) \\mathbb{P}(Z_i = k)$$\n",
    " $\\mathbb{P}(Z_i = k)$ représente la portion des individus utilisés pour estimer les paramètres de la classe $k$ , et  $\\mathbb{P}(X_i|Z_i = k) = f_k(X_i|\\theta_k)= \\mathcal{N}(X_i,\\mu_k,\\Sigma_k = \\sigma_k^2 I_p)$, et $\\hat{Z}$ représente la classe la plus probable pour l'individu $X_i$.\n",
    "\n",
    "### Question 2 ###"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 228,
   "id": "8e0c5b0d-606c-43ee-b541-6418a228aaa5",
   "metadata": {},
   "outputs": [],
   "source": [
    "def multivariate_normal(x, p, mean, sd):\n",
    "    x_m = x - mean\n",
    "    return ( (1. / (np.sqrt(2 * np.pi)**p * sd**p)) * np.exp(- x_m.dot(x_m)/(2*sd**2) ))\n",
    "    \n",
    "\n",
    "def param_estimation(X_train,y_train):\n",
    "    n, p = X_train.shape\n",
    "    mu1 = np.mean(X_train[y_train==-1,],axis=0)\n",
    "    mu2 = np.mean(X_train[y_train== 1,],axis=0)\n",
    "    \n",
    "    sigma1 = (np.sum( (X_train[y_train==-1,] - mu1)**2 )/(n*p))**0.5\n",
    "    sigma2 = (np.sum( (X_train[y_train== 1,] - mu2)**2 )/(n*p))**0.5\n",
    "    \n",
    "    pi1 = np.mean(y_train==-1)\n",
    "    pi2 = 1-pi1\n",
    "    return {'mu1':mu1,'mu2':mu2,'sigma1':sigma1,'sigma2':sigma2,'pi1':pi1 ,'pi2':pi2}\n",
    "\n",
    "def Naive_Bayes(X_test,parameters):\n",
    "    _,p = X_test.shape \n",
    "    mu1 = parameters['mu1']\n",
    "    mu2 = parameters['mu2']\n",
    "    \n",
    "    sigma1 = parameters['sigma1']\n",
    "    sigma2 = parameters['sigma2']\n",
    "    \n",
    "    pi1,pi2 = parameters['pi1'],parameters['pi2']\n",
    "    \n",
    "    prob_class1 = np.array([pi1 * multivariate_normal(row,p,mu1,sigma1) for row in X_test])\n",
    "    prob_class2 = np.array([pi2 * multivariate_normal(row,p,mu2,sigma2) for row in X_test])\n",
    "    \n",
    "    y_pred = prob_class1 < prob_class2\n",
    "    y_pred = np.where(y_pred==0,-1,1).astype('float32')\n",
    "    return y_pred\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 229,
   "id": "d93076ff-e922-4d54-8a5f-ae752776a86c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'mu1': array([3.43379684, 1.39853922, 2.70134151, 2.8445924 , 4.2646129 ,\n",
       "        1.6580854 , 1.57837985, 1.3536266 , 3.54980693, 0.62694133,\n",
       "        0.12219046, 3.2414309 , 2.58906729, 3.47442997, 3.25071783,\n",
       "        2.62496227, 1.7342745 , 3.55189402, 3.57585727, 3.58768637,\n",
       "        3.38713151, 2.06893318, 0.03960841, 1.16626074, 0.60862615,\n",
       "        0.600179  , 4.56821106]),\n",
       " 'mu2': array([4.09860046, 2.7545542 , 2.97647478, 3.24476847, 4.48651593,\n",
       "        2.82736209, 2.7393005 , 3.61775248, 4.00523323, 0.78838527,\n",
       "        1.91367111, 3.51747873, 3.02095262, 3.86870484, 4.06149283,\n",
       "        2.48775053, 0.45701799, 3.81093384, 3.91010096, 4.28065609,\n",
       "        3.05244012, 2.38093842, 2.67824179, 0.67521918, 2.77301047,\n",
       "        0.25260041, 5.16374081]),\n",
       " 'sigma1': 0.322941610105768,\n",
       " 'sigma2': 0.4145978063358746,\n",
       " 'pi1': 0.5,\n",
       " 'pi2': 0.5}"
      ]
     },
     "execution_count": 229,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "param_estimation(X, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8fb82947-3a53-4b30-98b9-14b52d53ffaa",
   "metadata": {},
   "source": [
    "# EXO 3  : Classifeur markovien\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0740f15a-c994-42b5-bbb5-b56de8fcc197",
   "metadata": {},
   "source": [
    "------------------------------------\n",
    "Supposons que : $f_k(x|\\theta_k) = \\mathcal{MC}(x,\\theta_k=(\\pi_k,A_k))$ \n",
    "\n",
    "------------------------------------\n",
    "\n",
    "### Question 1 : Estimation des paramètres des deux chaines :\n",
    "\n",
    "Pour cette question, on va faire l'estimation des paramètres d'une seule chaîne par le maximum de vraisemblance (MLE), et on appliquera le résultat pour les deux chaines.\n",
    "\n",
    "Pour une séquence $x_{1:T}$ des données d'entrainement on a la probabilité de réalisation d'une chaine est donnée par :\n",
    "\n",
    "\n",
    "$$p(x_{1:T}|\\theta) = \\pi(x_1) \\ A(x_1, x_2) \\ A(x_2, x_3) ... A(x_{T-1})A(x_{T})$$\n",
    "\n",
    "\n",
    "$$= \\pi(x_1) \\ \\prod_{t=2}^T  p(x_t|x_{t-1})$$\n",
    "\n",
    "\n",
    "$$= \\pi(x_1) \\ \\prod_{t=2}^T  A(x_t,x_{t-1}) $$\n",
    "\n",
    "\n",
    "$$= \\prod_{j=1}^K \\pi_j^{\\mathbb{1}_{\\{x_1=j\\}}} \\ \\prod_{t=2}^T \\prod_{j=1}^K \\prod_{k=1}^K (A_{jk})^{\\mathbb{1}_{\\{x_t=k,x_{t-1}=j\\}}}$$\n",
    "\n",
    "\n",
    " Or, pour un ensemble de séquences $\\mathcal{D} = \\{X_1,X_2,..,X_N\\}$ tel que $X_i=(x_{i,1},..,x_{i,T_i})$ \n",
    " \n",
    " nous avons :\n",
    "  $p_\\theta(\\mathcal{D}) = \\prod_{i=1}^N p_\\theta(X_i)$, le MLE pour $\\mathcal{D}$ est donné par:\n",
    "\n",
    "\n",
    "$$log \\ p_\\theta(\\mathcal{D}) = \\sum_{i=1}^N log (p_\\theta(X_i)) $$\n",
    "\n",
    "$$=  \\sum_{i=1}^N log \\big[\\prod_{j=1}^K \\pi_j^{\\mathbb{1}_{\\{x_1=j\\}}} \\ \\prod_{t=2}^{T_i} \\prod_{j=1}^K \\prod_{k=1}^K A_{jk}^{\\mathbb{1}_{\\{x_{i,t}=k,x_{i,t-1}=j\\}}} \\big]$$\n",
    "\n",
    "$$=\\sum_{j}N^{1}_{j}log(\\pi_j) + \\sum_{j}\\sum_{k}N_{jk}log(A_{jk})$$\n",
    "\n",
    "\n",
    "\n",
    "D'où l'expression du MLE: \n",
    "\n",
    "$$log \\ p_\\theta(\\mathcal{D}) = \\sum_{j} N_j log(\\pi_j) +\\sum_{j} \\sum_{k} N_{jk}log(A_{jk})$$ \n",
    "\n",
    "avec $N^1_j=\\sum_{i=1}^{N}\\mathbb{1}_{\\{x_{i1}=j\\}}$ et $N_{jk} = \\sum_{i=1}^{N}\\sum_{t=1}^{T_i-1}\\mathbb{1}_{\\{x_{i,t}=j, x_{i,t+1}=k\\}}$.\n",
    "\n",
    "\n",
    " On se réduit donc à un problème d'optimisation du log-vraisemblance ci-dessus, sous les contraintes suivantes :\n",
    " \n",
    " $$\\sum_{j=1}^K \\pi_j=1$$ et $$\\sum_{k=1}^K A_{jk}=1$$\n",
    " \n",
    "Il sera donc nécessaire de considérer le Lagrangien donné par : \n",
    "$$\\mathcal{L}(\\pi,A,\\lambda,\\nu)=log\\ p(\\mathcal{D}|\\theta)+\\lambda(\\sum_{j=1}^K \\pi_j-1)+\\nu(\\sum_{k=1}^K A_{jk}-1)$$\n",
    "\n",
    "\n",
    " \n",
    " \n",
    "L'estimateur de $\\pi_j$ : \n",
    "\n",
    "\n",
    "$$\\frac{\\partial}{\\partial \\pi_j } \\mathcal{L}(\\pi,A,\\lambda,\\nu) = 0$$ \n",
    "\n",
    "$$\\iff \\frac{\\partial}{\\partial \\pi_j } log\\ p(\\mathcal{D}|\\theta) + \\lambda = 0$$\n",
    "$$\\iff N_j\\frac{1}{\\pi_j}+ \\lambda = 0 $$\n",
    "$$\\iff N_j = -\\lambda \\pi_j$$\n",
    "$$\\iff \\sum_{j=1}^K N_j = -\\lambda \\sum_{j=1}^K \\pi_j$$\n",
    "\n",
    "\n",
    "Avec $\\sum_{j=1}^K \\pi_j=1$ on obtient : $\\lambda= -\\sum_{j=1}^K N_j$ donc $\\pi_j = \\frac{N_j}{\\sum_{j=1}^K N_j}$\n",
    "\n",
    "\n",
    "L'estimateur de $A_{jk}$ : \n",
    "\n",
    "$$\\frac{\\partial}{\\partial A_{jk} } \\mathcal{L}(\\pi,A,\\lambda,\\nu) = 0 $$\n",
    "$$\\iff \\frac{\\partial}{\\partial A_{jk} } log\\ p(\\mathcal{D}|\\theta) + \\lambda = 0 $$\n",
    "$$\\iff N_{jk}\\frac{1}{A_{jk}}+ \\nu = 0 $$\n",
    "$$\\iff N_{jk} = -\\nu A_{jk}$$\n",
    "$$\\iff \\sum_{k=1}^K N_{jk} = -\\nu \\sum_{k=1}^K A_{jk}$$\n",
    "\n",
    "\n",
    "Avec $\\sum_{k=1}^K A_{jk}=1$ on a : $\\nu = -\\sum_{k=1}^K A_{jk}$ donc $A_{jk}=\\frac{N_{jk}}{\\sum_{k=1}^K N_{jk}}$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 230,
   "id": "f77e6fa5-e429-44ec-9cc8-a311a7bb1e84",
   "metadata": {},
   "outputs": [],
   "source": [
    "alphabet = list('abcdefghijklmnopqrstuvwxyz ')\n",
    "\n",
    "def strip_accents(text):\n",
    "    import unicodedata\n",
    "    return ''.join(c for c in unicodedata.normalize('NFD', text)\n",
    "                  if unicodedata.category(c) != 'Mn')\n",
    "\n",
    "    \n",
    "    return(strip_accents(text.upper()))\n",
    "def cleaning_text(text):\n",
    "    special_char = {'~', ':', \"'\", '+', '[', '\\\\', '@', '^',\n",
    "                    '{', '%', '(', '-', '\"', '*', '|', ',', \n",
    "                    '&', '<', '`', '}', '.',\"’\", '_', '=', ']', \n",
    "                    '!', '>', ';', '?', '#', '$', ')', '/'}\n",
    "    for char in special_char:\n",
    "        text = text.replace(char, \"\")\n",
    "    \n",
    "    re_print = re.compile('[^%s]' % re.escape(string.printable))\n",
    "    table = str.maketrans('', '', string.punctuation)\n",
    "    text = normalize('NFD', text).encode('ascii', 'ignore')\n",
    "    text = text.decode('UTF-8')\n",
    "    text = text.split()\n",
    "    text = [word.lower() for word in text]\n",
    "\n",
    "    text = [word.translate(table) for word in text]\n",
    "    \n",
    "    text = [re_print.sub('', w) for w in text]\n",
    "\n",
    "    text = [word for word in text if word.isalpha()]\n",
    "\n",
    "    return ' '.join(text)\n",
    "\n",
    "\n",
    "\n",
    "os.chdir(r'/Users/mohamedennatiqi/Desktop/main dir/Bureaux/Unsupervised learning/projet/V2/data')\n",
    "myFiles = glob.glob('*.txt')\n",
    "\n",
    "text_data = pd.DataFrame(columns= [\"text\", \"language\"])\n",
    "texts = []\n",
    "\n",
    "def read_text_file(file_path):\n",
    "    with open(file_path, 'r') as f:\n",
    "        return(f.read())\n",
    "        \n",
    "        \n",
    "for file in myFiles:\n",
    "    if file.endswith(\".txt\"):\n",
    "        file_path = f\"{'/Users/mohamedennatiqi/Desktop/main dir/Bureaux/Unsupervised learning/projet/V2/data'}/{file}\"\n",
    "        texts.append(cleaning_text(read_text_file(file_path)))\n",
    "        \n",
    "text_data.text, text_data.language = texts, [1 if (\"english\" in file_name) else -1 for file_name in myFiles]\n",
    "\n",
    "alphabet = list('abcdefghijklmnopqrstuvwxyz ')\n",
    "N = len(alphabet)\n",
    "\n",
    "def to_sentences(doc):\n",
    "    return doc.strip().split('$')\n",
    "\n",
    "def nb_lettres(text):\n",
    "    freq = [text.count(x) for x in alphabet]\n",
    "    return freq\n",
    "\n",
    "nb_texts = text_data.shape[0]\n",
    "X = np.zeros((nb_texts, N))\n",
    "y = np.zeros(nb_texts)\n",
    "k = 0\n",
    "for text in text_data.text:\n",
    "    text_cleaned = cleaning_text(text)\n",
    "    X[k,:] = nb_lettres(text)\n",
    "    k=k+1\n",
    "    y = np.array(text_data.language.values)\n",
    "\n",
    "X = np.log(X+1)\n",
    "y = np.array(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 231,
   "id": "0b97539f-593f-4db0-8dc4-e3491c718a49",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_sequences = text_data.text.values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 232,
   "id": "0eb8128b-55e0-428a-8e72-ec24df6ef6c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "def MC_estimation_per_class(X):\n",
    "    def Ajk_compute(X,j,k):\n",
    "        Njk = 1\n",
    "        for i in range(len(X)):\n",
    "            Njk = Njk + ''.join(X[i]).count(j+k)\n",
    "        return Njk\n",
    "        \n",
    "    A = pd.DataFrame(np.zeros((N,N)),columns=alphabet,index=alphabet)\n",
    "    for j in alphabet:\n",
    "        for k in alphabet:\n",
    "            A[j][k] = Ajk_compute(X,j,k)\n",
    "            \n",
    "    A = A.div(A.sum(axis=1), axis=0)\n",
    "    \n",
    "    X_0  = np.array([sequence[0] for sequence in X])\n",
    "    \n",
    "    Njs  = np.array([sum(X_0==letter) for letter in alphabet]) \n",
    "    Pi = Njs/sum(Njs)\n",
    "    return A,Pi\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9fc848d4-9385-4bc7-885d-dec4483988f7",
   "metadata": {},
   "source": [
    "\n",
    "Lorsqu'on veut calculer les estimations numériquement on peut avoir des problèmes de division par zéros, donc il faut remplacer $N_{jk}$ par $N_{jk}+1$.\n",
    "\n",
    " On considère une variable Z qui prends des valeurs dans $\\{1,..,C\\}$ qui représente la classe des séquences, pour qu'on puisse déterminer la classe pour un nouvelle séquence, on utilise les paramètres $(\\pi_k,A_k)$ .\n",
    "\n",
    "\n",
    "Pour une nouvelle séquence $X_n=(x_{n,1},..,x_{n,T})$, on a :\n",
    "\n",
    "\\begin{align} \n",
    "\\hat{Z} &= \\underset{k \\in \\{1,..,C\\}}{Argmax} \\ \\ \\mathbb{P}(Z_t = k| X_n) \\\\\n",
    "&= \\underset{k}{Argmax} \\ \\ log(\\mathbb{P}(X_n|Z_t = k) \\mathbb{P}(Z_t = k))\n",
    "\\end{align}\n",
    "\n",
    "\n",
    "$\\mathbb{P}(Z_t = k)$ représente la portion des individus utilisés pour estimer les paramètres de la classe $k$ , et  $\\mathbb{P}(X_n|Z_t = k) = f_k(X_n|\\theta_k)= \\mathcal{MC}(X_n,\\theta_k=(\\pi_k,A_k))$, et $\\hat{Z}$ représente la classe la plus probable pour l'individu $X_n$.\n",
    "\n",
    "Puisque  $X_n=(x_{n,1},..,x_{n,T})$ on peut écrire $f_k(X_n|\\theta_k) = \\pi_k(x_{n,1}) \\ \\prod_{t=2}^T  A_k(x_{n,t},x_{n,t-1})$ donc : \n",
    "\n",
    "$$\\hat{Z} = \\underset{k \\in \\{1,..,C\\}}{Argmax} \\ \\ log \\big(\\pi_k(x_{n,1}) \\prod_{t=2}^T A_k(x_{n,t},x_{n,t-1}) \\big)$$\n",
    "$$= \\underset{k \\in \\{1,..,C\\}}{Argmax}\\ \\ log(\\pi_k(x_{n,1})) + \\sum_{t=2}^Tlog(A_k(x_{n,t},x_{n,t-1}))$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 233,
   "id": "82eb383f-1ebe-433c-aa57-d0ffbc6eeca5",
   "metadata": {},
   "outputs": [],
   "source": [
    "def MC_estimation(X_train,y_train):   \n",
    "    epsilon = 1e-20\n",
    "    A1,Pi1 = MC_estimation_per_class(X_train[y_train==-1,])\n",
    "    A2,Pi2 = MC_estimation_per_class(X_train[y_train==1,])\n",
    "    \n",
    "    Pi1 = np.where(Pi1<epsilon,epsilon,Pi1-epsilon)\n",
    "    Pi2 = np.where(Pi2<epsilon,epsilon,Pi2-epsilon)\n",
    "    \n",
    "    portion1 = np.mean(y_train==-1)\n",
    "    portion2 = 1-portion1\n",
    "    return {'A1':A1,'A2':A2,'Pi1':Pi1,'Pi2':Pi2,'portion1':portion1 ,'portion2':portion2}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 234,
   "id": "796ee5fe-79c3-4ba8-bac7-3c4e4f8f15b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "def MC_classifier(X_test,parameters):\n",
    "\n",
    "    # Récupération des paramètres \n",
    "    Pi1 = parameters['Pi1']\n",
    "    Pi2 = parameters['Pi2']\n",
    "    A1  = parameters['A1']\n",
    "    A2  = parameters['A2']\n",
    "    portion1   = parameters['portion1']\n",
    "    portion2   = parameters['portion2']\n",
    "    \n",
    "    X0_indexs = [alphabet.index(x[0]) for x in X_test]\n",
    "\n",
    "    prob_class1 = np.log(np.array([Pi1[i] for i in X0_indexs]))\n",
    "    prob_class2 = np.log(np.array([Pi2[i] for i in X0_indexs]))\n",
    "    \n",
    "    \n",
    "    for i in range(len(X_test)):\n",
    "        for t in range(1,len(X_test[i])):\n",
    "            j = X_test[i][t-1]\n",
    "            k = X_test[i][t]\n",
    "           \n",
    "            prob_class1[i]=prob_class1[i] + np.log(A1[j][k])\n",
    "            prob_class2[i]=prob_class2[i] + np.log(A2[j][k])\n",
    "            \n",
    "    y_pred = np.log(portion1)+prob_class1 < np.log(portion2)+ prob_class2\n",
    "    y_pred = np.where(y_pred==0,-1,1).astype('float32')\n",
    "    return y_pred\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 235,
   "id": "b015c00d-bee4-4577-875a-0e80c25493d0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]\n",
      "[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]\n",
      "[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]\n",
      "[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]\n",
      "[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]\n"
     ]
    }
   ],
   "source": [
    "kf = KFold(n_splits=5)\n",
    "\n",
    "for train_index, test_index in kf.split(X_sequences):\n",
    "    X_train, X_test = X_sequences[train_index], X_sequences[test_index]\n",
    "    y_train, y_test = y[train_index], y[test_index]\n",
    "    parameters = MC_estimation(X_train,y_train)  \n",
    "\n",
    "    y_pred = MC_classifier(X_test,parameters)\n",
    "    \n",
    "    print(abs(y_pred-y_test))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 236,
   "id": "c71e0679-c544-4ecb-bcc5-fc6837f27e7f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def MC_classifier(X_test,parameters):\n",
    "\n",
    "    # Récupération des paramètres \n",
    "    Pi1 = parameters['Pi1']\n",
    "    Pi2 = parameters['Pi2']\n",
    "    A1  = parameters['A1']\n",
    "    A2  = parameters['A2']\n",
    "    portion1   = parameters['portion1']\n",
    "    portion2   = parameters['portion2']\n",
    "    \n",
    "    X0_indexs = [alphabet.index(x[0]) for x in X_test]\n",
    "\n",
    "    prob_class1 = np.log(np.array([Pi1[i] for i in X0_indexs]))\n",
    "    prob_class2 = np.log(np.array([Pi2[i] for i in X0_indexs]))\n",
    "    \n",
    "    \n",
    "    for i in range(len(X_test)):\n",
    "        for t in range(1,len(X_test[i])):\n",
    "            j = X_test[i][t-1]\n",
    "            k = X_test[i][t]\n",
    "           \n",
    "            prob_class1[i]=prob_class1[i] + np.log(A1[j][k])\n",
    "            prob_class2[i]=prob_class2[i] + np.log(A2[j][k])\n",
    "            \n",
    "    y_pred = np.log(portion1)+prob_class1 < np.log(portion2)+ prob_class2\n",
    "    y_pred = np.where(y_pred==0,-1,1).astype('float32')\n",
    "    return y_pred\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 240,
   "id": "ec5f43e2-6fac-411e-a818-da7d9c8a7892",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]\n",
      "------------------------------------\n",
      "[ 1.  1. -1. -1.  1.  1.  1. -1. -1.  1.  1. -1. -1.  1.]\n",
      "[ 1  1 -1 -1  1  1  1 -1 -1  1  1 -1 -1  1]\n"
     ]
    }
   ],
   "source": [
    "kf = KFold(n_splits=5)\n",
    "\n",
    "for train_index, test_index in kf.split(X_sequences):\n",
    "    X_train, X_test = X_sequences[train_index], X_sequences[test_index]\n",
    "    y_train, y_test = y[train_index], y[test_index]\n",
    "    parameters = MC_estimation(X_train,y_train)  \n",
    "\n",
    "    y_pred = MC_classifier(X_test,parameters)\n",
    "    \n",
    "print(abs(y_pred-y_test))\n",
    "print(\"------------------------------------\")\n",
    "print(y_pred)\n",
    "print(y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1b84dbe-78d1-475d-b270-5991ab79e0d3",
   "metadata": {},
   "source": [
    "On remarque que l'algorithme a bien prédit les languages de l'ensemble des données test. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "78346e8d-3d5e-42c1-ad55-be974448b647",
   "metadata": {
    "tags": []
   },
   "source": [
    "# EXO 4 : algorithme de Viterbi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 209,
   "id": "b691cac0-0da3-4b68-999e-7b6f98e923c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_text_file(file_path):\n",
    "    with open(file_path, 'r') as f:\n",
    "        return(f.read())\n",
    "    \n",
    "path = '/Users/mohamedennatiqi/Desktop/main dir/Bureaux/Unsupervised learning/projet/V2/english.txt'\n",
    "def sample_sentences(path,nb_sentences):\n",
    "    txt = read_text_file(path)\n",
    "    sentences = to_sentences(txt)\n",
    "    tmp = []\n",
    "    for j in range(len(sentences)):\n",
    "        tmp.append(sentences[j])\n",
    "            \n",
    "\n",
    "    idx = sample(list(range(len(tmp))),nb_sentences)\n",
    "    sentences = [tmp[x] for x in idx]\n",
    "    return sentences\n",
    "\n",
    "\n",
    "direc = '/Users/mohamedennatiqi/Desktop/main dir/Bureaux/Unsupervised learning/projet/V2/english.txt'\n",
    "sentences_en = sample_sentences(direc,12)\n",
    "\n",
    "\n",
    "direc = '/Users/mohamedennatiqi/Desktop/main dir/Bureaux/Unsupervised learning/projet/V2/french.txt'\n",
    "sentences_fr = sample_sentences(direc,12)\n",
    "\n",
    "y_en = np.ones(len(sentences_en)) \n",
    "y_fr = -np.ones(len(sentences_fr)) \n",
    "\n",
    "sentences = sentences_en + sentences_fr\n",
    "y_s = np.concatenate((y_en, y_fr), axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 210,
   "id": "48d304db-4560-441b-95fc-714bce8b5ae1",
   "metadata": {},
   "outputs": [],
   "source": [
    "array_sentences = [cleaning_text(sentences[i]) for i in range(len(sentences))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 211,
   "id": "cf7a46c1-097f-4846-9d97-6aa93eb9362b",
   "metadata": {},
   "outputs": [],
   "source": [
    "array_sentences, y_s = shuffle(array_sentences, y_s, random_state=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 212,
   "id": "7b1e7dd3-27cc-4922-a0b7-96656da378ef",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "          fr        en\n",
      "fr  0.500000  0.500000\n",
      "en  0.454545  0.545455\n"
     ]
    }
   ],
   "source": [
    "epislon = 1e-20\n",
    "Ph = np.array([0,0])\n",
    "Ph[0]= y[0] ==-1\n",
    "Ph[1]= y[0] == 1\n",
    "\n",
    "\n",
    "Ph = np.where(Ph==0,epislon,Ph-epislon)\n",
    "\n",
    "\n",
    "Ah = pd.DataFrame(np.zeros((2,2)),columns=['fr','en'],index=['fr','en'])\n",
    "for j in range(len(y_s)-1):\n",
    "    k=j+1\n",
    "    if (y_s[j]==1 and y_s[k]==1):\n",
    "        Ah['en']['en'] = Ah['en']['en']+1\n",
    "    if (y_s[j]==-1 and y_s[k]==-1):\n",
    "        Ah['fr']['fr'] = Ah['fr']['fr']+1\n",
    "            \n",
    "    if (y_s[j]==1 and y_s[k]==-1):\n",
    "        Ah['en']['fr'] = Ah['en']['fr']+1\n",
    "    if (y_s[j]==-1 and y_s[k]==1):\n",
    "        Ah['fr']['en'] = Ah['fr']['en']+1\n",
    "        \n",
    "\n",
    "Ah = Ah.div(Ah.sum(axis=1), axis=0) \n",
    "print(Ah)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 213,
   "id": "df43759d-9a2b-4328-936c-5d2605e48dc7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def log_emission_proba(n,k):\n",
    "\n",
    "    Ak = parameters['A1'].copy()\n",
    "    Pi = parameters['Pi1'].copy()\n",
    "    if (k==1):\n",
    "        Ak = parameters['A2'].copy()\n",
    "        Pi = parameters['Pi2'].copy()\n",
    "        \n",
    "    Xn = array_sentences[n]\n",
    "    c = np.log(Pi[alphabet.index(Xn[0])])\n",
    "    for t in range(1,len(Xn)):\n",
    "        xt_1 = Xn[t-1]\n",
    "        xt = Xn[t]\n",
    "        c = c + np.log(Ak[xt_1][xt])\n",
    "    return c\n",
    "\n",
    "#"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d8fa0bc-06e2-43e1-8dd2-53f68b4578b8",
   "metadata": {},
   "source": [
    " On note $(X_n)$ les phrases de langue et $(Z_n)$ leur états cachées i.e (anglais ou français). \n",
    " $(Z_n)$ est une chaîne de Markov avec les paramètres $(A_h,\\pi_h)$ estimés au dessus.\n",
    " $(X_n)$ représente une phrase qui est sous forme d'une séquence de lettres alphabétiques qui a un comportement d'une chaîne de Markov de 27 états (les lettres d'alphabets et espace), Donc chaque phrase $(X_n)$ peut être présenté par une séquence $(x_n)_{\\{t=1,..,T\\}}$ , les paramètres $(A,\\pi)$ de cette chaîne de markov ont été estimés pendant l'exercice 3. Donc la probabilité d'observer une certaine phrase $(X_t^{obs})$ est :\n",
    " \n",
    "$$\\mathbb{P}(X_n^{obs}|Z_n = k) = \\mathbb{P}((x_n^{obs})_{\\{t=1,..,T\\}}|Z_n = k)$$\n",
    "$$= \\pi_k(x_1) \\ \\prod_{t=2}^T  A_k(x_t,x_{t-1})$$\n",
    "\n",
    " Par la suite, si on veut calculer la probabilité d'émission $\\mathbb{P}(X_n^{obs}|Z_n = k)$ on utilise la formule précédente. \n",
    " \n",
    " \n",
    " La fonction 'log_emission_proba' ci-dessous calcul le log des probabilités d'émission, puisque dans l'algorithme de Viterbi on fait introduire des logarithmes pour éviter des problèmes de multiplication par des réels très petits.\n",
    " \n",
    " \n",
    " La fonction 'log_emission_proba' ci-dessous, prends en paramètres (n,j) , n réfère à l'observation $X_n$ c'est à dire à la $n^{ème}$ phrase dans le mélange du texte et j réfère à sa classe (anglais ou français)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 246,
   "id": "1d36010c-66c8-4a6d-ac9c-14cff9520161",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "----------------------------------------------------------------------\n",
      "vrai Libellé : \n",
      "[ 1.  1. -1. -1. -1.  1. -1. -1. -1.  1.  1. -1.  1.  1.  1. -1.  1.  1.\n",
      " -1.  1.  1. -1. -1. -1.]\n",
      "----------------------------------------------------------------------\n",
      "le libellé de Viterbi\n",
      "[ 1.  1. -1. -1. -1.  1. -1. -1. -1.  1.  1. -1.  1.  1.  1. -1.  1.  1.\n",
      " -1.  1.  1. -1. -1. -1.]\n"
     ]
    }
   ],
   "source": [
    "def Viterbi(Pi,A,X):\n",
    "    K,_ = A.shape\n",
    "    T = len(X)\n",
    "    S = np.zeros((K,T))\n",
    "    logV = - np.ones((K,T)) * np.Inf\n",
    "    Zest = np.zeros(T)\n",
    "    for k in range(K):\n",
    "        logV[k,0] = log_emission_proba(0,k) + np.log(Pi[k])\n",
    "        S[k,0]= 0\n",
    "    \n",
    "    #forward\n",
    "    for t in range(1,T):\n",
    "        for k in range(K):\n",
    "            logV[k,t]= np.max(logV[:,t-1]+np.log(A.values[:,k]) + log_emission_proba(t,k))\n",
    "            S[k,t-1] = np.argmax(logV[:,t-1]+np.log(A.values[:,k]) + log_emission_proba(t,k))\n",
    "    \n",
    "    # backtracking       \n",
    "    Zest[T-1]= np.argmax(logV[:,T-1])\n",
    "    for t in range(T-2,-1,-1):\n",
    "        Zest[t] = S[round(Zest[t+1]),t]\n",
    "    Zest = np.where(Zest==0,-1,1)\n",
    "    return Zest.astype(np.float64)\n",
    "        \n",
    "y_hat = Viterbi(Ph,Ah,array_sentences)\n",
    "\n",
    "print(\"----------------------------------------------------------------------\")\n",
    "print(\"vrai Libellé : \")\n",
    "print(y_s)\n",
    "\n",
    "print(\"----------------------------------------------------------------------\")\n",
    "print('le libellé de Viterbi')\n",
    "print(y_hat)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f495041c-c229-4524-8fb6-4ebf363a01ab",
   "metadata": {},
   "source": [
    "Les parcours données par l'algorithme de Viterbi coincide bien avec celui donné par les vrai libéllés des textes. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1369cf8-1deb-4243-aecd-24ffdbed7238",
   "metadata": {},
   "source": [
    "# EXO 5(facultatif) : Baum Welch"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f41e5c7f-ae64-4a71-9959-de428e52df54",
   "metadata": {},
   "source": [
    "Le but de cet exercice est de coder l'algorithme de **Baum Welch**, qui consiste à trouver les paramètres du MLE, avec une approche **EM**. Comme d'habitude avec EM, nous devons veiller à initialiser les paramètres avec soin, afin de minimiser le risque de rester bloqué dans des optima locaux médiocres. Il existe plusieurs façons de le faire, comme par exemple : \n",
    "\n",
    "- Utilisez des données entièrement étiquetées pour initialiser les paramètres.\n",
    "- Ignorez initialement les dépendances de Markov et estimez les paramètres d'observation en utilisant les méthodes d'estimation standard des modèles de mélange, telles que **K-means** ou **EM**.\n",
    "- Initialiser les paramètres de façon **aléatoire**, utiliser plusieurs recommencements et choisir la meilleure solution.\n",
    "\n",
    "Mais en pratique d'après plusieurs travaux en littérature, l'initialisation des paramètres pour l'algorithme de Baum Welch se fait avec un apprentissage en utilisant l'algorithme de **Viterbi** vu dans l'exercice précedent."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2004b373-0199-4cb5-85df-4cfb95c2bd7b",
   "metadata": {},
   "source": [
    "### - l'étape E:\n",
    "Dans cette étape on calcule la quantité $Q(\\theta, \\theta_{old})$ donnée par :\n",
    "\n",
    "$$Q(\\theta, \\theta^{old}) = \\sum_{k=1}^{2}\\mathbb{E}(N^1_k) + \\sum_{k=1}^{2}\\sum_{j=1}^{2}\\mathbb{E}(N_{jk})log(A_{jk})+ \\sum_{i=1}^{N}\\sum_{t=1}^{T_i}\\sum_{k=1}^{2}p(z_t =k\\ / \\ X_i, \\theta^{old})\\ log \\ p(X_{it}, \\phi_k)$$\n",
    "\n",
    "Où les quantités $N^1_k$, $N_{jk}$, $N_{j}$ sont telque :\n",
    "$$\\mathbb{E}(N_{jk})=\\sum_{i=1}^{N}\\sum_{t=2}^{T_i}p(z_{i,t_i-1} =j, z_{i, t}=k\\ /\\ X_i, \\theta^{old})$$\n",
    "$$\\mathbb{E}(N^1_k)=\\sum_{i=1}^{N}p(z_{i,1} =k\\ /\\ X_i, \\theta^{old})$$\n",
    "$$\\mathbb{E}(N_{j})=\\sum_{i=1}^{N}\\sum_{t=1}^{T_i}p(z_{i,t_i} =j\\ /\\ X_i, \\theta^{old})$$\n",
    "\n",
    "Ces statistiques suffisantes attendues peuvent être calculées en exécutant des algorithmes **Backward** et **Forward** sur chaque séquence de texte : \n",
    "\n",
    "- Les étapes **Forward** : $\\alpha(X_k) = p(Y_{0:k}, X_k) = \\sum_{X_{k-1}}\\alpha(X_{k-1}) \\ p(X_k \\ / \\ X_{k-1}) \\ p(Y_k \\ / \\ X_{k}) $\n",
    "\n",
    "- Les étapes **Backward**$\\beta(X_k) = p(Y_{k+1:T}, X_k) = \\sum_{X_{k+1}}\\beta(X_{k+1}) \\ p(X_k \\ / \\ X_{k+1}) \\ p(Y_{k+1} \\ / \\ X_{k+1}) $\n",
    "\n",
    "\n",
    "### - l'étape M:\n",
    "\n",
    "Dans l'étape de maximisation on met à jour les quantités $a_{ij}^* = P(X_{k}=j\\ / \\ X_{k-1}=i)$ et $b_{ij}^* = P(Y_{k}=j\\ / \\ X_{k}=i)$ \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 217,
   "id": "63cfe09b-bf45-41be-9eb6-789402173cd9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.4 0.6]\n",
      " [0.4 0.6]]\n",
      "--------------------\n",
      "[[0.5 0.5]\n",
      " [0.5 0.5]]\n",
      "--------------------\n",
      "--------------------\n",
      "--------------------\n",
      "[[0.  0.5]\n",
      " [0.  0.5]]\n",
      "--------------------\n",
      "[[0.5 0.5]\n",
      " [0.5 0.5]]\n",
      "--------------------\n",
      "['1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1', '1']\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "def forward(V, a, b, initial_distribution):\n",
    "    alpha = np.zeros((V.shape[0], a.shape[0]))\n",
    "    alpha[0, :] = initial_distribution * b[:, V[0]]\n",
    " \n",
    "    for t in range(1, V.shape[0]):\n",
    "        for j in range(a.shape[0]):\n",
    "            # Matrix Computation Steps\n",
    "            #                  ((1x2) . (1x2))      *     (1)\n",
    "            #                        (1)            *     (1)\n",
    "            alpha[t, j] = alpha[t - 1].dot(a[:, j]) * b[j, V[t]]\n",
    " \n",
    "    return alpha\n",
    " \n",
    "def backward(V, a, b):\n",
    "    beta = np.zeros((V.shape[0], a.shape[0]))\n",
    " \n",
    "    # setting beta(T) = 1\n",
    "    beta[V.shape[0] - 1] = np.ones((a.shape[0]))\n",
    " \n",
    "    # Loop in backward way from T-1 to\n",
    "    # Due to python indexing the actual loop will be T-2 to 0\n",
    "    for t in range(V.shape[0] - 2, -1, -1):\n",
    "        for j in range(a.shape[0]):\n",
    "            beta[t, j] = (beta[t + 1] * b[:, V[t + 1]]).dot(a[j, :])\n",
    " \n",
    "    return beta\n",
    " \n",
    "\n",
    "def baum_welch(V, a, b, initial_distribution, n_iter=100):\n",
    "    M = a.shape[0]\n",
    "    T = len(V)\n",
    " \n",
    "    for n in range(n_iter):\n",
    "        alpha = forward(V, a, b, initial_distribution)\n",
    "        beta = backward(V, a, b)\n",
    " \n",
    "        xi = np.zeros((M, M, T - 1))\n",
    "        for t in range(T - 1):\n",
    "            denominator = np.dot(np.dot(alpha[t, :].T, a) * b[:, V[t + 1]].T, beta[t + 1, :])\n",
    "            for i in range(M):\n",
    "                numerator = alpha[t, i] * a[i, :] * b[:, V[t + 1]].T * beta[t + 1, :].T\n",
    "                xi[i, :, t] = numerator / denominator\n",
    " \n",
    "        gamma = np.sum(xi, axis=1)\n",
    "        a = np.sum(xi, 2) / np.sum(gamma, axis=1).reshape((-1, 1))\n",
    " \n",
    "        # Add additional T'th element in gamma\n",
    "        gamma = np.hstack((gamma, np.sum(xi[:, :, T - 2], axis=0).reshape((-1, 1))))\n",
    " \n",
    "        K = b.shape[1]\n",
    "        denominator = np.sum(gamma, axis=1)\n",
    "        for l in range(K):\n",
    "            b[:, l] = np.sum(gamma[:, V == l], axis=1)\n",
    " \n",
    "        b = np.divide(b, denominator.reshape((-1, 1)))\n",
    " \n",
    "    return (a, b)\n",
    " \n",
    "\n",
    "def viterbi(V, a, b, initial_distribution):\n",
    "    T = V.shape[0]\n",
    "    M = a.shape[0]\n",
    " \n",
    "    omega = np.zeros((T, M))\n",
    "    omega[0, :] = np.log(initial_distribution * b[:, V[0]])\n",
    " \n",
    "    prev = np.zeros((T - 1, M))\n",
    " \n",
    "    for t in range(1, T):\n",
    "        for j in range(M):\n",
    "            # Same as Forward Probability\n",
    "            probability = omega[t - 1] + np.log(a[:, j]) + np.log(b[j, V[t]])\n",
    " \n",
    "            # This is our most probable state given previous state at time t (1)\n",
    "            prev[t - 1, j] = np.argmax(probability)\n",
    " \n",
    "            # This is the probability of the most probable state (2)\n",
    "            omega[t, j] = np.max(probability)\n",
    " \n",
    "    # Path Array\n",
    "    S = np.zeros(T)\n",
    " \n",
    "    # Find the most probable last hidden state\n",
    "    last_state = np.argmax(omega[T - 1, :])\n",
    " \n",
    "    S[0] = last_state\n",
    " \n",
    "    backtrack_index = 1\n",
    "    for i in range(T - 2, -1, -1):\n",
    "        S[backtrack_index] = prev[i, int(last_state)]\n",
    "        last_state = prev[i, int(last_state)]\n",
    "        backtrack_index += 1\n",
    " \n",
    "    # Flip the path array since we were backtracking\n",
    "    S = np.flip(S, axis=0)\n",
    " \n",
    "    # Convert numeric values to actual hidden states\n",
    "    result = []\n",
    "    for s in S:\n",
    "        if s == 0:\n",
    "            result.append(\"1\")\n",
    "        else:\n",
    "            result.append(\"-1\")\n",
    " \n",
    "    return result\n",
    "  \n",
    "#V = text_data.sentences\n",
    "\n",
    "# Transition Probabilities\n",
    "A = np.ones((2, 2))\n",
    "A = A / np.sum(A, axis=1)\n",
    " \n",
    "# Emission Probabilities\n",
    "b = np.array(((2, 3), (4, 6)))\n",
    "b = b / np.sum(b, axis=1).reshape((-1, 1))\n",
    "print(b)\n",
    "print(\"--------------------\")\n",
    "print(A)\n",
    "print(\"--------------------\")\n",
    "\n",
    "# Equal Probabilities for the initial distribution\n",
    "initial_distribution = np.array((0.5, 0.5))\n",
    " \n",
    "A, b = baum_welch(V, A, b, initial_distribution, n_iter=1000)\n",
    "print(\"--------------------\")\n",
    "print(\"--------------------\")\n",
    "\n",
    "print(b)\n",
    "print(\"--------------------\")\n",
    "print(A)\n",
    "print(\"--------------------\")\n",
    "\n",
    "print(viterbi(V, A, b, initial_distribution))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1563120-060d-480a-b385-5598b14ec7ec",
   "metadata": {},
   "source": [
    "# EXO 6 (facultatif) : Blocs stochastiques : "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 247,
   "id": "4c545622-9ae4-49fa-ab5b-9f84657391a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def distance_pairwise(data):\n",
    "    from sklearn.feature_extraction.text import TfidfVectorizer\n",
    "    corpus = data.text                                                                                                                                                                                                 \n",
    "    vect = TfidfVectorizer(min_df=1, stop_words=\"english\")                                                                                                                                                                                                   \n",
    "    tfidf = vect.fit_transform(corpus)                                                                                                                                                                                                                       \n",
    "    pairwise_similarity = tfidf * tfidf.T\n",
    "\n",
    "    return pairwise_similarity.toarray()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 248,
   "id": "8fcef365-4386-471a-acba-27dd7c42e2b2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(70, 70)\n"
     ]
    }
   ],
   "source": [
    "distances = distance_pairwise(text_data)\n",
    "print(distances.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 249,
   "id": "bbc7398f-0077-4efc-bc13-42e0b1b8f4c6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Médiane des textes :  0.011940656648259648\n"
     ]
    }
   ],
   "source": [
    "median = np.median(distances)\n",
    "print(\"Médiane des textes : \", median)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 250,
   "id": "f85c2368-0a1a-441a-a1dd-406b5ea1fafe",
   "metadata": {},
   "outputs": [],
   "source": [
    "def binary_matrix(distance_matrix, median):\n",
    "    n = distance_matrix.shape[0]\n",
    "    matrix = np.zeros(distance_matrix.shape)\n",
    "    for i in range(n):\n",
    "        for j in range(n):\n",
    "            if distance_matrix[i][j] > median:\n",
    "                state = 0\n",
    "            else:\n",
    "                state = 1\n",
    "            matrix[i][j] = state\n",
    "    return(matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 251,
   "id": "7318fde9-a032-4746-bcfa-597ad74776da",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0. 1. 1. ... 1. 1. 0.]\n",
      " [1. 0. 0. ... 0. 0. 1.]\n",
      " [1. 0. 0. ... 0. 0. 1.]\n",
      " ...\n",
      " [1. 0. 0. ... 0. 0. 1.]\n",
      " [1. 0. 0. ... 0. 0. 1.]\n",
      " [0. 1. 1. ... 1. 1. 0.]]\n",
      "____________________________\n",
      "(70, 70)\n"
     ]
    }
   ],
   "source": [
    "K = binary_matrix(distances, median)\n",
    "print(K)\n",
    "print(\"____________________________\")\n",
    "print(K.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 252,
   "id": "5eb93cdf-82d6-4b03-8976-e05a350596b7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD7CAYAAACscuKmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAQa0lEQVR4nO3dUawcV33H8e+vTiJoCEpc30QWTmsiRYEINXbuyiQKqoKNkUsj6EOpiNTKqoL8EqoggSKnlSrxlqoSgocKyQqBSKSkUSBNlIeAlRK1DyhkL3HAiWMcqEssB3wNRFAeEAn/PuxcZ3256z27M3N2Z8/vI63u3fXM/Gf37t/nzMw581dEYGaL7w9mvQNmloeT3awQTnazQjjZzQrhZDcrhJPdrBC1kl3SPknHJb0s6WBTO2VmzdO019klbQJ+AOwFTgHPArdHxIvN7Z6ZNeWiGuvuAl6OiB8BSHoI+AgwMtm3bNkS27dvB2BlZWXDZZaXl2vs0myNek+pRr33utttI8bwdobXvdDfr8n3Mcqo/Rq1zLDU9zFqnZQYddcdt87Jkyc5e/asNlqmTsv+V8C+iPh49fxvgfdGxCdGrdPr9aLf76+tv+EyXR7RN+o9pRr13utut40Yw9sZXvdCf78m38coo/Zr1DLDUt/HqHVSYtRdd9w6vV6Pfr+/4UJ1jtk32uDv7aWkA5L6kvqrq6s1wplZHXW68aeAq4eebwNOr18oIg4BhwAkxdr/TDlasdzq9kpy9HaaijFNa56j15YSL+UzSP0e1vkeT7NunXh1WvZngWslvVPSJcDHgMdrbM/MWjR1yx4Rr0v6BPANYBNwf0S80NiemVmjpj5BN1Uw6VywRezGD5vmc3U3vj534yEiGj9BZ2Yd4mQ3K0Sds/ETW15eZp6vs6/fpzr7Ms3hyCJcZ592uaZi+zp7b+S6btnNCuFkNyuEk92sEE52s0I42c0K4WQ3K0TWS285TDqCqstTaqfR1HtPvYzWxojIRRllmZtbdrNCONnNCrFw3fiUrumoZUro0jc1om2WI+hsOm7ZzQrhZDcrxMJ14y3dInSxSzj0aopbdrNCONnNCtGpbvy8D4bpale4a1Lns9v53LKbFcLJblaITnXj2+6i1e0SLsLZ7S5YxM82x+HH2JZd0v2Szkg6OvTaZkmHJZ2ofl7R7m6aWV0p3fgvA/vWvXYQeCoirgWeqp6b2Rwbm+wR8V/Az9e9/BHgger3B4C/bHa3ZiMiznsMk7Tho0ltbXderP98N/qcR72+iIbfa46//bQn6K6KiFcBqp9XNrdLZtaG1s/Gu2Sz2XyYNtl/KmkrQPXzzKgFI+JQRPQiore0tDRluIFZdnNTuqBNxlhEKYdCKX/jefqc6uxL7u/ztMn+OLC/+n0/8Fgzu2NmbUm59PZV4NvAdZJOSboDuBfYK+kEsLd6bmZzbOygmoi4fcQ/7Wl4X8aa90E1lsc8Xa2Yp30Zx8NlzQrhZDcrxNyNjU/tFk1aEneae8hPOqW2bpcuxxTeprqddbeTozz1KE1+LyaNN42U9VOWcctuVggnu1khnOxmhZi7Y/Ymj1UnLRgx6/nsuS/11dnfuu91lpes6hQSWf9vKe+j7vmJlHhry/R6vZHbcctuVggnu1kh5q4b35aUyy1dHjE3Te01K4tbdrNCONnNClFMN37RdWlChs2GW3azQjjZzQpRTDfeZ6GtdG7ZzQrhZDcrRNZu/MrKyrmzxm3VRaszJ7yE21LlPmvvqwTzwy27WSGc7GaFmIuz8Tm69LnXn7ZoQFPbajtG6uFSjlttTRov5TOoe3u0tqa+1omXct/4qyV9S9IxSS9Iuqt63WWbzTokpRv/OvCpiHg3cBNwp6Trcdlms05JKdn8akR8t/r9V8Ax4B0saNlms3HmqdbcJCY6QSdpO7ATeAaXbTbrlOQTdJLeBnwN+GRE/DL1BIakA8CB6XbPzJqS1LJLuphBoj8YEV+vXk4q2zxcsrmJHTabtVmWDq8j5Wy8gC8CxyLis0P/5LLNZh2icScZJL0P+G/g+8Dvqpf/gcFx+8PAHwM/Bj4aET8fs61zwWZZ/icHX2efbLmm+Do7RMSGGxib7E1ysl+Yk70+J/voZPdwWbNCONnNCuFkNytE1okwy8vL9Pt9IM/x6bBpjuXq7EuTtd6aPI/RVIxpjm+beh+p9fkm/Y5Nc36hzvd4mnVd683MxnKymxViLuaz51C3TO8iynGZs417FSzK5dnc3LKbFcLJblYIj6BriUfQTbZcUzyCziPozIrnZDcrhJPdrBBOdrNCONnNCuFkNyuEk92sEE52s0IUMza+rrbq0Znl4pbdrBBOdrNCuBufyF1367qUIhFvkfQdSc9XJZs/U73uks1mHZLSjf8NsDsibgB2APsk3YRLNpt1SkrJ5oiI/6ueXlw9ApdstgRdLW+8iFILO26SdIRB8cbDEeGSzWYdk5TsEfFGROwAtgG7JL0nNYCkA5L6kvpT7qOZNWCiS28R8RrwNLAPl2w265Sxl94kLQG/jYjXJL0V+ADwz7xZsvleEks21y0S0fYtjqa5HdCo9btcJGLSQgupn9ss72Y7/Hrdz6BOIYpR8VLXrVMkIuU6+1bgAUmbGPQEHo6IJyR9G3hY0h1UJZsTtmVmMzI22SPie8DODV7/GbCnjZ0ys+Z1agRdShetrXi5ecRes9r6vnRpgpTHxpsVwsluVggXiWiJi0RMtlxTXCTCRSLMiudkNytEp87GW32ekFIut+xmhXCymxXC3fjC5Djjb/PJLbtZIZzsZoVwspsVwsluVggnu1khsp6N79qdaurE6PKdaibdTuq6s7xTzahlhk3z/erSnWrcspsVwsluVghPcW2Jp7hOtlxTPMXVU1zNiudkNytE1mRfXl523S+buVLrzyUne1Xv7TlJT1TPXbLZrEMmadnvAo4NPXfJZrMOSa3iug34C+C+oZddstmsQ1JH0H0OuBu4bOi180o2S5qoZPM8jqCrq26tt2FNFR9Y/zm1ccmza7Xe6iy/3ryMoEsxtmWXdBtwJiJWpgkwXLJ5dXV1mk2YWQNSuvG3AB+WdBJ4CNgt6StMUbJ5aWmpod02s0lNNIJO0q3ApyPiNkn/AvwsIu6VdBDYHBF3j1l/ZiPoZtmdTOURdPVNOoKuyVFsqeu3EbvtEXT3AnslnQD2Vs/NbE4VMzbeLXuzMdyyp+/HpNtpq2X33WWHrP8gF2VSjo1W9wz8PF4VGsVj480K4WQ3K0Qx3Xh3ya10btnNCuFkNytEMd34lEtv7urbmtyXauvy3WXN7Bwnu1khiunGd6ErZvm1NYKujXUvJOX77ZbdrBBOdrNCFNONN9tI7rHxsxxX75bdrBBOdrNCzF03vu7N9iYt2TtPZ+nndb+sOSld97a6927ZzQrhZDcrxNx146fpvk7a/Z3XLvK87pc1J8ctsUZxy25WCCe7WSGc7GaFyHrMvry8TL/fB9q7bfK83O2z7rmHuttqO0bqeZJ5qfU2aplhXa71ljKfPSnZq9JPvwLeAF6PiJ6kzcC/A9uBk8BfR8QvUrZnZvlN0o1/f0TsiIi1/zpcn92sQ5IqwlQtey8izg69dhy4tSrXvBV4OiKuG7OdmVWEyc0VYSZbril1Lm1Nc8hR53s8r7XeAvimpBVJB6rXzqvPDmxYn324ZHNiLDNrQeoJulsi4rSkK4HDkl5KDRARh4BDcH7LbmZ5JbXsEXG6+nkGeBTYRWJ9dpsvEXHusYjxbLSxyS7pUkmXrf0OfBA4CjwO7K8W2w881tZOmll9Kd34q4BHqxMAFwH/FhFPSnoWeFjSHcCPgY+2t5tmVlfW+uy9Xi/aHlRTx/p9qrMvpQ6qqXs/gro8qKZHv9+vdTbezDrOyW5WiLmbzz5LPmNc34UOhSYdtLLot+ny3WXNrBVOdrNCZD0b77HxF+ax8fV5bHz9sfFm1nFOdrNCLNydaurwoJrpttPknWrqnoH3oJrRd6pxy25WCCe7WSE8qGZICYNq5v1qx7zvX5e5ZTcrhJPdrBBOdrNC+Jh9SJOX3kZtZ9bHpLMcueiJMOfzRBgza4WT3awQRU6EydE97MJEmKZGCHoiTHqMuut6IoyZjeVkNytEkWfjF/HM7jRyd7FttpJadkmXS3pE0kuSjkm6WdJmSYclnah+XtH2zprZ9FK78Z8HnoyIdwE3AMdwyWazThl7Nl7S24HngWtiaOFpSjbPe5GIJnk++2TLNcXz2esVibgGWAW+JOk5SfdVNd8mLtm8urqaEM7M2pCS7BcBNwJfiIidwK+ZoMseEYciohcRvaWlpSl308zqSkn2U8CpiHimev4Ig+RfuJLNks57mC2SsckeET8BXpG0djy+B3gRl2w265TU6+x/Dzwo6RLgR8DfMfiPwiWbzToiKdkj4giw0W0r9zS6N2ZDUspC+3ArnYfLmhXCyW5WiGKmuM5ycEcq13qrb9Iprk1OM01dv43YnuJqZuc42c0KsdBTXH3W1uxNbtnNCuFkNytEMSWbc0+d9RTXyZZrKranuLpks1nxnOxmhXCymxVioS+9TaqtWm9m88Atu1khnOxmhShmIkxunggz2XJNca03T4QxK56T3awQTnazQjjZzQrhZDcrRDETYVI0OajGE2Gmi1F34ownwngijFnxnOxmhcg9qGYV+F9gC3A2W+DzObZjL3LsP4mIDSuoZk32c0GlfkSMPrhwbMd27Ma5G29WCCe7WSFmleyHZhTXsR27lNi/ZybH7GaWn7vxZoXImuyS9kk6LullSQdbjnW/pDOSjg69tlnSYUknqp9XtBT7aknfknRM0guS7soVX9JbJH1H0vNV7M/kij20D5skPSfpiZyxJZ2U9H1JRyT1M8e+XNIjkl6q/u435/zMU2RLdkmbgH8F/hy4Hrhd0vUthvwysG/daweBpyLiWuCp6nkbXgc+FRHvBm4C7qzea474vwF2R8QNwA5gn6SbMsVecxdwbOh5ztjvj4gdQ5e8csX+PPBkRLwLuIHB+8/5vseLiCwP4GbgG0PP7wHuaTnmduDo0PPjwNbq963A8Uzv/TFgb+74wB8C3wXemys2sI3BF3s38ETOzx04CWxZ91rrsYG3A/9DdQ5s1t+3UY+c3fh3AK8MPT9VvZbTVRHxKkD188q2A0raDuwEnskVv+pGHwHOAIcjIlts4HPA3cDvhl7LFTuAb0pakXQgY+xrgFXgS9Xhy32SLs0UO1nOZN9ous5CXwqQ9Dbga8AnI+KXueJGxBsRsYNBK7tL0ntyxJV0G3AmIlZyxNvALRFxI4NDxTsl/VmmuBcBNwJfiIidwK+ZdZd9AzmT/RRw9dDzbcDpjPEBfippK0D180xbgSRdzCDRH4yIr+eODxARrwFPMzh3kSP2LcCHJZ0EHgJ2S/pKpthExOnq5xngUWBXptingFNVDwrgEQbJn/XvPU7OZH8WuFbSOyVdAnwMeDxjfKp4+6vf9zM4lm6cBpOOvwgci4jP5owvaUnS5dXvbwU+ALyUI3ZE3BMR2yJiO4O/739GxN/kiC3pUkmXrf0OfBA4miN2RPwEeEXSddVLe4AXc8SeSM4TBMCHgB8APwT+seVYXwVeBX7L4H/eO4A/YnDy6ET1c3NLsd/H4BDle8CR6vGhHPGBPwWeq2IfBf6pej3Lex/aj1t58wRdjvd9DfB89Xhh7fuV8W++A+hXn/t/AFfk/szHPTyCzqwQHkFnVggnu1khnOxmhXCymxXCyW5WCCe7WSGc7GaFcLKbFeL/ATgZXwkge389AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(K, cmap='Greys',  interpolation='nearest')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "158684f1-7168-4fb3-af29-09e68c86af47",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "247acb1a-ab53-4e00-9cc0-88a86af9a70d",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "py36",
   "language": "python",
   "name": "py36"
  },
  "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}