2-4.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "uBLrtV3fB-Br"
   },
   "source": [
    "# Chapter 5 Information Criteria"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "aLlQv0y4B-Bt"
   },
   "source": [
    "## 5.1 Information Criteria"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 69
    },
    "colab_type": "code",
    "id": "T3QLDaDDB-Bu",
    "outputId": "6b44f716-7b51-4e92-aa9c-9022f604567d"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import scipy\n",
    "from scipy import stats\n",
    "from numpy.random import randn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "d3rg5wT-B-B1"
   },
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LinearRegression\n",
    "import itertools  # enumerate the combination"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "HoT1Xo-4B-B6"
   },
   "outputs": [],
   "source": [
    "res = LinearRegression()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "V3-MD0pfB-CA"
   },
   "outputs": [],
   "source": [
    "def RSS_min(X,y,T):\n",
    "    S_min=np.inf\n",
    "    m=len(T)\n",
    "    for j in range(m):\n",
    "        q=T[j]\n",
    "        res.fit(X[:,q],y)\n",
    "        y_hat=res.predict(X[:,q])\n",
    "        S=np.linalg.norm(y_hat-y)**2\n",
    "        if S<S_min:\n",
    "            S_min=S\n",
    "            set_q=q\n",
    "    return(S_min, set_q)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "mEotWhgcB-CF"
   },
   "outputs": [],
   "source": [
    "from sklearn.datasets import load_boston"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "E0c404m5B-CP"
   },
   "outputs": [],
   "source": [
    "boston = load_boston()\n",
    "X=boston.data[:,[0,2,4,5,6,7,9,10,11,12]]\n",
    "y=boston.target"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 35
    },
    "colab_type": "code",
    "id": "OXlWMiHAB-CT",
    "outputId": "309a8a07-6f81-4d68-bfce-4b594fad89c9"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1619.7876085566147 (0, 2, 3, 5, 7, 8, 9)\n"
     ]
    }
   ],
   "source": [
    "n,p=X.shape\n",
    "AIC_min=np.inf\n",
    "for k in range(1,p+1,1):\n",
    "    T=list(itertools.combinations(range(p), k))\n",
    "    # each column consists of the combinations (k out of p)\n",
    "    S_min, set_q=RSS_min(X,y,T)\n",
    "    AIC=n*np.log(S_min/n)+2*k  ##\n",
    "    if AIC<AIC_min:\n",
    "        AIC_min=AIC\n",
    "        set_min=set_q\n",
    "print(AIC_min,set_min)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "4iRrC8TxB-Cc"
   },
   "outputs": [],
   "source": [
    "y_bar=np.mean(y)\n",
    "TSS=np.linalg.norm(y-y_bar)**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 35
    },
    "colab_type": "code",
    "id": "H0hXfapUB-Ch",
    "outputId": "50ad5572-5e54-417e-c19b-40dcf7793292"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7130213433266738 (0, 2, 3, 5, 7, 8, 9)\n"
     ]
    }
   ],
   "source": [
    "D_max=-np.inf\n",
    "for k in range(1,p+1,1):\n",
    "    T=list(itertools.combinations(range(p), k))\n",
    "    S_min, set_q=RSS_min(X,y,T)\n",
    "    D=1-(S_min/(n-k-1))/(TSS/(n-1))\n",
    "    if D>D_max:\n",
    "        D_max=D\n",
    "        set_max=set_q\n",
    "print(D_max,set_max)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "oEhowEYmB-Cl"
   },
   "outputs": [],
   "source": [
    "def IC(X,y,k):\n",
    "    n,p=X.shape\n",
    "    T=list(itertools.combinations(range(p), k))\n",
    "    S, set_q=RSS_min(X,y,T)\n",
    "    AIC=n*np.log(S/n)+2*k\n",
    "    BIC=n*np.log(S/n)+k*np.log(n)\n",
    "    return {'AIC':AIC, 'BIC':BIC}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "_U8mIKRoB-Cp"
   },
   "outputs": [],
   "source": [
    "AIC_seq=[]; BIC_seq=[]\n",
    "for k in range(1,p+1,1):\n",
    "    AIC_seq.append(IC(X,y,k)['AIC'])\n",
    "    BIC_seq.append(IC(X,y,k)['BIC']) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "aDqmjW_AB-Cu"
   },
   "outputs": [],
   "source": [
    "x_seq=np.arange(1,p+1,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "mrr4XViHB-Cy",
    "outputId": "15077e53-54f7-4ae0-97b0-278ac82eb044"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x20e93b8a460>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEWCAYAAABBvWFzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA8xklEQVR4nO3dd3iUVfr/8fcnCRCqUmKhCKiAYkPFwtrAir2iiUbR1XXdte/X9aeuXXd1dV277uragYAiIlZQsSsiKIoKKCIK0kFApCf374/zhExCJpmEmUzK/bquuTJ56j31nlOec2RmOOecc5sqI90BOOecqx88oTjnnEsKTyjOOeeSwhOKc865pPCE4pxzLik8oTjnnEsKTyjOOeeSwhNKHSdpG0krJGUmsG1fSbMrWP+kpFuTG2FiJL0j6bw0nftGSYPSce5NkebXS5KekPSLpPE1cL7XJA1McNuZkg6Ns67Cz0A1YztA0rRkHrOuqtMJJfoS+kVSkzLLS33QJDWOvjS+k/Rb9IZ7XFKXMtssktQiOu7q6It6maT3JO0Ss22pL6Dow3WJpK+i48+W9FzsPqliZj+ZWQszK0z1ueqDVHyhNFD7A4cBHc1s71SfzMyONLOnUn2e6jCz982sR7rjqA3qbEKJksEBgAHHVbL58Gib04HNgN2AicAhMdscCEwysxXR/xeZWQugDfAO8EwFx78XuBS4JNq+OzASODrRx1MdkrJSeXzXcCRSwi2jMzDTzH5LRTzFoh9rdfZ7qrrq6uOucwHHOAsYBzwJxC0KR0Xfw4DjzexTM1tvZsvM7EEzeyxm06OAV8vuH/3yHwr0jHP8bsCFQJ6ZjTWzNWa20swGm9nt5Wx/mqQJZZZdLmlUdP9oSZ9LWi5plqQbY7brIskknSvpJ2BszLKsaJtzJE2R9KukGZL+WE4M10SlsZmSzqjguTtG0iRJSyV9JGnXmHX/T9LP0XmmSTqknP27RvtmRP8/KmlBzPpnJF0Ws0tnSR9GxxwjqV3MtvtGMSyV9IWkvjHr3pF0S7x9Y7ZrDrwGtI9KnysktY9WN5b0dLT/15J6x+zXXtLzkhZK+kHSJRU8Z09KelDSK9GxPpG0XbSu1GsVE/t50f2zo8dwd/Q4Z0j6XbR8lqQF2rjap52kN6JzvSupc8yxd4jWLYleo1PLxPmwpFcl/Qb0K+extJc0Ktp/uqQ/RMvPBf4H9Imew5vK7Nckin/nmGU5klZJ2kJSa0kvR8/nL9H9jmWek79L+hBYCWxb5nnaTtJYSYuj9/FgSZuXCX8vSd9Ex39CUnac1yvuaytpb0kTFD6L8yX9O84xSpV6FT5XV0j6UqGGY1js+SUdr/C5Wi7pe0n9K3jcFb2GFX1XZEsaFD1HSyV9KmnLaN1mkh6TNFfhM3yroh8UkraP3kfLoud2WHmPOS4zq5M3YDrwZ2BPYB2wZcy6J4Fbo/u3A+8mcLypQI/o/jvAedH9xsDfgfditr0RGBTdvwD4sQpxNwN+BbrFLPsUyI3u9wV2IST7XYH5wAnRui6EEtnTQHOgacyyrGibo4HtAAEHRW/MPWKOvR74N9AkWv9bzOOOfd52BxYA+wCZhKQ9M9qvBzALaB8T13ZxHu9PwJ7R/WnADGDHmHW7xzzn3xNKd02j/2+P1nUAFhOSfgbhB8JiIKeyfcuJpy8wu8yyG4HV0fEzgduAcdG6DEJp9vrovbBt9BiOiHP8J6PY9gaygMHA0DKvX1bM9u9Q8l47O3p9zoniuDV6jh6MnvfDCe+dFjHn+pVQum5CKCl/EK1rHr1G50Rx7A4sAnrG7LsM2C96jNnlPJb3gIeAbKAXsBA4OCbWDyp4nz8O/D3m/wuB16P7bYGTCZ+FlsBzwMgyz8lPwE5R7I3KPE/bR++BJkBOFOc9MfvPBL4COhFqDD6k5H294fWv7LUFPgbOjO63APZN5D0VnX880D46/xTggmjd3tHzflh0/g7ADnEe92aVvIZ9if9d8Ufgpeg5ziR8T7aK1r0A/JfwHtkiivWP0boC4G/F7wlg/6p8L9fJEoqk/QlF7mfNbCLhy+T0OJu3BeZWcrztCB/y2Ia1+yQtJXxgLwJuKm/fRI4fy8xWAi8CedG5uwE7AKOi9e+Y2WQzKzKzLwkv8EFlDnOjmf1mZqvKOf4rZva9Be8CYwhVg7Gus1CSehd4BTi17HGA84H/mtknZlZoof56DbAvUEj4MPeU1MjMZprZ93Ee8rvAQZK2iv4fHv3fFWgFfBGz7RNm9m30uJ4lfIkB5AOvmtmr0fPyBjCBkAAq2zdRH0THLyRUb+4WLd+LkLhuNrO1ZjYDeBTIreBYL5jZeDNbT0goVYnlBzN7IopjGOFL8ebo9RoDrCV8oRZ7xczeM7M1hC+CPpI6AccQqqSesFAq/xx4HhgQs++LZvZh9Jyujg0iOsZ+wP8zs9VmNolQKjkrwccxhNLP0enRMsxssZk9b6Ek/yvhB1vZ9/iTZvZ1FPu62BVmNt3M3oiek4WEH0hl93/AzGaZ2ZLo+HnlxFjZa7sO2F5SOzNbYWbjEnzsAPeZ2Zzo/C9R8h44F3g8ir/IzH42s6nlPW6gPxW8hpV8V6wjfDdtH31+J5rZ8qiUchRwWfQdsgC4u8xj7kz4sbjazD6owmOumwmF8Gt5jJktiv4fQvxqr8XA1pUc7yhCVUisS8xsc8Iv3mOA4Yqp8qni8csaQskb/HTCr7OVAJL2kfR2VARfRigBla2+mRXvwJKOlDQuKiIvJTy22P1/sdL13j8SfkmV1Rn4v6i4vDQ6VifCG206cBnhl/0CSUNVUnVU1ruEX1IHEn5JvkN40x8EvG9mRTHbzou5v5Lwq7A4lgFlYtmf0s97vH0TVXb/bIWqqc6EKrLYc18DbFmFY1Ullvkx91cBmFnZZbHH2/BesND+t4TwenYG9ikT9xnAVuXtW472wJLoC7/Yj4Rf1Il4G2gWvZ+7EL5QXwCQ1EzSfyX9KGk54X2xuUq341T0Ht8yes/9HO0/iIo/IxW9xyt6bc8llHqnRlVGxyT0yIN474FOhB/A8cTGXeFrWMl3xTPAaGCopDmS7pDUKDpmI2BuzDH/SyipAFxJqN0Yr1D1+/sqPOa6l1AkNSX8oj5I0jxJ84DLgd0k7VbOLm8CeyumjrYc5bafAETZ/31CFdvh5WzyFtBRMXXuCXgDyJHUi5BYhsSsG0IorXQys82A/xBe4FJhlXdQhd5uzwP/IlQBbk54XLH7t1ZoSyi2DTCnnMPNIlRZbB5za2ZmBQBmNsTMikuKBvwzzmN9l1BC6hvd/4Dwy/eg6P9EzAKeKRNLcyunjSoBVZ2vYRah1BB77pZmdlSle26sOJE3i1m2VXkbVkGn4juSijuRzCHE/W6ZuFuY2Z9i9q3ouZgDtJHUMmbZNsDPiQQVlbCeJby/84CXY5LT/xGqTfcxs1aEHxtQ+n1aUWz/iNbvEu2fz8afkU4x9yt6j8d9bc3sOzPLI3zZ/pPwo7J5OcepilmEKul4Yh93Za9h3O8KM1tnZjeZWU/gd4QfxWdFx1wDtIs5Zisz2ynab56Z/cHM2hOqzR6SFFsirlCdSyjACYQql56EXz29gB2B9ymnOG5mbxK+wF+QtKekLEktJV0g6feSmhHqNd+Od0JJfaLzfV3O8b8j1DMXKDTONY4axHIlXVXe8aIi/HPAnYQvgDdiVrck/DJcLWlv4lfllacxoSpqIbBe0pGUnwRviuI8gPBGe66cbR4FLoh+BUlS86gRsKWkHpIOjhLYasKv5qJyjlH8/KwifOjfNbPlhF/hJ5N4QhkEHCvpCEmZ0fPbt5IfCfHMB9pK2izB7ccDvyp0QmganX9nSXtV9cRR9czPQH50nN9T8ZdLIo6StL+kxsAthLafWcDLQHdJZ0pqFN32krRjgrHOAj4Cboue710Jv9ircr3OEOA0wq/q2B9NLQnviaWS2gA3VOGYxfuvAJZJ6gD8tZxtLpTUMTr+3wjVh2VV+NpKypeUE5Wil0b7lPs+r4LHgHMkHSIpQ1IHSTvE2bay1zDud4WkfpJ2iUp9ywlVWUVmNpdQDX6XpFZRDNtJOijab0DM5+oXQoJL+DHXxYQykFBf/lOUTeeZ2TzgAeAMld+V9hTCL/VhhAaxr4DehNLLwcDHZeuQgQcU9QQiFB+vNbOy1WLFLonO/yDhjfc9cCKh7jSeIcChwHNRfWmxPwM3S/qV0Fj4bAXHKCX6BXhJtM8vhDfYqDKbzYvWzSHU719Qpg63+FgTgD9Ej+sXQgnt7Gh1E0Jnh0XR8bYArq4gtHeBxdGXVPH/Aj5L8HHNAo4nVEcsJPzK+ivVeP9Gj7UAmBEV+eNV1RVvX0hIur2AHwiP+X+EBtPq+AMh9sWExtePqnmcYkMIX8hLCA2v+bDhvXA4oW58DuF1+ifhtUtUHqEjwRxCddUN0Q+0hJjZJ4RSWXtKVynfQ6hKXkToqfl6FWKC0J65B+Gz/AowopxthhC+OGcQPo8bXQCawGvbH/g6+g64l9BxZqN2y6ows/GERva7o/jfJZTyy9u2stewou+KrQjtlcsJnQLepeTSh7MIPz6/IXy2h1NSfbwX8En0mEcBl0ZtS0RVYHF7hQLIrGHP2CjpIeArM3so3bE451xd5hfGwSQqLkk455xLQIMvoTjnnEuOutiG4pxzrhaqt1Ve7dq1sy5duqQ7DOecqzMmTpy4yMxyqrt/vU0oXbp0YcKECZVv6JxzDgBJP27K/l7l5ZxzLik8oTjnnEsKTyjOOeeSot62oTjnXDKtW7eO2bNns3p12UE16p7s7Gw6duxIo0aNknpcTyjOOZeA2bNn07JlS7p06YJUdizKusPMWLx4MbNnz6Zr165JPbZXeTnnXAJWr15N27Zt63QyAZBE27ZtU1LS8oTinHMJquvJpFiqHocnlBirVsFdd8HYsemOxDnn6p6UJRRJj0taIOmrmGW9FGYTnCRpQjSGP9HcFsui5ZMkXR+zT39J0yRNjze/SLI0agR33gkPPJDKszjnXPWNHDkSSUydGmadmDlzJjvvvPOG9ePHj+fAAw+kR48e7L777px33nmsXLmyRmJLZQnlScJ8ArHuAG4ys16E8fvviFn3vpn1im43A0STwzwIHEmY4CpPUs9UBZyVBXl58MorsGRJqs7inHPVV1BQwP77709BQcFG6+bPn8+AAQP45z//ybRp0/j888/p378/v/76azlHSr6UJRQze48w6U+pxUCr6P5mlD8tZ6y9gelmNsPM1gJDCRMtpUz+GcbatTB8eCrP4pxzVbdixQo++OADHnvsMYYOHbrR+gcffJCBAwfSp0+fDctOOeUUttxyyxqJr6a7DV8GjJb0L0Iy+13Muj6SviAkmSvM7GugA2F2vmKzgX3iHVzS+cD5ANtss03Vo/vtN/YYsAs75Ixn0KB2nH9+1Q/hnGsALrsMJk1K7jF79YJ77qlwkxdffJH+/fvTvXt32rZty8SJE2nbtu2G9V999RUDBw5MblxVUNON8n8CLjezTsDlhPmVIUwF29nMdgPuB0ZW5+Bm9oiZ9Taz3jk51Rgws3lz1H5r8jMKeP99mDmzOlE451xqFBQUkJubC0Bubm651V7pVNMllIHApdH95wjzN2Nmy4s3MLNXJT0kqR3wM9ApZv+O0bLUycvj9Ivv4louZsgQuOaalJ7NOVcXVVKSSIUlS5YwduxYJk+ejCQKCwuRxIUXXrhhm5122omJEydy/PEpbRmIq6ZLKHOAg6L7BwPfAUjaSlHH6KjnVwawGPgU6Capq6TGQC4wKqURDhhA14yf2L/TTJ55BnxCS+dcbTB8+HDOPPNMfvzxR2bOnMmsWbPo2rUrs2aVtApcdNFFPPXUU3zyyScblo0YMYL58+fXSIyp7DZcAHwM9JA0W9K5wB+Au6K2kn8QtXcApwBfRcvvA3ItWA9cBIwGpgDPRm0rqbPllnDwweSv+h9Tp8Lnn6f0bM45l5CCggJOPPHEUstOPvlkbrvttg3/b7nllgwdOpQrrriCHj16sOOOOzJ69GhatmxZIzHW2znle/fubdWeYOvxx1ly7hVslbWIiy7O4N//Tm5szrm6Z8qUKey4447pDiNpyns8kiaaWe/qHtOvlC/PiSfSptEKju78NQUFsH59ugNyzrnazxNKeVq3hiOPJP+X+5k3z4dicc65RHhCiSc3l6OXPM1mLdYzaFC6g3HOudrPE0o8xx1HdrNMBnQcxwsvwG+/pTsg55yr3TyhxNO8ORx7LPlz72TFChiV2s7KzjlX53lCqUheHgcse4lOOau82ss55yrhCaUi/fuTsVkrztjyLUaPhgUL0h2Qc64hy8zMpFevXuy2227ssccefPTRR0DtGcLeE0pFmjSBk04if8bNFBbCsGHpDsg515A1bdqUSZMm8cUXX3Dbbbdx9dVXb7RNOoew94RSmbw8dlr5Kb26LvVqL+dcrbF8+XJat2690fJ0DmFf04ND1j39+sEWW5DfchRXjD+Lb7+F7t3THZRzLp3SNHo9q1atolevXqxevZq5c+cytpyL5NI5hL2XUCqTlQUDBpA37UYkY/DgdAfknGuoiqu8pk6dyuuvv85ZZ51FbRo+y0soicjNpf2DD3LwTvMYNGhrbrwRwtjIzrmGKA2j12+kT58+LFq0iIULF5Zans4h7L2Ekojf/Q46dSI/cygzZsC4cekOyDnX0E2dOpXCwsJSMzZCeoew9xJKIjIyIDeXk/59C3/KvoxBg0RMe5dzztWI4jYUADPjqaeeIjMzs9Q2sUPYL1iwgIyMDA488ED69++f8vh8+PpEffYZ7Lknub2n8+YP2zFnDjRunLzDO+dqNx++vnJe5ZWo3XeH7t3JX/MYixfD6NHpDsg552oXTyiJkiA3lyMm/4u2rYv8mhTnnCvDE0pV5ObSiHXk7jSZUaNg+fJ0B+Scq0n1pYkgVY/DE0pV7Lgj7LYb+b/cz+rVMGJEugNyztWU7OxsFi9eXOeTipmxePFisrOzk35s7+VVVXl57HPVVWzX+WEGDWrE2WenOyDnXE3o2LEjs2fP3ui6j7ooOzubjh07Jv24nlCq6rTT0FVXkb/dOG4eewA//wwdOqQ7KOdcqjVq1IiuXbumO4xazau8qqpLF+jThzN+vgMzKChId0DOOVc7eEKpjrw8uk17mX12Xem9vZxzLuIJpToGDICMDPLbj+WLL2Dy5HQH5Jxz6ecJpTq22gr69eO0abeQmekjEDvnHHhCqb7cXHJ+GE//PssYPBiKitIdkHPOpZcnlOo66SRo1Ij8zV5m9mx47710B+Scc+nlCaW62rSBI47guM9vokUL88Z551yD5wllU+Tl0WzOdE7efwHPPQerV6c7IOecSx9PKJviuOOgaVPyGw1l+XJ4+eV0B+Scc+njCWVTtGgBxx5Lv4/+wdZbe7WXc65h84SyqfLyyFy8gLw+P/Lqq7B4cboDcs659PCEsqn694dWrchf9wTr1sFzz6U7IOecS4+UJhRJj0taIOmrmGW9JI2TNEnSBEl7R8sl6T5J0yV9KWmPmH0GSvouug1MZcxVlp0NJ51Er3fuoeeOPvGWc67hSnUJ5Umgf5lldwA3mVkv4Prof4AjgW7R7XzgYQBJbYAbgH2AvYEbJLVOcdxVk5uLfl1O/p5T+PBD+OGHdAfknHM1L6UJxczeA5aUXQy0iu5vBsyJ7h8PPG3BOGBzSVsDRwBvmNkSM/sFeIONk1R6HXII5ORw+i8PATBkSJrjcc65NEhHG8plwJ2SZgH/Aq6OlncAZsVsNztaFm/5RiSdH1WjTajRSXCysmDAADqPfYID9ytk0CCo45O6OedclaUjofwJuNzMOgGXA48l68Bm9oiZ9Taz3jk5Ock6bGJyc2HVKvJ3mMDUqfDZZzV7euecS7d0JJSBQPFs7M8R2kUAfgY6xWzXMVoWb3ntst9+0LEjp8y6m8aN8cZ551yDk46EMgc4KLp/MPBddH8UcFbU22tfYJmZzQVGA4dLah01xh8eLatdMjLgtNNo/fYIjjl8LQUFsH59uoNyzrmak+puwwXAx0APSbMlnQv8AbhL0hfAPwg9ugBeBWYA04FHgT8DmNkS4Bbg0+h2c7Ss9snLg3XryO/0LvPnw1tvpTsg55yrObJ62nrcu3dvmzBhQs2e1Ay6d2dNp+3Z6vPXOOYYeOaZmg3BOeeqS9JEM+td3f39SvlkkiAvjybvjuHUY1YyYgSsWJHuoJxzrmZ4Qkm23FwoKiK/zausXAkvvpjugJxzrmZ4Qkm2nj1h113Zb/zddO7svb2ccw2HJ5RUyMsjY9xHnHH0UsaMgfnz0x2Qc86lnieUVDjtNADOaDKcoiIYOjTN8TjnXA3whJIKXbvCvvvSc+wD7LGHV3s55xoGTyipkpsLX3xB/mHzmTABpk1Ld0DOOZdanlBS5dRTISOD3LVPk5EBgwenOyDnnEstTyipsvXW0LcvW7/8KIceaj4CsXOu3vOEkkq5ufDdd5yx34/88AN8/HG6A3LOudTxhJJKJ58MWVmcuPh/NG3qjfPOufrNE0oqtWkDRxxByxee5oTjjWHDYO3adAflnHOp4Qkl1fLyYNYs8vf4hiVL4PXX0x2Qc86lhieUVDvuOMjO5rAZ/yUnx6u9nHP1V6UJRdKlklpFE189JukzSYfXRHD1QsuWcOyxNBoxjNxTixg1CpYtS3dQzjmXfImUUH5vZssJMyW2Bs4Ebk9pVPVNXh4sWEB+9/GsWQPPP5/ugJxzLvkSSSiK/h4FPGNmX8csc4k48kho1Yq9Pn+Ebt282ss5Vz8lklAmShpDSCijJbUEilIbVj2TnQ0nnoheGEF+7nreeQdmzUp3UM45l1yJJJRzgauAvcxsJdAYOCelUdVHubmwbBlntH8bMygoSHdAzjmXXIkkFAN6ApdE/zcHslMWUX11yCHQrh3bvfMYffp4tZdzrv5JJKE8BPQB8qL/fwUeTFlE9VWjRnDKKfDSS+QPWMPkyfDll+kOyjnnkieRhLKPmV0IrAYws18I1V6uqvLyYOVKTm32MllZXkpxztUviSSUdZIyCVVfSMrBG+WrZ//9oUMH2r36NEceCUOGQGFhuoNyzrnkSCSh3Ae8AGwh6e/AB8A/UhpVfZWREaYHfu018k9Ywc8/w3vvpTso55xLjkoTipkNBq4EbgPmAieY2XOpDqzeysuDdes4ds1wWrb0ai/nXP2RyNAr2wArgZeAUcBv0TJXHXvuCdttR9MRgznlFBg+HFatSndQzjm36RKp8noFeDn6+xYwA3gtlUHVa1IopYwdS/5RS1i+HF5+Od1BOefcpkukymsXM9s1+tsN2BvwuQc3RW4uFBVx0JwCOnTwai/nXP1Q5eHrzewzYJ8UxNJw7LQT7LILmcOGcPrp8OqrsGhRuoNyzrlNk0gbyl9ibldIGgLMqYHY6re8PPjoI87oN4f16+E57+bgnKvjEimhtIy5NSG0pRyfyqAahNNOA2DXLwex885e7eWcq/uyKtvAzG6qiUAanG23hX32QUMLyM+/kquughkzwmLnnKuL4pZQJL0kaVS8W00GWW/l5sKkSZy+93QABg9OczzOObcJZGblr5AOqmhHM3u3wgNLjwPHAAvMbOdo2TCgR7TJ5sBSM+slqQswBZgWrRtnZhdE++wJPAk0BV4FLrV4Qcfo3bu3TZgwobLN0mvOHOjYEa6/nn7v3sicOTB1auhZ7JxzNU3SRDPrXd3941Z5VZYwEvAk8ADwdMwxTyu+L+kuIHZ29e/NrFc5x3kY+APwCSGh9Ke+XAfTvj0cdBAMHUr+FTdw3h/EhAmw117pDsw556oukV5e3SQNl/SNpBnFt8r2M7P3gCVxjingVKDCaaYkbQ20MrNxUankaeCEys5dp+TlwbRpnNx9Mk2aeOO8c67uSqSX1xOEUsJ6oB/hS31Tv/YOAOab2Xcxy7pK+lzSu5IOiJZ1AGbHbDM7WlZ/nHwyZGWx+SuDOfbYMJPjunXpDso556oukYTS1MzeIrS3/GhmNwJHb+J58yhdOpkLbGNmuwN/AYZIalXVg0o6X9IESRMWLly4iSHWkLZt4fDDQ7XX6UUsXAhvvpnuoJxzruoSSShrJGUA30m6SNKJQIvqnlBSFnASMKx4mZmtMbPF0f2JwPdAd+BnoGPM7h2jZeUys0fMrLeZ9c7JyaluiDUvLw9++okjW4+jTRuv9nLO1U2JJJRLgWaEOeX3BPKBgZtwzkOBqWa2oSpLUk40iReStgW6ATPMbC6wXNK+UbvLWcCLm3Du2un44yE7m8bPF3DqqTByJKxYke6gnHOuahJJKIVmtsLMZpvZOWZ2spmNq2wnSQWEQSR7SJot6dxoVS4bN8YfCHwpaRIwHLjAzIob9P8M/A+YTii51I8eXrFatoRjjoFnnyU/dz0rV4ak4pxzdUnc61A2bCC9DWxF+KIfZmZf1URgm6pOXIcS6/nn4ZRTsDFvsO35h9KjB7z+erqDcs41JJt6HUoiw9f3I/TuWgj8V9JkSddW94QujqOOgpYto6FY4I03YN68dAflnHOJS2j4ejObZ2b3ARcAk4DrUxlUg9S0KZx4IowYQf6pazGDq69Od1DOOZe4RC5s3FHSjZImA/cDH1G655VLltxcWLqUHj+8zrXXwpNPwlNPpTso55xLTCIllMeBX4AjzKyvmT1sZgtSHFfDdOih4bqUggJuuAH69oU//xm++SbdgTnnXOUSaUPpY2b3mplPqpVqjRrBKafAqFFkrv6NIUOgRQsYMAB++y3dwTnnXMWqPAWwS7G8PFi5El56ia23Dhc5TpkCF1+c7sCcc65inlBqm/33D6MQDx0KwGGHwbXXwhNPeHuKc652q2iCrWxJG41fEl3Vnp3asBqwzMwwPfBrr8HSpQDenuKcqxMqKqHcRxgVuKz9gbtTE44D4PTTYe1auP9+IOQYb09xztV2FSWUPc1sRNmFZvYCYagUlyq9e4cuxDffDJ99BuDtKc65Wq+ihNKsmvu5ZHjwQdhyS8jPh1WrAG9Pcc7VbhUlhgWS9i67UNJehGFYXCq1aRMyx5QppS6ZL25P+dOf4Ouv0xeec86VVVFC+SvwbHSV/LHR7Sbg2WidS7XDDoNLLoF7790w61Zxe0rLlnDqqd6e4pyrPeImFDMbD+wDCDg7ugnYx8w+qYngHHD77bDDDnD22fDLL0BoTxk8OBReLrooveE551yxCttCzGy+md0QzYFyspld78Ou1LCmTUNr/Pz5cOGFGxYfeihcd10Y7+vJJ9MWnXPObVDRdSiTJX1Zzm2ypC9rMsgGb8894cYboaAg3CLXX19yfYq3pzjn0i3uBFuSOle0o5n9mJKIkqTOTbBVmfXr4cADQz3X5MnQMQz4PHcu9OoF7drB+PHQvHl6w3TO1V0pm2DLzH4sewN+A36q7cmkXsrKgqefhnXrQntKURHg7SnOudqjoiqvfSW9I2mEpN0lfQV8BcyX1L/mQnQbbL893H03vPUWPPDAhsXenuKcqw0qqvKaAFwDbAY8AhxpZuMk7QAUmNnuNRdm1dW7Kq9iZnDccaEb8cSJ0LMnAIWFoZfxuHHw6aew005pjtM5V+ekck75LDMbY2bPAfPMbByAmU2t7slcEkjw6KNhYK8zzwxjflFyfUqrVj7el3MuPSpKKEUx91eVWVd+scbVjK22gkceCeN83XxzqcWDB8PUqaV6GDvnXI2oKKHsJmm5pF+BXaP7xf/vUkPxuXhOPBHOOQduuw0++mjD4kMOCd2Jn3rK21OcczUrbhtKXVdv21BiLV8Ou+0GGRnwxRehGgxvT3HOVU/K2lAkTZR0r6T+PqFWLdWqVehK/MMP8Je/bFjs7SnOuXSoqMprH+AFoC/wrqRXJV0qqXuNROYSc8ABcOWVoaH+pZc2LPb2FOdcTUu4yktSe6B/dNseGGdmf05hbJukQVR5FVuzBvbZJ1w2P3kybLHFhlU33gg33RRGwj/77LRF6JyrA1JZ5VVqTnkzm2NmjwMXEaYBHlzdk7oka9IkDCC5dCmcf364ViVy3XXQr18Y7+urr9IXonOu/qvOnPL7AXeZ2YepCclVy847hx5fL74YiiORsu0pK1akMUbnXL3mc8rXJ5ddFoojl14KM2ZsWLzVViGpTJsWSir1tGOfcy7NfE75+iQjI1x8kpEBZ50V+g9HDj44TB/8zDN+fYpzLjV8Tvn6Zptt4MEH4cMP4c47S6269tqQWC680NtTnHPJV9HgkHsT5o9/EpgYLe4NnAXk1vZpgBtUL6+yzOC002DkSPjkE9i9ZBzPefPC/CmtW4eLHqNrIZ1zLqXzofic8nWVBA8/HGbdys+H1as3rPL2FOdcqqRsTnlJj0taEM2jUrxsmKRJ0W2mpEkx666WNF3SNElHxCzvHy2bLumqajzGhqlt29Db65tv4JprSq3y9hTnXCpUVOU1mfJHFRZgZrZrhQeWDgRWAE+b2c7lrL8LWGZmN0vqCRQAewPtgTeB4ivyvwUOA2YDnwJ5ZvZNZQ+sQVd5xbrootCm8uabYeTISGEhHH44fPxxmDp4541eIedcQ7OpVV5ZFaw7proHBTCz9yR1KW+dJAGnAgdHi44HhprZGuAHSdMJyQVgupnNiPYbGm1baUJxkTvugDfeCJfJT54Mm28OhOtTBg8O7SkDBnh7inNu01VpTvloLvlOwJWbeN4DgPlm9l30fwdgVsz62dGyeMvLJel8SRMkTVi40DuiAdCsWbiKfu7cjSadL25P+fZb+NOfvD3FObdpErqeJJpT/k5JM4FbgE2dtTGPUMWVVGb2iJn1NrPeOTk5le/QUOy1V5gkZfBgGDas1Kri9pRBg0pdYO+cc1VW0Vhe3SXdIGkqcD/wE6HNpZ+Z3V/dE0rKAk4CYr/ZfiaUfIp1jJbFW+6q6pprYO+9Q1Hk59JP4d/+FppXLrww1Io551x1VFRCmUpo4zjGzPaPkkhhBdsn6lBgqpnNjlk2CsiV1ERSV6AbMJ7QCN9NUldJjYHcaFtXVVlZoVvXmjVhpseikhmei9tTNt8cTj3Vx/tyzlVPRQnlJGAu8LakRyUdQujhlRBJBcDHQA9JsyWdG63KpUx1l5l9TbiI8hvgdeBCMys0s/WE0Y1HA1OAZ6NtXXV07w533RUa6R96qNSqLbf09hTn3KapdD4USc0JPavyCCWWp4EXzGxM6sOrPu82HIcZHH00vP02fP457LBDqdU33xzaVB57DH7/+zTF6JxLi03tNlylOeUltQYGAKeZ2SGVbZ9OnlAqMHcu7LILdOkSLkRp1GjDqsJCOOKIMBTY+PFhM+dcw5CyoVfKY2a/RD2panUycZXYemt45BGYOBFuuaXUKm9Pcc5Vlw9D31CddBIMHAh//zuMG1dqVWx7yh//WGooMOeci8sTSkN2773QqROceeZGRZF+/cJ89EOGhDEmBwwIJZelS9MSqXOuDvCE0pBtthk89RR8/z1cccVGq6+9FkaPDvnmww/DwMU5OWEMsIce2uhyFudcA1elRvm6xBvlq+Cvf4V//Qtefjn0ACtHUVFopB85El54IVSHQbgI/4QT4MQTQ4cxJdyx3DlX29RoL6+6xBNKFaxZEzLDggVhKsd27SrdZerUkFhGjgyJBsJlLiecEG777BNmInbO1R2eUOLwhFJFX34ZksrRR8Pzz1epqPHzzzBqVEguY8fC+vVh4Mnjjgsll379oEmT1IXunEuOGu027OqxXXeFW28NxY6nn67Srh06hKvrR4+GhQtD4/0BB4QG/SOPDO0uubkwdCgsW5ai+J1zaeclFFeisDAMP/z556HE0qXLJh1u9epQYhk5El58MdSoNWoUTnHiiaEEs/XWSYncOZcEXuUVhyeUapo5M5RWevUKw7NkZiblsIWF4XKX4kb9778Py/fdt6TdpUePpJzKOVdNXuXlkqtLF7j/fnj/ffj3v5N22MxM2G8/uPNO+O670PZ/662wbh1cdVXoIdazZxhlf/z4UoMhO+fqCC+huI2ZwSmnwEsvhbmBd9stpaebNStUiY0cCe+8E0oz7dvD8ceHkkvfvtC4cUpDcM7hVV5xeULZRIsWwc47hy7Ezz0HO+5YI6ddsgReeSUkl9dfh5Urw/WXJ50EN90ULux3zqWGV3m51GjXLlxFP316qIs65BAYMSL0CU6hNm3ClfnPPx9y2qhRIZkUFIRqsdtuC5fNOOdqH08oLr4jjgj1UbfdFhLLySdD166h8WP+/JSfvmlTOPZYePxxmDIlhHPNNaHg9OqrKT+9c66KPKG4iuXkhFbzGTNCPdSOO8J114W6pzPOgI8+qpHpHbt0CQWk0aNDA//RR4duxzNmpPzUzrkEeUJxicnMDK3kY8aEcVf+/Ocw9td++8Eee4QpHleuTHkYhx8eLpG5445wjUvPnnD99TVyaudcJTyhuKrr0QPuuSeMufKf/4RuWeedBx07wv/9X6geS6HGjcN4ltOmhVq4W24JBacRI2qksOSci8MTiqu+Fi3CDFxffAHvvQeHHQb33QfdusFRR4XuWoWFKTt9hw5hmJd33w09wU4+ObSzTJ2aslM65yrgCcVtOikM3jVsGPz4Y5iZa9IkOOaYkFzuvBMWL07Z6Q88ED77LOSy8eNhl13gyivh119TdkrnXDk8objkat8ebrghJJZhw0Lj/ZVXhuqw3/8+zGOfAllZcPHFYZ6WgQNDDuvRI5RgvBrMNQRr1sA334TKgnTxCxtd6k2eDA8+CM88E1rP990XLrwwzCuconHtP/kELroIJkwIhaf770/5Bf/OpVxRUWi6/Pbb0IYY+3fmzLA+JycMxFodfqV8HJ5QaqFly8LFkg8+GD4BOTmhMf+CC2CbbZJ+uqKicA3LVVfBL7+Ejmk33wytWyf9VM4l1bJlGyeMadPCOHixPRqbNQsT2/XoUfK3Rw/Yc8/qzZ7qCSUOTyi1WFERvPVWSCwvvRSWHXdcKLUcckjS5xFesiR0LX744XAl/u23wznn+IySLr3Wrg3XUZVNGtOmlS5hZGSE64nLJo7u3UPHlGR+XDyhxOEJpY748cfQ9fh//wtjreywQyhKDBwIrVol9VSTJoVqsA8/DJNTPvAA7L13Uk/hXClmMHdu+aWNH34o3QkyJ6d0KaP4/rbb1tyMp55Q4vCEUsesXh0GoXzggdBVq3lzOOusUGrZaaekncYsNNT/9a8wbx6ce24YWSYnJ2mncA3QsmXh8quyiePbb2HFipLtsrNDoiivtFEbqmI9ocThCaUOmzAhVIcVFISuKwcdBJdfHqrFklS+X748tKfce2+4nOaWW0JTTlZWUg7v6pn162H27FBFNWNGmCCu+P6MGaFatZgEnTtvnDB69AidHWtzVasnlDg8odQDixaFVvWHHw5dWE49taQhJEmmTAndjd96K0xU+cADoVeYa3iWL4+fMGbOLD3QdlZWSBrbbReqpIpvPXqEZU2bpu1hbBJPKHF4QqlHCgvD4F3XXw9bbBF6ih16aNIObxaGbbn88jC48hlnhNO1b5+0U7haoLAwdLktL2F8//3G1962aVM6WcQmj44d62dp1hNKHJ5Q6qHPPoP8/FCsuPTS0PiRxJ+CK1eGQ95xRxgv7IYb4JJLfLbIuuTXX0Njd3kJY+bMMOV0sczM8ksZxbfNN0/Xo0gfTyhxeEKpp1atCheW3HdfGGp48GDo1Supp5g+PZRWXn45dDq7774wTJlLPzNYuDC8RtOnh0RR/HfGjLAuVuvW8UsZnTrVz1LGpvCEEocnlHpuzBg4++zQznLLLXDFFeEnZxK98kooCH3/fZg18t//Dr9oXWoVd7UtThplb7FjtGVkhGtit9uu/JJGbeg5VZfU2oQi6XHgGGCBme0cs/xi4EKgEHjFzK6U1AWYAkyLNhtnZhdE2+8JPAk0BV4FLrUEgvaE0gAsXhy6Zg0fHlrSn346zMSVRKtXh0Ry663hi+6cc0KBqGfPcEti/4AGpago9JoqL2F8/33pq8GzssKFfdtvv/GtSxevkkym2pxQDgRWAE8XJxRJ/YC/AUeb2RpJW5jZgiihvBybeGKOMx64BPiEkFDuM7PXKju/J5QGwiyMEXbRReH/++8P168k+Wr7n36C//f/woX9v/1WsnyrrUqSS8+e4ZKZnj2hXbuknr5OWr8+PG/lJY0ZM0KP8GKNG4cSRnlJY5ttvGqqptTahAJQNlFIehZ4xMzerGi7mOVbA2+b2Q7R/3lAXzP7Y2Xn9oTSwMycGRLJ+++HiVH++19o2zbppykqCj3Bvv46jOwae4utisnJKZ1gim9bbJH0XJdWa9eGp768pPHDD6W72jZtWn7C2H77MIRIkmssXTVsakKp6bzfHThA0t+B1cAVZvZptK6rpM+B5cC1ZvY+0AGYHbP/7GhZuSSdD5wPsE0KBht0tViXLvD223DXXXDttWGu+yeeCDNuJVFGRmhH6dw5zCFWzCxU4cQmmK+/Dn0Gli0r2a5t2/JLNFttVXsSzapVoXF7wYKSv7H3Y5fNnh2SbLGWLUOC6NULTjmldNLYeuva8xhdatR0QskC2gD7AnsBz0raFpgLbGNmi6M2k5GSqjzehpk9AjwCoYSSvLBdnZCZGeZeOfzw0L24f/9w1eI//5nyK82k0GuoU6fSOay4gbk4wRQnm2efDSMgF9t88/JLNMkY/G/NmpAEKkoOsX9jhwqJ1aRJKGHl5IS/O+4Y8nhs0sjJ8aTRkNV0QpkNjIga1cdLKgLamdlCYA2AmU2U9D2hNPMz0DFm/47RMufi69UrDN9y9dVwzz3wxhuhqLDHHjUeihQukGzfvvS1mGYwf/7G1WYvvACPPlqyXatW5ZdoGjWqvARR/Hf58vJja9SoJDnk5ISEEJswyv5t0cKThatYTSeUkUA/4G1J3YHGwCJJOcASMyuMSizdgBlmtkTSckn7EhrlzwLur+GYXV2UnQ133w1HHx26F++zD9x0U2hZrwWV9VKo5tpqKzj44NLrFi4sXW32zTehC/Pjj1d8zMzM0klgr73iJ4ecHNhsM08QLrlS2curAOgLtAPmAzcAzwCPA72AtYQ2lLGSTgZuBtYBRcANZvZSdJzelHQbfg242LsNuypZsgT+9KdQz7TffqFXWNeu6Y6qyhYvDoMEfPNNaLcomyA237x2Dzzoar9a3csrnTyhuFLMYMiQMBx+YWG4/P3ss/0nunMxNjWh+O8Z1zBIYdTHL7+E3r3h978P3YsXLUp3ZM7VG55QXMOyzTZhrPo77wwNE7vsAq9Vep2scy4BnlBcw5OREcb++vTTcEn7UUeFqrDY8T6cc1XmCcU1XLvuGpLKX/4CDz0UuhV7u5tz1eYJxTVs2dnh6vq33gqDdPXpE0aCjB0zxDmXEE8ozkG4GOTLL2HAALjuOjjwwDDsrXMuYZ5QnCvWunXoWjxkSLjYo1cveOyx0OXYOVcpTyjOlZWXB5Mnh0vNzzsPTjxx46kAnXMb8YTiXHk6dYI33wztK6+9FroXv/JKuqNyrlbzhOJcPBkZoQfYhAlhfJNjjgljsj/1FMyZk+7onKt1PKE4V5lddgndi6+6KkzgdfbZYVz5nXeGyy8PJZjYaRyda6B8LC/nqqKoKPQGe+ONcHvvvTDhSOPGYeDJww4Ltz328JEaXZ3jg0PG4QnF1YhVq+CDD0JyGTMGvvgiLG/bFg45pCTBdO6c3jidS4AnlDg8obi0mD8/XCQ5ZkxIMsVtLd27lySXfv3CzFnO1TKeUOLwhOLSzixMYFKcXN55J4wXlpkJ++4bksvhh4fuyVk1PdedcxvzhBKHJxRX66xZAx9/XNL+MmFCSDqbbRau1C8uwWy3nc/T4kozg9WrYcWK0AGk+G/s/eK/WVlw0UXVOo0nlDg8obhab/FiGDu2pP3lxx/D8q5dS5LLwQdDmzbpjdMlrrAw/hd9ecsSXffbb6FDSCLatav2hbieUOLwhOLqFDOYPr2k9DJ2LCxfHnqK9e5dkmD69Ak9ylyJwsJQ+iu+rV5d+v/Klidrn9Wrw60qmjWD5s2hRYvwN/Z+ZX8rWteoUbWeSk8ocXhCcXXa+vUwfnxJ6eWTT8IXZ/Pm0LcvdOsWEkujRqX/lresonWVbZ+ZWbW4i4rCl+qqVSV/K7slul28bdetS85znpEBTZqU3LKzS/9f9lZ2fXZ21b74mzWrdV3LPaHE4QnF1SvLloVG/TFjwpAw8+bB2rXhy7SwMHXnlSpORGvXlv5yX7Om+udq1AiaNg237OyS+5XdsrPL//KvakLwjhGbnFD8GXSuLthsMzj++HArq7AwJJZ160qSzNq1pe8nsq462zdunNgXfiLb+Bd6neevoHN1XWZmuGVnpzsS18DVrgo855xzdZYnFOecc0nhCcU551xSeEJxzjmXFJ5QnHPOJYUnFOecc0nhCcU551xSeEJxzjmXFPV26BVJC4Ef0x3HJmoHLEp3ELWEPxel+fNRmj8fJTbluehsZjnVPXG9TSj1gaQJmzKuTn3iz0Vp/nyU5s9HiXQ+F17l5ZxzLik8oTjnnEsKTyi12yPpDqAW8eeiNH8+SvPno0TangtvQ3HOOZcUXkJxzjmXFJ5QnHPOJYUnlFpGUidJb0v6RtLXki5Nd0y1gaRMSZ9LejndsaSTpM0lDZc0VdIUSX3SHVM6Sbo8+px8JalAUoOaZUzS45IWSPoqZlkbSW9I+i7627qm4vGEUvusB/7PzHoC+wIXSuqZ5phqg0uBKekOoha4F3jdzHYAdqMBPyeSOgCXAL3NbGcgE8hNb1Q17kmgf5llVwFvmVk34K3o/xrhCaWWMbO5ZvZZdP9XwhdGh/RGlV6SOgJHA/9LdyzpJGkz4EDgMQAzW2tmS9MaVPplAU0lZQHNgDlpjqdGmdl7wJIyi48HnoruPwWcUFPxeEKpxSR1AXYHPklzKOl2D3AlUJTmONKtK7AQeCKq/vufpObpDipdzOxn4F/AT8BcYJmZjUlvVLXClmY2N7o/D9iypk7sCaWWktQCeB64zMyWpzuedJF0DLDAzCamO5ZaIAvYA3jYzHYHfqMGqzNqm6ht4HhCom0PNJeUn96oahcL14XU2LUhnlBqIUmNCMlksJmNSHc8abYfcJykmcBQ4GBJg9IbUtrMBmabWXGJdTghwTRUhwI/mNlCM1sHjAB+l+aYaoP5krYGiP4uqKkTe0KpZSSJUEc+xcz+ne540s3MrjazjmbWhdDgOtbMGuSvUDObB8yS1CNadAjwTRpDSrefgH0lNYs+N4fQgDspxBgFDIzuDwRerKkTe0KpffYDziT8Ep8U3Y5Kd1Cu1rgYGCzpS6AX8I/0hpM+UUltOPAZMJnwfdaghmCRVAB8DPSQNFvSucDtwGGSviOU4m6vsXh86BXnnHPJ4CUU55xzSeEJxTnnXFJ4QnHOOZcUnlCcc84lhScU55xzSeEJxaWVpLYx3aPnSfo5ur9UUq27xkJSl9iRXVN4niaS3oyei9OSfOwLJJ1VyTZPSjqlnOV9G/qIzy6+rHQH4Bo2M1tMuJ4CSTcCK8zsX9E4ZvXui0tSlpmtT2DT3QHMrFcKzv+fZB7TuWJeQnG1WaakR6P5LsZIagogaTtJr0uaKOl9STuU3VHSjdFcEe9ImiHpkmh5qRKGpCuiREa07d2SJkRzjewlaUQ0r8StMYfPkjQ42ma4pGbR/ntKejeKa3TM8BfvSLpH0gTCMPyxcbaRNFLSl5LGSdpV0hbAIGCvqISyXcz2O0gaH/N/F0mTo/vXS/o0mhvkkejq8Y3OHz03V0Tr/hDt84Wk54sfS+TQ6Ln4NhpTrexz3Dx6jsdHg1UeHy3fKVo2KXpc3Sp9pV294AnF1WbdgAfNbCdgKXBytPwR4GIz2xO4Angozv47AEcAewM3RGOkVWatmfUG/kMYsuJCYGfgbElto216AA+Z2Y7AcuDP0bHvB06J4noc+HvMcRubWW8zu6vM+W4CPjezXYFrgKfNbAFwHvC+mfUys++LNzazqUBjSV2jRacBw6L7D5jZXtHcIE2B2CQQ7/wjon2K51Y5N2ZdF8JzdzTwH208edXfCEPh7A30A+5UGP34AuDeqHTVmzAGmWsAvMrL1WY/mNmk6P5EoIvCKMy/A56LfoADNImz/ytmtgZYI2kBiQ3jPSr6Oxn4ungYcEkzgE6ExDbLzD6MthtEmOTpdULieSOKK5MwpHqxYZRvf6JEaWZjozalVpXE+Cwhkdwe/S1uY+kn6UrCvCBtgK+Blyo5/85R6WtzoAUwOvY8ZlYEfBc9/rIlwcMJA3deEf2fDWxDGArkbwrz2Iwws+8qeTyunvCE4mqzNTH3Cwm/ujOApQm2LZTdP4swI2Zsybzsr+7ifYrK7F9Eyeel7HhFBoiQgOJNyftbAvEmahghoY4gjFD+XVR6eIgwe+GsqBov9rHFO/+TwAlm9oWks4G+MevKe5yxBJxsZtPKLJ8i6RNCyeZVSX80s7GJPTRXl3mVl6tTorlhfpA0AMLozJJ2q8Ih5gNbRCWBJpSuFkrUNiqZy/104ANgGpBTvFxSI0k7JXCs94Ezon36Aosqm/8mqgIrBK6jpORRnDwWRaW4jXpoxdESmBtV2Z1RZt0ASRlRG862hMcYazRwcUxbze7R322BGWZ2H6HacNcEY3F1nCcUVxedAZwr6QtCtc7xie4YzZtxMzAeeAOYWo3zTwMulDQFaE2Y8Got4Uv8n1Fck0hsbo4bgT0VRg++nZJhxyszDMgnVH8RTQX8KPAV4Yv+0wSPcx1hRtAP2fi5+InwPL0GXGBmq8usvwVoBHwp6evof4BTga8kTSJUAz6dYCyujvPRhp1zziWFl1Ccc84lhScU55xzSeEJxTnnXFJ4QnHOOZcUnlCcc84lhScU55xzSeEJxTnnXFL8fzo0IvFhNCNpAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x_seq,AIC_seq,c=\"red\",label=\"AIC\")\n",
    "plt.plot(x_seq,BIC_seq,c=\"blue\",label=\"BIC\")\n",
    "plt.xlabel(\"The number of variables\")\n",
    "plt.ylabel(\"AIC/BIC values\")\n",
    "plt.title(\"AIC/BIC variables when the number of variables increases.\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "f78qcjRVB-C1"
   },
   "source": [
    "## 5.2  Efficient Estimator and Fisher Information Matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "fwAxJje9B-C2"
   },
   "source": [
    "## 5.3 Kullback-Leibler Divergence"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Lr4o5w7oB-C3"
   },
   "source": [
    "## 5.4 Derivation of Akaike's Information Criterion (AIC)"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "name": "2-4.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}