{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "dfe0f828",
   "metadata": {},
   "source": [
    "### Exercise 32"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f663470",
   "metadata": {},
   "source": [
    "#### Import"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "id": "933c4e44",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy.linalg import inv\n",
    "from scipy.linalg import solve\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42b30c1e",
   "metadata": {},
   "source": [
    "#### Data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "52be1e89",
   "metadata": {},
   "source": [
    "$$f(x) = \\sin(2\\pi x)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "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": "734f8bbb",
   "metadata": {},
   "source": [
    "#### Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "id": "652a38b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "mplot = 200\n",
    "epsilon = .0001\n",
    "space = np.linspace(0, 1, mplot)\n",
    "\n",
    "plt.plot(space, np.sin(2*np.pi*space), linestyle='dashed', alpha=.5, color='black')\n",
    "plt.scatter(x, y, color='black')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "plt.savefig('sinusoidal_figure.pdf', format='pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "161d69ec",
   "metadata": {},
   "source": [
    "#### Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "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"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b5c2919a",
   "metadata": {},
   "source": [
    "#### Frequentist solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "id": "40904e80",
   "metadata": {},
   "outputs": [],
   "source": [
    "ps = [2, 3, 4, n]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "id": "f182977b",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_1880556/400215955.py:8: LinAlgWarning: Ill-conditioned matrix (rcond=2.34821e-17): result may not be accurate.\n",
      "  wMLE = solve(Phi, y)\n"
     ]
    }
   ],
   "source": [
    "fig, axs = plt.subplots(2, 2, figsize=((10, 8)))\n",
    "counter = 0\n",
    "\n",
    "for p in ps:\n",
    "    Phi = create_Phi_poly(x, p, n)\n",
    "    wMLE = inv(Phi.T @ Phi) @ Phi.T @ y\n",
    "    if p == n:\n",
    "        wMLE = solve(Phi, y)\n",
    "        \n",
    "    curve = fitted_curve_poly(space, wMLE, mplot, p)\n",
    "    axs[(counter//2, counter%2)].plot(space, curve, color='black')\n",
    "    axs[(counter//2, counter%2)].scatter(x, y, color='black')\n",
    "    axs[(counter//2, counter%2)].set_ylim(np.min(y)-.2, np.max(y)+.2)\n",
    "    axs[(counter//2, counter%2)].set_title(f'degree = {p-1}')\n",
    "    \n",
    "    counter += 1\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c681fa82",
   "metadata": {},
   "source": [
    "#### Compare MSE and $\\lVert w\\rVert^2 / p$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "id": "312b4bc1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p = 2:\n",
      "MSE = 0.204\n",
      "norm squared = 2.26\n",
      "\n",
      "p = 3:\n",
      "MSE = 0.204\n",
      "norm squared = 1.697\n",
      "\n",
      "p = 4:\n",
      "MSE = 0.025\n",
      "norm squared = 329.972\n",
      "\n",
      "p = 20:\n",
      "MSE = 0.0\n",
      "norm squared = 2.949344409000133e+26\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_1880556/3514712903.py:5: LinAlgWarning: Ill-conditioned matrix (rcond=2.34821e-17): result may not be accurate.\n",
      "  wMLE = solve(Phi, y)\n"
     ]
    }
   ],
   "source": [
    "for p in ps:\n",
    "    Phi = create_Phi_poly(x, p, n)\n",
    "    wMLE = inv(Phi.T @ Phi) @ Phi.T @ y\n",
    "    if p == n:\n",
    "        wMLE = solve(Phi, y)\n",
    "        #for j in range(p):\n",
    "        #    wMLE[j] = wMLE_mp[j]\n",
    "    norm = wMLE.T @ wMLE\n",
    "    curve = fitted_curve_poly(x, wMLE, n, p)\n",
    "    MSE = np.mean((curve - y)**2)\n",
    "    \n",
    "    print(f'p = {p}:')\n",
    "    print(f'MSE = {np.round(MSE, 3)}')\n",
    "    print(f'norm squared = {np.round(norm/p, 3)}')\n",
    "    print('')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc6dfc59",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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
}