{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "dfe0f828",
   "metadata": {},
   "source": [
    "### Exercise 35"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f663470",
   "metadata": {},
   "source": [
    "#### Import"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "933c4e44",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy.linalg import inv, det\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42b30c1e",
   "metadata": {},
   "source": [
    "#### Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f3dad044",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv('sinusoidal.csv')\n",
    "\n",
    "x = df[['x']].to_numpy()\n",
    "y = df[['y']].to_numpy()\n",
    "n = len(x)\n",
    "\n",
    "x = x.reshape((n,))\n",
    "y = y.reshape((n,))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "161d69ec",
   "metadata": {},
   "source": [
    "#### Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "0491c1a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_Phi_poly(x, p, n):\n",
    "    X = np.zeros((n, p))\n",
    "    \n",
    "    for i in range(n):\n",
    "        for j in range(p):\n",
    "            X[i][j] = x[i]**j\n",
    "    return X\n",
    "\n",
    "\n",
    "def create_Phi_gauss(x, p, n):\n",
    "    X = np.zeros((n, p))\n",
    "    centres = np.linspace(0, 1, p)\n",
    "    \n",
    "    for i in range(n):\n",
    "        for j in range(p):\n",
    "            X[i][j] = np.exp(-50*(x[i] - centres[j])**2)\n",
    "    return X\n",
    "\n",
    "\n",
    "def fitted_curve_poly(space, w, mplot, p):\n",
    "    result = np.zeros(mplot)\n",
    "    for j in range(p):\n",
    "        result += w[j] * space**j\n",
    "    return result\n",
    "\n",
    "\n",
    "def fitted_curve_gauss(space, w, mplot, p):\n",
    "    result = np.zeros(mplot)\n",
    "    centres = np.linspace(0, 1, p)\n",
    "    for j in range(p):\n",
    "        result += w[j] * np.exp(-50*(space - centres[j])**2)\n",
    "    return result\n",
    "\n",
    "\n",
    "def get_log_ML(alpha, beta, x, y, n, p):\n",
    "    Phi = create_Phi_poly(x, p, n)\n",
    "    A = alpha*np.eye(p) + beta*Phi.T @ Phi\n",
    "    detA = det(A)\n",
    "    S = inv(A)\n",
    "    m = beta*S @ Phi.T @ y\n",
    "    Em = .5*beta*(y - Phi @ m).T @ (y - Phi @ m) + .5*alpha* m.T @ m\n",
    "    return (p/2)*np.log(alpha) + (n/2)*np.log(beta) - Em - (1/2)*np.log(detA) - (p/2)*np.log(2*np.pi)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20ca04e7",
   "metadata": {},
   "source": [
    "#### Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "c70f08fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "ps = np.arange(1, 10+1)\n",
    "alpha = .0001\n",
    "beta = 10"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "398cba48",
   "metadata": {},
   "source": [
    "#### Calculations and plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "9f1e4de2",
   "metadata": {},
   "outputs": [],
   "source": [
    "log_MLs = np.zeros(len(ps))\n",
    "for j in range(len(ps)):\n",
    "    p = ps[j]\n",
    "    log_MLs[j] = get_log_ML(alpha, beta, x, y, n, p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "8dedad52",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEGCAYAAACO8lkDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAAsTAAALEwEAmpwYAAAm2klEQVR4nO3dd3hVZd7u8e8DSBGQohQREAvSIkQSQwBpwgFBlAFRQTLWGRseYZw5zqCgvAcUh1esCBxwrEF0NsKLDgw4gAgJHSQUASkSepMmQQIhz/kjOzEhENreedbOuj/XtS+SvXfWut2XcGf9VjPWWkRERHIr5jqAiIh4j8pBRETyUTmIiEg+KgcREclH5SAiIvmUcB0gFK666ipbp04d1zFERCLKsmXL9ltrq5zptSJRDnXq1GHp0qWuY0iIDBs2DIABAwY4TiJStBljUs/2WpEoByla6tev7zqCiO+pHMRzunfv7jqCiO9ph7SIiOSjchDPGTJkCEOGDHEdQ8TXNFYSz2nSpInrCCK+p3IQz7n77rtdRxDxPY2VREQkH5WDeM7gwYMZPHiw6xgivqaxknhOXFyc6wgivqdyEM/p0qWL6wgivqexkoiI5KNyEM8ZOHAgAwcOdB1DxNc0VhLPadWqlesIIr6nchDP6dSpk+sIIr6nsZKIiOSjchDPGTBggO7lIOKYxkriOR06dHAdQcT3PFsOxpg7gLeB4sD71trXHEeSQtK+fXvXEUR8z5NjJWNMceA9oDPQEOhtjGnoNpUUloyMDDIyMlzHEPE1T5YDEAdstNZuttaeAD4HujnOJIVk0KBBDBo0yHUMEV/z6ljpGmBbru+3A81yv8EY8zjwOEDt2rULL5mEXceOHV1HEPE9r5bDOVlrxwJjAWJjY63jOBJC7dq1cx1BxPe8OlbaAdTK9X3N4HPiA+np6aSnp7uOIeJrXi2HJUBdY8x1xpiSQC/gK8eZpJDofg4i7nlyrGStzTDGPAPMIOtQ1g+stWscx5JCokt2i7jnyXIAsNZOA6a5ziGFTxfeE3HPq2Ml8bG0tDTS0tJcxxDxNc9uOYh/DR06FIBhw4Y5TiLiXyoH8Zy77rrLdQQR31M5iOe0aNHCdQQR39M+B/GcI0eOcOTIEdcxRHxN5SCeM2zYMO1vEHFMYyXxnO7du7uOIOJ7KgfxnLi4ONcRRHxPYyXxnIMHD3Lw4EHXMUR8TeUgnjN8+HCGDx/uOoaIr2msJJ7Ts2dP1xFEfE/lIJ4TExPjOoKI72msJJ6zf/9+9u/f7zqGiK+pHMRTMjIyiI2NpXHjxrz77rvs3bvXdSQRX1I5iKeMHDmS1NRULrvsMp599llq1KjBnXfeyYQJEzh27JjreCK+oXIQz9i5cycvvfQSnTt3ZsuWLaxcuZK//OUvrFy5kgceeIBq1arx0EMP8Z///IdTp065jitSpBlrresMlyw2NtYuXbrUdQy5RL1792by5MnMmTOHOnXqUL16dQAyMzOZN28eiYmJBAIBDh8+zNVXX03v3r1JSEggOjoaY4zj9CKRxxizzFobe8bXVA7iBbNmzaJDhw68/PLLpKenA2e+n8Px48eZOnUqiYmJTJ06lZMnT9KwYUP69OnDAw88QJ06dQo5uUjkUjmIp504cYImTZpw4sQJVq9ezaZNmwCIiooq8OcOHDhAIBAgMTGRpKQkIOsWowkJCdx7771UqlQp7NlFIllB5aB9DuLcG2+8wbp163j33XcpU6YMUVFR5ywGgMqVK/PEE08wb948Nm/ezNChQ9m3bx9PPPEE1atXp0ePHkyaNClnS0REzp+2HMSprVu30qBBAzp16sSkSZMA2LFjBwDXXHPNBS/PWsvy5ctJTExkwoQJ7Nmzh4oVK3LvvfeSkJDAbbfdRrFi+p1IBDRWEg/r0aMHM2bMYO3atdSuXRuAAQMGAJd+D+mMjAxmzZpFYmIikydPJi0tjdq1a9OnTx8SEhJo2LDhJecXiWQqB/GkadOmceeddzJs2DD+9re/5Ty/du1aABo0aBCydaWlpTFlyhQSExP55ptvOHXqFLfccgsJCQn06tWLGjVqhGxdIpFC5SCe8+uvvxIVFUXJkiVJSUmhZMmShbbuPXv28MUXX5CYmMiSJUsoVqwYt99+OwkJCfTo0YPy5csXWhYRl7RDWjzn73//O5s3b+a9997LVwypqamkpqaGbd3VqlXj2WefZfHixaxbt44XX3yRTZs28fDDD1OtWjV69+6dc5isiF9py0EK3aZNm2jUqBE9evTgs88+y/d6qPY5XAhrLQsWLCAxMZEvvviCAwcOcNVVV9GxY0fatm1LmzZtqFu3rk62kyJFYyXxDGstXbp0ITk5mfXr13P11Vfne8+GDRsAqFu3bmHHA7LOu5gxYwYTJkxg9uzZ7NmzB4Dq1avTpk2bnEeDBg1UFhLRVA7iGZMmTeKee+7hzTffpH///q7jnJO1lh9//JHvvvuOOXPm8N1337Fz504AqlSpQuvWrXPKIioqSofJSkRROYgnpKWl0aBBAypVqsSyZcsoUeLM95ravHkzANdff31hxjsv1lo2bdrEd999l/PYunUrkHVSXqtWrXLKokmTJhQvXtxxYpGzK6gcdCc4KTRDhgxh27ZtTJgw4azFADBu3DigcPc5nC9jDDfeeCM33ngjjz32GABbtmzJUxZTpkwBoEKFCtx22205ZdG0adMC/7tFvERbDlIo1q5dS+PGjUlISODDDz8s8L1e3nI4H9u3b89TFj/++CMA5cqVo2XLljllERsbW6iH8IqcTmMlccpaS/v27VmxYgXr16+nSpUqriMVql27djF37tycsvjhhx8AuPzyy2nevHlOWTRr1oxSpUo5Tit+4rlyMMbcCwwGGgBx1tqluV4bADwGnAKetdbOONfyVA7eNmHCBB544AFGjx7Nk08+ec73uz5aKdz27duXpyxWrlwJQKlSpYiPj885dDY+Pp4yZco4TitFmRfLoQGQCfw/4C/Z5WCMaQhMAOKAGsBM4CZrbYG3/VI5eNfhw4epX78+NWvWZOHChee1g9bFeQ4u/fzzz8ybNy+nLFasWIG1lpIlSxIXF0fr1q1p1aoVzZs3p0KFCq7jShHiuXLIWbkxc8hbDgMArLXDgt/PAAZbaxcUtByVg3f179+fd955h8WLFxMbe8b/B/PJPjv62muvDWc0zzp06BBJSUk5ZbF8+XJOnTqFMYbGjRtz22235Txq1qzpOq5EsEg6WukaYGGu77cHn5MIlJKSwrvvvssTTzxx3sUA/i2FbBUrVqRr16507doVgF9++YVFixaRlJREcnIyH330Ee+99x6Q9VnlLouGDRvqXAsJibCVgzFmJlD9DC+9aK2dEoLlPw48DuRc6lm8IzMzk6effprKlSvzyiuvXNDPhuOqrJGsfPnydOjQgQ4dOgBZlyJPSUkhKSmJpKQkZs2axfjx44GsYmnZsmVOWcTGxlK6dGmX8SVCha0crLUdLuLHdgC1cn1fM/jcmZY/FhgLWWOli1iXhNHHH3/M/Pnz+fDDD6lcufIF/ewnn3wC+Gefw4UqUaIEMTExxMTE0K9fP6y1bN68OacskpKSmDp1KgAlS5YkNjY2pyxatGjBlVde6fi/QCKB1/Y5NAI+47cd0rOAutohHVkOHDhAvXr1qFevHnPnzr3gMcel3AlOsuzfv5/58+fnlMXSpUtzrjLbsGHDPKOoOnXq6BpRPuW5HdLGmO7Au0AV4BCwwlrbKfjai8CjQAbQ31r773MtT+XgLU8++STvv/8+y5cvp3Hjxq7jCFn3z1iyZElOWcyfP5/Dhw8DUKNGDW677baccVTjxo11JrdPeK4cQk3l4B1LliyhWbNm9OvXjzfffPOilrF69WoAoqKiQhlNcsnMzGTNmjV5RlHZ14gqV64czZs3z9myaNasGWXLlnWcWMJB5SCF4tSpUzRr1oydO3eybt06rrjiiotajt/Oc/CKrVu3kpycnFMWq1atwlpL8eLFadq0KS1btsx5nOlS6xJ5VA5SKEaPHs3TTz/NhAkT6NWr10UvZ/fu3UDW/RPEnUOHDrFgwYKcwli0aBHHjx8H4LrrrstTFo0aNdIhtBFI5SBht3fvXurVq0fTpk2ZOXOmdnAWQSdOnOD7778nOTk555F9I6QKFSrQvHnznLKIi4vTKCoCqBwk7B555BHGjx/PypUrqV+//iUta8WKFQBER0dfejAJm+xDaHOXxZo1awAoXrw40dHRebYudPSZ96gcPG7NmjVMmjSJ/v37U758eddxLlhSUhKtWrXib3/7W0j2E2ifQ+Q6ePAgCxYsYP78+SQnJ7No0SJ+/fVXIOts7txlERUVpZshOaZy8LgePXowefJkrrvuOj799FNatmzpOtJ5y8jIoGnTphw+fJgffvghJKOE/fv3A3DVVVdd8rLErZMnT7JixYo8Wxe7du0C4IorriA+Pp4WLVrQsmVLmjVrFpG/HEUylYOHHT16lCpVqtCmTRt+/PFHUlNTef755/mv//qviLgRzJtvvslzzz3HpEmT6N69u+s44nHWWrZs2ZKnLFavXo21lmLFitGkSZM8Wxe1atU690LloqkcPOzzzz+nd+/ezJ07l+joaP70pz/xj3/8g+joaD799FNPH+u/c+dO6tevT6tWrfjXv/4Vsp3Qy5YtAyAmJiYkyxNvO3ToEAsXLiQ5OZn58+ezcOFCjh07BkCtWrXylMXNN9+sE/RCSOXgYT179iQ5OZnt27fnzF+/+uor/vCHP3D48GGGDRtG//79PXmYYO/evZk8eTJr1qzhhhtuCNlytc/B37IvLJh76yL7kirlypWjWbNmxMfHEx8fT7NmzXx3Z8FQKqgcsNZG/CMmJsZGoqNHj9oyZcrYvn375nttz549tlu3bhawbdu2tVu2bHGQ8OxmzpxpATt48OCQL/vAgQP2wIEDIV+uRKbMzEy7ZcsWO378ePv000/bW265xRYvXtwCFrA33HCDfeCBB+w777xjFy9ebNPT011HjhjAUnuWf1e15eDQP//5T+6//37mzJlDmzZt8r1ureXDDz+kX79+FCtWjHfffZff//73zs8hSE9Pp0mTJmRkZLB69WpdEloK3bFjx1i2bBkLFy5k4cKFLFiwIGdHd6lSpYiJicnZsoiPj6dWrVrO/954kcZKHnXfffcxd+5cduzYUeAhfT/99BMPPvggSUlJ3HPPPYwZM8bpkTzDhg3jhRdeYNq0aXTu3Dnky1+8eDEAcXFxIV+2FE3WWrZv386iRYtyCmPZsmU5Z3RfffXVOaOo+Ph4YmJidJIeKgdPOnbsGFWqVOGhhx5i1KhR53z/qVOnGDFiBAMHDuTKK6/kgw8+CMs/zOeSmppKgwYN6Ny5M19++WVY1qF9DhIKJ06cYOXKlXkKY+PGjUDWSXo333xznsKoW7euJ/fthZPKwYMmTpzIvffey+zZs2nXrt15/1xKSgoJCQmsXr2aJ598ktdff71QfwPq3r0733zzDevWrQvbYYZHjhwBuOgL94mczf79+1m0aFFOYSxatCjn/7dKlSoRFxeXUxZxcXEXfKOqSKNy8KD777+fb7/9lp07d17woXnHjx9n0KBBjBgxghtuuIFPP/2U+Pj4MCX9zdSpU+natSuvvfYaf/3rX8O+PpFwy8zMZN26dTlbFosWLWL16tVkZmYCUK9evTxHRhW1Q2lVDh5z7NgxqlatSkJCAmPGjLno5Xz33Xc89NBDbNu2jRdeeIGXXnqJyy67LIRJf/Prr7/SqFEjSpcuzYoVK8J6gt78+fMBaNGiRdjWIXI2v/zyC0uXLs0pjIULF7J3714ALr/8cmJjY3O2LG699daI3tldUDkUnQqMINOnTyctLY177733kpbTpk0bVq5cSb9+/Rg6dCjTpk0jMTGRBg0ahCjpb1577TV++uknZs+eHfYzt7/++mtA5SBulC9fnnbt2uWMe621pKam5imLt956ixMnTgBQtWpVbr311jyPonDuhbYcHOjduzczZ85k165dIdtEnTx5Mo8//jhHjx7l73//O88880zIdq5t3LiRqKgoevTowWeffRaSZRYkLS0NQEeTiGelp6eTkpLCkiVLch5r164l+9/T2rVr5ymLmJgYKlSo4Dh1fhorecivv/5K1apV6d27N2PHjg3psnfv3s0f/vAHpk6dSvv27fnoo4+oWbPmJS3TWkuXLl1ITk5m/fr1ugOYyFkcPXqU5cuX5ymMzZs357xer169PIURHR1NmTJlHCbWWMlTZsyYwdGjRy95pHQm1atX5+uvv2bcuHE899xz3HzzzYwaNYrevXtf9DInT57M9OnTeeuttwqtGObNmwdAq1atCmV9IqFQrlw5WrduTevWrXOe+/nnn1m6dGlOWcyaNYvExEQASpQoQVRUVE5ZxMbGEhUVFbb9hhdKWw6FrE+fPkyfPp3du3eH9X+CjRs38uCDD7JgwQLuv/9+Ro0adcGH5R09epSGDRtSuXJlli5dWmhHaeg8BynKduzYkVMW2cVx8OBBAEqXLk10dHSeLYybbropbOdfaKzkEcePH6dq1arcd999vP/++2FfX0ZGBsOHD+fll1+matWqfPjhh3Ts2PG8f/6vf/0rw4cPJzk5uVB3DqenpwNZl0EQKeps8I56ucdRy5cvz9n3dsUVVxATE5OnMGrXrh2SI6R04T2PmDJligXs9OnTC3W9y5Ytsw0bNrSA7du3r01LSzvnz6xZs8aWKFHCPvLII4WQUERyy8jIsKtWrbIffPCBfeqpp2xsbKy97LLLci42WKVKFdu5c2f70ksv2YULF170eijgwnvO/2EPxSNSyiEhIcFWqlTJnjhxotDXfezYMdu/f38L2JtuuskuXrz4rO/NzMy0bdu2tZUqVbJ79+4txJRZZs+ebWfPnl3o6xXxsuPHj9vFixfb9957zz788MO2UaNGtlixYnbgwIEXvcyCykE7pAtJeno6X331Fffcc4+THU5lypThzTffpGvXrjz88MM0b96cQYMG8cILL+TLM2HCBObMmcOYMWOcHK/9zTffAFzQZUVEirpSpUrljJWyHT16NOd8i5A7W2sU9AD6X8zPhesRCVsOX3/9tQXstGnTXEexBw8etAkJCRawcXFxdv369TmvHTp0yFavXt3eeuutNiMjw0m+kydP2pMnTzpZt4ifUMCWw8XuAn8uVOXkF4FAgIoVK9K+fXvXUahYsSKffvopX3zxBRs2bCA6OppRo0ZhreWll15iz549jBo1qsDLiIdTiRIlitT1a0Qi0cWWQ2ReSMSREydOMGXKFLp16xb2S09ciPvuu4/Vq1fTunVr+vbtS9u2bRk5ciRPPvkksbFnPoChMMyaNYtZs2Y5W7+IXHw5RP7xr4Vo5syZHD58OCwnvl2qGjVq8O9//5v33nuPJUuWcOWVV/LKK684zTRz5kxmzpzpNIOI3531PAdjzC+cuQQMUMZa65ntfq+f5/DII48wefJk9uzZ4+lj97du3UpGRgbXX3+96ygiUggu6vIZ1try4YvkHydOnOB//ud/uPvuuz1dDJB1sTARESigHIwxBV5rwVp7IPRxip7Zs2dz6NAhT46UvGrGjBkAdOrUyXESEf8qaDS0H9gOZAS/z70T2gKaPZyHQCBA+fLlL+iyFX6XfeE9lYOIOwWVwztAOyAZmAAk2bPtoLhAxpj/Bu4CTgCbgEestYeCrw0AHgNOAc9aa2eEYp0unDx5MmJGSl4ydOhQ1xFEfO+sRytZa/sD0UAA+D3wvTFmuDHmuhCs9z9AlLW2MfAjMADAGNMQ6AU0Au4ARhlj3BxsHwLffvstBw4c0EhJRCJOgYeyBk+i+xZ4HhgDPAJ0uNSVWmu/sdZmj6sWAtl3pOkGfG6tTbfW/gRsBOIudX2uBAIBypUrp/HIBZo2bRrTpk1zHUPE185aDsaYssaYB4wxU4BpQDkgxlo7LsQZHgX+Hfz6GmBbrte2B587U77HjTFLjTFL9+3bF+JIly4jI4PJkydz1113Ubp0addxIsrixYtZvHix6xgivlbQPoe9wAbg8+CfFog1xsQCWGsnFbRgY8xMoPoZXnrRWjsl+J4XydrhPf5Cg1trxwJjIes8hwv9+XCbM2cOP//8s0ZKF2Hw4MGuI4j4XkHlECCrEOoFH7lZoMBysNYWOH4yxjwMdAXa59rRvQOolettNYPPRZxAIEDZsmW54447XEcREblgBZ0E93C4VmqMuYOs/RhtrLXHcr30FfCZMeYNoAZQF4i4+UL2SKlr167ObyAeib766isA7r77bsdJRPwrPDcmPbeRQHngP8aYFcaYMQDW2jXAP4EfgOlAX2vtKUcZL9rcuXPZt2+fRkoXKSUlhZSUFNcxRHxN95AOg6eeeopPPvmEffv2cfnll7uOIyJyRgVdW8nVlkORderUKSZNmsSdd96pYhCRiHXOK6saY3qc4enDwCpr7d7QR4ps8+bNY+/evRopXYLJkycD0L17d8dJRPzrfC67/RjQHPg2+H1bYBlwnTHm/1prPw1TtogUCAQoU6YMXbp0cR0lYq1bt851BBHfO59yKAE0sNbuATDGVAM+AZoBcwGVQ1D2SKlLly6ULVvWdZyINWDAANcRRHzvfPY51MouhqC9wecOACfDEysyJScns3v3bo2URCTinc+WwxxjzL/IOikOoGfwubLAoXAFi0SBQIDSpUtz5513uo4S0SZOnAhAz549HScR8a/zKYe+QA/gtuD3HwNfBs9qbheuYJEmMzOTL7/8ks6dO1OuXDnXcSLa5s2bXUcQ8b1zloO11hpjksi694IFFofqvg5Fyfz589m1a5dGSiHw/PPPu44g4nvn3OdgjLmPrEtY9ATuAxYZY7S9f5pAIECpUqXo2rWr6ygiIpfsfMZKLwK3Zp/TYIypAswEJoYzWCTJzMxk4sSJ3HHHHZQvX951nIj3+eefA9CrVy/HSUT863zKodhpJ7v9jM6szmPhwoXs3LlTI6UQ2bEjIi/EK1KknE85TDfGzCDrPtIA95N18x8JCgQClCxZkrvuust1lCLhz3/+s+sIIr53Pjuk/48x5h6gZfCpsdbayeGNFTmyR0qdOnXiiiuucB1HRCQkzmfLAWvtl8CXYc4SkRYvXsz27dt59dVXXUcpMsaPz7oxYJ8+fRwnEfGvs5aDMeYXsg5dzfcSWUe46tdkskZKl112mW5ME0JevCe4iN8UdCc4HXZzDtZaJk6cSMeOHalQoYLrOEVG//79XUcQ8T0ddXQJlixZwtatW3WUkogUOSqHS6CRUnh8/PHHfPzxx65jiPjaee2QlvystQQCATp06EClSpVcxylSfvnlF9cRRHxP5XCRli1bRmpqKi+//LLrKEXOM8884zqCiO9prHSRAoEAJUqUoFu3bq6jiIiEnMrhImSPlNq3b0/lypVdxylyPvjgAz744APXMUR8TeVwEb7//nt++uknHaUUJunp6aSnp7uOIeJr2udwEQKBAMWLF+d3v/ud6yhF0lNPPeU6gojvacvhAmWPlG6//XauvPJK13FERMJC5XCBUlJS2LRpk0ZKYTRu3DjGjRvnOoaIr6kcLlD2SKl79+6uo4iIhI32OVyA7JFS27Ztueqqq1zHKbL++Mc/uo4g4nvacrgAq1atYsOGDRopiUiRp3K4AIFAgGLFimmkFGajR49m9OjRrmOI+JrGSucpe6TUpk0bqlat6jpOkVaqVCnXEUR8T+VwntasWcP69evp16+f6yhF3qOPPuo6gojvORkrGWOGGGNWGmNWGGO+McbUCD5vjDHvGGM2Bl9v6iLfmQQCAYwx9OjRw3UUEZGwc7XP4b+ttY2ttdHAv4CXgs93BuoGH48Dnhk8BwIBWrduTbVq1VxHKfJGjhzJyJEjXccQ8TUn5WCtPZLr27L8dq/qbsAnNstCoKIx5upCD3iaH374gbVr1+oopUJSvnx5ypfXXWpFXHK2z8EY8wrwIHAYaBd8+hpgW663bQ8+t+sMP/84WVsX1K5dO6xZs0dK99xzT1jXI1keeugh1xFEfC9sWw7GmJnGmNVneHQDsNa+aK2tBYwHLvjuLtbasdbaWGttbJUqVUIdP49AIECrVq2oXr16WNcjIuIVYdtysNZ2OM+3jgemAS8DO4BauV6rGXzOmbVr17JmzRreeecdlzF85a233gKgf//+TnOI+Jmro5Xq5vq2G7Au+PVXwIPBo5bigcPW2nwjpcI0ceJEjZQKWZUqVQj31qCIFMzVPofXjDH1gEwgFXgy+Pw0oAuwETgGPOIm3m8CgQAtW7akRo0arqP4Rp8+fVxHEPE9J+VgrT3jr+HWWgv0LeQ4Z7V+/XpWrVqVM+YQEfELXVupABMnTgTQSKmQjRgxghEjRriOIeJrunxGAQKBAC1atKBmzZquo/jKNddc4zqCiO+pHM5iw4YNpKSk8MYbb7iO4ju9evVyHUHE9zRWOovskVLPnj0dJxERKXwqh7MIBALEx8dTq1atc79ZQmr48OEMHz7cdQwRX9NY6Qw2bdrE999/z+uvv+46ii9df/31riOI+J7K4Qw0UnJLn7uIexornUEgECAuLo5rr73WdRQRESdUDqfZvHkzy5Yt02+vDg0bNoxhw4a5jiHiaxornebLL78ENNpwqX79+q4jiPieybpiRWSLjY21S5cuDcmy4uLisNayZMmSkCxPRMSrjDHLrLWxZ3pNY6VctmzZwpIlS7TVICK+p3LIJXukpNuBujVkyBCGDBniOoaIr2mfQy6BQICmTZvqOHvHmjRp4jqCiO+pHIK2bt3KokWLePXVV11H8b27777bdQQR39NYKUgjJRGR36gcggKBANHR0dx4442uo/je4MGDGTx4sOsYIr6msRKwbds2FixYwNChQ11HEbIOJxYRt1QOwKRJkwCNlLyiS5curiOI+J7GSmSNlBo3bsxNN93kOoqIiCf4vhx27NhBcnKyTnzzkIEDBzJw4EDXMUR8zfdjJY2UvKdVq1auI4j4nu/LIRAIEBUVpYu9eUinTp1cRxDxPV+PlXbt2kVSUpJGSiIip/F1OUyfPh1rrUZKHjNgwAAGDBjgOoaIr/l6rPTwww8TFxdHw4YNXUeRXDp06OA6gojv+bocjDE0atTIdQw5Tfv27V1HEPE9X4+VxJsyMjLIyMhwHUPE13y95SDeNGjQIADdR1rEIZWDeE7Hjh1dRxDxPZWDeE67du1cRxDxPe1zEM9JT08nPT3ddQwRX9OWg3hO9r0ctM9BxB2nWw7GmD8bY6wx5qrg98YY844xZqMxZqUxpqnLfOJGly5ddNluEcecbTkYY2oBHYGtuZ7uDNQNPpoBo4N/io/ownsi7rnccngTeB6wuZ7rBnxisywEKhpjrnaSTpxJS0sjLS3NdQwRX3NSDsaYbsAOa23KaS9dA2zL9f324HNnWsbjxpilxpil+/btC1NScWHo0KG6ZauIY2EbKxljZgLVz/DSi8ALZI2ULpq1diwwFiA2Ntae4+0SQe666y7XEUR8L2zlYK0949XTjDE3A9cBKcYYgJrAcmNMHLADqJXr7TWDz4mPtGjRwnUEEd8r9LGStXaVtbaqtbaOtbYOWaOjptba3cBXwIPBo5bigcPW2l2FnVHcOnLkCEeOHHEdQ8TXvHaewzSgC7AROAY84jaOuJB9foPOcxBxx3k5BLcesr+2QF93acQLunfv7jqCiO85LweR08XFxbmOIOJ7uraSeM7Bgwc5ePCg6xgivqZyEM8ZPnw4w4cPdx1DxNc0VhLP6dmzp+sIIr6nchDPiYmJcR1BxPc0VhLP2b9/P/v373cdQ8TXVA7iOSNGjGDEiBGuY4j4msZK4jn333+/6wgivqdyEM+Jjo52HUHE9zRWEs/ZvXs3u3fvdh1DxNdUDuI5b7/9Nm+//bbrGCK+prGSeE6fPn1cRxDxPZWDeE5UVJTrCCK+p7GSeM6OHTvYsUP3eBJxSeUgnjNy5EhGjhzpOoaIr2msJJ7z4IMPuo4g4nsqB/GcBg0auI4g4nsaK4nnpKamkpqa6jqGiK+pHMRzxowZw5gxY1zHEPE1jZXEcx599FHXEUR8T+UgnlO3bl3XEUR8T2Ml8ZzNmzezefNm1zFEfE3lIJ4zbtw4xo0b5zqGiK9prCSe88c//tF1BBHfM9Za1xkumTFmHxDpxz5eBejemL/R55GXPo/f6LPI61I+j2uttVXO9EKRKIeiwBiz1Fob6zqHV+jzyEufx2/0WeQVrs9D+xxERCQflYOIiOSjcvCOsa4DeIw+j7z0efxGn0VeYfk8tM9BRETy0ZaDiIjko3IQEZF8VA6OGWNqGWO+Ncb8YIxZY4zp5zqTa8aY4saY740x/3KdxTVjTEVjzERjzDpjzFpjTHPXmVwyxvwp+PdktTFmgjGmtOtMhckY84ExZq8xZnWu5yobY/5jjNkQ/LNSKNalcnAvA/iztbYhEA/0NcY0dJzJtX7AWtchPOJtYLq1tj7QBB9/LsaYa4BngVhrbRRQHOjlNlWh+wi447Tn/gbMstbWBWYFv79kKgfHrLW7rLXLg1//QtZf/mvcpnLHGFMTuBN433UW14wxFYDWwD8ArLUnrLWHnIZyrwRQxhhTArgc2Ok4T6Gy1s4FDpz2dDfg4+DXHwO/C8W6VA4eYoypA9wCLHIcxaW3gOeBTMc5vOA6YB/wYXDM9r4xpqzrUK5Ya3cArwNbgV3AYWvtN25TeUI1a+2u4Ne7gWqhWKjKwSOMMeWAL4H+1tojrvO4YIzpCuy11i5zncUjSgBNgdHW2luANEI0MohEwVl6N7JKswZQ1hiT4DaVt9iscxNCcn6CysEDjDGXkVUM4621k1zncaglcLcxZgvwOXC7MSbRbSSntgPbrbXZW5ITySoLv+oA/GSt3WetPQlMAlo4zuQFe4wxVwME/9wbioWqHBwzxhiyZsprrbVvuM7jkrV2gLW2prW2Dlk7Gmdba337m6G1djewzRhTL/hUe+AHh5Fc2wrEG2MuD/69aY+Pd9Dn8hXwUPDrh4ApoVioysG9lsDvyfoteUXw0cV1KPGM/w2MN8asBKKBV93GcSe4BTURWA6sIuvfL19dSsMYMwFYANQzxmw3xjwGvAb8L2PMBrK2rl4Lybp0+QwRETmdthxERCQflYOIiOSjchARkXxUDiIiko/KQURE8lE5iIhIPioHERHJR+UgEgbGmDrBezCMD96HYaIx5nLXuUTOl8pBJHzqAaOstQ2AI8DTjvOInDeVg0j4bLPWJge/TgRucxlG5EKoHETC5/Rr0+haNRIxVA4i4VM71z2fHwCSXIYRuRAqB5HwWU/WPcHXApWA0Y7ziJy3Eq4DiBRhGX6+H4VENm05iIhIPrqfg4iI5KMtBxERyUflICIi+agcREQkH5WDiIjko3IQEZF8/j+ARKyIjErt/AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(ps, log_MLs,color='black')\n",
    "plt.xlabel('p')\n",
    "plt.vlines(1 + np.argmax(log_MLs), -100, 100, color='black', linestyle='dotted', alpha=.7)\n",
    "plt.ylim(np.min(log_MLs) - 10, np.max(log_MLs) + 10)\n",
    "plt.ylabel('log ML')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}