HW4_Classification.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4eeb8c94",
   "metadata": {},
   "source": [
    "# Mathematical exercises\n",
    "\n",
    "**Exercise 1** (3 points): Problem 4.10 from the PRML textbook.\n",
    "\n",
    "**Exercise 2** (3 points): Problem 4.16 from the PRML textbook.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c0f2e6a7",
   "metadata": {},
   "source": [
    "# Coding exercises\n",
    "\n",
    "In the exercises below, we will practice using logistic regression for classification, starting with some low-dimensional synthetic data, then moving on to image classification with real data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "b7de1fec",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "00d8270e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'x2')"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAREAAAEGCAYAAABCR6GtAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAfJ0lEQVR4nO2dfZAcdZnHv8/Mbl6QCHshZTDJZuGChBBBsmuMxaWON8+gCHeAp3CoHHIp77BO6zg8EQ897yjxKCwtoU5zSsFJAN+CUBzvkDvjyZLshAQTXiTGLBuMAnFBLEI2M/PcHz2z2zPpnunuX3f/ume+n6rU7uzsTD+98PvO8/Z7fqKqIISQqBRsG0AIyTcUEUKIERQRQogRFBFCiBEUEUKIET22DQjDEUccoQMDA7bNIKTrKJVKL6vqHK/nciUiAwMDGBkZsW0GIV2HiIz6PcdwhhBiBEWEEGIERYQQYgRFhBBiBEWEEGIERYQQYgRFpEsojY7jxvU7UBodt20K6TBy1SdColEaHcdffXsYE+UqpvUUsPbSFRhc2GfbLNIh0BPpAoZ37sVEuYqqAgfKVQzv3GvbJNJBUES6gBVHz8a0ngKKAvT2FLDi6Nm2TSIdhLVwRkRmAPgJgOk1O36oql+wZU8nM7iwD2svXYHhnXux4ujZDGVIrNjMiewHcJqq/kFEegH8VETuU9VhizZ1LIML+ygeJBGsiYg6w13/UHvYW/vHga+E5AyrORERKYrIFgAvAnhIVR/3+J3VIjIiIiMvvfRS6jYSQlpjVURUtaKq7wAwH8ByEVnq8TtrVHVIVYfmzPEcZ0AIsUgmqjOq+gqA9QBWWTaFEBISayIiInNE5PDa9zMBvAfAM7bsIYREw2Z15kgAt4hIEY6YfV9V77FoDyEkAjarM08COMnW9bNKaXSc/RwkV3DvTIbgHheSRzKRWCUO3ONC8ghFJENwjwvJIwxnMgT3uJA8QhHJGNzjQvIGwxlCiBEUEUKIERQRQogRFBFCiBEUEUKIERQRQogRFBFCiBEUEUKIERQRQogRFBFCiBEUEUKIERQRQogRFBFCiBEUEUKIERQRQogRFBFCiBEUEUKIERQRQogRFBFCiBEUEUKIETbP4l0gIutF5CkR2S4in7JlCzmY0ug4bly/A6XRcdumkIxjc9p7GcDlqrpZRGYBKInIQ6r6lEWbCHgSHwmHNU9EVfeo6uba968BeBrAPFv2kCl4Eh8JQyZyIiIyAOdw78c9nlstIiMiMvLSSy+lbls3wpP4SBhEVe0aIHIogP8FcI2qrmv1u0NDQzoyMpKOYV1OaXScJ/GRSUSkpKpDXs9ZPQFPRHoB/AjA2nYCQtKFJ/GRoNiszgiA7wB4WlW/assOQogZNnMiJwP4CIDTRGRL7d/7LNpDCImAtXBGVX8KQGxdnxASD5mozpAuZ2wjsOF65yvJHVYTq4RgbCNwy9lAZQIoTgM+djewYLltq0gI6IkQu+za4AiIVpyvuzbYtoiEhCJC7DKw0vFApOh8HVhp2yISEoYz3czYRueTf2ClvRBiwXInhLFtB4kMRaRbyVIuYsFyikeOYTjTrTAXQWKCItKt2M5FsKzbMTCc6VZs5iKyFEoRYygi3YytXIRXKEURyS0UkU5hbCNe2PIgHqsswVEnnZrtHbj1UKruibCsm2soIp3A2EZUb/4A3lKewPvRg78ufR5XXPrR7AoJy7odBUWkE6iFBz1SBbSMQd2O4Z17sysiAMu6HQSrMzmmPpH9mRknAsVpKGsBB9CDkhzfdqQhp7mTuKAnklPcE9m/0VPAj89ei1m/GcZjlSW4ok1OhNPcSZxQRHJK80T2R/4wgMs+cAbOj/DazIc+JNMwnMkpJhPZOc2dxIn1ae9h4LT3RkwmsnOaOwlDq2nvFJEuh2JCgpDZIyOIXcIkWCk2xA+KSBcTNMFaGh3HBf85jAPlKnp7Crj9b6bEhuJCKCJdTD3BWhcHvwTrus27MVGuAgAmylWs27wbgwv7WComACgiXc3gwj6svXRFW0+iOWtWfxzGk6G30rlQRLqcIMdlnrdsPn44MoYDFUVvUXDesvkAgnky9FY6H9tn8d4E4CwAL6rqUpu2RCILM0pTYHBhH25f/e6DvIkgngwb2zof257IzQBuAPBflu0ITwcN1vEKN257/Hnct20Pzlx6JC58V7+vx9LOkwmadyH5xaqIqOpPRGTApg2RSXOwToIej1e48exvXsPn7vw5AGDDcy8DAC58V3+k9w+adyH5xbYn0hYRWQ1gNQD090f7HzkRTAfrBBWGiB5P0GSmV7gxvHNvw+/ct20Pjp07K7IQBMm7TOL6u5Sqx2B45170HTIN469PUIQySuZFRFXXAFgDOB2rls2ZwmSwThhhiODxhElmeoUbfYdMm/RAAOD4I98cW3K0pbi5/i7VQi+um/gcNpYXoarOye/Te5mYzSKZF5FME3WwThhhiODxhElmeoUb9d/93qbn8ZY3z8Dv95cjJ0fdogGgtRg1/F3gDFfSRQCcsjITs9mEImJKlHxFGGGI4PEELb16CUedY+fOwrO/fQ0/f+FV9BQL6CkIKlVFsVjAC6/sQ2l0vO1ibvaIzl02v7UYuf8uhV6UKsejIEBVne3mTMxmE9sl3tsBnALgCBHZDeALqvodmzaFImqFJqwwhPR42iUzg4Q7bm+mUqniw8udfNQPRsZwx8bnsW7zbs/XucWp2SMSwFPcpl5zDAZrf5fCwEpcwZxILrBdnbnA5vWNManQtBKGGKoxrZKZQcKdZm/m3GXzMbxzL8pV9X1dszhdfdbxB71H/X3qguApaCudex6s3QfJNgxnTEji6IMU+k+C9m6cu2w+pPYVAF54ZR96igVUKt6vaxan8dcnPD0itzDE1YzG1np7UERMSOLogxT6T8KGO8e/9TB86Z7tmChXUSgIls47DB96Z39b78Uv39LuNWFha71dKCKm1Bd4/UBs0wUfwLuJ41M3TLhz37Y9k4+rFcWTu1/Fs7/djmPnzmp4jyiNZWFf43XvbK23C0XElLjDjzbejeenbuG5WL2hZu/gzKVHYtOu32H/gSoUrcutoRrLQr7Gz+Nga71dKCKmJBF+NCddXYnW4Z1/1PCp+6sn1mNw29/FmkPx8g6OnTsL6zbvxg9GxlCpaqTFaupB+XkcbK23C0XElKTPlR25Gbj3ckCrQHE6Tn/vrfiG61P33cWnEsmhNHsH9cf16srph+7C4udvAgrBvJ848hatPI4oHhCJB4qIKUkkV+uex8zZjoBUy87PK/ux+I2tWHvpJZOfuvMKs4GtN6R2OPbgwj4nfLrlosn29HVv/4+pQ8R9ytNR8xbN3gs9juxBEYmDOM+VdedYRIBqZeo5KQADKzG4wP2pm87h2A2L+fmpEK5aVvxq5H58vjQTPz67F4sfuMgztIqSt/DzXige2YIikjXcORYUgEIRqFaBQgF43/XeIpHw4djNi/nHZ5+IxcVpqJYncABFPFY5DgdQxfhTj/qGVkG8iGavg1WXfEARyRrNOZZV1wL79lqdnuZ1ZOfij92NPVsexOUbZ2ErFqG3p4C+JacBYzdFCq28vI6w3kvQxC0b0+KFIpIkUdrXk8ixRLWlhudiXrAI8xYsxxUnTC3IxQv7gLnetrdLrHp5HZeduihwDsT9/j3FAs4fnI/zls2PtG+IhIMikhQm/SNxhydtbGn3ydwqFDkoR+Fje7vQxM/rCJoDcb//RLmK2x/33iTIECl+KCJJkeb4RANbgn4yu7tD3Y+D0i40Me1crb9/u4Y4NqbFD0UkKbz6R2xNh2/RyxLm7BiTMCCISJh2rq69dEXbhjiWieOHIpIUzbkNwN50+BZ5lqCfzJHCgCbRbCcSJnNh6+/tbogzFSsSDIpIkrjzAxuutxve+OQqgn4yhw4DQuaE/LwLL2EJEhpRJNKDIpIWcbbHxxwWBVl0ocOAkDkhL+8C8J7JypAkW1BE0iKu0m1cu4YjCFGoT/iQounlXbQKoehtZIeWIiIibwYwR1V/2fTzE1T1yUQt60TiKN26P+HLbwBbbw//nmmc3hdSNP28C1ZSso+viIjIXwL4GoAXRaQXwMWquqn29M0AliVuXbcQxisYWOnsodEKAAWeuBU48YJwIuARatQPimobHoSxNcKAadMhRyR9WnkinwMwqKp7RGQ5gO+KyJWqeiecs4Q6Dxsl2ChegVanvq+Uwydpm0KNZ2acONXtWRB8cGgBzq11ezYkNgvPJebB+FVmGLZkn1YiUlTVPQCgqhtF5FQA94jIAji9PJ2FrQO6wzal7dqAhj9/oRA+SesKNZ6ZcSKueXIW9h94GQpgoqJY+/jz+N6mMVz6J0fh5sd2TSY2H1lewrwEKkxsRc83hRbPvSYif1x/UBOUUwCcA+D4hO1KH6/FnAZ1r0CKwao2AyuB4nQnpJEi8LZV0a67YDlK/Zfgz+8+gP/b8fJBnwrlqmLNhp3Yf2AqsflYZUk4WwPiV5kh+aCVJ/K3AAoiskRVnwIAVX1NRFYB+HAq1qVJ0hPK/IhykNXH7nYSqk/cCjx7H7DjkUiek3vxFgR4+7zDsO3Xv0el6kiKKlAsCFSd7s+jTjoVGIx/cyBb0fONr4io6lYAEJFtIvJdAP8OYEbt6xCA76ZiYVrEWYKNsnM3yO+63/uw+c7AIp/9MEGSkc2L9+oPHI9nf/Marr5rG6qqkwdQNZ4+F//skiQTqNz2nzxB+kTeBeArAH4GYBaAtQBOjuPiNa/m6wCKAL6tqtfG8b6RMS3BJplXaX7vVdd6ek5h8gtei7c+lDntheeVQDUVAOZa0iGIiBwAsA/ATDieyK9U3eWBaIhIEcCNAN4DYDeATSJydz10yiW7NgCV/U71pLI/3tb25pzNvr2enlPYPS5eizeuioiJCMQhAH5/C3on8RJERDYBuAvAOwEcAeCbInKeqn7Q8NrLAexQ1Z0AICJ3wEna5ldEZs6eKr9q1XkcF145Gw/PKSv5BVMRiGPuh9ffgt5J/AQRkY+r6kjt+z0AzhGRj8Rw7XkAxlyPd8MJnRoQkdUAVgNAf39/DJdNkH174RS8qs7XfTFWGQLmbLLSoGUqAnGIodff4sb1OziUKGbaiohLQNw/Sy2pqqprAKwBgKGhoWz3pwysBHqmJ1fhCZizSatBq1VYYCoCQcUwyFQ2DiVKFlG1sy5F5N0Avqiq7609vhIAVPXLfq8ZGhrSkZGDNC1bxNX1amuAUUCChAV+CzyunETU0IQ5kfCISElVh7yes7mLdxOAY0TkKAAvwOk9udCiPfEQxyY7W92zITA5jMokJ+EWgHWbd0+OQwxjA1vp48WaiKhqWUQ+CeABOCXem1R1uy17MkWW5rP60C4s8BMLk1xJ80T3arU62WlbLDI0sYXVeSKqei+Ae23akElsdc+GoF3Owk8sTHISze9ZRwCcP3jw8RAkHTiUKIt4zWfdcH24/EgKOZVWYUGrIyCiVo/c71ksFgDVyYHM5y2bH8s9kfBYS6xGIReJ1bjxyo8ArQUi7elnPr+XRALT/Z4AmCBNiawmVkk7xjYC//NlVxfshLPxbsvtrQUijpxKUCFq8XtJJzAHF/Y5M0523QUUslnF6gYoIlmlvjjL+wFUna3/xWkAtL1AxJFTCSpECSeBmz2PxoPFe7H4gYsyXcXqBigiWaW+OOvdr0efApxypfPcljtaC0QcO5KDClGCSeDmCs+5y+Y3JFbHn3rUWMDYM2IORSRr1PMLM2c3Ls5TrpxaIEEEwrRfxS1EM2dPDWlqfs8YDyBvXtDN1RhB4+DmviWnAWM3RRYw7qOJB4pIlvDa7r9v78GLs/6938KOi/r7tsuNeAhWXRD6DpnWNI/EG68F3XfINBREADgVmHOXzW843W7xwj5gbnQB4+He8UARsUlzVcNru//Ky71fZ1J9aVV1aWdTgJDBLQhVdfo4pve2/qRvXtA/2rwb6zbvRqWqKBYEV591fMOZM5OE9Ljc3g730cQDRcQWXkLgl19wL2zg4IpNmFxAKwHyem7mbEAE0ELgkMEtCAACtaU394Bsf+FVTJSdjlRVxfjrE8HurwVe3k4WdjznHYqILbw+4VdefnB+wb2wCz1wqjNlNFRswuQCWnkWzc/Vy8ladabKr7o2kFjVBWFyfivQ9pO+3oS2bvNu/GBkDFt3vwrA8WLi8hK8wpfLTl1E8TCEImILP6+j2T1vWNj1Vm9FQ8XG4MyZBgFqfm6ynFwFIIHno7i7UoPmROqvG965FwcqUw2QImgIZcLC8CV5KCK2CFrVcC/suidSrRxcsYnjul7t9u3KyT5EaTQrjY7j16/sQ0EAl45EDmUYvqQDRcQmbq/DL9nptbBNy6mtkpHNz8VUvm2He8EXCoKiOrmQaQYeA8OXdKCIJE2Q/Sftqi3NYpMmBtWPqDNVpar40PJ+zDt8ZqDScFLT1UgwKCJJErQUG7SMGsfGupGbgafvAo47Bxi6OMpd+WLSvNW84M9b1n5rf7vrZWXebKfT6hhNYkrQozmDHqXpfr/yG8DW28LZM3IzcM+ngF8+6nwduTnc69sQ5jjM0ug4bly/A6XRcQBTC/4f/uzYwOIT5HqDC/tw2amLAKDheiQ+6IkkSatKSHOYEzTJWigClQoABZ64DTjxwvZ5lTpP33XwYxNvpOl6QcOH2x5/vuGUvbpoBE3Gurthg1yP7e3JQhFJEj9x8AtL6n0hrQYQzX078ELJ+b5angp9goQ6x53jeCHux1HxuN7gwuWN4UPhOWBD472XRsdx9V3bUK51ok0YjEj0PubzYNjeniwUEVPaffp7JSb9ciCBukn3196kqdGsVV7FbeNZX48nJ+JzvUlvwudehnfunTwwHAAKIkYjEsdfn5gMV/xggjVZKCImRE10+oU5gbpJa52qzY1mrVrmm21sFo8ooxTbjQDwuZcVR8/G9N4CJg44pdwvnbM08UOtmGBNFoqICVEH8viFOWG6SZsbzfzes52NUYWwXR7H515MF3TU17fKt3CmiBkUERNMBvJ4hTlhukmDhk4RPYbI9xDAXtOxiXGOXWTS1RyKiAkxDuRpeM84Bw1F9BhiwXQwUgow6WoORcSUJBeKX64ibA4joseQJs0hRVohBpOu5lgRERH5IIAvAjgOwHKvQ8O7Hr9cRRJHbFr2GLzKtl+6Z3sqZ+wy6WqOLU9kG4BzAXzL0vWzR9CJYjk4YjMszSHFfdv2hA4xTHIbPJvXDCsioqpPA4CI2Lh89ggz5SwHR2yGpTmkOHPpkdi063eBQ4zS6Di+9vAvmNuwROZzIiKyGsBqAOjv77dsTUK0mnK29XYAruEaXqMBwh6xmTG8Qopj584KFGLUPZD9B5xRigWZmoQWNa/Ckm84EhMREXkYwFyPp65S1bs8fu6Jqq4BsAZwjtGMybxsUfcuyvudUV4zXZ+89dPuttzR2B6fVH7Ej4TP9m0OKYKGGPVQqDbrDScvOgKfPuNtABApvGHJNzyJ7eJV1TNUdanHv8ACEhv1/ShhZ3FEfV1YFix35pcWCk5H6v2fnVq0rXYBB90lbHovdbF69Brna9ozTVpQD4WKAkzrLeDTZ7zN88yaVjuK3UR9XTeT+XDGmKif1lEO0jZh315AtXGCe7v8R9D8iKnHkuFkrl91JWrpliXf8Ngq8f4FgG8AmAPgv0Vki6q+N5GLRV0AB00+v61x3mjcoYOXILTr4Qja42EqAhlP5nqFPibt8Sz5hsNWdeZOAHemcrGoC+CgyeeS7KexnyC06+EI0uNhKgIZaUgLS9TSLUu+4RDV/OQqh4aGdGQkQl9a1KRg86FRaSUxkyDhxCjpbESkpKpDns91hYjEBReiNTqp7JrHe2klIp2XWE1yoedgQ1kn0kll1066lzqdNajZtBQZpAyaVtmXTNJJZddOupc6neWJmFQhgpRB02zuyhoWQ7lOKrt20r3U6SwRMalCBBGguPol6gty5mynPyTrORbL4tlJZddOupc6nSUiJqXIIAIUR7+Ee+CyVgEUgJ7p2fZqMtBsFlfZNQtJzU4rIXeWiADRk59Bxw+a9ku4By4DAKqZ6wI9iIw3mwWlE5OaWaDzRMSEIAJkWqGZXJAuTyTrCzOnzWbNJDUKMQvejU0oImnjXpB5yYkAHVHeTiKpSe+GIhKOuCoUHbAg80i7pGYUj4KDnikiwenm8m4H4ZfUjOpRrDh6NnoKggMVRbEQ7jS/TqGzms2SJMzsjjRh81tLSqPjuHH9DpRGx1v+3KgJrD7ms0vHfdITAYKFKVmsUNA7aomfd+H186j5kuGde1GuOJPVKhWGM91J0IWYxQpFBvo3soxfvsLr55eduihSE1gndqCGhSISZiFmLSGaRe8oQ/gtcL+fR2kC68QO1LBwFEDeQwKOJ2iJX8Wl23s7wsJ5In7EtYeFC5l0ON01TyQozR7IqmunKi5hp5/F6clQkEjO6F4RaciF7AfuvdyZth5ECNwLPc7kZt5DK9KVdG+fSD0pKUVACkC1Gvz8Fvfgo5mzp97HNLnp14vCXhCSYbrXE2new3L/Z4NVOZoX+r695qVfd26mudpC74RknO4VEaCxZPuWJcGEwO98mLhCmFXXNiZ5N1zPXpCIsAKTDt0tIm6CCkHcTWdens3Ky6eeZy9IJLK4u7ZTRY0iEoU4m87aiUQWO2VzQNZ212ZR1OLCSmJVRK4TkWdE5EkRuVNEDrdhRyaoi8RpV7Vuua97J0ywBsJ90HcW2tE7ccp7HVueyEMArlTVsoh8BcCVAP7Jki32aefZjG0Ett4OPHErUK0wwRqArLWjd/IeG1tn8T7oejgM4HwbduSCeuK1/AaAWncxE6yByNJA5LRFLc38SxZyIpcA+J7fkyKyGsBqAOjv70/LpuxQT7zWBQTCBGtOSUvU0s6/JJYTEZGHRWSbx79zXL9zFYAygLV+76Oqa1R1SFWH5syZk5S52cXdFFecDgxdzFCGtCTt/EtinoiqntHqeRG5GMBZAE7XPO0CTJs8Vme4/8cqaedfrIQzIrIKwGcA/Kmqvm7DhlyRtTkmrWCHrXXSzr/YyoncAGA6gIfEmUs5rKqfsGQLiRNOW8sEaSaVbVVnFtm4bmbxc//zGBaww7bryEJ1prvxc//zGhbkMYdDjKCI2MbP/c9zWJCnHA4xpnvniWSFhhKuy/33+zkhGYOeSFjizlP4uf8MC0hOoIiEIak8hZ/7z7CA5ACGM2FI4ihNjj4kOYeeSBjiLl/mtQJDiAuKSBiSnmqWpwoMITUoImFJc6oZIQkS17gAiohNWIEhCdPqGNG4xgVQRGzDCgxJiFZCEecMWlZnCOlQWs0ViXMGLT0RQjqUVnNF4hwXIHmaBzQ0NKQjIyO2zSAkN8SVPBWRkqoOeT1HT4SQDiaNuSLMiWQFdq6SnEJPJAuwc5XkGHoiWWDrbc65MnHuySEkJSgithnbCDxxGybPlSkU2blKcgVFxDa7NgDVcu2BACddxFCG5AqKiG3cE8x6ZgAnXmDbIkJCwcSqbbh/huQcikgW4P4ZkmMYzhBCjLAiIiLyryLypIhsEZEHReStNuwghJhjyxO5TlVPUNV3ALgHwNWW7CCEGGJFRFT1966Hb8JkkwQhJG9YS6yKyDUAPgrgVQCntvi91QBWA0B/f386xhFCApPYKAAReRjAXI+nrlLVu1y/dyWAGar6hXbvyVEAhNih1SgA6/NERKQfwL2qujTA774EYDR5qyJzBICXbRuRAJ16X0Dn3lvc97VQVed4PWElnBGRY1T1udrDcwA8E+R1fjeRFURkxE+t80yn3hfQufeW5n3ZyolcKyLHAqjC8Sw+YckOQoghVkREVc+zcV1CSPywYzVe1tg2ICE69b6Azr231O7LemKVEJJv6IkQQoygiBBCjKCIxIyIXCciz9Q2GN4pIofbtskEEVklIs+KyA4R+axte+JARBaIyHoReUpEtovIp2zbFDciUhSRJ0TknqSvRRGJn4cALFXVEwD8AsCVlu2JjIgUAdwI4EwASwBcICJL7FoVC2UAl6vqEgArAFzWIffl5lMAnk7jQhSRmFHVB1W1PjR1GMB8m/YYshzADlXdqaoTAO6A0xyYa1R1j6purn3/GpzFNs+uVfEhIvMBvB/At9O4HkUkWS4BcJ9tIwyYB2DM9Xg3OmixAYCIDAA4CcDjlk2Jk68B+AycZs7E4XjECATZXCgiV8Fxm9emaRsJjogcCuBHAD7dNJ4it4jIWQBeVNWSiJySxjUpIhFQ1TNaPS8iFwM4C8Dpmu9GnBcALHA9nl/7We4RkV44ArJWVdfZtidGTgZwtoi8D8AMAG8WkVtV9aKkLshms5gRkVUAvgrgT1X1Jdv2mCAiPXCSw6fDEY9NAC5U1e1WDTNERATALQB+p6qftmxOYtQ8kX9U1bOSvA5zIvFzA4BZAB6qzZD9pm2DolJLEH8SwANwko/fz7uA1DgZwEcAnFb7b7Sl9slNIkBPhBBiBD0RQogRFBFCiBEUEUKIERQRQogRFBFCiBEUEZIaInK/iLySxs5Skh4UEZIm18HpzyAdBEWExI6IvLM2T2WGiLypNrNjqao+AuA12/aReOHeGRI7qrpJRO4G8G8AZgK4VVW3WTaLJARFhCTFl+DstXkDwN9btoUkCMMZkhSzARwKZx/RDMu2kAShiJCk+BaAf4YzT+Urlm0hCcJwhsSOiHwUwAFVva02p/VnInIagH8BsBjAoSKyG8DHVfUBm7YSc7iLlxBiBMMZQogRFBFCiBEUEUKIERQRQogRFBFCiBEUEUKIERQRQogR/w+iaW1EebNLWwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 288x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Make up some two-dimensional data:\n",
    "nA = 100\n",
    "muA = np.array([1., 0])\n",
    "sigA = np.array([[1,-0.8],[-0.8,1]])\n",
    "xA = np.random.multivariate_normal(muA, sigA, nA)\n",
    "nB = 100\n",
    "muB = np.array([-1., -0.5])\n",
    "sigB = np.array([[1,0.3],[0.3,1]])\n",
    "xB = np.random.multivariate_normal(muB, sigB, nB)\n",
    "\n",
    "# Reformat the data:\n",
    "X = np.vstack((xA, xB))  # X should be n_samples by n_variables, where n_samples=nA+nB\n",
    "y = np.append(np.zeros(nA), np.ones(nB))  # y should have length n_samples\n",
    "\n",
    "# Plot the data:\n",
    "plt.figure(figsize=(4,4))\n",
    "plt.plot(xA[:,0], xA[:,1], '.')\n",
    "plt.plot(xB[:,0], xB[:,1], '.')\n",
    "plt.xlabel('x1')\n",
    "plt.ylabel('x2')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "679bd9ef",
   "metadata": {},
   "source": [
    "## Logistic regression with simulated data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4344d195",
   "metadata": {},
   "source": [
    "**Exercise 1** (5 points): Fit a logistic regression classifier to the synthetic data from above using `LogisticRegression` from the `sklearn` library. Alternatively, you can use `LogisticRegressionCV`, which includes L2 regularization of the weights, with built-in cross-validation to determine the regularization parameter. Report the classification accuracy and plot the result as in the previous exercise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "90e108ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LogisticRegression\n",
    "\n",
    "## Solution ##\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0df4d1c",
   "metadata": {},
   "source": [
    "## Logistic regression on MNIST\n",
    "\n",
    "The MNIST dataset is an extremely common dataset for testing algorithms in machine learning. The dataset consists of a large number of handwritten images, where the images are represented by one-dimensional arrays, where each component of the vector represents the grayscale intensity of one image pixel. The images are 28-by-28 pixels, flattened into arrays of length $784=28^2$. The image labels are integers from 0 to 9, representing the identity of each handwritten digit. Below, we will train a classifier with logistic regression to take one of the images as input and predict the identity of the that image.\n",
    "\n",
    "Note that, because of the large dimension $M=784$ of the data, logistic regression is a better choice than the generative classification approaches implemented above (LDA and QDA). While logistic regression requires us to fit $~M$ parameters, the generative approaches would require $~M^2$ parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "44b56a32",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[5, 0, 4, 1, 9, 2, 1, 3, 1, 4, 3, 5, 3, 6, 1, 7, 2, 8, 6, 9]\n"
     ]
    }
   ],
   "source": [
    "from sklearn.datasets import fetch_openml\n",
    "\n",
    "# Import the data:\n",
    "images0, labels0 = fetch_openml(\"mnist_784\", version=1, return_X_y=True, as_frame=False)\n",
    "labels0 = [int(label) for label in labels0]  # Convert labels from strings to integers\n",
    "\n",
    "# Look at the labels:\n",
    "print(labels0[:20])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "65fe6387",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(70000, 784) (70000,)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x12068dac0>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD4CAYAAAAq5pAIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAOZ0lEQVR4nO3dbYxc5XnG8euKbezamMQbB9chLjjgFAg0Jl0ZEBZQoVCCKgGqArGiyKG0ThOchNaVoLQqtKKVWyVElFIkU1xMxUsgAeEPNAm1ECRqcFlcY2wIb8Y0NmaNWYENIX5Z3/2w42iBnWeXmTMv3vv/k1Yzc+45c24NXD5nznNmHkeEAIx/H+p0AwDag7ADSRB2IAnCDiRB2IEkJrZzY4d5ckzRtHZuEkjlV3pbe2OPR6o1FXbb50m6QdIESf8WEctLz5+iaTrV5zSzSQAFa2NN3VrDh/G2J0i6SdLnJZ0oaZHtExt9PQCt1cxn9gWSXoiIzRGxV9Ldki6opi0AVWsm7EdJ+sWwx1try97F9hLbfbb79mlPE5sD0IyWn42PiBUR0RsRvZM0udWbA1BHM2HfJmnOsMefqC0D0IWaCfvjkubZnmv7MElflLS6mrYAVK3hobeI2G97qaQfaWjobWVEbKqsMwCVamqcPSIelPRgRb0AaCEulwWSIOxAEoQdSIKwA0kQdiAJwg4kQdiBJAg7kARhB5Ig7EAShB1IgrADSRB2IAnCDiRB2IEkCDuQBGEHkiDsQBKEHUiCsANJEHYgCcIOJEHYgSQIO5AEYQeSIOxAEoQdSIKwA0kQdiCJpmZxRffzxPJ/4gkfm9nS7T/7F8fUrQ1OPVBc9+hjdxTrU7/uYv3V6w+rW1vX+73iujsH3y7WT713WbF+3J8/Vqx3QlNht71F0m5Jg5L2R0RvFU0BqF4Ve/bfi4idFbwOgBbiMzuQRLNhD0k/tv2E7SUjPcH2Ett9tvv2aU+TmwPQqGYP4xdGxDbbR0p6yPbPI+LR4U+IiBWSVkjSEe6JJrcHoEFN7dkjYlvtdoek+yUtqKIpANVrOOy2p9mefvC+pHMlbayqMQDVauYwfpak+20ffJ07I+KHlXQ1zkw4YV6xHpMnFeuvnPWRYv2d0+qPCfd8uDxe/JPPlMebO+k/fzm9WP/HfzmvWF978p11ay/te6e47vL+zxXrH//JofeJtOGwR8RmSZ+psBcALcTQG5AEYQeSIOxAEoQdSIKwA0nwFdcKDJ792WL9+ttuKtY/Nan+VzHHs30xWKz/zY1fKdYnvl0e/jr93qV1a9O37S+uO3lneWhuat/aYr0bsWcHkiDsQBKEHUiCsANJEHYgCcIOJEHYgSQYZ6/A5GdfKdaf+NWcYv1Tk/qrbKdSy7afVqxvfqv8U9S3Hfv9urU3D5THyWf9838X66106H2BdXTs2YEkCDuQBGEHkiDsQBKEHUiCsANJEHYgCUe0b0TxCPfEqT6nbdvrFgOXnl6s7zqv/HPPEzYcXqw/+fUbP3BPB12383eK9cfPKo+jD77xZrEep9f/AeIt3yyuqrmLniw/Ae+zNtZoVwyMOJc1e3YgCcIOJEHYgSQIO5AEYQeSIOxAEoQdSIJx9i4wYeZHi/XB1weK9ZfurD9WvunMlcV1F/zDN4r1I2/q3HfK8cE1Nc5ue6XtHbY3DlvWY/sh28/XbmdU2TCA6o3lMP42Se+d9f4qSWsiYp6kNbXHALrYqGGPiEclvfc48gJJq2r3V0m6sNq2AFSt0d+gmxUR22v3X5U0q94TbS+RtESSpmhqg5sD0Kymz8bH0Bm+umf5ImJFRPRGRO8kTW52cwAa1GjY+23PlqTa7Y7qWgLQCo2GfbWkxbX7iyU9UE07AFpl1M/stu+SdLakmba3SrpG0nJJ99i+TNLLki5uZZPj3eDO15taf9+uxud3//SXni7WX7t5QvkFDpTnWEf3GDXsEbGoTomrY4BDCJfLAkkQdiAJwg4kQdiBJAg7kARTNo8DJ1z5XN3apSeXB03+/eg1xfpZX7i8WJ/+vceKdXQP9uxAEoQdSIKwA0kQdiAJwg4kQdiBJAg7kATj7ONAadrk1792QnHd/1v9TrF+1XW3F+t/efFFxXr874fr1ub8/c+K66qNP3OeAXt2IAnCDiRB2IEkCDuQBGEHkiDsQBKEHUiCKZuTG/ij04v1O675drE+d+KUhrf96duXFuvzbtlerO/fvKXhbY9XTU3ZDGB8IOxAEoQdSIKwA0kQdiAJwg4kQdiBJBhnR1GcMb9YP2L51mL9rk/+qOFtH//wHxfrv/239b/HL0mDz29ueNuHqqbG2W2vtL3D9sZhy661vc32+trf+VU2DKB6YzmMv03SeSMs/25EzK/9PVhtWwCqNmrYI+JRSQNt6AVACzVzgm6p7Q21w/wZ9Z5ke4ntPtt9+7Snic0BaEajYb9Z0rGS5kvaLuk79Z4YESsiojcieidpcoObA9CshsIeEf0RMRgRByTdImlBtW0BqFpDYbc9e9jDiyRtrPdcAN1h1HF223dJOlvSTEn9kq6pPZ4vKSRtkfTViCh/+ViMs49HE2YdWay/cslxdWtrr7yhuO6HRtkXfemlc4v1Nxe+XqyPR6Vx9lEniYiIRSMsvrXprgC0FZfLAkkQdiAJwg4kQdiBJAg7kARfcUXH3LO1PGXzVB9WrP8y9hbrf/CNK+q/9v1ri+seqvgpaQCEHciCsANJEHYgCcIOJEHYgSQIO5DEqN96Q24HFs4v1l/8QnnK5pPmb6lbG20cfTQ3DpxSrE99oK+p1x9v2LMDSRB2IAnCDiRB2IEkCDuQBGEHkiDsQBKMs49z7j2pWH/um+Wx7lvOWFWsnzml/J3yZuyJfcX6YwNzyy9wYNRfN0+FPTuQBGEHkiDsQBKEHUiCsANJEHYgCcIOJME4+yFg4tyji/UXL/143dq1l9xdXPcPD9/ZUE9VuLq/t1h/5IbTivUZq8q/O493G3XPbnuO7YdtP217k+1v1Zb32H7I9vO12xmtbxdAo8ZyGL9f0rKIOFHSaZIut32ipKskrYmIeZLW1B4D6FKjhj0itkfEutr93ZKekXSUpAskHbyWcpWkC1vUI4AKfKDP7LaPkXSKpLWSZkXEwYuPX5U0q846SyQtkaQpmtpwowCaM+az8bYPl/QDSVdExK7htRiaHXLEGSIjYkVE9EZE7yRNbqpZAI0bU9htT9JQ0O+IiPtqi/ttz67VZ0va0ZoWAVRh1MN425Z0q6RnIuL6YaXVkhZLWl67faAlHY4DE4/5rWL9zd+dXaxf8nc/LNb/9CP3FeuttGx7eXjsZ/9af3it57b/Ka474wBDa1Uay2f2MyR9WdJTttfXll2toZDfY/sySS9LurglHQKoxKhhj4ifShpxcndJ51TbDoBW4XJZIAnCDiRB2IEkCDuQBGEHkuArrmM0cfZv1q0NrJxWXPdrcx8p1hdN72+opyos3bawWF938/xifeb3NxbrPbsZK+8W7NmBJAg7kARhB5Ig7EAShB1IgrADSRB2IIk04+x7f7/8s8V7/2ygWL/6uAfr1s79jbcb6qkq/YPv1K2duXpZcd3j//rnxXrPG+Vx8gPFKroJe3YgCcIOJEHYgSQIO5AEYQeSIOxAEoQdSCLNOPuWC8v/rj138r0t2/ZNbxxbrN/wyLnFugfr/bjvkOOve6lubV7/2uK6g8UqxhP27EAShB1IgrADSRB2IAnCDiRB2IEkCDuQhCOi/AR7jqTbJc2SFJJWRMQNtq+V9CeSXqs99eqIqP+lb0lHuCdONRO/Aq2yNtZoVwyMeGHGWC6q2S9pWUSssz1d0hO2H6rVvhsR366qUQCtM5b52bdL2l67v9v2M5KOanVjAKr1gT6z2z5G0imSDl6DudT2Btsrbc+os84S2322+/ZpT3PdAmjYmMNu+3BJP5B0RUTsknSzpGMlzdfQnv87I60XESsiojcieidpcvMdA2jImMJue5KGgn5HRNwnSRHRHxGDEXFA0i2SFrSuTQDNGjXsti3pVknPRMT1w5bPHva0iySVp/ME0FFjORt/hqQvS3rK9vrasqslLbI9X0PDcVskfbUF/QGoyFjOxv9U0kjjdsUxdQDdhSvogCQIO5AEYQeSIOxAEoQdSIKwA0kQdiAJwg4kQdiBJAg7kARhB5Ig7EAShB1IgrADSYz6U9KVbsx+TdLLwxbNlLSzbQ18MN3aW7f2JdFbo6rs7eiI+NhIhbaG/X0bt/siordjDRR0a2/d2pdEb41qV28cxgNJEHYgiU6HfUWHt1/Srb11a18SvTWqLb119DM7gPbp9J4dQJsQdiCJjoTd9nm2n7X9gu2rOtFDPba32H7K9nrbfR3uZaXtHbY3DlvWY/sh28/XbkecY69DvV1re1vtvVtv+/wO9TbH9sO2n7a9yfa3ass7+t4V+mrL+9b2z+y2J0h6TtLnJG2V9LikRRHxdFsbqcP2Fkm9EdHxCzBsnynpLUm3R8RJtWX/JGkgIpbX/qGcERFXdklv10p6q9PTeNdmK5o9fJpxSRdK+oo6+N4V+rpYbXjfOrFnXyDphYjYHBF7Jd0t6YIO9NH1IuJRSQPvWXyBpFW1+6s09D9L29XprStExPaIWFe7v1vSwWnGO/reFfpqi06E/ShJvxj2eKu6a773kPRj20/YXtLpZkYwKyK21+6/KmlWJ5sZwajTeLfTe6YZ75r3rpHpz5vFCbr3WxgRn5X0eUmX1w5Xu1IMfQbrprHTMU3j3S4jTDP+a5187xqd/rxZnQj7Nklzhj3+RG1ZV4iIbbXbHZLuV/dNRd1/cAbd2u2ODvfza900jfdI04yrC967Tk5/3omwPy5pnu25tg+T9EVJqzvQx/vYnlY7cSLb0ySdq+6binq1pMW1+4slPdDBXt6lW6bxrjfNuDr83nV8+vOIaPufpPM1dEb+RUl/1Yke6vT1SUlP1v42dbo3SXdp6LBun4bObVwm6aOS1kh6XtJ/Serpot7+Q9JTkjZoKFizO9TbQg0dom+QtL72d36n37tCX21537hcFkiCE3RAEoQdSIKwA0kQdiAJwg4kQdiBJAg7kMT/A65XcTMQuIbWAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# The size of the full dataset:\n",
    "print(np.shape(images0), np.shape(labels0)) \n",
    "\n",
    "# To speed things up, we'll just use a fraction of the data (feel free to use more if you like):\n",
    "n_samples = 5000  \n",
    "images = images0[:n_samples,:]\n",
    "labels = labels0[:n_samples]\n",
    "\n",
    "# Keep some separate data for testing:\n",
    "images_test = images0[n_samples:2*n_samples,:]\n",
    "labels_test = labels0[n_samples:2*n_samples]\n",
    "\n",
    "# To show an example image, reshape the length 784 vector of pixel intensities into a 28x28 matrix:\n",
    "plt.imshow(np.reshape(images[0,:], (28,28)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43100093",
   "metadata": {},
   "source": [
    "In class, we discussed logistic regression with $K\\geq 2$ classes. The probability that an input $\\vec{x}$ belongs to class $C_k$ is given by the softmax function:\n",
    "$$\n",
    "p(C_k | \\vec{x}) = y_k \n",
    "= \\frac{\\exp(\\vec{w}_k \\cdot \\vec{\\phi}(\\vec{x}))}{\\sum_j \\exp(\\vec{w}_j \\cdot \\vec{\\phi}(\\vec{x}))},\n",
    "$$\n",
    "and the target label $\\vec{t}_n$ for input $\\vec{x}_n$ is given by a one-hot vector, i.e. a vector with 1 in the $k$th entry (if the image belongs to class $C_k$) and zeros elsewhere. We showed that the \"cross-entropy\" error to be minimized is given by\n",
    "$$\n",
    "E(\\vec{w}_1, \\ldots, \\vec{w}_K) = -\\sum_{n=1}^N \\sum_{k=1}^K t_{nk} \\ln y_k(\\vec{x}_n).\n",
    "$$\n",
    "Taking the derivative of this error with respect to $\\vec{w}_i$ gave\n",
    "$$\n",
    "\\frac{\\partial E}{\\partial \\vec{w}_i} = \n",
    "-\\sum_n \\sum_k t_{nk} [\\delta_{ki} - y_i(\\vec{x}_n)] \\vec{\\phi}(\\vec{x}_n).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3afa8d6a",
   "metadata": {},
   "source": [
    "**Exercise 2** (5 points): In this exercise we will use the above expressions to implement logistic regression from scratch using only basic Numpy functions. In this exercise, no built-in functions performing logistic regression may be used. First, create one-hot target vectors from the image labels. Next, initialize the weights $\\vec{w}_k$ randomly and use the above expressions to update the weights with gradient descent (or, if you prefer, stochastic gradient descent). You will have to use trial and error to choose a learning rate that is small enough to avoid diverging $\\vec{w}_k$ but large enough to successfully train in a reasonable number of steps. At each iteration of gradient descent, record the error, then make a plot that shows the error decreasing during training. Once the error has decreased to a small value, freeze the weights, then compute and report the percentage of correct classifications for the images in the testing set. You should be able to obtain testing accuracy of around 70-90%.\n",
    "\n",
    "Note (Optional): If you're feeling fancy, feel free to try implementing the iterative reweighted least squares algorithm described in Section 4.3.3 of PRML. This algorithm uses second-derivative information in addition to the first derivative shown above. This means that the algorithm should converge in a smaller number of steps than gradient descent, though each step will be more expensive to compute. Another (also optional) thing to try if you're feeling fancy: Include an L2 regularization term for the weights in the error defined above, and see whether this leads to improved performance compared with the model without regularization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "04ef8fa9",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Solution ##\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cdbc201c",
   "metadata": {},
   "source": [
    "**Exercise 3** (3 points): Repeat the exercise above, now doing things the easy way. Use `LogisticRegression` from the `sklearn` library to perform logistic regression on the MNIST data, then report the percentage of correct classifications using the held-out testing data. If you're feeling fancy (optional) try using L1 and/or L2 regularization or any other tricks you can think of to try to boost the final testing performance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "801d0eae",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LogisticRegression\n",
    "\n",
    "## Solution ##\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "229e8b39",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}