notebook/Unit2-2-Understanding_Experimental_Data.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# UNDERSTANDING EXPERIMENTAL DATA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This chapter is about **understanding experimental data**. \n",
    "\n",
    "We will\n",
    "\n",
    "* make extensive use of `plotting to visualize the data`, and \n",
    "\n",
    "* show how to use `linear regression` to build a model of experimental data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## The Behavior of Projectiles\n",
    "\n",
    "We decided to use one of our springs to build `a device capable of launching a projectile`. \n",
    "\n",
    "We used the device **four** times to fire a projectile at **a target** 1080 inches from the **launching point**. \n",
    "\n",
    "Each time, we measured `the height of the projectile` at various `distances from the launching point`.\n",
    "\n",
    "* The `launching point` and the `target` were at `the same height`, which we treated as `0.0` in our measurements.\n",
    "\n",
    "![projectile](./img/projectile.jpg)\n",
    "\n",
    "The data was stored in a file **./data/launcherData.csv**\n",
    "\n",
    "* The `first column` contains `distances of the projectile from the target`. \n",
    "\n",
    "* The `other columns` contain `the height of the projectile at that distance` for each of the four trials. \n",
    "\n",
    "All of the measurements are in **inches**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load ./data/launcherData.csv\n",
    "Distance,Trial1,Trial2,Trial3,Trial4\n",
    "1080,0.0,0.0,0.0,0.0\n",
    "1044,2.25,3.25,4.5,6.5\n",
    "1008,5.25,6.5,6.5,8.75\n",
    "972,7.5,7.75,8.25,9.25\n",
    "936,8.75,9.25,9.5,10.5\n",
    "900,12.0,12.25,12.5,14.75\n",
    "864,13.75,16.0,16.0,16.5\n",
    "828,14.75,15.25,15.5,17.5\n",
    "792,15.5,16.0,16.6,16.75\n",
    "756,17.0,17.0,17.5,19.25\n",
    "720,17.5,18.5,18.5,19.0\n",
    "540,19.5,20.0,20.25,20.5\n",
    "360,18.5,18.5,19.0,19.0\n",
    "180,13.0,13.0,13.0,13.0\n",
    "0,0.0,0.0,0.0,0.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following code was used to plot :\n",
    "\n",
    "* the **mean altitude of the projectile(Y)** against the **distance from the point of launch(X)**.\n",
    "\n",
    "* the best `linear` and `quadratic` fits to the points.\n",
    "\n",
    "```python\n",
    "d2h={'d':None,'h':[[],[],[],[]]}\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import csv\n",
    "\n",
    "def get_trajectory_data(file_name):\n",
    "    csvfile = open(file_name, 'r')\n",
    "    csvdata = csv.DictReader(csvfile)\n",
    "    # fields\n",
    "    fields = csvdata.fieldnames\n",
    "    \n",
    "    trials = len(fields)-1\n",
    "    d2hs = {'d': [], 'h': [[] for i in range(trials)]}\n",
    "    \n",
    "    # values\n",
    "    rows=0\n",
    "    for row in csvdata:\n",
    "        # the distance in first column\n",
    "        d2hs['d'].append(float(row[fields[0]]))\n",
    "        for i in range(trials):\n",
    "            cur_height = float(row[fields[i+1]])\n",
    "            d2hs['h'][i].append(cur_height)\n",
    "        rows+=1    \n",
    "    \n",
    "    csvfile.close()\n",
    "    \n",
    "    return d2hs, trials,rows\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "file_name='./data/launcherData.csv'\n",
    "\n",
    "d2hs,trials,rows = get_trajectory_data(file_name)\n",
    "\n",
    "print(f\"trials：{trials}, rows={rows}\")\n",
    "\n",
    "for key in d2hs.keys():\n",
    "    print(key,end=\"\\t\")\n",
    "print(\"\")\n",
    "\n",
    "for i in range(len(d2hs['d'])):\n",
    "    print(d2hs['d'][i],end=\"\\t\")\n",
    "    \n",
    "    for j in range(len(d2hs['h'])):\n",
    "        cur_column=d2hs['h'][j]\n",
    "        print(cur_column[i],end=\"\\t\")\n",
    "    print(\"\")   \n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1 Linear Regression to Find a Fit\n",
    "\n",
    "Whenever we fit any curve (including a line) to data we need some way to decide `which curve is the best fit for the data`. \n",
    "\n",
    "This means that we need to define <b style=\"color:blue\">an objective function</b> that provides `a quantitative assessment of how well the curve fits the data`. \n",
    "\n",
    "Once we have such `a function`, finding the best fit can be formulated as\n",
    "\n",
    "* finding a curve that <b style=\"color:blue\">minimizes (or maximizes)</b> the value of that function, i.e., as **an optimization problem** \n",
    "\n",
    "The most commonly used objective function is called <b style=\"color:blue\">least squares(最小二乘)</b>,\n",
    "\n",
    "The objective function is then defined as\n",
    "\n",
    "$$\\sum_{i=0}^{len(observed)-1}(observed[i]-predicted[i])^2$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Numpy.polyfit**\n",
    "\n",
    "http://docs.scipy.org/doc/numpy/reference/generated/numpy.polyfit.html\n",
    "\n",
    "`Numpy` provides a function, `polyfit`, that finds the best Least squares polynomial fit.\n",
    "```python\n",
    "numpy.polyfit(x, y, deg, rcond=None, full=False, w=None, cov=False)\n",
    "```\n",
    "The algorithm used by `polyfit` is called <b style=\"color:blue\">linear regression</b>.\n",
    "\n",
    "Fit a polynomial of degree `deg` to points (x, y). \n",
    "\n",
    "$$p(x) = p[0] * x^{deg} + ... + p[deg]$$\n",
    "\n",
    "Returns a vector of **coefficients**  <b style=\"color:blue\">$p$</b> that minimises the squared error.\n",
    "\n",
    "\n",
    "\n",
    "```python\n",
    "numpy.polyfit(observedXVals, observedYVals, n)\n",
    "numpy.polyfit(observedXVals, observedYVals, 1) # y = ax + b\n",
    "numpy.polyfit(observedXVals, observedYVals, 2) # y = ax^2 + bx+c\n",
    "```\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "x = np.array([0.0, 1.0, 2.0, 3.0,  4.0,  5.0])\n",
    "y = np.array([0.0, 0.8, 0.9, 0.1, -0.8, -1.0])\n",
    "fit = np.polyfit(x, y, 1)\n",
    "fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fit = np.polyfit(x, y, 2)\n",
    "fit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**numpy.polyval(p, x)**\n",
    "\n",
    "Evaluate a polynomial at specific values.\n",
    "\n",
    "This function takes two arguments:\n",
    "* <b style=\"color:blue\">$p$</b>: a sequence of polynomial coefficients  \n",
    "* <b style=\"color:blue\">$x$</b>: a sequence of values at which the polynomial is to be evaluated.\n",
    "> https://numpy.org/doc/stable/reference/generated/numpy.polyval.html\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v=np.array([1.1,2.3])\n",
    "z= np.polyval(fit,v)\n",
    "z"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Elementwise operations** on arrays \n",
    "\n",
    "linear fit : `a,b = np.polyfit(x,y, 1)`\n",
    "  \n",
    "quadratic fit: `a,b,c = np.polyfit(x, y, 2)`\n",
    "\n",
    "```python\n",
    " y = a*x^2 + b*x+c\n",
    "```\n",
    "x,y : the array\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a,b,c=np.polyfit(x, y, 2)\n",
    "v=np.array([1.1,2.3])\n",
    "#z=[]\n",
    "#for item in v:\n",
    "#   z.append(a*(item**2)+b*item+c\n",
    "z=a*(v**2)+b*v+c\n",
    "z"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2 Using numpy.polyfit to fit the experimental data\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def get_mean_height(d2hs,trials,rows):\n",
    "    distances = np.array(d2hs['d'])\n",
    "    \n",
    "    mean_heights = np.empty(rows)\n",
    "    \n",
    "    for i in range(rows):\n",
    "        mean_heights[i] = d2hs['h'][0][i]\n",
    "        for j in range(trials-1):\n",
    "            mean_heights[i] += d2hs['h'][j+1][i]\n",
    "            mean_heights[i] /= trials\n",
    "   \n",
    "    return distances,mean_heights \n",
    "\n",
    "def process_trajectories(distances,  mean_heights):\n",
    "    plt.title(f'Trajectory of Projectile (Mean Height of {trials} Trials)')\n",
    "    plt.xlabel('Inches from Launch Point')\n",
    "    plt.ylabel('Mean Height of the projectile')\n",
    "\n",
    "    # the experimental data\n",
    "    plt.plot(distances, mean_heights, 'ro')\n",
    "\n",
    "    # Linear Fit\n",
    "    fit = np.polyfit(distances, mean_heights, 1)\n",
    "    altitudes=np.polyval(fit,distances)\n",
    "\n",
    "    plt.plot(distances, altitudes, 'g-*', label=\"Linear Fit\")\n",
    "\n",
    "    # Quadratic Fit\n",
    "    fit = np.polyfit(distances, mean_heights, 2)\n",
    "    altitudes=np.polyval(fit,distances)\n",
    "  \n",
    "    plt.plot(distances, altitudes, 'b:x', label=\"Quadratic Fit\")\n",
    "   \n",
    "    plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "distances,mean_heights =get_mean_height(d2hs,trials,rows)\n",
    "process_trajectories(distances,mean_heights)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A quick look at the plot on the right makes it quite clear that\n",
    "\n",
    "* a **quadratic** fit is far better than a linear one\n",
    "\n",
    "<b style=\"color:blue;font-size:120%\">how bad a fit is the line and how good is the quadratic fit?</b>\n",
    "\n",
    "* **Coefficient of Determination(决定系数)**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "####  something about plotting lines with matplotlib.pyplot\n",
    "\n",
    "* [Unit2-1-Matplotlib](./Unit2-1-Matplotlib.ipynb)\n",
    "\n",
    "The function `process_trajectories` plot lines with a legend(图例说明)  \n",
    "\n",
    "**line styles**\n",
    "\n",
    "* the linear fit for the data.:  \n",
    "```python\n",
    "plt.plot(distances, altitudes, 'g-*', label = 'Linear Fit')\n",
    "```\n",
    "\n",
    "* the quadratic fit for the data.: \n",
    "```python\n",
    "plt.plot(distances, altitudes, 'b:x', label = 'Quadratic Fit')\n",
    "```\n",
    "\n",
    "**Legend of matplotlib.pyplot**\n",
    "\n",
    "* **label**：a legend of line      \n",
    "       \n",
    "* **Place a legend on the fig**\n",
    "\n",
    "  *  `plt.legend(loc ='best')` \n",
    "\n",
    "     * Location String： **`'best'`**   \n",
    "\n",
    "> reference: [matplotlib.pyplot.legend](https://matplotlib.org/devdocs/api/_as_gen/matplotlib.pyplot.legend.html): \n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def r_squared(measured, predicted):\n",
    "    \"\"\"Assumes measured a one-dimensional array of measured values\n",
    "               predicted a one-dimensional array of predicted values\n",
    "       Returns coefficient of determination\"\"\"\n",
    "    # RSS: residual sum of squares\n",
    "    estimate_error = ((predicted - measured)**2).sum()\n",
    "    \n",
    "    # the mean of the observed values.\n",
    "    mean_measured = measured.mean()\n",
    "    # TSS: total sum of squares\n",
    "    variability = ((measured - mean_measured)**2).sum()\n",
    "    \n",
    "    return 1.0 - estimate_error/variability"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def process_trajectories( distances , mean_heights):\n",
    "    plt.title(f'Trajectory of Projectile (Mean Height of {trials} Trials)')\n",
    "    plt.xlabel('Inches from Launch Point')\n",
    "    plt.ylabel('Mean Height of the projectile')\n",
    "\n",
    "    # the experimental data\n",
    "    plt.plot(distances, mean_heights, 'ro')\n",
    "\n",
    "    # Linear Fit\n",
    "    fit= np.polyfit(distances, mean_heights, 1)\n",
    "    altitudes =np.polyval(fit,distances)\n",
    "\n",
    "    r_line = r_squared(mean_heights, altitudes)\n",
    "    plt.plot(distances, altitudes, 'g-*', label=f\"Linear Fit, $R^2$= {r_line:.2f}\")\n",
    "\n",
    "    # Quadratic Fit\n",
    "    fit = np.polyfit(distances, mean_heights, 2)\n",
    "    altitudes =np.polyval(fit,distances)\n",
    "\n",
    "    r_quad = r_squared(mean_heights, altitudes)\n",
    "    plt.plot(distances, altitudes, 'b:x',label=f\"Quadratic Fit, $R^2$= {r_quad:.2f}\")\n",
    "    plt.legend()\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "process_trajectories(distances , mean_heights)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This tells us that \n",
    "\n",
    "* 1 less than 2% of the variation in the measured data can be explained by the linear model,\n",
    "\n",
    "* 2 more than **98%** of the  variation can be explained by the quadratic model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Matplotlib.pyplot: Formatting text: LaTeX\n",
    "```python\n",
    "   plt.plot(distances, altitudes, 'g-*',label=f\"Linear Fit, $R^2$= {r_line:.2f}\")\n",
    "```\n",
    "Matplotlib has great support for LaTeX. All we need to do is to use dollar signs encapsulate LaTeX in any text (legend, title, label, etc.).\n",
    "\n",
    "* For example, `$R^2$`. $R^2$"
   ]
  }
 ],
 "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.5"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "165px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  },
  "widgets": {
   "state": {},
   "version": "1.1.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}