Clean up.
This commit is contained in:
Danilo Enoque Ferreira de Lima committed
1 parent
65893bbabf
commit
1cd365f487
16 files changed
-1254
No files matched your search
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
@@ -1,284 +0,0 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "1bba0128",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"# Bayesian fits\n",
|
||||
"\n",
|
||||
"This is a notebook showing a visualization of a Byesian interpretation for the line fit.\n",
|
||||
"The idea is that we generate all possible line fits and reject them if they are not compatible with the data.\n",
|
||||
"Naturally, this is not how most fits are done, but it illustrates the logic behind the Bayes Theorem when applied for fitting lines.\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
"id": "23feddde",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"from typing import Tuple\n",
|
||||
"\n",
|
||||
"import pandas as pd\n",
|
||||
"import numpy as np\n",
|
||||
"import matplotlib.pyplot as plt\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "bb1286f0",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"We start by generating some fake dataset, which is simple enough that we can visualize the results easily. For this reason, the dataset will contain only two variables.\n",
|
||||
"\n",
|
||||
"The simulated example data will be $f(x) = 2 x + 1 + \\epsilon$, where $\\epsilon \\sim \\mathcal{N}(\\mu=0, \\sigma=0.5)$.\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 8,
|
||||
"id": "aef54435-78af-4b13-aee2-2bc576c773cd",
|
||||
"metadata": {
|
||||
"tags": []
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"true_alpha = 0.3\n",
|
||||
"true_beta = 0.5\n",
|
||||
"true_eps = 0.5\n",
|
||||
"def generate_data(N: int) -> np.ndarray:\n",
|
||||
" x = 2*np.random.randn(N, 1)\n",
|
||||
" epsilon = true_eps*np.random.randn(N, 1)\n",
|
||||
" z = true_beta*x + true_alpha + epsilon\n",
|
||||
" return np.concatenate((x, z), axis=1).astype(np.float32)\n",
|
||||
"\n",
|
||||
"data = generate_data(N=100)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "04432613-39cd-4209-8833-7af585e7eed9",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Given some parameters $\\alpha$ and $\\beta$, we don't know the true $f(x)$, but we assume it is something like $f(x) = \\beta x + \\alpha$."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 9,
|
||||
"id": "a795239a-03eb-483a-8b07-635e3d9dbc55",
|
||||
"metadata": {
|
||||
"tags": []
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def f(x: np.ndarray, alpha: float, beta: float) -> np.ndarray:\n",
|
||||
" return beta*x + alpha"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "a2fba887-dc1a-4400-ab40-59fec452dc83",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Let's plot it for reference:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 48,
|
||||
"id": "fead461d-4d3a-48f2-a525-c518be4f1cf8",
|
||||
"metadata": {
|
||||
"tags": []
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAp4AAAH5CAYAAAAobz7sAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB1XElEQVR4nO3dd3hUZdoG8HsmvQ4EEyZAgFAEQugCUtSAIAFBkKKyIgvroiC6gusnoGICFtauWLCLK7YVRKpBEBBBECSGXgRCCwkEEiYhpM6c74/jhEwy5ZzJnDPt/l1Xrs9M3pl5M7ifN295Ho0gCAKIiIiIiBSmdfcEiIiIiMg/MHgSERERkSoYPImIiIhIFQyeRERERKQKBk8iIiIiUgWDJxERERGpgsGTiIiIiFQR6O4J2GMymXDu3DlERUVBo9G4ezpEREREVIsgCCguLkaTJk2g1dpf0/To4Hnu3DkkJCS4expERERE5MCZM2fQrFkzu2M8OnhGRUUBEH+R6OhoN8+GiIiIiGorKipCQkJCdW6zx6ODp3l7PTo6msGTiIiIyINJORbJy0VEREREpAoGTyIiIiJSBYMnEREREamCwZOIiIiIVMHgSURERESqYPAkIiIiIlUweBIRERGRKhg8iYiIiEgVDJ5EREREpAoGTyIiIiJSBYMnEREREamCwZOIiIiIVMHgCcBoErD9+CWsyMrB9uOXYDQJir/npEmToNFooNFoEBQUhMaNG2Pw4MH45JNPYDKZZL1WSkoKZsyYocxEiYiIiFwk0N0TcLeM/bmYt+ogcg1l1Y/F60KRNiIJqcnxir53amoqPv30UxiNRpw/fx4ZGRl49NFHsXTpUqxcuRKBger+8ezatQsLFy7Eb7/9hnPnzkGv16Nz5864fPkyGjRogIkTJ2LUqFGqzomIiHyP0SRgZ3YBLhSXIS4qFL0SYxCg1bh7WqQCv17xzNifi2lLMi1CJwDkGcowbUkmMvbnKvr+ISEh0Ov1aNq0Kbp3744nn3wSK1aswA8//IDFixeLc8zIQP/+/dGgQQM0atQIw4cPx/Hjx6tfY9KkSfj555/x5ptvVq+gnjx50uHzatu5cyeOHDmC999/HzNmzEBJSQnWrl2L7777Dhs3bsSoUaNw55134u2331b0MyEiIt+WsT8X/V/ciPEf7sCjX2dh/Ic70P/FjYr/N5c8g0YQBOX3lZ1UVFQEnU4Hg8GA6Ohol7620SSg/4sb64ROMw0AvS4UW2cNVORvYZMmTcLly5fx/fff1/lZ165d0aRJE6xduxbLli2DRqNB586dceXKFTzzzDM4efIksrKyoNVqYTAYMHToUCQnJ2P+/PkAgNjYWHz//fd2n2fPvffeiz179mD//v0Wj/ft2xelpaX4448/XPY5EBGR/zAv+NQOHub/yi6a0F3x3UZyPTl5zW+32ndmF9gMnQAgAMg1lGFndgH6tG6k3sQAtG/fHnv37gUAjBkzxuJnn3zyCWJjY3Hw4EEkJydDp9MhODgY4eHh0Ov11eMcPc+WyspKrFmzBtOnT6/zs+bNm2PXrl31+dWIiMhPGU0C5q06WCd0AuJ/czUA5q06iMFJem67+zC/3Wq/UGw7dDozzpUEQYBGI/6P7s8//8T48ePRqlUrREdHo2XLlgCA06dP230NZ5+3adMmGAwGq2c5T58+jVatWsn+fYiIiOQs+JDv8tsVz7ioUJeOc6VDhw4hMTERADBixAi0aNECH374IZo0aQKTyYTk5GRUVFTYfQ1nn7dixQo0a9YMN9xwg8XjJ0+exG+//YaFCxfW75cjIiK/5MkLPqQevw2evRJjEK8LRZ6hzOqyv/mMZ6/EGFXntXHjRuzbtw8zZ87EpUuXcOTIEXz44Ye46aabAABbt26t85zg4GAYjcbq76U+rzZBELBy5Urccccd1SuuAFBRUYEpU6agQ4cO+Oc//1nfX5GIiPyQJy/4kHr8NngGaDVIG5GEaUsyoQEswqc5cqWNSFL0nEl5eTny8vIsyiktWLAAw4cPx8SJE6HRaNCoUSN88MEHiI+Px+nTpzF79uw6r9OyZUv89ttvOHnyJCIjIxETEyPpebVt374dZ8+exYgRIwAAJSUlWL16NebOnYuQkBCsWbMGISEhLv8ciIjI93nqgg+py2/PeAJAanI8Fk3oDr3O8m9Xel2oKjfrMjIyEB8fj5YtWyI1NRWbNm3CwoULsWLFCgQEBECr1eLrr7/G7t27kZycjJkzZ+Lll1+u8zqPP/44AgICkJSUhNjYWJw+fVrS82o6ePAgxo8fDwB44IEHkJycjKSkJHz22WeYM2cOMjMz0aJFC0U+ByIi8n3mBR/g2gKPmVoLPuR+fltOqSYWsiUiIlKHOxu3+IXSUuDAAUCjAa6/HoiKUvwt5eQ1Bk8iIiJSFRd8FPK//8F04AByCktRUlGFsAbRaDZjGgL0jRV9W9bxJCIiIo8VoNWoXiPbp129Crz0Eo5dKMbmI/m4Ul711w/ykJvzPgbMmuIxq8mKnvFctGgROnfujOjoaERHR6NPnz744YcflHxLIiIiIo9hNAnYfvwSVmTlYPvxSzCaXLzRvH9/dehcvTe3RugUVVwuUqUNuFSKrng2a9YM//nPf9C2bVsIgoDPPvsMI0eOxB9//IGOHTsq+dZEREREbqXoeVZBAF5+Gbh6FSZBwOYj+VaHnWwgvo+ndIVSdMVzxIgRGDZsGNq2bYvrr78ezz//PCIjI7Fjxw4l35aIiIhIcfZWM8196Wt3a8ozlNV/BfLkSWDePHGLHUBOYWmdlU6zP5q086iuUKqd8TQajfj2229RUlKCPn36WB1TXl6O8vLy6u+LiorUmh4RERGRZPZWMwcn6ZXrS//aa0BREUyCUH2JqOBK3a6EJxs2wfcdB1g85gldoRQPnvv27UOfPn1QVlaGyMhILF++HElJSVbHLliwAPPmzVN6SkREREROM69m1g6W5tXMGYOul9yXXvIlq6oq4LnnAMDKJSJL3yel4GRM0zqPe0JXKMXLKVVUVOD06dMwGAxYunQpPvroI/z8889Ww6e1Fc+EhASWUyIiIiKPYDQJ6P/iRpvBUgNAFx6Ey1crHb7Wm/d0xciudQNiHVlZwPffA0D1JSJb3r3xLlQEBtWZk14Xiq2zBipyxtOjyikFBwejTZs2AIAePXpg165dePPNN/H+++/XGRsSEsKWjEREROSxdmYXOFzNlBI6AYkrkOnp1f9o7xLR6QZ6fJd8a53HPa0rlOp1PE0mk8WqJhEREZG3kHpOskFYEAyllc73pS8pEW+t12DrEtHG1j2xN/56qy+j97CuUIoGzzlz5mDo0KFo3rw5iouL8eWXX2Lz5s1Yt26dkm9LHm7x4sV46aWXcPLkSbRo0QKvvPIKbr/9dndPi4iIyCGp5yQn92uJNzb8CQ1gET4lrUBu2ABs3Vrn4ZKKuqHz7T53oyrAMs49PKAN2jaO9MiuUIqWU7pw4QImTpyIdu3a4dZbb8WuXbuwbt06DB48WMm39XgajcbuV3qNZXVfs2zZMjz88MOYO3cu9u/fjyFDhmDq1KnunhYREZEkvRJjEK8Lha0op4F4u/3hgW2xaEJ36HWWQVWvC8WiCd1tr0Cmp1sNnQAQEXwtYJYHBuON/vfWCZ0A0K/NdRjZtSn6tG7kUaETYK92t8jLy6v+52+++QbPPPMMjhw5Uv1YZGQkIiMjLZ5TUVGB4OBg1eaolH79+mHQoEHV1QvWr1+PcePG4fLly/V63ZSUFEyaNAmTJk2q/ySJiIjsMN9qB6yvZtYMlpL70ufnA++8Y/d9TYKAT7ZmY0mbm3CyYZM6P1f6EpEtHnW5SDUrVwIXLrh3DnFxwB13OBym1+ur/1mn00Gj0Vg8BohBKjk5GYGBgViyZAk6deqETZs2oWXLlpgxYwZmzJhRPbZr164YNWoU0tPTYTKZ8OKLL+KDDz5AXl4err/+esydOxdjx461OZ9mzZrhySefxEMPPVT92K+//opBgwbh0KFDaNGihYwPwbbi4mLs2LEDr732WvVj69atQ7du3aq//+qrr/CPf/wDJ06cQHy8+D/ayZMnY/fu3fjll1+g0+lcMpedO3fiiSeewG+//YYWLVpgyZIlyMzMxOrVq7Fy5UqXvAcREfmm1OR4LJrQvU4dT2vnKSX1pf/yS+DoUYfvq9Vo0OyN/+DUF384t4XvAXwneF64AJw96+5ZuNRnn32GadOmYdu2bZKfs2DBAixZsgTvvfce2rZtiy1btmDChAmIjY3FLbfcYvU5vXv3xq5du6q/FwQBM2bMwMyZM+uEzhdeeAEvvPCC3TkcPHgQzZs3r/P4nj17oNVq0aVLF1y9ehVffvklFi5ciOXLl1ePueeee/Cf//wHL7zwAt566y2kpaVhw4YN2LFjh8tC544dOzBgwADMnz8fH374IZ544gnMnz8fBw4cwNKlS13yHkRE5NtSk+MxOEkvbTXTFkEQOxBJ0aoVMHEiUgEsmqCRFHo9ke8ETx/Utm1bvPTSS5LHl5eX44UXXsCGDRuqu0O1atUKW7duxfvvv28zeN5444347LPPqr///PPPcebMGcyZM6fO2KlTp+Kuu+6yO48mTeou/wNAVlYW2rdvj927d6N///4AgNGjR2Po0KHVYzQaDZ5//nmMHTsWer0eb731Fn755Rc0bSqhzplEjz32GMaNG4f/+7//AwCMHz8e48ePx8iRIy1WX4mIiOyRtJppy4kTwH//K23sQw+Ju6p/cUnodRMGTw/Wo0cPWeOPHTuGq1ev1rm8VVFRYTdQ3XjjjZg9ezauXLkCjUaDJ598Es8991ydc6YAEBMTg5gYO+Uf7MjKykL37t3RqVMn/Pbbb9i2bRuefvppzJ8/3+JC1fDhw5GUlIT58+fjxx9/RMeOHeu8Vu2V19LSUuzYsQMPP/xw9WPWVl7Pnj2L7du345VXXql+LDAwEIIgsGsWERGLine truncated
|
||||
"text/plain": [
|
||||
"<Figure size 800x600 with 1 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"fig = plt.figure(figsize=(8, 6))\n",
|
||||
"plt.scatter(data[:,0], data[:,1], label=r\"Data $\\mathcal{D}$\")\n",
|
||||
"plt.plot(data[:,0], f(data[:,0], true_alpha, true_beta), c='r', ls='-', alpha=0.5, lw=4, label=r\"True $y=\\beta x + \\alpha$\")\n",
|
||||
"xs = np.linspace(-2, 2, 10)\n",
|
||||
"plt.plot(xs, f(xs, true_alpha, true_beta) - 2*true_eps, c='r', ls='--', alpha=0.5, lw=4)\n",
|
||||
"plt.plot(xs, f(xs, true_alpha, true_beta) + 2*true_eps, c='r', ls='--', alpha=0.5, lw=4)\n",
|
||||
"xa = 2\n",
|
||||
"ya = f(xa, true_alpha, true_beta)\n",
|
||||
"plt.annotate(xy=(xa, ya + 2*true_eps), xytext=(xa, ya - 2*true_eps), arrowprops=dict(arrowstyle='|-|', lw=3, alpha=0.5, color=\"green\"), text=r\"\")\n",
|
||||
"plt.text(xa+0.5, ya, s=r\"$\\epsilon$ noise\", size=16)\n",
|
||||
"plt.legend(frameon=False)\n",
|
||||
"plt.savefig(\"data.png\")"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "00e098e8-0082-4451-927a-073cf242787e",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Let's show several random lines and reject the ones that do not agree with the likelihood:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 49,
|
||||
"id": "08d9daaf-25c5-4d06-b4a6-6cecaaa35651",
|
||||
"metadata": {
|
||||
"tags": []
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"import matplotlib.pylab as pylab\n",
|
||||
"params = {'legend.fontsize': 'x-large',\n",
|
||||
" 'figure.figsize': (15, 5),\n",
|
||||
" 'axes.labelsize': 'x-large',\n",
|
||||
" 'axes.titlesize':'x-large',\n",
|
||||
" 'xtick.labelsize':'x-large',\n",
|
||||
" 'ytick.labelsize':'x-large'}\n",
|
||||
"pylab.rcParams.update(params)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 59,
|
||||
"id": "4c7002b2-5cb5-4929-8e01-964144ee60b6",
|
||||
"metadata": {
|
||||
"tags": []
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def get_log_prior(A, B):\n",
|
||||
" #return 0.25*A/A/1/5\n",
|
||||
" return -0.5*A**2/1**2 - 0.5*B**2/1**2 #- (2*np.pi*p_width**2)\n",
|
||||
"def sample_prior() -> Tuple[float, float]:\n",
|
||||
" #return (2*np.random.rand() - 1)*1, (2*np.random.rand() - 1)*5\n",
|
||||
" return 1*np.random.randn(), 1*np.random.randn()\n",
|
||||
"def get_log_likelihood(x, y, A, B, sigma):\n",
|
||||
" error = f(x, A, B) - y\n",
|
||||
" return np.sum(-0.5*error**2/sigma**2, axis=0) # - (2*np.pi*sigma**2)\n",
|
||||
"def plot(I: int):\n",
|
||||
" # how many lines to show in last column\n",
|
||||
" K = 20\n",
|
||||
" # assumed noise std. dev.\n",
|
||||
" sigma=0.5\n",
|
||||
" # region of alpha and beta to plot\n",
|
||||
" R = 4\n",
|
||||
" # take a subset of data points\n",
|
||||
" D = data[:I,:]\n",
|
||||
" # make plot\n",
|
||||
" fig, ax = plt.subplots(figsize=(4*5, 4), ncols=4)\n",
|
||||
" # create alpha and beta grid\n",
|
||||
" A = np.linspace(-R, R, 50)\n",
|
||||
" B = np.linspace(-R, R, 50)\n",
|
||||
" A, B = np.meshgrid(A, B)\n",
|
||||
" \n",
|
||||
" # calculate the log prior\n",
|
||||
" log_prior = get_log_prior(A, B)\n",
|
||||
" # calculate the log likelihood\n",
|
||||
" log_likelihood = get_log_likelihood(D[:, 0, np.newaxis, np.newaxis],\n",
|
||||
" D[:, 1, np.newaxis, np.newaxis],\n",
|
||||
" A[np.newaxis, ...],\n",
|
||||
" B[np.newaxis, ...],\n",
|
||||
" sigma=sigma)\n",
|
||||
" # the log posterior is, except for a constant factor, just the sum of those\n",
|
||||
" log_posterior = log_prior + log_likelihood\n",
|
||||
"\n",
|
||||
" m=ax[0].imshow(np.exp(log_prior), extent=(-R, R, -R, R), aspect=\"equal\")\n",
|
||||
" #plt.colorbar(m,ax=ax[0], orientation=\"horizontal\")\n",
|
||||
" ax[0].set(xlabel=r\"$\\alpha$\", ylabel=r\"$\\beta$\", title=r\"$p(\\alpha, \\beta)$\")\n",
|
||||
" \n",
|
||||
" m=ax[1].imshow(np.exp(log_likelihood), extent=(-R, R, -R, R), aspect=\"equal\")\n",
|
||||
" #plt.colorbar(m,ax=ax[1], orientation=\"horizontal\")\n",
|
||||
" \n",
|
||||
" ax[1].set(xlabel=r\"$\\alpha$\", ylabel=r\"$\\beta$\", title=r\"$p(\\mathcal{D}|\\alpha,\\beta)$\")\n",
|
||||
" m=ax[2].imshow(np.exp(log_posterior), extent=(-R, R, -R, R), aspect=\"equal\")\n",
|
||||
" #plt.colorbar(m,ax=ax[2], orientation=\"horizontal\")\n",
|
||||
" ax[2].set(xlabel=r\"$\\alpha$\", ylabel=r\"$\\beta$\", title=r\"$p(\\alpha, \\beta|\\mathcal{D})$\")\n",
|
||||
"\n",
|
||||
" if I > 0:\n",
|
||||
" ax[3].scatter(D[:, 0], D[:,1], label=\"Data\")\n",
|
||||
" ax[3].set(xlabel=\"x\", ylabel=\"y\", title=\"Data\", xlim=(-5, 5), ylim=(-5, 5))\n",
|
||||
"\n",
|
||||
" # sample alpha and beta\n",
|
||||
" def get_line():\n",
|
||||
" a, b = sample_prior()\n",
|
||||
" if len(D) == 0:\n",
|
||||
" return a, b\n",
|
||||
" for _ in range(1000):\n",
|
||||
" ap, bp = a + 0.1*np.random.randn(), b + 0.1*np.random.randn()\n",
|
||||
" A = min(1, np.exp(get_log_likelihood(D[:, 0], D[:, 1], ap, bp, sigma=sigma) - get_log_likelihood(D[:, 0], D[:, 1], a, b, sigma=sigma) ))\n",
|
||||
" if np.random.rand() <= A:\n",
|
||||
" a = ap\n",
|
||||
" b = bp\n",
|
||||
" return a, b\n",
|
||||
" for _ in range(K):\n",
|
||||
" a, b = get_line()\n",
|
||||
" ax[3].axline((0, a), slope=b, lw=2, ls='-', c='r', alpha=0.5)\n",
|
||||
" #plt.legend(frameon=False)\n",
|
||||
" plt.tight_layout()\n",
|
||||
" plt.savefig(f\"bayes_reg_{I}.png\")"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 67,
|
||||
"id": "693827d3-c5ac-4b6d-b019-233799cd5d78",
|
||||
"metadata": {
|
||||
"tags": []
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAB5EAAAGGCAYAAABSR3Y3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAADkBElEQVR4nOzdd3hc5Z02/vtM16h3yZbkXmTcK9hgDNgOSZaQhc1uNmWBd8kCSVgCbPKDZJeyyQtJNtkFNnnT2DRIJYEAyQI22KYYcK/YuMpdvZeZ0ZTz++Pm6GisM0YybpLuz3XNJfuZU54zI/vYuuf7fQzTNE2IiIiIiIiIiIiIiIiIiIgAcJ3vCYiIiIiIiIiIiIiIiIiIyIVDIbKIiIiIiIiIiIiIiIiIiPRQiCwiIiIiIiIiIiIiIiIiIj0UIouIiIiIiIiIiIiIiIiISA+FyCIiIiIiIiIiIiIiIiIi0kMhsoiIiIiIiIiIiIiIiIiI9FCILCIiIiIiIiIiIiIiIiIiPRQii4iIiIiIiIiIiIiIiIhID4XIIiIiIiIiIiIiIiIiIiLSQyGyiIiIiIiIiIiIiIiIiIj0UIgsIiIiIiIiIiIiIiIiIiI9FCKLfAA33HADioqK0NnZec7PHYlE8PDDD2PatGlIS0tDWloaZs6cif/5n/9J2m7Tpk0wDAOPP/74OZ+jiIiIiIiIiIiIiIiIDD4KkUVO04YNG/DEE0/gnnvuQXp6+jk9d3d3N5YvX46vfvWr8Hg8uPXWW3HDDTegqqoKN998M1566aWebefMmYOPf/zj+Ld/+zd0dHSc03mKiIicS2frw10PPPAADMPAoUOHzuhx+0sfHBMRkaFsMHw4Gzh7/x7Q/VtERIYSwzCSHn6/H4WFhZg9ezZuvvlmvPDCC4jH42fkXD//+c9hGAZ+/vOfn5HjiUhfCpFFTtPXvvY1ZGVl4bbbbjvn53700Ufx2muv4ZZbbsHmzZvxX//1X/jhD3+IH/zgBwCAN954I2n7e++9FzU1NXjsscfO+VxFRETOhfP54a6zSR8cExGRoWywfDj7bNL9W0REhqL7778f999/P77yla/gk5/8JHJycvDEE0/gIx/5CC6++GLs3bv3fE9RRPrBc74nIDIY7d27Fy+//DJuvvlmpKWlnfPz//CHP0QwGMR3v/tdGIbRM+7x8I90fn5+0vbz58/H5MmT8aMf/Qj33HMPXC59fkRERIaWU3246/jx4/jJT36CV199Ffv370djYyMikQjGjh2LzMxMhEIhjBw5EgsWLMAnP/lJTJs27TxcgbPeHxz7wQ9+0HPfX7x4MT796U/jjTfewIc+9KGe7e+9914sWLAAjz32GL761a+er2mLiIj0y4Xy4ez+3GNP5YP+W0P3bxERGWoeeOCBPmO1tbW4/fbb8dRTT2Hp0qXYuHEjioqKzv3kRKTflCSJAHjrrbdgGAZuv/12/OY3v8Fll12G7Oxs+P1+zJ8/H6tXr07a/qc//SlM08Tf/d3fpTxmNBrFI488gpkzZyItLQ1lZWW488470d3dja6uLhQXF+PTn/70gOd6+PBhHDx4EFdddVWfT2r//ve/BwBceeWVffb75Cc/iSNHjmDlypUDPqeIiMiFzPpw19/+7d8mfbjLNE185zvfwbe//W0sWbIEjz76KNLT0xEKhXD99ddj37592Lx5M3bs2IF/+Zd/wTPPPIPp06fj7rvvPo9Xk+yDfHAskUic07mKiIgMRKr797ky0HuskzP1bw3dv0VEZDgoLi7Gb3/7WyxZsgRHjx7FQw89lPT8pk2bcMcdd2DGjBnIy8tDIBDAhAkTcPfdd6O5uTlp2yVLluCmm24CANx0001JLbStpSdOnDiBf//3f8eiRYtQUlICn8+HESNG4FOf+hR27dp1Tq5ZZLBTiCwC3qAA4OWXX8ZnP/tZ5OXl4dZbb8XixYuxYcMGfOQjH8GRI0d6tn/55Zfhdrtx8cUXOx6vqakJixYtwp133ony8nL88z//M7Kzs/HII4/gscceww9/+EM0NTXhwQcfHPBcN27cCABYsGBBz5hpmnj00Ufxxz/+EUuXLsX06dP77Ldo0SIAUIgsIiIXvDP14a6amhp86EMfwqOPPoolS5Zg+vTpPT8Y/sxnPtOzncfjwdVXX43Vq1cjOzsb//mf/zmg+6U+OCYiIjI8Ppx9sjP5bw3dv0VEZDhwuVz413/9VwDAb37zG5im2fPcT37yE/z2t7/FpEmTcNNNN+G2225DaWkp/vM//xOLFi1Ce3t7z7Y33ngjrr32WgDAtdde29M++/7770dOTg4A4LXXXsM3v/lN5OTk4Prrr8edd96Jiy++GH/4wx8wf/58bNu27dxduMhgZYqIedNNN5kAzKysLPP1119Peu722283AZh33323aZqm2dHRYbrdbnPq1Kkpj7d06VITgPnoo4/2jNXX15t+v9+88sorzZKSEvPmm28+rbnee++9JgDzhRdeMFetWmXefPPN5qRJk0wA5owZM8yamhrH/VpaWkwA5rx5807rvCIiIufKf//3f5sAzMmTJ5tut9v82Mc+Zn7lK1/pub8GAgHz8OHDPdvPmTPHdLvdZkdHxymPe+DAgZ79Ozs7Hbf5x3/8RxOA+bnPfa5n7P777zcBmFVVVX22b2xsNOfNm2cCMP/qr/7K/MpXvmJOmTLFBGD+x3/8h/nd737X9Hg85r59+wb8OvzhD38wAZjf+MY3esYSiYT5yCOPmADMpUuXOu63cuXKpH+7iIiInAtn+v59od1jT/XvAcvp/lvDNHX/FhGRoQGA+X6xUzgcNj0ejwnAPHjwYM/4oUOHzFgs1mf7xx9/3ARgfvOb30wa/9nPfmYCMH/2s585nqe2ttZsa2vrM75161YzPT3dvPrqq/txRSLDmyqRRQBs3rwZAPDf//3fuPTSS5Oeu/nmmwGgp8XF8ePHEY/HUVpa6nisl19+GS+//DIuu+wy3H777T3jBQUFGD16NFatWoXm5mbcd999pzVXq2p67ty5+NGPfoTHH38ce/bsAQBMnjwZ8Xjccb/s7GwEAoGkimoREZELkXVfPnHiBNasWYNnn30W3/rWt7By5UrcfvvtCIfDeOyxxwAAnZ2d2Lp1KyorK/tUEp3sueeeAwBcddVVCAaDjtuMHTsWAHD06NF+zfXv/u7vsGHDBjz66KN4/vnn8a1vfQuvvvoq/H4/XnjhBfzHf/wHbrzxRowfP75fx+vNuufPmTMHq1evxuc+9zlUVlbiS1/6EmbMmIEnn3zScb958+YB4KeuRUREzpUzff++EO+x7+eD/FtD928RERku/H5/z7IR9fX1PeOjRo2C2+3us/3/+T//B1lZWXjppZcGdJ6ioiJkZmb2GZ8xYwauvPJKrF69GtFodICzFxleFCLLsBeJRLBr1y6Ul5cntZuyWDc064bS2NgIAMjNzXU83hNPPAEA+NKXvpS0rhIABAIBAMAtt9yC8vLy05rvpk2bMGrUKBQUFOBXv/oVGhoasGbNGlx//fX43e9+h6uvvjrlvnl5eWhoaDit84qIiJwrZ/LDXb09//zzAIC/+qu/SrlNR0cHAMDn873v8fTBMREREdtw+HD2+/kg/9bQ/VtERIYT87021r1/fh6NRvG9730Pl156KfLy8uB2u2EYBlwuF9ra2nD8+PEBn+cvf/kLrrnmGpSWlsLr9fasm/z8888jEonoZ+Ui78Nzvicgcr7t2LED0WgUH/rQh+By9f1cxaFDhwAAFRUVAIC0tDQAQDgcdjzeq6++Cq/XmzLMDQaD+OpXv3pacz18+DAaGxuxZMkSAIDb7UZ+fj4uv/xyXH755Zg5cya2bduGgwcP9ny6ubdQKNQzfxERkQvRmf5wl6WtrQ2vv/46gFP/YPfw4cMAgIsuuuh953quPzj2/e9/Hzt37sR///d/43e/+x127dqF7du3O+6bl5eH2tra0zqviIjIQA32D2cP5B6bypn4t4bu3yIiMhyEw2E0NTUBAAoLC3vG/+7v/g7PPPMMxo4di2uvvRYlJSXw+/0AgEceeQSRSGRA53n00UfxpS99Cbm5uVi2bBkqKioQDAZhGAb+9Kc/Ydu2bQM+pshwoxBZhj3r09KjR492fP4vf/kLAGD58uUA2AYDsP/T21soFMKRI0cwbty4Pq2rDh48iHfffRcLFy5EcXHxac2196elnVj/AXdq05FIJNDS0oIxY8ac1rlFRETOhTP94S7Liy++iGg0ilmzZqGsrMxxm3g83tMe66Mf/ej7zlUfHBMREaHh9OHsVM7EvzV0/xYRkeHgjTfeQCwWQ3Fxcc/P5Ddu3IhnnnkGS5cuxQsvvACPx46uEokEvv3tbw/oHLFYDA888ABKSkqwefPmPt1P3nrrrQ98HSLDgdpZy7BnhcgtLS19nmtqasKPf/xjlJeX49prrwUAlJaWorCwsKfVVW+hUAimaTr+p/nOO+9EJBJJugEOVO91m5zmunbtWkybNi3pE1yWPXv2wDRNzJw587TPLyIicradyQ939fbnP/8ZAHDNNdek3OZnP/sZGhsbsWDBAixatOiUx7M+ODZq1KiUHxxbsGDBef3gmPXaiIiInG1n48PZF+I99lQ+6L81dP8WEZHhIJFI4P/+3/8LAPjUpz7VM75//34AwMc+9rE+Pz9fv349QqFQn2NZ6yc7LUPR0NCAlpYWLFy4sE+A3NHR0fNvFxE5NYXIMuxZ/4H8wx/+gM7Ozp7xjo4OfOpTn0JjYyMeeeSRnpZZhmFg8eLFaGho6Lm5WXJzc5GRkYH9+/cntb76wQ9+gOeeew6Ac1gNADfeeCMMw8DPf/7z953r7373u551IwCgu7sbt9xyC6LRKO6++27Hfd9++20AwBVXXJHy+CIiIufbmfxwlyWRSOB///d/AaT+we7+/ftx7733wu/34wc/+MH7zlMfHBMREbENlw9np3Im/q2h+7eIiAx1dXVLine truncated
|
||||
"text/plain": [
|
||||
"<Figure size 2000x400 with 4 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"plot(I=7)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"id": "a1fbd00b-e7f2-4eff-90db-e1d44d104fba",
|
||||
"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.6"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 5
|
||||
}
|
||||
@@ -1,667 +0,0 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "6a764d92",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"# Gaussian Processes\n",
|
||||
"\n",
|
||||
"If one is ready to assume the data available is such that the variable being fit is Gaussian for any pair of samples taken, one can use fit the data and obtain a reliable uncertainty with a very small set of assumptions.\n",
|
||||
"\n",
|
||||
"For this example, we will use the dataset collected in Mauna Loa showing the amount of CO2 present in the atmosphere. This data may be downloaded from https://gml.noaa.gov/ccgg/trends/data.html and it should already be available in this repository.\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "c8131588",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"We can take a quick look at the file to see what is inside:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 1,
|
||||
"id": "07ecc0d3",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"# --------------------------------------------------------------------\n",
|
||||
"# USE OF NOAA GML DATA\n",
|
||||
"# \n",
|
||||
"# These data are made freely available to the public and the\n",
|
||||
"# scientific community in the belief that their wide dissemination\n",
|
||||
"# will lead to greater understanding and new scientific insights.\n",
|
||||
"# The availability of these data does not constitute publication\n",
|
||||
"# of the data. NOAA relies on the ethics and integrity of the user to\n",
|
||||
"# ensure that GML receives fair credit for their work. If the data \n",
|
||||
"# are obtained for potential use in a publication or presentation, \n",
|
||||
"# GML should be informed at the outset of the nature of this work. \n",
|
||||
"# If the GML data are essential to the work, or if an important \n",
|
||||
"# result or conclusion depends on the GML data, co-authorship\n",
|
||||
"# may be appropriate. This should be discussed at an early stage in\n",
|
||||
"# the work. Manuscripts using the GML data should be sent to GML\n",
|
||||
"# for review before they are submitted for publication so we can\n",
|
||||
"# ensure that the quality and limitations of the data are accurately\n",
|
||||
"# represented.\n",
|
||||
"# \n",
|
||||
"# Contact: Pieter Tans (303 497 6678; pieter.tans@noaa.gov)\n",
|
||||
"# \n",
|
||||
"# File Creation: Fri Jan 21 05:00:12 2022\n",
|
||||
"# \n",
|
||||
"# RECIPROCITY\n",
|
||||
"# \n",
|
||||
"# Use of these data implies an agreement to reciprocate.\n",
|
||||
"# Laboratories making similar measurements agree to make their\n",
|
||||
"# own data available to the general public and to the scientific\n",
|
||||
"# community in an equally complete and easily accessible form.\n",
|
||||
"# Modelers are encouraged to make available to the community,\n",
|
||||
"# upon request, their own tools used in the interpretation\n",
|
||||
"# of the GML data, namely well documented model code, transport\n",
|
||||
"# fields, and additional information necessary for other\n",
|
||||
"# scientists to repeat the work and to run modified versions.\n",
|
||||
"# Model availability includes collaborative support for new\n",
|
||||
"# users of the models.\n",
|
||||
"# --------------------------------------------------------------------\n",
|
||||
"# \n",
|
||||
"# \n",
|
||||
"# See www.esrl.noaa.gov/gmd/ccgg/trends/ for additional details.\n",
|
||||
"# \n",
|
||||
"# NOTE: DATA FOR THE LAST SEVERAL MONTHS ARE PRELIMINARY, ARE STILL SUBJECT\n",
|
||||
"# TO QUALITY CONTROL PROCEDURES.\n",
|
||||
"# NOTE: The week \"1 yr ago\" is exactly 365 days ago, and thus does not run from\n",
|
||||
"# Sunday through Saturday. 365 also ignores the possibility of a leap year.\n",
|
||||
"# The week \"10 yr ago\" is exactly 10*365 days +3 days (for leap years) ago.\n",
|
||||
"# \n",
|
||||
"# Start of week CO2 molfrac (-999.99 = no data) increase\n",
|
||||
"# (yr, mon, day, decimal) (ppm) #days 1 yr ago 10 yr ago since 1800\n",
|
||||
" 1974 5 19 1974.3795 333.37 5 -999.99 -999.99 50.40\n",
|
||||
" 1974 5 26 1974.3986 332.95 6 -999.99 -999.99 50.06\n",
|
||||
" 1974 6 2 1974.4178 332.35 5 -999.99 -999.99 49.60\n",
|
||||
" 1974 6 9 1974.4370 332.20 7 -999.99 -999.99 49.65\n",
|
||||
" 1974 6 16 1974.4562 332.37 7 -999.99 -999.99 50.06\n",
|
||||
" 1974 6 23 1974.4753 331.73 5 -999.99 -999.99 49.72\n",
|
||||
" 1974 6 30 1974.4945 331.68 6 -999.99 -999.99 50.02\n",
|
||||
" 1974 7 7 1974.5137 331.46 6 -999.99 -999.99 50.20\n",
|
||||
" 1974 7 14 1974.5329 330.83 5 -999.99 -999.99 50.01\n",
|
||||
" 1974 7 21 1974.5521 330.76 7 -999.99 -999.99 50.41\n",
|
||||
" 1974 7 28 1974.5712 329.80 4 -999.99 -999.99 49.97\n",
|
||||
" 1974 8 4 1974.5904 329.85 5 -999.99 -999.99 50.54\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"with open(\"co2_weekly_mlo.txt\", \"r\") as f:\n",
|
||||
" for i, line in enumerate(f.readlines()):\n",
|
||||
" if i > 60: # let's print only the first 60 lines only\n",
|
||||
" break\n",
|
||||
" print(line[:-1])"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "fce4d8e8",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"We start by loading the necessary Python modules."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
"id": "300cf8d3",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"import matplotlib.pyplot as plt\n",
|
||||
"\n",
|
||||
"import pandas as pd\n",
|
||||
"import numpy as np\n",
|
||||
"from sklearn.gaussian_process import GaussianProcessRegressor\n",
|
||||
"from sklearn.gaussian_process.kernels import RBF, ExpSineSquared, WhiteKernel, RationalQuadratic"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "0ecd6a69",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Let us now load the data downloaded using pandas. If you are not familiar with Pandas, take a look at the documentation for a quick introduction at https://pandas.pydata.org/docs/user_guide/10min.html\n",
|
||||
"\n",
|
||||
"We can easily parse the data file using Pandas' read_csv function. It would require some coding to do it manually. As Pandas does not detect the column names, we have to specify it manually in the names parameter below."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 3,
|
||||
"id": "4959a292",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"data = pd.read_csv(\"co2_weekly_mlo.txt\",\n",
|
||||
" delim_whitespace=True,\n",
|
||||
" comment='#',\n",
|
||||
" header=None,\n",
|
||||
" names=[\"year\", \"month\", \"day\", \"decimal\",\n",
|
||||
" \"co2\",\n",
|
||||
" \"ndays\",\n",
|
||||
" \"last_year\", \"last_decade\",\n",
|
||||
" \"increase\"])"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "d8295e8a",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Let's print out the dataset read first."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 4,
|
||||
"id": "024fb65a",
|
||||
"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>year</th>\n",
|
||||
" <th>month</th>\n",
|
||||
" <th>day</th>\n",
|
||||
" <th>decimal</th>\n",
|
||||
" <th>co2</th>\n",
|
||||
" <th>ndays</th>\n",
|
||||
" <th>last_year</th>\n",
|
||||
" <th>last_decade</th>\n",
|
||||
" <th>increase</th>\n",
|
||||
" </tr>\n",
|
||||
" </thead>\n",
|
||||
" <tbody>\n",
|
||||
" <tr>\n",
|
||||
" <th>0</th>\n",
|
||||
" <td>1974</td>\n",
|
||||
" <td>5</td>\n",
|
||||
" <td>19</td>\n",
|
||||
" <td>1974.3795</td>\n",
|
||||
" <td>333.37</td>\n",
|
||||
" <td>5</td>\n",
|
||||
" <td>-999.99</td>\n",
|
||||
" <td>-999.99</td>\n",
|
||||
" <td>50.40</td>\n",
|
||||
" </tr>\n",
|
||||
" <tr>\n",
|
||||
" <th>1</th>\n",
|
||||
" <td>1974</td>\n",
|
||||
" <td>5</td>\n",
|
||||
" <td>26</td>\n",
|
||||
" <td>1974.3986</td>\n",
|
||||
" <td>332.95</td>\n",
|
||||
" <td>6</td>\n",
|
||||
" <td>-999.99</td>\n",
|
||||
" <td>-999.99</td>\n",
|
||||
" <td>50.06</td>\n",
|
||||
" </tr>\n",
|
||||
" <tr>\n",
|
||||
" <th>2</th>\n",
|
||||
" <td>1974</td>\n",
|
||||
" <td>6</td>\n",
|
||||
" <td>2</td>\n",
|
||||
" <td>1974.4178</td>\n",
|
||||
" <td>332.35</td>\n",
|
||||
" <td>5</td>\n",
|
||||
" <td>-999.99</td>\n",
|
||||
" <td>-999.99</td>\n",
|
||||
" <td>49.60</td>\n",
|
||||
" </tr>\n",
|
||||
" <tr>\n",
|
||||
" <th>3</th>\n",
|
||||
" <td>1974</td>\n",
|
||||
" <td>6</td>\n",
|
||||
" <td>9</td>\n",
|
||||
" <td>1974.4370</td>\n",
|
||||
" <td>332.20</td>\n",
|
||||
" <td>7</td>\n",
|
||||
" <td>-999.99</td>\n",
|
||||
" <td>-999.99</td>\n",
|
||||
" <td>49.65</td>\n",
|
||||
" </tr>\n",
|
||||
" <tr>\n",
|
||||
" <th>4</th>\n",
|
||||
" <td>1974</td>\n",
|
||||
" <td>6</td>\n",
|
||||
" <td>16</td>\n",
|
||||
" <td>1974.4562</td>\n",
|
||||
" <td>332.37</td>\n",
|
||||
" <td>7</td>\n",
|
||||
" <td>-999.99</td>\n",
|
||||
" <td>-999.99</td>\n",
|
||||
" <td>50.06</td>\n",
|
||||
" </tr>\n",
|
||||
" <tr>\n",
|
||||
" <th>...</th>\n",
|
||||
" <td>...</td>\n",
|
||||
" <td>...</td>\n",
|
||||
" <td>...</td>\n",
|
||||
" <td>...</td>\n",
|
||||
" <td>...</td>\n",
|
||||
" <td>...</td>\n",
|
||||
" <td>...</td>\n",
|
||||
" <td>...</td>\n",
|
||||
" <td>...</td>\n",
|
||||
" </tr>\n",
|
||||
" <tr>\n",
|
||||
" <th>2482</th>\n",
|
||||
" <td>2021</td>\n",
|
||||
" <td>12</td>\n",
|
||||
" <td>12</td>\n",
|
||||
" <td>2021.9466</td>\n",
|
||||
" <td>416.45</td>\n",
|
||||
" <td>7</td>\n",
|
||||
" <td>414.00</td>\n",
|
||||
" <td>392.19</td>\n",
|
||||
" <td>137.16</td>\n",
|
||||
" </tr>\n",
|
||||
" <tr>\n",
|
||||
" <th>2483</th>\n",
|
||||
" <td>2021</td>\n",
|
||||
" <td>12</td>\n",
|
||||
" <td>19</td>\n",
|
||||
" <td>2021.9658</td>\n",
|
||||
" <td>417.49</td>\n",
|
||||
" <td>6</td>\n",
|
||||
" <td>414.86</td>\n",
|
||||
" <td>392.54</td>\n",
|
||||
" <td>137.93</td>\n",
|
||||
" </tr>\n",
|
||||
" <tr>\n",
|
||||
" <th>2484</th>\n",
|
||||
" <td>2021</td>\n",
|
||||
" <td>12</td>\n",
|
||||
" <td>26</td>\n",
|
||||
" <td>2021.9849</td>\n",
|
||||
" <td>417.46</td>\n",
|
||||
" <td>7</td>\n",
|
||||
" <td>415.32</td>\n",
|
||||
" <td>392.43</td>\n",
|
||||
" <td>137.66</td>\n",
|
||||
" </tr>\n",
|
||||
" <tr>\n",
|
||||
" <th>2485</th>\n",
|
||||
" <td>2022</td>\n",
|
||||
" <td>1</td>\n",
|
||||
" <td>2</td>\n",
|
||||
" <td>2022.0041</td>\n",
|
||||
" <td>417.43</td>\n",
|
||||
" <td>6</td>\n",
|
||||
" <td>415.44</td>\n",
|
||||
" <td>393.18</td>\n",
|
||||
" <td>137.40</td>\n",
|
||||
" </tr>\n",
|
||||
" <tr>\n",
|
||||
" <th>2486</th>\n",
|
||||
" <td>2022</td>\n",
|
||||
" <td>1</td>\n",
|
||||
" <td>9</td>\n",
|
||||
" <td>2022.0233</td>\n",
|
||||
" <td>418.33</td>\n",
|
||||
" <td>7</td>\n",
|
||||
" <td>415.06</td>\n",
|
||||
" <td>393.99</td>\n",
|
||||
" <td>138.10</td>\n",
|
||||
" </tr>\n",
|
||||
" </tbody>\n",
|
||||
"</table>\n",
|
||||
"<p>2487 rows × 9 columns</p>\n",
|
||||
"</div>"
|
||||
],
|
||||
"text/plain": [
|
||||
" year month day decimal co2 ndays last_year last_decade \n",
|
||||
"0 1974 5 19 1974.3795 333.37 5 -999.99 -999.99 \\\n",
|
||||
"1 1974 5 26 1974.3986 332.95 6 -999.99 -999.99 \n",
|
||||
"2 1974 6 2 1974.4178 332.35 5 -999.99 -999.99 \n",
|
||||
"3 1974 6 9 1974.4370 332.20 7 -999.99 -999.99 \n",
|
||||
"4 1974 6 16 1974.4562 332.37 7 -999.99 -999.99 \n",
|
||||
"... ... ... ... ... ... ... ... ... \n",
|
||||
"2482 2021 12 12 2021.9466 416.45 7 414.00 392.19 \n",
|
||||
"2483 2021 12 19 2021.9658 417.49 6 414.86 392.54 \n",
|
||||
"2484 2021 12 26 2021.9849 417.46 7 415.32 392.43 \n",
|
||||
"2485 2022 1 2 2022.0041 417.43 6 415.44 393.18 \n",
|
||||
"2486 2022 1 9 2022.0233 418.33 7 415.06 393.99 \n",
|
||||
"\n",
|
||||
" increase \n",
|
||||
"0 50.40 \n",
|
||||
"1 50.06 \n",
|
||||
"2 49.60 \n",
|
||||
"3 49.65 \n",
|
||||
"4 50.06 \n",
|
||||
"... ... \n",
|
||||
"2482 137.16 \n",
|
||||
"2483 137.93 \n",
|
||||
"2484 137.66 \n",
|
||||
"2485 137.40 \n",
|
||||
"2486 138.10 \n",
|
||||
"\n",
|
||||
"[2487 rows x 9 columns]"
|
||||
]
|
||||
},
|
||||
"execution_count": 4,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"data"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "1c178424",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"We can plot this fairly easily using Matplotlib."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 5,
|
||||
"id": "e63b38c5",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAtAAAAK/CAYAAABA9foTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAByMklEQVR4nO3dd3gU1eLG8XfTE1IIpFFC6L2DQkAUFBO9gKKCikoTERBUwIvIVQH1KioqogK2K3AVRFHxWhAITcFEei+hF4GEnkCA1Pn9kV/WrJsNGVI2JN/P8+QhO3N25syeZPNy9sw5FsMwDAEAAAAoEBdnVwAAAAC4nhCgAQAAABMI0AAAAIAJBGgAAADABAI0AAAAYAIBGgAAADCBAA0AAACYQIAGAAAATCBAAwAAACYQoAEAAAATCNAAAACACQRoAEVi4sSJslgsOn36dImcb9asWbJYLDp06FCJnA9lkzN/jvgZBq5fBGigmEyfPl0Wi0Xt2rXLc39sbKwmTpyo8+fPl2zF4ND+/fs1ZMgQ1a5dW15eXvL391fHjh01depUXb582absjh079Mgjj6hatWry9PRU1apV9fDDD2vHjh1Oqv31rbh/H8rL71tOKLdYLFq9erXdfsMwFB4eLovFou7duzuhhkUv55rXr1/v7KqgHCFAA8Vkzpw5qlmzptauXat9+/bZ7Y+NjdVLL71U5v+gF5e+ffvq8uXLioiIKJLj/fzzz2rWrJm+/vpr9ejRQ++//74mTZqkGjVqaMyYMXr66aetZb/77ju1bt1ay5Yt08CBAzV9+nQNGjRIK1asUOvWrbVgwYIiqVN5Uty/D46OX9Q/R2YU57m9vLw0d+5cu+2//vqr/vzzT3l6ehb5OYHyxM3ZFQDKooMHDyo2NlbfffedhgwZojlz5mjChAnOrlaZkJKSogoVKsjV1VWurq5FcsyDBw/qwQcfVEREhJYvX64qVapY9w0fPlz79u3Tzz//LCm7l7pv376qXbu2fvvtNwUHB1vLPv300+rUqZP69u2rrVu3qnbt2kVSP9jL+TkorKL8OSpN5/7HP/6h+fPn67333pOb219/6ufOnas2bdqU2FAroKyiBxooBnPmzFFgYKC6deumXr16ac6cOTb7J06cqDFjxkiSatWqZf3I9dChQ9axxHv27NEjjzyigIAABQcH68UXX5RhGDp69Kjuvvtu+fv7KywsTG+//Xaeddi0aZPuvPNO+fv7y9fXV7fddpv++OMPmzIXLlzQyJEjVbNmTXl6eiokJES33367Nm7caK2nxWLR7t27df/998vf31+VK1fW008/rStXruR53vPnz2vAgAGqWLGiAgICNHDgQF26dMmu3LFjx/Too48qNDRUnp6eatKkiT777DO718lisWjnzp166KGHFBgYqJtuuklS3uNHjx07pkGDBqlq1ary9PRUrVq1NGzYMKWlpeXTWtKbb76pixcv6j//+Y9NeM5Rt25daw/05MmTdenSJX388cc24VmSgoKC9NFHHyklJUVvvvlmvucsqIJcU0HaOue13LdvX4Hb52rnNdOG+Z03v9+H3MfI6+fg8OHDeuKJJ9SgQQN5e3urcuXK6t27t83PRX7HdzQOuThe07/L69yFPWaOPn366MyZM4qJibFuS0tL0zfffKOHHnrIrnxBXkdJGjBggGrWrGn3/Jx657XtatdS0HMXVkHatKTqgusfPdBAMZgzZ47uvfdeeXh4qE+fPpoxY4bWrVunG264QZJ07733as+ePfryyy81ZcoUBQUFSZJNIHvggQfUqFEjvf766/r555/173//W5UqVdJHH32kW2+9VW+88YbmzJmjf/7zn7rhhht08803W5+7Y8cOderUSf7+/nr22Wfl7u6ujz76SJ07d9avv/5qHZc9dOhQffPNNxoxYoQaN26sM2fOaPXq1dq1a5dat25tPd7999+vmjVratKkSfrjjz/03nvv6dy5c/rvf/9rd+3333+/atWqpUmTJmnjxo369NNPFRISojfeeMNaJjExUe3bt5fFYtGIESMUHBysX375RYMGDVJycrJGjhxpc8zevXurXr16eu2112QYRp6v+fHjx3XjjTfq/Pnzevzxx9WwYUMdO3ZM33zzjS5duiQPDw+H7fXjjz+qdu3a6tChg8MyucvWrFlTnTp1ynP/zTffrJo1a1p7rAujINdU0LbOUZD2Kch5zbZhfuctyO+DlPfPwbp16xQbG6sHH3xQ1atX16FDhzRjxgx17txZO3fulI+PT4GPn6M4XlOzCnvMmjVrKjIyUl9++aXuvPNOSdIvv/yipKQkPfjgg3rvvfdsyhfkdSyuaynOc+coaJuWRF1QRhgAitT69esNSUZMTIxhGIaRlZVlVK9e3Xj66adtyk2ePNmQZBw8eNBm+4QJEwxJxuOPP27dlpGRYVSvXt2wWCzG66+/bt1+7tw5w9vb2+jfv7/NMXr27Gl4eHgY+/fvt247fvy44efnZ9x8883WbQEBAcbw4cMdXktOXe666y6b7U888YQhydiyZYtd2UcffdSm7D333GNUrlzZZtugQYOMKlWqGKdPn7bZ/uCDDxoBAQHGpUuXbI7Zp08fu7rNnDnT5vXr16+f4eLiYqxbt86ubFZWlsNrTEpKMiQZd999t8MyOc6fP1+gsnfddZchyUhOTr7qMfNTkGsqaFubaZ+CnNdsG17tvI5+H3IfI6+fg5zz5BYXF2dIMv773/9e9fh//zkyjOJ5TfOS17mL6pjr1q0zPvjgA8PPz8/6GvXu3dvo0qWLYRiGERERYXTr1s36vIK+jv379zciIiLsyubUO69tV7uWgp67INfsSEHbtLB1QfnBEA6giM2ZM0ehoaHq0qWLJMliseiBBx7QvHnzlJmZWeDjPPbYY9bvXV1d1bZtWxmGoUGDBlm3V6xYUQ0aNNCBAwes2zIzM7VkyRL17NnTZgxulSpV9NBDD2n16tVKTk62Pn/NmjU6fvx4vnUZPny4zeMnn3xSkrRw4UK7skOHDrV53KlTJ505c8Z6TsMw9O2336pHjx4yDEOnT5+2fkVHRyspKck6hMTRMf8uKytL33//vXr06KG2bdva7f/7R8u55dTLz88v33NI2UNeClI2Z3/Osa9FQa7JTFvnuFr7FOS8RdGGfz9vQeT1c+Dt7W39Pj09XWfOnFHdunVVsWJFuzoURHG8pteiKI55//336/Lly/rpp5904cIF/fTTT3kO35CK/nXM7WrXUpznlsy1aXHXBWUHARooQpmZmZo3b566dOmigwcPat++fdq3b5/atWunxMRELVu2rMDHqlGjhs3jgIAAeXl5WT9+zr393Llz1senTp3SpUuX1KBBA7tjNmrUSFlZWTp69Kik7LG/27dvV3h4uG688UZNnDjRJoznqFevns3jOnXqyMXFJc9xgX+vd2BgoCRZ63jq1CmdP3/eOoY499fAgQMlSSdPnrQ5Rq1atezOk9upU6eUnJyspk2b5lsuL/7+/pL+Csf5yQnGVyubO2inpqbq0UcfVY0aNeTv76/27dsrLi7uqucqyDWZaescBWmfgpzXbBte7bwFkdfPweXLlzV+/HiFh4fL09NTQUFBCg4O1vnz55WUlFTgY+cojtf0WhTFMYODg9W1a1fNnTtX3333nTIzM9WrV688yxb165jb1a6lOM8tmWvT4q4Lyg7GQANFaPny5Tpx4oTmzZunefPm2e2fM2eOoqKiCnSsvO7Od3THvuFgXPDV3H///erUqZMWLFigJUuWaPLkyXrjjTf03XffWcdN5iW/Ht2r1TErK0uS9Mgjj6h///55lm3evLnN49y9QkXN399fVatW1fbt269aNiAgQFWqVNHWrVvzLbd161ZVq1ZN/v7+SklJUc2aNbV69WpVr17dOk3eoUOH5OvrW1SXUWBF8TN0LW1YFOfN6+fgySef1MyZMzVy5EhFRkYqICBAFotFDz74oLWexa2ofy+L8pgPPfSQBg8erISEBN15552qWLFinuUK+jo6+t3P79O1q11LaWjDHKWpLijdCNBAEZozZ45CQkI0bdo0u33fffedFixYoA8//FDe3t75htDCCA4Olo+Pj+Lj4+327d69Wy4uLgoPD7duq1Klip544gk98cQTOnnypFq3bq1XX33VJkDv3bvXpvdv3759ysrKyvNu/ILUz8/PT5mZmeratavp5zs6pr+/f4FCcF66d++ujz/+WHFxcYqMjLxq2U8++USrV6+2zgSR26pVq3To0CENGTJEklShQgWNHz/euv/BBx/U6NGjFR8frzZt2hTqmsy2dUEU9LxF3YbX+vvwzTffqH///jaz0Vy5csVuvueCHr84XlNnuueeezRkyBD98ccf+uqrrxyWK+jrGBgYmOdc3YcPH77mOhb03NfKTJsWd11QdjCEAygily9f1nfffafu3burV69edl8jRozQhQsX9MMPP0iSdQ7bon5jdnV1VVRUlP73v//ZDLFITEzU3LlzddNNN8nf31+ZmZl2H0mGhISoatWqSk1Ntdn+9/8QvP/++5KUby91fvW777779O233+YZ0k6dOmX6mC4uLurZs6d+/PHHPFcju1qv3bPPPqsKFSroscceU2Jiot3+/fv3a+rUqZKkMWPGyNvbW0OGDNGZM2dsyp09e1ZDhw6Vj4+Pddq0v9u7d6/Onj2runXrFvqaCtrWZhT0vEXdhtf6++Dq6mrXvu+//75dj2hBj18cr6kz+fr6asaMGZo4caJ69OjhsFxBX8c6deooKSnJ5lOYEydOFGrxoIKeuzDHL2ibFnddUHbQAw0UkR9++EEXLlzQXXfdlef+9u3bKzg4WHPmzNEDDzxg7X18/vnn9eCDD8rd3T3fP3Bm/Pvf/1ZMTIxuuukmPfHEE3Jzc9NHH32k1NRU6/zEFy5cUPXq1dWrVy+1aNFCvr6+Wrp0qdatW2c3t/TBgwd111136Y477lBcXJyLine truncated
|
||||
"text/plain": [
|
||||
"<Figure size 800x800 with 1 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"fig, ax = plt.subplots(figsize=(8, 8))\n",
|
||||
"data.plot(x=\"decimal\", y=\"co2\", legend=False, ax=ax)\n",
|
||||
"ax.set(xlabel=\"Time [years]\", ylabel=r\"CO$_2$ [ppm]\", title=\"Atmospheric CO$_2$ concentration in Mauna Loa\")\n",
|
||||
"plt.show()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "efd249fe",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"The spikes are related to missing data. We can clean those up and plot again."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 6,
|
||||
"id": "0837b3ff",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"data = data[data.co2 > -999]"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 7,
|
||||
"id": "333581b5",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAArsAAAK/CAYAAABp8gWUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAADHGUlEQVR4nOzdd3gc5fU24Gf7SiutenGRe8cFbINx6JhuSggldEgI1RDKF5LwI4WQAiEJCb0kBFIMJBAgQOi9uGBswL13q9eVtpf5/ph5Z2d2Z1YrW7JX0nNfly9mdt+dGVkuh+PznmORJEkCEREREdEAZD3QD0BERERE1FcY7BIRERHRgMVgl4iIiIgGLAa7RERERDRgMdglIiIiogGLwS4RERERDVgMdomIiIhowGKwS0REREQDFoNdIiIiIhqwGOwSERER0YDFYJeIiIiIBiwGu0SD1J133gmLxYLm5ub9cr+nn34aFosF27dv3y/3o4HpQP464q9hov6JwS6R4pFHHoHFYsGcOXMM31+0aBHuvPNOtLe3798HI1NbtmzBNddcgzFjxsDtdsPr9eKII47A/fffj2AwqFu7Zs0aXHLJJRg2bBhcLheGDh2Kiy++GGvWrDlAT9+/9fXvh8Hy+00E0BaLBZ9++mna+5IkoaamBhaLBaeffvoBeMLeJ77mL7744kA/Cg0SDHaJFAsXLsSoUaPw+eefY/PmzWnvL1q0CL/4xS8G/F++feXSSy9FMBjEyJEje+V6//vf/zBt2jT8+9//xhlnnIEHH3wQd999N0aMGIHbbrsNN910k7r2xRdfxMyZM/Hee+/hO9/5Dh555BFceeWV+OCDDzBz5ky89NJLvfJMg0lf/34wu35v/zrqib68t9vtxjPPPJP2+kcffYTdu3fD5XL1+j2JBgv7gX4Aolywbds2LFq0CC+++CKuueYaLFy4ED//+c8P9GMNCH6/Hx6PBzabDTabrVeuuW3bNlxwwQUYOXIk3n//fQwZMkR9b8GCBdi8eTP+97//AZCzv5deeinGjBmDjz/+GBUVFeram266CUcddRQuvfRSrFy5EmPGjOmV56N04tfBvurNX0e5dO/TTjsNzz//PB544AHY7cm/mp955hnMmjVrv5UbEQ1EzOwSQc7qlpSUYP78+Tj33HOxcOFC3ft33nknbrvtNgDA6NGj1X923L59u1r7unHjRlxyySUoKipCRUUFfvrTn0KSJOzatQtnnXUWvF4vqqur8Yc//MHwGb788kuceuqp8Hq9KCgowLx587BkyRLdms7OTtx8880YNWoUXC4XKisrceKJJ2LFihXqc1osFqxfvx7nn38+vF4vysrKcNNNNyEUChnet729HVdccQWKi4tRVFSE73znOwgEAmnr9uzZg+9+97uoqqqCy+XCQQcdhL/+9a9pP08WiwVr167FRRddhJKSEhx55JEAjOsd9+zZgyuvvBJDhw6Fy+XC6NGjcd111yESiWT4bgH33nsvurq68OSTT+oCXWHcuHFqZvd3v/sdAoEAnnjiCV2gCwDl5eV4/PHH4ff7ce+992a8Z7ay+Zqy+V6Ln8vNmzdn/f3p7r49+R5mum+m3w/aaxj9OtixYweuv/56TJw4EXl5eSgrK8N5552n+3WR6fpmdbN98XOayuje+3pN4cILL0RLSwveeecd9bVIJIIXXngBF110Udr6bH4eAeCKK67AqFGj0j4vntvote6+lmzvva+y+Z7ur2eh/o2ZXSLIwe63vvUtOJ1OXHjhhXj00UexbNkyHHrooQCAb33rW9i4cSOeffZZ/PGPf0R5eTkA6IKnb3/725g8eTLuuece/O9//8OvfvUrlJaW4vHHH8fxxx+P3/72t1i4cCF+8IMf4NBDD8XRRx+tfnbNmjU46qij4PV68cMf/hAOhwOPP/44jj32WHz00UdqHfG1116LF154ATfccAOmTJmClpYWfPrpp1i3bh1mzpypXu/888/HqFGjcPfdd2PJkiV44IEH0NbWhr///e9pX/v555+P0aNH4+6778aKFSvwl7/8BZWVlfjtb3+rrmloaMDhhx8Oi8WCG264ARUVFXjjjTdw5ZVXwufz4eabb9Zd87zzzsP48ePxm9/8BpIkGf6c19bW4rDDDkN7ezuuvvpqTJo0CXv27MELL7yAQCAAp9Np+v169dVXMWbMGHzjG98wXaNdO2rUKBx11FGG7x999NEYNWqUmgneF9l8Tdl+r4Vsvj/Z3Len38NM983m9wNg/Otg2bJlWLRoES644AIMHz4c27dvx6OPPopjjz0Wa9euRX5+ftbXF/ri57Sn9vWao0aNwty5c/Hss8/i1FNPBQC88cYb6OjowAUXXIAHHnhAtz6bn8e++lr68t5Ctt/T/fEsNABIRIPcF198IQGQ3nnnHUmSJCmRSEjDhw+XbrrpJt263/3udxIAadu2bbrXf/7zn0sApKuvvlp9LRaLScOHD5csFot0zz33qK+3tbVJeXl50uWXX667xje/+U3J6XRKW7ZsUV+rra2VCgsLpaOPPlp9raioSFqwYIHp1yKe5cwzz9S9fv3110sApK+//jpt7Xe/+13d2rPPPlsqKyvTvXbllVdKQ4YMkZqbm3WvX3DBBVJRUZEUCAR017zwwgvTnu2pp57S/fxddtllktVqlZYtW5a2NpFImH6NHR0dEgDprLPOMl0jtLe3Z7X2zDPPlABIPp+v22tmks3XlO33uiffn2zu29PvYXf3Nfv9oL2G0a8DcR+txYsXSwCkv//9791eP/XXkST1zc+pEaN799Y1ly1bJj300ENSYWGh+nN03nnnSccdd5wkSZI0cuRIaf78+ernsv15vPzyy6WRI0emrRXPbfRad19LtvfO5ms2k+33dF+fhQYHljHQoLdw4UJUVVXhuOOOAwBYLBZ8+9vfxnPPPYd4PJ71db73ve+pxzabDbNnz4YkSbjyyivV14uLizFx4kRs3bpVfS0ej+Ptt9/GN7/5TV3N6JAhQ3DRRRfh008/hc/nUz+/dOlS1NbWZnyWBQsW6M5vvPFGAMDrr7+etvbaa6/VnR911FFoaWlR7ylJEv7zn//gjDPOgCRJaG5uVn+cfPLJ6OjoUMsozK6ZKpFI4OWXX8YZZ5yB2bNnp72f+s+rWuK5CgsLM94DkMs+slkr3hfX3hvZfE09+V4L3X1/srlvb3wPU++bDaNfB3l5eepxNBpFS0sLxo0bh+Li4rRnyEZf/Jzujd645vnnn49gMIjXXnsNnZ2deO211wxLGIDe/3nU6u5r6ct7Az37nvb1s9DAwGCXBrV4PI7nnnsOxx13HLZt24bNmzdj8+bNmDNnDhoaGvDee+9lfa0RI0bozouKiuB2u9V/gtW+3tbWpp43NTUhEAhg4sSJadecPHkyEokEdu3aBUCuVV29ejVqampw2GGH4c4779QFzsL48eN152PHjoXVajWsY0t97pKSEgBQn7GpqQnt7e1qzav2x3e+8x0AQGNjo+4ao0ePTruPVlNTE3w+H6ZOnZpxnRGv1wsgGchmIoLY7tZqg+JwOIzvfve7GDFiBLxeLw4//HAsXry423tl8zX15HstZPP9yea+Pf0ednffbBj9OggGg/jZz36GmpoauFwulJeXo6KiAu3t7ejo6Mj62kJf/Jzujd64ZkVFBU444QQ888wzePHFFxGPx3Huuecaru3tn0et7r6Wvrw30LPvaV8/Cw0MrNmlQe39999HXV0dnnvuOTz33HNp7y9cuBAnnXRSVtcy2qVttnNbMqlj7c7555+Po446Ci+99BLefvtt/O53v8Nvf/tbvPjii2qdn5FMmdLunjGRSAAALrnkElx++eWGa6dPn64712ZbepvX68XQoUOxevXqbtcWFRVhyJAhWLlyZcZ1K1euxLBhw+D1euH3+zFq1Ch8+umnGD58uNrabPv27SgoKOitLyNrvfFraG++h71xX6NfBzfeeCOeeuop3HzzzZg7dy6KiopgsVhwwQUXqM/Z13r792VvXvOiiy7CVVddhfr6epx66qkoLi42XJftz6PZ7/1M/2rV3deSC99DIZeehXIXg10a1BYuXIjKyko8/PDDae+9+OKLeOmll/DYY48hLy8vY8C4LyoqKpCfn48NGzakvbd+/XpYrVbU1NSorw0ZMgTXX389rr/+ejQ2NmLmzJn49a9/rQt2N23apMuqbd68GYlEwnBXdjbPV1hYiHg8jhNOOKHHnze7ptfrzSpgNXL66afjiSeewOLFizF37txu1/75z3/Gp59+qnYE0Prkk0+wfft2XHPNNQAAj8eDn/3sZ+r7F1xwAW699VZs2LABs2bN2qevqaff62xke9/e/h7u7e+HF154AZdffrmuK0koFErrp5vt9fvi5/RAOvvss3HNNddgyZIl+Ne//mW6Ltufx5KSEsNeyDt27NjrZ8z23nurJ9/Tvn4WGhhYxkCDVjAYxIsvvojTTz8d5557btqPG264AZ2dnXjllVcAQO0R2tt/iNpsNpx00kn473//qyszaGhowDPPPIMjjzwSXq8X8Xg87Z/lKisrMXToUITDYd3rqcH7gw8+CAAZs7+Znu+cc87Bf/7zH8OAqqmpqcfXtFqt+OY3v4lXX33VcIpSd9mwH/7wh/B4PPje976HhoaGtPe3bNmC+++/HwBw2223IS8vD9dccw1aWlp061pbW3HttdciPz9fbXWVatOmTWhtbcW4ceP2+WvK9nvdE9net7e/h3v7+8Fms6V9fx988MG0TGO21++Ln9MDqaCgAI8++ijuvPNOnHHGGabrsv15HDt2LDo6OnT/ulFXV7dPg1Syvfe+XD/b72lfPwsNDMzs0qD1yiuvoLOzE2eeeabh+4cffjgqKiqwcOFCfPvb31azenfccQcuuOACOByOjH8Z9cSvfvUrvPPOOzjyyCNx/fXXw2634/HHH0c4HFb7v3Z2dmL48OE499xzMWPGDBQLine truncated
|
||||
"text/plain": [
|
||||
"<Figure size 800x800 with 1 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"fig, ax = plt.subplots(figsize=(8, 8))\n",
|
||||
"data.plot(x=\"decimal\", y=\"co2\", legend=False, ax=ax)\n",
|
||||
"ax.set(xlabel=\"Time [years]\", ylabel=r\"CO$_2$ [ppm]\", title=\"Atmospheric CO$_2$ concentration in Mauna Loa\")\n",
|
||||
"plt.show()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "c62fd7f3",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"That's better.\n",
|
||||
"\n",
|
||||
"We can now try to model this behaviour with Gaussian Processes. The key parameter we must choose is the covariance matrix analytic form, which establishes how correlated two points in the time axis are.\n",
|
||||
"\n",
|
||||
"One simple assumption would be that two values of the CO2 concentration are highly correlated they happened at close by times. We use a Gaussian-based model (named here \"Radial Basis Function\", as this does not refer to a probability distribution) to establish that if two points are one year apart they are likely to be highly correlated.\n",
|
||||
"\n",
|
||||
"Additionally, we assume that there is a correlation between points within certain period, since we know that data from the same season in an year is expected to be correlated to data if that same season in the year (ie: CO$_2$ emissions have a rough periodicity of one year due to the season changes).\n",
|
||||
"\n",
|
||||
"Another term of the covariance matrix must be related to the random noise appearing in the data.\n",
|
||||
"\n",
|
||||
"All those assumptions can be motivated by an analysis of the covariance matrix of the data. That is, we can take the covariance matrix of all the data points and analyse it to establish whether our assumptions here are valid or not. A very long and detailed explanation on this point can be found in Bishop (2006), which has a very pedagogical chapter on Gaussian Processes. Online resources can be found on: http://www.gaussianprocess.org/\n",
|
||||
"\n",
|
||||
"Note that all of those are very strong assumptions, but they are made directly on the covariance matrix of the fit function. The actual numbers used below are taken from an example, but they will be optimized in the fit below (unless explicitly asked not to).\n",
|
||||
"\n",
|
||||
"Note: this notebook is a simplified example with extra explanations, with lots of material taken from https://scikit-learn.org/stable/auto_examples/gaussian_process/plot_gpr_co2.html -- take a look at the original example for more details.\n",
|
||||
"\n",
|
||||
"We also assume the data within an year from each other is highly correlated and should show a periodicity of 1 year."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 8,
|
||||
"id": "a31d7eca",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Kernel with parameters given in GPML book\n",
|
||||
"\n",
|
||||
"# these hyper parameters are optimised in the GP fit\n",
|
||||
"# we provide here only a starting value\n",
|
||||
"\n",
|
||||
"# Kernel with optimized parameters\n",
|
||||
"k1 = 50.0 ** 2 * RBF(length_scale=50.0) # long term smooth rising trend\n",
|
||||
"k2 = (\n",
|
||||
" 2.0 ** 2\n",
|
||||
" * RBF(length_scale=100.0)\n",
|
||||
" * ExpSineSquared(length_scale=1.0, periodicity=1.0, periodicity_bounds=\"fixed\")\n",
|
||||
") # seasonal component\n",
|
||||
"# medium term irregularities\n",
|
||||
"k3 = 0.5 ** 2 * RationalQuadratic(length_scale=1.0, alpha=1.0)\n",
|
||||
"k4 = 0.1 ** 2 * RBF(length_scale=0.1) + WhiteKernel(\n",
|
||||
" noise_level=0.1 ** 2, noise_level_bounds=(1e-5, np.inf)\n",
|
||||
") # noise terms\n",
|
||||
"kernel = k1 + k2 + k3 + k4"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "0ad8ddf0",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"The kernel hyper-parameters are:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 9,
|
||||
"id": "032ada0d",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"50**2 * RBF(length_scale=50) + 2**2 * RBF(length_scale=100) * ExpSineSquared(length_scale=1, periodicity=1) + 0.5**2 * RationalQuadratic(alpha=1, length_scale=1) + 0.1**2 * RBF(length_scale=0.1) + WhiteKernel(noise_level=0.01)\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"print(kernel)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 10,
|
||||
"id": "4412a9a7",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"gp = GaussianProcessRegressor(kernel=kernel, normalize_y=True)\n",
|
||||
"time = data.decimal.to_numpy()[:, np.newaxis]\n",
|
||||
"co2 = data.co2.to_numpy()[:, np.newaxis]\n",
|
||||
"gp = gp.fit(time, co2)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "eb69bb80",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"These are the hyper-parameters after the fit and the log-likelihood given the data:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 11,
|
||||
"id": "667db862",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"GPML kernel: 3.19**2 * RBF(length_scale=42.7) + 0.114**2 * RBF(length_scale=174) * ExpSineSquared(length_scale=1.31, periodicity=1) + 0.295**2 * RationalQuadratic(alpha=0.000879, length_scale=6.84) + 0.0101**2 * RBF(length_scale=0.0188) + WhiteKernel(noise_level=0.00014)\n",
|
||||
"Log-marginal-likelihood: 6660.284\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"print(\"GPML kernel: %s\" % gp.kernel_)\n",
|
||||
"print(\"Log-marginal-likelihood: %.3f\" % gp.log_marginal_likelihood(gp.kernel_.theta))"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "4968439e",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Now we can use the fit model to predict how the CO2 concentration will evolve at other times."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 12,
|
||||
"id": "ef6b7762",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"new_time = np.linspace(time.min(), time.max()+20, 1000)[:, np.newaxis]\n",
|
||||
"co2_pred, co2_std = gp.predict(new_time, return_std=True)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 14,
|
||||
"id": "5b263358",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAArsAAAK/CAYAAABp8gWUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAD3pUlEQVR4nOzdeXxcddn//9fMZCYzmWSyL23TlZallIKsFnEDZCtFRUCgQNEKqAVBvz9UbjdcQeTWWxbZREEtoiAoVGSTnRYsIFIKQgvd0jb7MllmSTLn98dnzpyZZNKmbZJmeT8fjzySOefMmXMakLefXue6XJZlWYiIiIiIjEPuvX0BIiIiIiLDRWFXRERERMYthV0RERERGbcUdkVERERk3FLYFREREZFxS2FXRERERMYthV0RERERGbcUdkVERERk3FLYFREREZFxS2FXRERERMYthV0RERERGbcUdkUmqKuvvhqXy0VjY+OIfN5dd92Fy+Vi48aNI/J5Mj7tzX+O9M+wyNiksCuS9Ktf/QqXy8VRRx2Vdf/KlSu5+uqraW1tHdkLkwG99957XHLJJcyaNQu/308oFOJDH/oQv/zlL4lEIhnHrl27lvPOO48pU6aQm5vL5MmTWbx4MWvXrt1LVz+2Dfe/DxPl3zc7QLtcLl544YV++y3LYurUqbhcLk499dS9cIVDz77nV155ZW9fikwQCrsiScuXL2fGjBn861//Yv369f32r1y5ku9///vj/j++w+X8888nEokwffr0ITnf3//+dw466CD+/Oc/s2jRIm688UauueYapk2bxpVXXsnll1+eOvaBBx7g0EMP5Z///Cef+9zn+NWvfsXSpUt5+umnOfTQQ3nwwQeH5JomkuH+92Gg8w/1P0e7Yjg/2+/3c8899/Tb/uyzz1JTU0Nubu6Qf6bIRJGzty9AZDTYsGEDK1eu5IEHHuCSSy5h+fLlfO9739vblzUudHZ2EgwG8Xg8eDyeITnnhg0bOPvss5k+fTpPPfUUkyZNSu1btmwZ69ev5+9//ztgVn/PP/98Zs2axXPPPUd5eXnq2Msvv5wPf/jDnH/++bzxxhvMmjVrSK5P+rP/OdhTQ/nP0Wj67FNOOYX77ruPG264gZwc5z/N99xzD4cddtiIlRuJjEda2RXBrOoWFxezcOFCzjjjDJYvX56x/+qrr+bKK68EYObMmam/dty4cWOq9vXdd9/lvPPOo7CwkPLycr7zne9gWRZbtmzhk5/8JKFQiKqqKv73f/836zX8+9//5uSTTyYUCpGfn89xxx3HSy+9lHFMe3s7V1xxBTNmzCA3N5eKigo+8YlP8Nprr6Wu0+Vy8d///pezzjqLUChEaWkpl19+OdFoNOvntra2cuGFF1JUVERhYSGf+9zn6Orq6nfc1q1b+fznP09lZSW5ubkceOCB/OY3v+n35+RyuXjrrbc499xzKS4u5phjjgGy1ztu3bqVpUuXMnnyZHJzc5k5cyZf+tKXiMfjO/htwXXXXUdHRwd33nlnRtC1zZ49O7Wy+7Of/Yyuri5uv/32jKALUFZWxm233UZnZyfXXXfdDj9zsAZzT4P5Xdt/luvXrx/072dnn7srv8Mdfe6O/n1IP0e2fw42bdrEl7/8Zfbbbz8CgQClpaWceeaZGf9c7Oj8A9XNDsefaV/ZPntPz2k755xzaGpq4oknnkhti8fj3H///Zx77rn9jh/MnyPAhRdeyIwZM/q9377ubNt2di+D/ew9NZjf6Uhdi4xtWtkVwYTd008/HZ/PxznnnMMtt9zC6tWrOeKIIwA4/fTTeffdd/njH//IL37xC8rKygAywtNnP/tZDjjgAK699lr+/ve/86Mf/YiSkhJuu+02jj32WH7605+yfPly/r//7//jiCOO4CMf+UjqvWvXruXDH/4woVCIr3/963i9Xm677TY+9rGP8eyzz6bqiL/4xS9y//33c+mllzJ37lyampp44YUXePvttzn00ENT5zvrrLOYMWMG11xzDS+99BI33HADLS0t/O53v+t372eddRYzZ87kmmuu4bXXXuPXv/41FRUV/PSnP00dU1dXxwc/+EFcLheXXnop5eXl/OMf/2Dp0qWEw2GuuOKKjHOeeeaZzJkzh5/85CdYlpX1z3zbtm0ceeSRtLa2cvHFF7P//vuzdetW7r//frq6uvD5fAP+vh5++GFmzZrF0UcfPeAx6cfOmDGDD3/4w1n3f+QjH2HGjBmpleA9MZh7Guzv2jaY389gPndXf4c7+tzB/PsA2f85WL16NStXruTss8+murqajRs3csstt/Cxj32Mt956i7y8vEGf3zYcf6a7ak/POWPGDBYsWMAf//hHTj75ZAD+8Y9/0NbWxtlnn80NN9yQcfxg/hyH616G87Ntg/2djsS1yDhgiUxwr7zyigVYTzzxhGVZlpVIJKzq6mrr8ssvzzjuZz/7mQVYGzZsyNj+ve99zwKsiy++OLWtp6fHqq6utlwul3Xttdemtre0tFiBQMBasmRJxjk+9alPWT6fz3rvvfdS27Zt22YVFBRYH/nIR1LbCgsLrWXLlg14L/a1nHbaaRnbv/zlL1uA9Z///KffsZ///Oczjv30pz9tlZaWZmxbunSpNWnSJKuxsTFj+9lnn20VFhZaXV1dGec855xz+l3bb3/724w/vwsuuMByu93W6tWr+x2bSCQGvMe2tjYLsD75yU8OeIyttbV1UMeedtppFmCFw+GdnnNHBnNPg/1d78rvZzCfu6u/w5197kD/PqSfI9s/B/bnpFu1apUFWL/73e92ev6+/xxZ1vD8mWaT7bOH6pyrV6+2brrpJqugoCD1Z3TmmWdaH//4xy3Lsqzp06dbCxcuTL1vsH+OS5YssaZPn97vWPu6s23b2b0M9rMHc88DGezvdE+vRSYGlTHIhLd8+XIqKyv5+Mc/DoDL5eKzn/0s9957L729vYM+zxe+8IXUzx6Ph8MPPxzLsli6dGlqe1FREfvttx/vv/9+altvby+PP/44n/rUpzJqRidNmsS5557LCy+8QDgcTr3/5ZdfZtu2bTu8lmXLlmW8vuyyywB45JFH+h37xS9+MeP1hz/8YZqamlKfaVkWf/nLX1i0aBGWZdHY2Jj6OvHEE2lra0uVUQx0zr4SiQR//etfWbRoEYcffni//X3/ejWdfV0FBQU7/AwwZR+DOdbeb597dwzmnnbld23b2e9nMJ87FL/Dvp87GNn+OQgEAqmfu7u7aWpqYvbs2RQVFfW7hsEYjj/T3TEU5zzrrLOIRCKsWLGC9vZ2VqxYkbWEAYb+zzHdzu5lOD8bdu13OtzXIuODwq5MaL29vdx77718/OMfZ8OGDaxfv57169dz1FFHUVdXxz//+c9Bn2vatGkZrwsLC/H7/am/gk3f3tLSknrd0NBAV1cX++23X79zHnDAASQSCbZs2QKYWtU333yTqVOncuSRR3L11VdnBGfbnDlzMl7vs88+uN3urHVsfa+7uLgYIHWNDQ0NtLa2pmpe078+97nPAVBfX59xjpkzZ/b7nHQNDQ2Ew2HmzZu3w+OyCYVCgBNkd8QOsTs7Nj0Ux2IxPv/5zzNt2jRCoRAf/OAHWbVq1U4/azD3tCu/a9tgfj+D+dxd/R3u7HMHI9s/B5FIhO9+97tMnTqV3NxcysrKKC8vp7W1lba2tkGf2zYcf6a7YyjOWV5ezvHHH88999zDAw88QG9vL2eccUbWY4f6zzHdzu5lOD8bdu13OtzXIuODanZlQnvqqafYvn079957L/fee2+//cuXL+eEE04Y1LmyPaU90JPb1gB1rDtz1lln8eEPf5gHH3yQxx9/nJ/97Gf89Kc/5YEHHkjV+WWzo5XSnV1jIpEA4LzzzmPJkiVZj50/f37G6/TVlqEWCoWYPHkyb7755k6PLSwsZNKkSbzxxhs7PO6NN95gypQphEIhOjs7mTFjBi+88ALV1dWp1mYbN24kPz9/qG5j0Ibin6Hd+R0Oxedm++fgsssu47e//S1XXHEFCxYsoLCwEJfLxdlnn526zuE21P9eDuU5zz33XC666CJqa2s5+eSTKSoqynrcYP8cB/p3f0d/a7WzexkNv0PbaLoWGb0UdmVCW758ORUVFdx888399j3wwAM8+OCD3HrrrQQCgR0Gxj1RXl5OXl4e77zzTr99//3vf3G73UydOjW1bdKkSXz5y1/my1/+MvX19Rx66KH8+Mc/zgi769aty1hVW79+PYlEIutT2YO5voKCAnp7ezn++ON3+f0DnTMUCg0qsGZz6qmncvvtt7Nq1SoWLFiw02PvuOMOXnjhhVRHgHTPP/88Gzdu5JJLLgEgGAzy3e9+N7X/7LPP5mtf+xrvvPMOhx122B7d067+rgdjsJ871L/D3f334f7772fJkiUZXUmi0Wi/frqDPf9w/JnuTZ/+9Ke55JJLeOmll/jTn/404HGD/XMsLi7O2gt506ZNu32Ng/3s3bUrv9PhvhYZH1TGIBNWJBLhgQce4NRTT+WMM87o93XppZfS3t7OQw89BJDqETrU/yPq8Xg44YQT+Nvf/pZRZlBXV8c999zDMcccQygUore3t99fy1VUVDB58mRisVjG9r7h/cYbbwTY4ervjq7vM5/5DH/5y1+yBqqGhoZdPqfb7eZTn/oUDz/8cNYpSjtbDfv6179OMBjkC1/4AnV1df32v/fee/zyl78E4MorryQQCHDJJZfQ1NSUcVxzczNf/OIXycvLS7W66mvdunU0Nzcze/bsPb6nwf6ud8VgP3eof4e7+++Dx+Pp9/u98cYb+600Dvb8w/Fnujfl5+dzyy23cPXVV7No0aIBjxvsn+M+++xDW1tbxt9ubN++fY8GqQz2s/fk/IP9nQ73tcj4oJVdmbAeeugh2tvbOe2007Lu/+AHP0h5eTnLly/ns5/9bGpV71vf+hZnn302Xq93h/8x2hU/+tGPeOKJJzjmmGP48peLine truncated
|
||||
"text/plain": [
|
||||
"<Figure size 800x800 with 1 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"fig, ax = plt.subplots(figsize=(8, 8))\n",
|
||||
"ax.scatter(time, co2, color='red', marker='.', label=\"Data\")\n",
|
||||
"ax.fill_between(new_time[:,0],\n",
|
||||
" co2_pred - co2_std,\n",
|
||||
" co2_pred + co2_std,\n",
|
||||
" color=\"b\",\n",
|
||||
" alpha=0.4,\n",
|
||||
" label=\"Predictions\")\n",
|
||||
"ax.set(xlabel=\"Time [years]\", ylabel=r\"CO$_2$ [ppm]\", title=\"Atmospheric CO$_2$ concentration in Mauna Loa\")\n",
|
||||
"ax.legend(frameon=False, loc=\"lower right\")\n",
|
||||
"plt.show()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "b625b695",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Many more details on the kernel selection as well as further details in this example can be seen in the following references:\n",
|
||||
"\n",
|
||||
"https://scikit-learn.org/stable/auto_examples/gaussian_process/plot_gpr_co2.html\n",
|
||||
"\n",
|
||||
"https://scikit-learn.org/stable/auto_examples/gaussian_process/plot_gpr_noisy_targets.html#sphx-glr-auto-examples-gaussian-process-plot-gpr-noisy-targets-py\n",
|
||||
"\n",
|
||||
"To read more on Gaussian Processes:\n",
|
||||
"\n",
|
||||
"http://www.gaussianprocess.org/gpml/chapters/RW.pdf"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"id": "6d44d78a",
|
||||
"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.6"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 5
|
||||
}
|
||||
@@ -1,199 +0,0 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "1bba0128",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"# Simple fits\n",
|
||||
"\n",
|
||||
"Here we are fitting a line from scratch.\n",
|
||||
"In the next notebook, we will do fancier fits with neural networks, but let's start with a basic problem and complicate it as we go along.\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 1,
|
||||
"id": "23feddde",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"from typing import Tuple\n",
|
||||
"\n",
|
||||
"import pandas as pd\n",
|
||||
"import numpy as np\n",
|
||||
"import matplotlib.pyplot as plt\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "bb1286f0",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"We start by generating some fake dataset, which is simple enough that we can visualize the results easily. For this reason, the dataset will contain only two variables.\n",
|
||||
"\n",
|
||||
"The simulated example data will be $f(x) = 3 x + \\epsilon$, where $\\epsilon \\sim \\mathcal{N}(\\mu=0, \\sigma=0.5)$.\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
"id": "5d457cd8",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def generate_data(N: int) -> np.ndarray:\n",
|
||||
" x = 2*np.random.randn(N, 1)\n",
|
||||
" epsilon = 0.5*np.random.randn(N, 1)\n",
|
||||
" z = 3*x + epsilon\n",
|
||||
" return np.concatenate((x, z), axis=1).astype(np.float32)\n",
|
||||
"\n",
|
||||
"data = generate_data(N=1000)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "48433f6f",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"We can fit this line from scratch, assuming $y = f(x) = \\beta x + \\alpha + \\epsilon$, where $\\epsilon$ is a zero-mean Gaussian noise."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "6a170344-34ff-4af6-9a6c-daa44b95408e",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"We can do this in two ways. In both cases, we start assuming the knowledge of $f(x)$ and use the Bayes Theorem:\n",
|
||||
"\n",
|
||||
"$ p(\\alpha, \\beta|\\mathcal{D}) \\, \\propto p(\\mathcal{D}|\\alpha, \\beta) \\, p(\\alpha) \\, p(\\beta)$\n",
|
||||
"\n",
|
||||
"We assume $p(\\alpha)$ and $p(\\beta)$ are constants, so we just need to maximize the likelihood $p(\\mathcal{D}|\\alpha, \\beta)$.\n",
|
||||
"That likelihood is a Gaussian distribution, because we assumed $y$ is exactly equals to a line plus some Gaussian noise coming from $\\epsilon$.\n",
|
||||
"\n",
|
||||
"$p(\\alpha,\\beta|\\mathcal{D}) \\propto \\prod \\mathcal{N}(y_i|\\text{mean}=f(x_i), \\text{std. dev.}=\\sigma_\\epsilon)$\n",
|
||||
"\n",
|
||||
"This leads to minimizing $L = \\sum_i (y_i - \\beta x_i - \\alpha)^2$, as seen in the lecture (simply take $- \\log p$ and to maximize $p$, one has to minimize $-\\log p$).\n",
|
||||
"\n",
|
||||
"We can try finding out $\\beta$ and $\\alpha$ by taking the derivative of $L$ and setting it to zero:\n",
|
||||
"\n",
|
||||
"$\\frac{d L}{d\\beta} = 2 \\sum_i (y_i - \\beta x_i - \\alpha) x_i = 0$\n",
|
||||
"\n",
|
||||
"$\\frac{d L}{d\\alpha} = -2 \\sum_i (y_i - \\beta x_i - \\alpha) = 0$\n",
|
||||
"\n",
|
||||
"Re-arranging:\n",
|
||||
"\n",
|
||||
"$2 \\sum_i (y_i x_i - \\beta x_i^2 - \\alpha x_i) = 0$\n",
|
||||
"\n",
|
||||
"$2 \\sum_i (y_i - \\beta x_i - \\alpha) = 0$\n",
|
||||
"\n",
|
||||
"This leads to:\n",
|
||||
"\n",
|
||||
"$ \\beta \\sum_i x_i^2 = \\sum_i (y_i - \\alpha) x_i$\n",
|
||||
"\n",
|
||||
"$ \\alpha = 1/N \\sum_i y_i - \\beta x_i = \\mathbb{E}[y - \\beta x]$\n",
|
||||
"\n",
|
||||
"We can then substitute $\\alpha$ in $\\beta$:\n",
|
||||
"\n",
|
||||
"$ \\beta \\sum_i x_i^2 = \\sum_i y_i x_i - 1/N \\sum_i \\sum_j (y_j - \\beta x_j) x_i$\n",
|
||||
"\n",
|
||||
"And expand it:\n",
|
||||
"\n",
|
||||
"$ \\beta (\\sum_i x_i^2 - 1/N \\sum_i \\sum_j x_i x_j) = \\sum_i y_i x_i - 1/N \\sum_i \\sum_j y_j x_i$\n",
|
||||
"\n",
|
||||
"Taking $x_j$ from a sum that does not depend on $j$, we see that the second term is just the square of the sum over $x_i$:\n",
|
||||
"\n",
|
||||
"$ \\beta (\\sum_i x_i^2 - 1/N(\\sum_{i} x_i)^2 ) = \\sum_i y_i x_i - 1/N \\sum_i \\sum_j y_j x_i$\n",
|
||||
"\n",
|
||||
"Dividing by the number of points $N$, several sums reduce to averages, which make this much easier to interpret if we use some standard identities from statistics.\n",
|
||||
"For example $1/N \\sum_i x_i^2 - 1/N^2 (\\sum_i x_i)^2$ is a known expression for the variance of $X$. \n",
|
||||
"\n",
|
||||
"$ \\beta \\text{var}[X] = \\mathbb{E}[XY] - \\mathbb{E}[X]\\mathbb{E}[Y]$\n",
|
||||
"\n",
|
||||
"And the right-hand-side is a known expression for the covariance between $X$ and $Y$:\n",
|
||||
"\n",
|
||||
"$ \\beta \\text{var}[X] = \\text{cov}[X,Y]$\n",
|
||||
"\n",
|
||||
"This leads to:\n",
|
||||
"\n",
|
||||
"$\\beta = \\frac{\\text{cov}[X,Y]}{\\text{var}[X]}$\n",
|
||||
"\n",
|
||||
"This, together with the expression from above for $\\alpha$, allows us to easily fit a line in a few lines of code:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 3,
|
||||
"id": "30205402",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"cov = np.cov(data[:,0], data[:,1])[0,1]\n",
|
||||
"varX = np.std(data[:,0])**2\n",
|
||||
"beta = cov/varX\n",
|
||||
"alpha = np.mean(data[:,1] - beta*data[:,0])"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 4,
|
||||
"id": "cabdfc12-3033-4a2f-96ed-4ab1ea8312ab",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"<matplotlib.lines._AxLine at 0x7f5db8662e30>"
|
||||
]
|
||||
},
|
||||
"execution_count": 4,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAioAAAGdCAYAAAA8F1jjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABZHElEQVR4nO3de1zT1f8H8Nc2YFyEKXJPVLyHeENTUSvvooaYZZlpaX2tFO+WguWFsvCSt7x21czMrCzFjPJuGUZqqIiaKN5geEM2FBmwfb5/EIPJgG0MNtjr+Xjw+LWzz+ezI7++7tW5vI9IEAQBRERERFZIbOkOEBEREZWFQYWIiIisFoMKERERWS0GFSIiIrJaDCpERERktRhUiIiIyGoxqBAREZHVYlAhIiIiq2Vn6Q5UlkajQXp6OlxdXSESiSzdHSIiIjKAIAjIzs6Gn58fxOKyx01qfFBJT0+Hv7+/pbtBREREJrh27RoaNGhQ5vtVGlRiYmKwfft2nDt3Dk5OTujWrRsWLVqEli1baq/Jzc3FjBkzsHXrVqhUKgwYMABr166Ft7e3QZ/h6uoKoPAP6ubmViV/DiIiIjIvpVIJf39/7fd4WURVedZPaGgoRowYgcceewwFBQWYPXs2kpKSkJycDBcXFwDA+PHj8fPPP2Pjxo2QyWSYOHEixGIxjhw5YtBnKJVKyGQyKBQKBhUiIqIawtDv7yoNKg+7desWvLy8cOjQITzxxBNQKBTw9PTEli1b8OyzzwIAzp07h0cffRTx8fHo2rVrhc9kUCEiIqph1Goob96EzM+vwu/val2jolAoAADu7u4AgOPHjyM/Px99+/bVXtOqVSs0bNiwzKCiUqmgUqm0r5VKZRX3moiIiMzmxg3gxRcBZ2eDLq+27ckajQZTp05F9+7dERQUBADIyMiAg4MD6tatq3Ott7c3MjIy9D4nJiYGMplM+8OFtERERDXImjXAvn1AbKxBl1dbUImIiEBSUhK2bt1aqedERUVBoVBof65du2amHhIREVGVe+cdoHNnwMvLoMurZepn4sSJ2LVrFw4fPqyzBcnHxwd5eXnIysrSGVW5ceMGfHx89D5LKpVCKpVWdZeJiIjIHAoKALsSccPBAfj+e0ClApo3r/D2Kh1REQQBEydOxI8//oj9+/cjICBA5/2OHTvC3t4e+/bt07adP38eV69eRUhISFV2jYiIiKra4cNAq1bAqVO67f7+1jGiEhERgS1btmDHjh1wdXXVrjuRyWRwcnKCTCbDq6++iunTp8Pd3R1ubm6YNGkSQkJCDNrxQ0RERFZIowEWLSqc5tFogOHDgWPHgApqpuhTpUFl3bp1AICePXvqtG/YsAFjxowBACxfvhxisRjPPPOMTsE3IiIiqoFu3wZGjwbi4orbGjQAcnNNCirVWkelKrCOChERkfmpNQISUjNxMzsXXq6O6BzgDom4gjP1jhwBnn8eSEsrfC0SAXPnAnPmABKJzqWGfn/X+LN+iIiIyLzikuSIjk2GXJGrbfOVOWJeWCBCg3xL36DRAB9+CMyeDajVhW1eXsDXXwMlaqWZotq2JxMREZH1i0uSY/zmEzohBQAyFLkYv/kE4pLkujfcuQMMGQLMmlUcUnr2BBITKx1SAAYVIiIi+o9aIyA6Nhn61oQUtUXHJkOtKXHFv/8Wr0cRiQqnefbuBXz1jLyYgEGFiIiIAAAJqZmlRlJKEgDIFblISM0sbgwJAT74APD0LAws775baj1KZTCoEBEREQDgZnbZIaWIW+493My6r9v45ptAUhLQv7/Z+8SgQkRERAAAL1fHct/vkHYOuzdMQofN63TfEIsNLuBmLAYVIiIiAgB0DnCHr8wRpTYhCwJeTfgR27bMQgPlLfh/tBg4cKBa+sSgQkRERAAAiViEeWGBOm1uuffwyY/vY86Bz2GvKdzVI+rWzaBzesyBQYWIiIh0yJztAQDt0s9j94bJ6H/haPGbkZGFoyklDhmuSiz4RkRERACKa6gIgoCxx3ci6sAGOGgKAAB3HV1xceladJowqlr7xKBCRERE2hoqLqocfLh7OUL/jde+d9yvFSaFz4Kg8MUfGqHiUvpmxKkfIiIi0tZQKRCLEZCZpm1f33kYnh+5EOlunqVrqFQDBhUiIiLS1lDJtXfEhPAoXHfzxCvPzMXCXq+gQGJX6rrqwqkfIiIiW6ZQAAoFvFxdtE0XPfzR67VPkC+xL3V5RbVWzI0jKkRERLbqxAkgOBgID0dnX2edGioPhxQRCk9Q7hzgXq1dZFAhIiKyNYIArF1beE7PpUtAYiIk77ytraHy8FLZotfzwgKrdSEtwKBCRERkW5RKYMQIICICyMsrbOvcGZg8GaFBvlg3Khg+Mt3pHR+ZI9aNCkZokHlORDYG16gQERHZisREYPhwICWluG3qVGDRIsDBAQAQGuSLfoE+SEjNxM3sXHi5Fk73VPdIShEGFSIiotpOEIBPPgGmTAFUqsI2mQzYsAF4+ulSl0vEIoQ0rV/NndSPQYWIiKg2EwRg9Gjg66+L2zp1ArZtAwICLNcvA3GNChERUS2l1giIv5SJ804lRkcmTwb++KNGhBSAIypERES1UlySHNGxyZArciGp1wfrmifgwGP98eS4NxAqlVq6ewYTCYIgWLoTlaFUKiGTyaBQKODm5mbp7hAREVnWvXs49nUshqe64eEv+KLlsJbawVOSod/fnPohIiKqLU6fhtCpE9pFvIS26edLvV0UXKJjk6HW1IxxCgYVIiKiGkatERB/8Q52JKYh/uIdqNUa4PPPgc6dITp/HvbqAiyMW1W4kPYhAmCRwwVNxTUqRERENUjJtScA4JSXi2UH12PgP3u115zxaoIJQyMBUdm1T6r7cEFTMagQERHVEHFJcozffEI7hdP81hWs3bEQze9c016TMXIshvmGQWXnUO6zqvtwQVNx6oeIiKgGUGsERMcma0PKs6f3Yuem6dqQcs/BCXOefxvuX34G9/pupc7rKWKpwwVNxaBCRERUAySkZmqne2Yc/gof7l4Bp4LCKrNnPRtjyEvL8VXjEBy/ctcqDxc0FYMKERFRDVByTcmvLUKgkhSu3tjSbgCGjl6KS/UbaK+zxsMFTcU1KkRERDVAyTUlST7NMK/vG3hgL8WO1r30XmdthwuaikGFiIjImuXkAKtWofO06fCVOSJDkQsBwNb2oTqXiVA4YlJy7Yk1HS5oKk79EBERWZmiOin7fjyEnA6dgMhISObNrVVrTwxVpUHl8OHDCAsLg5+fH0QiEX766Sed98eMGQORSKTzExoaqv9hRERENiAuSY4ei/bjm2kL0fX5UDj/exYAULDyI4R6iGrN2hNDVenUz/3799GuXTu88sorGDZsmN5rQkNDsWHDBu1raQ06KImIiMic4pLkmLohHnP3fYKRJ3/Vtv9bvyEihkZixm2h1qw9MVSVBpWBAwdi4MCB5V4jlUrh4+NTld0gIiKyemqNgC82/Ibtm6MReDNV2/59UB/M6TceuQ6OiI5NRr9An1qx9sRQFl9Me/DgQXh5eaFevXro3bs3FixYgPr1y/7lq1QqqFQq7WulUlkd3SQiIqpSF1d/ji9WT0advAcAgAd2UszpPx7ft+mrvabojB5bCSmAhRfThoaGYtOmTdi3bx8WLVqEQ4cOYeDAgVCr1WXeExMTA5lMpv3x9/evxh4TERFVgW3b0GLKOG1ISXFvgPCXluqElCI15YwecxEJgp6jFavig0Qi/Pjjjxg6dGiZ11y6dAlNmzbF3r170adPH73X6BtR8ff3h0KhgJubm7m7TUREVPVUKtx7rCvqnE7ED617YU7/CchxcNJ76TfjutaKERWlUgmZTFbh97fFp35KatKkCTw8PJCSklJmUJFKpVxwS0RENYZaI2gXvnq4SAERcPueSncRrFQKpx9/wIJZ6/B50ych6Dn1WF+dFFtgVUHl+vXruHPnDnx9a9/2KiIisj1xSXJExyZrz+gBAIeCfMw6tBGr2/bHvWYtMS8sEKFBvpA0bYJO86fj880nIAJQcrqjNtdJqUiVrlG5d+8eEhMTkZiYCABITU1FYmIirl69inv37uGtt97C0aNHcfnyZezbtw/h4eFo1qwZBgwYUJXdIiIiqnJxSXKM33xCJ6Q0vCvHD5vfxKvHdmDtjoVQ3MrC+M0nEJckB4BadUaPuVTpGpWDBw+iV69epdpffvllrFu3DkOHDsU///yDrKws+Pn5oX///njvvffg7e1t8GcYOsdFRERUlXSmeOpIMWNbIjKUxWsqQ88fweLdK+GWlwMAUEnsMWZ4NI42agsfmSP+mNVbO1pS8lm1tU6Kod/f1baYtqowqBARkaXpm+Ip4lCQj6iDX2Ds8Vht20X3RzAxfBbOejXRttWWRbKGqpGLaYmIiGqaoikeff/V3yArA6t3LkJ7+QVt285Hn0DUgIm4L3XWudbWth0bikGFiIjICGqNgKOX7iD+4h0IgoAtCVf1hpT+/8Zjye4VkKnuAyic6onu+xq2tAsF9Ozq8XJ1LNVGDCpEREQGi0uSI3L7aWTl5Jd7XYOsDKz9KQZ2ggYAkFrPFxPDI3HGu2mpa21127GhLFqZloiIqKaIS5Ljjc0nKgwpAHC9rg+WPT4KALCrZQ+EvbyyzJAC2Oa2Y0NxRIWIiKgCao2A+TuTjbpnXddnccGjIfY066J3qgcoHEkpqqNC+jGoEBERVeDopTvIUOpf7GqnLsCsQxuR6SzDuq7Dte2CSIw9zbvqvaeukz3WvBiMrk3qcySlAgwqRERE/9FXv2RPcgYifzit93o/5U2s3rEIwennoRaJcfyRR5HgH1Th52Q9yIdYJGJIMQCDChEREfTXQqnrbF/mmpTeKQlY9vMy1M29BwBQi8RomJVhUFABuB3ZUAwqRERk88qqhaIvpNipC/DW4U14PWG7tu2qzBsR4ZE47dvc4M/kdmTDMKgQEZFNU2sERMcm662Line truncated
|
||||
"text/plain": [
|
||||
"<Figure size 640x480 with 1 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"plt.scatter(data[:,0], data[:,1])\n",
|
||||
"plt.axline((0, alpha), slope=beta, c='r', lw=2, ls='--', label=\"Line fit\")"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"id": "a1fbd00b-e7f2-4eff-90db-e1d44d104fba",
|
||||
"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.6"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 5
|
||||
}
|
||||
-104
@@ -1,104 +0,0 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "1bba0128",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"# Simple fits\n",
|
||||
"\n",
|
||||
"Here we are fitting a line from scratch.\n",
|
||||
"In the next notebook, we will do fancier fits with neural networks, but let's start with a basic problem and complicate it as we go along.\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 1,
|
||||
"id": "23feddde",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"from typing import Tuple\n",
|
||||
"\n",
|
||||
"import pandas as pd\n",
|
||||
"import numpy as np\n",
|
||||
"import matplotlib.pyplot as plt\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "bb1286f0",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"We start by generating some fake dataset, which is simple enough that we can visualize the results easily. For this reason, the dataset will contain only two variables.\n",
|
||||
"\n",
|
||||
"The simulated example data will be $f(x) = 3 x + \\epsilon$, where $\\epsilon \\sim \\mathcal{N}(\\mu=0, \\sigma=0.5)$.\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
"id": "5d457cd8",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def generate_data(N: int) -> np.ndarray:\n",
|
||||
" x = 2*np.random.randn(N, 1)\n",
|
||||
" epsilon = 0.5*np.random.randn(N, 1)\n",
|
||||
" z = 3*x + epsilon\n",
|
||||
" return np.concatenate((x, z), axis=1).astype(np.float32)\n",
|
||||
"\n",
|
||||
"data = generate_data(N=1000)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "48433f6f",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"We can fit this line from scratch, assuming $y = f(x) = \\beta x + \\alpha + \\epsilon$, where $\\epsilon$ is a zero-mean Gaussian noise.\n",
|
||||
"\n",
|
||||
"How would you do it? Feel free to use standard Python modules. Look at the solution for a simple mathematical expression for this fit with a full derivation.\n",
|
||||
"\n",
|
||||
"Tip: Look for the documentation for `numpy.linalg.lstsq`."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"id": "85040eef-3558-4531-9bce-4f29d520f86b",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": []
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"id": "5e9170af-ca77-4526-b3d2-e19e77fef44e",
|
||||
"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.6"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 5
|
||||
}
|
||||
Binary file not shown.
|
Before Width: | Height: | Size: 77 KiB |
Binary file not shown.
|
Before Width: | Height: | Size: 78 KiB |
Binary file not shown.
|
Before Width: | Height: | Size: 69 KiB |
Binary file not shown.
|
Before Width: | Height: | Size: 69 KiB |
Reference in new issue
Block a user