ScipyExercises.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "toc": true
   },
   "source": [
    "<h1>Table of Contents<span class=\"tocSkip\"></span></h1>\n",
    "<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#Exercises-for-Numpy,-Scipy-and-Matplotlib\" data-toc-modified-id=\"Exercises-for-Numpy,-Scipy-and-Matplotlib-1\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span>Exercises for Numpy, Scipy and Matplotlib</a></span><ul class=\"toc-item\"><li><span><a href=\"#Notebook-examples\" data-toc-modified-id=\"Notebook-examples-1.1\"><span class=\"toc-item-num\">1.1&nbsp;&nbsp;</span>Notebook examples</a></span></li><li><span><a href=\"#Data-analysis-1\" data-toc-modified-id=\"Data-analysis-1-1.2\"><span class=\"toc-item-num\">1.2&nbsp;&nbsp;</span>Data-analysis-1</a></span></li><li><span><a href=\"#Data-analysis-2\" data-toc-modified-id=\"Data-analysis-2-1.3\"><span class=\"toc-item-num\">1.3&nbsp;&nbsp;</span>Data-analysis-2</a></span></li><li><span><a href=\"#Solving-differential-equation\" data-toc-modified-id=\"Solving-differential-equation-1.4\"><span class=\"toc-item-num\">1.4&nbsp;&nbsp;</span>Solving differential equation</a></span></li><li><span><a href=\"#Data-fitting---half-time-of-radioactive-source\" data-toc-modified-id=\"Data-fitting---half-time-of-radioactive-source-1.5\"><span class=\"toc-item-num\">1.5&nbsp;&nbsp;</span>Data fitting - half time of radioactive source</a></span></li><li><span><a href=\"#Polynom-fit----which-polynom-degree-is-needed?\" data-toc-modified-id=\"Polynom-fit----which-polynom-degree-is-needed?-1.6\"><span class=\"toc-item-num\">1.6&nbsp;&nbsp;</span>Polynom fit -- which polynom degree is needed?</a></span></li></ul></li></ul></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises for Numpy, Scipy and Matplotlib\n",
    "\n",
    "### Notebook examples\n",
    "* Check in detail the examples presented\n",
    "* Test variations (data creation, plot options, ...)\n",
    "\n",
    "### Data-analysis-1\n",
    "The file [rohr1.dat](http://www.etp.physik.uni-muenchen.de/kurs/comp14/uebungen/source/rohr1.dat) contains a list of measurements of wire-positions of drift tubes used in the ATLAS\n",
    "Muon-Chamber Sytem   \n",
    " Read the numbers with:\n",
    " ```Ipython\n",
    "data = numpy.loadtxt('rohr1.dat')\n",
    "```\n",
    "1. Determine mean and standard-deviation *(Hint: numpy-Functions)*\n",
    "\n",
    "1. Fill the values in a histogram and plot it.\n",
    "\n",
    "In a similar way read the (x,y) coordinates of file \n",
    "[rohr2.dat](http://www.etp.physik.uni-muenchen.de/kurs/comp14/uebungen/source/rohr2.dat) using\n",
    "```Ipython\n",
    "x,y = numpy.loadtxt('rohr2.dat',unpack=True)\n",
    "```\n",
    "\n",
    "Determine for both x and y mean and standard-deviation as well as the correlation.\n",
    "\n",
    "Visualize the data:\n",
    "1. 1D histogram of both x and y\n",
    "1. (x,y) point plot\n",
    "1. 2D histogram\n",
    " \n",
    " \n",
    " \n",
    "### Data-analysis-2\n",
    "The file [faithful.csv](https://people.sc.fsu.edu/~jburkardt/data/csv/faithful.csv) contains measurement data of the **Old Faithful Geysir** , i.e. duration of the eruption and time since the last eruption.\n",
    "\n",
    "Download the file ( `wget http://...`) and read the data with numpy:\n",
    "```Ipython\n",
    "data=numpy.loadtxt('faithful.csv',delimiter=',',skiprows=1)\n",
    "```\n",
    "(Why  the options in in loadtxt(...)?)\n",
    "\n",
    "1. Determine again mean and standard-deviation for duration and wait-time\n",
    "1. Fill histograms for both values and plot it. Is it compatible with Gaussian distribution?\n",
    "1. Are there correlations between wait-time and duration or the duration of sub-sequent eruptions?\n",
    "(Make (x,y) plots for both)\n",
    "\n",
    "\n",
    "### Solving differential equation\n",
    "\n",
    "(a) Extend the example for free fall and introduce a term for air-friction:\n",
    "\n",
    "$$F_R = \\frac{1}{2} c_W \\rho_{Luft} A v^2 $$\n",
    "\n",
    "This results in a non-leinear DE which can only be solved numerically: \n",
    "\n",
    "$$ y'' = - g + \\frac{1}{2} c_W \\rho_{Luft} A / m_{fb} \\times y'^{ 2}$$\n",
    "\n",
    "Take numbers for a football:\n",
    "\n",
    "        # Luftwiderstandskraft Fussball\n",
    "        rhol = 1.184 # Luftdichte (kg/m^3)\n",
    "        rfb  =  0.11 # Fussballradius (m)\n",
    "        cwfb =  0.45 # Kugel CW-Wert\n",
    "        mfb  =  0.43 # Masse Fussball (kg)\n",
    "        g = 9.81 # Erdbeschleunigung\n",
    "        v = y[1] # Geschwindigkeit\n",
    "        fR = 0.5 * cwfb * rhol * math.pi * rfb**2 * v**2\n",
    "\n",
    "\n",
    "\n",
    "(b) The next step would be the extension to a trajectory of a thrown object with air resistance (\"Schiefer Wurf mit Luftreibung\"), the combined movement in horizontal and vertical direction.\n",
    "\n",
    "How far will a football fly which is shot under 45 degree angle?\n",
    "\n",
    "(The aerodynamic of a football is rather complex, see an interesting discussion at \n",
    "http://www.weltderphysik.de/gebiet/leben/fussball/coca-cola-formel/  )\n",
    "\n",
    "\n",
    "### Data fitting - half time of radioactive source\n",
    "In an experiment you measure the activity of a radioactive source with a counter. The counter takes readings every 30 seconds and shows the integrated number of counts during that period.\n",
    "The list of numbers below is the result of such a measurement. \n",
    "What is the lifetime of the radioactive source? *Hint: Fit an exponential* \n",
    "\n",
    "```IPython\n",
    "counts = [728, 612, 581, 460, 381, 324, 305, 211, 203, 166, 139, 110, 120, 92, 85, 77, 48, 53, 47, 38, 30, 20, 22, 17, 17, 16, 18, 11, 9]\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Polynom fit -- which polynom degree is needed?\n",
    "\n",
    "Below is notebook with test-data and polynom fit, check which degree of polynom is sensible"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-2.07592094 19.23496798]\n",
      "Chi2, n, chi2/ndf: 105.59783156674631 1 7.039855437783087\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEICAYAAABS0fM3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3XmcFNW5//HPw76DyiICMqggKioygPDTK6hBjPv1aoJ3VFQi7luMwVySG4zhusS4JVGjUVAhoCExel3iAozEKEkENxRRLgKyKItCBkaYhef3R1UzPTM9CzM909VT3/frVa/uOqeWp7tn6ulzquq0uTsiIhJfzTIdgIiIZJYSgYhIzCkRiIjEnBKBiEjMKRGIiMScEoGISMwpEUiVzGyKmc3IdBxNmZntb2bbzKx5pmPJBP2NRYMSQQyY2Uoz+yY84HxpZtPMrEOm40qH8EBSHL62LWb2ppmNzHRcteXuq929g7uXpnvbZuZmtj18bxLTD9O9H8l+SgTxcbq7dwCGAMOAH2c4nnR6KnxtXYH5wB8aYidm1qIhttvAjgwTTWK6M9MBSfQoEcSMu68FXgIGAZjZfmb2nJl9ZWbLzezSVOuZ2Qtmdk2FsvfN7KzwuZvZ5Wb2qZl9bWa/MTML65qZ2Y/NbJWZbTCzJ8ysc1iXE657sZl9Hq57uZkNC7e/xcx+XcvXVgLMBHqZWbekOE8zs3eTWgxHJNUNMbN3zKzAzP5gZk+Z2c/DutFmtsbMJpnZF8C0WmxvkpmtDbe3zMxODMuHm9nbZvavsFV2d4XX36KmzyNs/Twdvn8FZvahmQ2tzXtTkZm9aGa/TJp/ysweC58faGbzzGyzmW0ys5lm1iVp2ZVmdlP4+Ww3s0fNrIeZvRTG9ZqZ7VXh9U00s3Vmtt7MbqwmrhHhe7rFzN4zs9F1eX2yh9xdUxOfgJXAt8LnfYAPgVvD+deBB4A2wGBgI3BiWDcFmBE+/w7w96RtHglsBlqF8w48D3QB9g+3c3JYdwmwHDgA6AD8CXgyrMsJ130ojOEkYAfwZ6A70AvYAIyq4rUlx9gKuB3YBLQIy4aE6x8NNAfGh+9H63D5VcB1QEvgbKAI+Hm47migBLgjXL5tDds7GPgc2C/ptR0YPn8LuCB83gEYUeH1t6jl57EDOCXc923Awmo+dwcOqqJu3/B1nADkASuAjmHdQcCY8DV1AxYA91b4e1oI9Ej6fBYDR4XrzAN+WuH1zQLaA4eHr+lbKT6/XgR/U6cQfEkdE853y/T/UFOfMh6Apkb4kIN/3G3AlvDA90B4UOsDlCYOAOGytwHTw+fJ/6Stga+A/uH8XcADSes5cGzS/NPAzeHzucCVSXUHA8VAi6QDRa+k+s3Ad5Pm/whcX8Vrm0Jw8N4SvpbNwOik+gcJk15S2TJgFHAcsBawpLo3KJ8IioA2tdzeQeFB8VtAywrLLABuAbpWKE+8/ha1/DxeS6o7FPimms/dgX+F701iGptUfzZB4tqU/Nml2M5ZwDsV/p7yKnw+DybNXwP8ucLrG5hUfyfwaIq/sUmEXxCSln0ZGJ/p/6GmPqlrKD7Ocvcu7t7X3a9092+A/YCv3L0gablVBN/MynH3nQQH9/PNrBlwHvBkhcW+SHpeSPDNl3A/qyrsowXBN8qEL5Oef5NivrqT20+7e5dwe0uA3KS6vsCNYVfDFjPbQnDA3S+c1np4xAl9XmHbG919R2225+7LgesJDm4bzGy2me0XrjcBGAB8bGb/NLPTUryO2nweFd/jNjWcuxgSfu6J6eWkuucJWhbL3P2NRKGZdQ9jX2tm/wJmEJx/Sbann1fy+7qK4LVW1Bc4t8J7eyzQs5rXJ2mgRBBv64C9zaxjUtn+BN+SU3mcoBvhRKDQ3d/ag/30rbCPEsofPOrN3TcBlwFTzCxx8PgcmFrhYNjO3WcB6wnOJ1jSZvpU3GyF+eq2h7v/3t2PJXi9TtCthLt/6u7nEXR33QHMMbP2Fba9p59HfU0FlgI9zey8pPLbwtiPcPdOwPmApVh/TyS/r/sTvNaKPidoESS/t+3d/fZ67ltqoEQQY+7+OfAmcJuZtQlPek4gOOGaavm3gF3AL6ncGqjOLOAGM+tnwWWr/0NwpU9JvV5A6hg/JuhOSFwm+QhwuZkdbYH2ZnZqeLB9i6Ar5moza2FmZwLDa9hFldszs4PN7AQza03Ql/9NuH3M7Hwz6+buuwi6aEjUJcW+R59HfZjZccDFwIXh9CszS7Q8OhJ2JYZlN6Vhlz8xs3Zmdli436dSLDMDON3MxppZ8/A9GG1mvdOwf6mGEoGcR9CPuw54huAk36vVLP8EwQm/PbkJ6DGCxLEA+IzgIHlNtWvUzy+AiWbW3d3fBi4Ffg18TXDS+iIAdy8i6CefQHBwPp+gu2RnVRuubnsE51ESJ6u/IPj2/19h3cnAh2a2DbgPGFehyylhTz+Pmrxn5e8juNfMOhF8jle7+9qwW+hRYFrYOrqF4KT4VuAFgpP79fU6wXs1F7jL3V+puECYCM8keM82ErQQbkLHqQZn5btHRapnZhcCE8PujybHzP4OPOTu0zIdS1NgZjkEyb9lQ7QAJT2UaaXWzKwdcCXwcKZjSRczG2Vm+4ZdQ+OBI4C/ZDoukcakRCC1YmZjCZrrXwK/z3A46XQw8B5BN8iNwDnuvj6zIYk0LnUNiYjEnFoEIiIxlxWDaHXt2tVzcnIqlW/fvp327Steih0dUY8Poh9j1OOD6McY9fhAMaZDqvgWLVq0yd27VbFKmUzf2lybKTc311OZP39+yvKoiHp87tGPMerxuUc/xqjH564Y0yFVfMDbriEmRESkJkoEIiIxp0QgIhJzWXGyWESyX3FxMWvWrGHHjsoja3Tu3JmlS5dmIKrai3KMbdq0ofzYiXtGiUBEGsWaNWvo2LEjOTk5lQ5aBQUFdOzYsYo1oyGqMbo7mzdvrtcVTeoaEpFGsWPHDvbZZ596fXOVysyMffbZh+bNm9d5G0oEItJolAQaRn3fVyUCEZGYUyIQEcmg//7v/+a1115LWbdsWTA1NJ0sFpFY2n1XbbO6fR8uKSmhRYv6H0J/9rOf1Xsb9aUWgYjExsqVKznkkEO48sorGTJkCE8++SQjR45kyJAhnHvuuWzbtg2AF198kYEDB3Lsscdy7bXXctpppwEwZcoUJk6cyEknncSFF15IaWkpN910E8OGDeOII47gt7/9LQDr16/nuOOOY/DgwQwaNIi//vWvlJaWctFFFzFo0CAOP/xw7rnnHgAuuugi5syZA8DcuXM56qijOPzww7nkkksoKgp+LC8nJ4ef/vSnDBkyhMMPP5yPP/44re+LWgQikhmjR+9+2ra0FJo3h+98B668EgoL4ZRTKq9z0UXBtGkTnHNO+br8/FrtdtmyZUybNo2f/exnnH322bz22mu0b9+eO+64g7vvvpsf/vCHXHbZZSxYsIB+/fpx3nnnlVt/0aJFvPHGG7Rt25aHH36Yzp07889//pOdO3dyzDHHcNJJJ/GnP/2JsWPHMnnyZEpLSyksLOTdd99l7dq1LFmyBIAtW7aU2+6OHTu46KKLmDt3LgMGDODCCy9k1qwHGT/+egC6du3K4sWLeeCBB7jrrrv43e9+V6vXWxtqEYhIrPTt25cRI0awcOFCPvroI4455hgGDx7M448/zqpVq/j444854IAD6NevH0ClRHDGGWfQtm1bAF555RWeeOIJBg8ezNFHH83mzZv59NNPGTZsGNOmTWPKlCl88MEHdOzYkQMOOIAVK1ZwzTXX8Je//IVOnTqV2+6yZcvo168fAwYMAGD8+PG8/faC3fVnn302ALm5uaxcuTKt74laBCKSGUnf4L+peLNWu3bVf8Pv2rXWLYCKEjdeuTtjxoxh1qxZ5erfeeedWq2f2MavfvUrxo4dW2m5BQsW8MILL3DBBRdw0003ceGFF/Lee+/x8ssv85vf/Iann36axx57rNy2qtO6dWsAmjdvTklJen/+WS0CEYmlESNG8Le//Y3ly5cDUFhYyCeffMLAgQNZsWLF7m/dTz31VJXbGDt2LA8++CDFxcUAfPLJJ2zfvp1Vq1bRvXt3Lr30UiZMmMDixYvZtGkTu3bt4j/+4z+49dZbWbx4cbltDRw4kJUrV+6O58knn2TYsFEN8MorU4tARGKpW7duTJ8+nfPOO4+dO4OTsj//+c8ZMGAADzzwACeffDJdu3Zl+PDhVW7je9/7HitXrmTIkCG4O926dePPf/4z+fn5/OIXv6Bly5Z06NCBJ554grVr13LxxReza9cuAG677bZy22rTpg3Tpk3j3HPPpaSkhGHDhjFu3OUN9wYkq82PFmR60g/TNJyoxxj1+NyjH2NU4vvoo4+qrPvXv/7ViJHUrKCgwN3dd+3a5VdccYXffffdGYnx44+DqTYWL15cqQz9MI2ISN088sgjDB48mMMOO4ytW7dy2WWXZTqkBqWuIRGRCm644QZuuOGGcmUFBQUZiqbhqUUgIhJzSgQiIjGnRCAiEnNKBCISWaNHlxuJQhqIEoGISAodOnTIdAiNRolARCJp5kxYuBBefx1ycoJ5aRhpSQRm1sXM5pjZx2a21MxGmtneZvaqmX0aPu4VLmtmdr+ZLTez981sSDpiEJGmY+ZMmDgRwht+WbUqmK9PMpg0aRIPPPDA7vkpU6Zwyy23cOKJJ+4e3vnZZ5+ttF5+fv7uYagBrr76aqZPnw4EI5GOGjWK3Nxcxo4dy/r16wG4//77OfTQQzniiCMYN25c3YNuJOlqEdwH/MXdBwJHAkuBm4G57t4fmBvOA3wb6B9OE4EH0xSDiDQRkycHI1EnKywMyutq3Lhx5cYNevrpp7n44ot55plnWLx4MfPnz+fGG2+scfC3hOLiYq655hrmzJnDokWLuOSSS5gcBnj77bfzzjvv8P777/PQQw/VPehGUu8bysysE3AccBGAuxcBRWZ2JjA6XOxxIB+YBJwJPBHe/rwwbE30dPf19Y1FRJqG1av3rLw2jjrqKDZs2MC6devYuHEje+21Fz179uSGG25gwYIFNGvWjLVr1/Lll1+y77771ri9ZcuWsWTJEsaMGQNAaWkpPXv2BOCII44gLy+Ps846i7POOqvuQTeSdNxZfACwEZhmZkcCi4DrgB6Jg7u7rzez7uHyvYDPk9ZfE5aVSwRmNpGgxUCPHj3ITzHk7LZt21KWR0XU44Poxxj1+CD6MUYlvs6dO1d5d25paWm5ut692/P555U7LHr33kVBwfY6x3D66aczY8YMNmzYwFlnncWjjz7K+vXryc/Pp2XLlgwaNIhNmzbtHmq6oKCAoqIiioqKdsdYUFDAjh072LZtGwMHDmTu3Lnl9lFQUMDs2bP529/+xosvvsgtt9zCP/7xjzr9rGVJSdtwm9/UuKy71/1zrs2ARNVNwFCgBDg6nL8PuBXYUmG5r8PHF4Bjk8rnArnV7UODzjWcqMcY9fjcox9jVOLbk0HnZsxwb9fOHcqmdu2C8vpYsmSJjxw50vv37+/r1q3ze++916+++mp3d583b54D/tlnn7m7e/v27d3dffXq1d63b1/fuHGjb9myxXNycnzatGm+c+dOP/DAA/3NN990d/eioiJfsmSJl5aW7t5GUVGRd+/e3b/++us6xdtYg86lo0WwBljj7n8P5+cQnA/4MtHlY2Y9gQ1Jy/dJWr83sC4NcYhIE5GXFzxOmBCcMO7bF6ZOLSuvq8MOO4yCggJ69epFz549ycvL4/TTT2fo0KEMHjyYgQMHVlqnT58+fOc732HkyJEcfPDBHHXUUQC0atWKOXPmcO2117J161ZKSkq4/vrrGTBgAOeffz5bt27F3bnhhhvo0qVL/QJvYPVOBO7+hZl9bmYHu/sy4ETgo3AaD9wePiZOxz8HXG1ms4Gjga2u8wMiUkFeHjzySPA8nT1bH3zwwe7nXbt25a233kq5XOKH7AHuvPNOfvKTn5T/FTVg8ODBLFiwoOKqvPHGG2mKtnGka/TRa4CZZtYKWAFcTHBF0tNmNgFYDZwbLvsicAqwHCgMlxURqSQCpzZiIS2JwN3fJThXUNGJKZZ14Kp07FdEROpPdxaLSKPxWl6jL3umvu+rEoGINIo2bdqwefNmJYM0c3c2b95MaWlpnbehXygTkUbRu3dv1qxZw8aNGyvV7dixgzZt2mQgqtrLRIxffBE8hr93X6U2bdqwfXvd769QIhCRRtGyZUv69euXsi4/P3/3ZZlRlYkYr7gise+al121alWd96OuIRGRmFMiEBGJOSUCEZGYUyIQEYk5JQIRkQhqzF9oUyIQEYmYhviFtuooEYiIRExD/EJbdZQIREQipiF+oa06SgQiIhGz//57Vl5fSgQiIhEzdSq0a1e+rF27oLwhKBGIiERMXh48/DC0bh3M9+0bzNf3F9qqorGGREQiqKF+oS0VtQhERGJOiUBEJOaUCEREYk6JQEQk5pQIRERiTolARCTmlAhERGJOiUBEJOaUCEREYk6JQEQk5pQIRERiTolARCTmlAhERGJOiUBEJOaUCEREYk6JQEQk5pQIRERiTolARCTmlAhERGJOiUBEJOaUCEREYi5ticDMmpvZO2b2fDjfz8z+bmafmtlTZtYqLG8dzi8P63PSFYOIiOy5dLYIrgOWJs3fAdzj7v2Br4EJYfkE4Gt3Pwi4J1xOREQyJC2JwMx6A6cCvwvnDTgBmBMu8jhwVvj8zHCesP7EcHkREckAc/f6b8RsDnAb0BH4AXARsDD81o+Z9QFecvdBZrYEONnd14R1/wcc7e6bKmxzIjARoEePHrmzZ8+utN9t27bRoUOHesffUKIeH0Q/xqjHB9GPMerxgWKsyvXXDwbg3nvfrXHZVPEdf/zxi9x9aI0ru3u9JuA04IHw+WjgeaAbsDxpmT7AB+HzD4HeSXX/B+xT3T5yc3M9lfnz56csj4qox+ce/RijHp979GOMenzuirEqo0YFU22kig9422txHG9Ri6RUk2OAM8zsFKAN0Am4F+hiZi3cvQToDawLl18TJoY1ZtYC6Ax8lYY4RESkDup9jsDdf+Tuvd09BxgHzHP3PGA+cE642Hjg2fD5c+E8Yf28MHOJiEgGNOR9BJOA75vZcmAf4NGw/FFgn7D8+8DNDRiDiIjUIB1dQ7u5ez6QHz5fAQxPscwO4Nx07ldEROpOdxaLiMScEoGISMyltWtIRETSJz+/cfajFoGISMwpEYiIxJwSgYhIzCkRiIjEnBKBiEjMKRGIiMScEoGISMwpEYiIxJwSgYhIzCkRiIjEnBKBiEjMKRGIiKTJ6NHBlG2UCEREYk6JQEQk5pQIRERiTolARCTmlAhERGJOiUBEJOaafCLI1su5REQaS5NPBCIiUj0lAhGRmFMiEBGJOSUCEZGYUyIQEYk5JQIRkZhTIhARSYOZM2HhQnj9dcjJCeazhRKBiEg9zZwJEyfCzp3B/KpVwXy2JAMlAhGRepo8GQoLy5cVFgbl2UCJQESknlav3rPyqGnSiSCb++xEJHvsv/+elUdNk00E2d5nJyLZY+pUaNeufFm7dkF5NmiyiSDb++xEJHvk5cHDD0Pr1sF8377BfF5eZuOqrRaZDqChZHufnYhkl7w8eOSR4Hl+fkZD2WNNtkWQ7X12IiKNpd6JwMz6mNl8M1tqZh+a2XVh+d5m9qqZfRo+7hWWm5ndb2bLzex9MxtS3xhSyfY+OxGRxpKOFkEJcKO7HwKMAK4ys0OBm4G57t4fmBvOA3wb6B9OE4EH0xBDJdneZyci0ljqfY7A3dcD68PnBWa2FOgFnAmMDhd7HMgHJoXlT7i7AwvNrIuZ9Qy3k1bZ3GcnItJYLDgep2ljZjnAAmAQsNrduyTVfe3ue5nZ88Dt7v5GWD4XmOTub1fY1kSCFgM9evTInT17dqX9bdu2jQ4dOlQb0/XXDwbg3nvfrfPrqqvaxJdpUY8x6vFB9GOMenzQdGKM2vHm+OOPX+TuQ2tc2d3TMgEdgEXA2eH8lgr1X4ePLwDHJpXPBXKr23Zubq6nMn/+/JTlyUaNCqZ02ZPt1Sa+TIt6jFGPzz36MUY9PvemE2O6jzd7IlV8wNtei+N3Wq4aMrOWwB+Bme7+p7D4SzPrGdb3BDaE5WuAPkmr9wbWpSMOERHZc+m4asiAR4Gl7n53UtVzwPjw+Xjg2aTyC8Orh0YAW70Bzg+IiEjtpOOGsmOAC4APzCzRMfZfwO3A02Y2AVgNnBvWvQicAiwHCoGL0xCDiIjUUTquGnoDsCqqT0yxvANX1Xe/IiKSHk32zmIREakdJYI9oGGtRaQparKDzqVbVcNag+5WFpFAtt64qhZBLWlYaxFpqpp8iyBdGVrDWotIU6UWQS1pWGsRaaqUCGpJw1qLSFOlRFBLGtZaRJqqJn+OIJ00rLWINEVqEYiIxJwSgYhIzCkRiIjEnBJBho0eHUwiIpmiRCAiEnO6amgP6WohEWlq1CIQEYk5JYIMiuuw1jovIhItSgQZUtWw1nFJBiISHUoEGRLXYa3j2goSiTIlggyJ47DWagWJRJMSQYbEcVjruLaCRKJOiSBDsmlY63Sd3I1jK0gkGygRZEgch7WOYytIJBsoEWRQXh6MGAGjRsHKlU07CUB2tYJE4kSJQKqVzqt84tgKEskGGmJCqlTVVT5Q94N3Q/y4T+L8hYb/EKkbtQgyLD8/ugewbLjKR/cliNSfWgRSpYa6yiddia8hWiwicaQWgVQp6lf5ZEOLRSQbKBFIlaJ+lY/uSxBJDyUCqVLUr/KJeotFJFsoETRB6RzmOcr3OkS9xSKSLZQIJGtFvcUiki2a/lVDzzwD558PrVqVn+bMgaOOguefh9tvLytv2TJ4vOeeoI9h/nz44x8rr3/NNdC5MyxeHEwV6086Kdj/6tWwcWP5bbdqBb16gRkUF0OzZtC8eWbfp2pE9fJWCA76OvCL1E/TTwQHHABXXAFFRcFBt6gomDp1CuqbN4c2bYKywsKy+uLioP6TT2D27LLyxLWKF18cJIL//V+YMqXyfrdsCR5/9Su4667K9SUlwb6vvRYeeihIBokk0blz2RnPm26Cl14Kylu3Dh67d4c//CGov+sueO+93evO/GwkH7xxJl+VdiYnB6ae+iZ5B7xVPkn16AGnnAJApw8/hNLS8tvv3Bn69Qu2v3Fj+dhatgzmRbLYzJnB1WWrVwff96ZOjfcXiqafCI48Mpiq8u1vB1NVLrssmBLcgwNn4hv8978Pl1xSPskUFUGHDkH9JZfAcceVrysqKlv/jDOC1kEiyRQVlT/Q9u4NBx9cef2ETz+FN9+EoiJmbj2ViQXnUEjQcb5qFUx8aAjs+jV5zCpbZ+jQ3Ymg/333BdtINmpUWTPgmGMq1595Jvz5z8HzI4+Er78un2hOPbWso/7f/70s0SSmE06ACy8M3ssf/ahya2rYsGC/xcV0mzev8vb79Qv+e4uLYcWKyuu3axckLJEUEvefJC491v0ncUgE6WYGLZLeto4dg6kqhxwSTFWpKRFdd10wVeW3v939dHIOFBaUry7c1YbJfWaQ986vUiaajydNYtiAAeWTTJcuZRv46U9h8+by9f37l9WPGQNffVW+fu+9y+o3bYLt28vX77dfUFdSEnTBJSc2gEmTgkRQUMBht95a+TXfeiv8+Mewfj0MHFi5/p574PrrYenSoPsv0dJJTHfeCeeeC0uWwKWXVkokM4+8g8nTB7B6tbN/+6+YeswL5B3+QVmr6T//Ew46KEhCr77KvitWwLp1Zds49tjgPdiwITjDXrFbsGfP4LG4OEiGLVsGf1fSKKq7/0SJoJGZ2cnAfUBz4HfufnumYmkqqryufk0z2GeflHXbDzwQ/u3fqt5oTf8Zqbq9kv31r1XXtWwZJCf3ICkkuuQS3+Y7deIf06czfPDg8okkJyeo33tvmDWrLMElpuOOC+r32itIohVbUz16BPXNmgVdhIluwS1bmPnFCUx86UAKdwIYq7btw8SXz4F588jbNSNo3YwcGSSCxYvh8suplIrefDNY5oUXghZhRe+/D4cfHnQJXntt2XuRSBTvvhu0eB56CO6/v3KL55lngu673/8+2EfyuolE16IFvPIKvPcevVatCpJiq1ZBN2jiM33/ffjyy/Lrtm0Lhx4a1H/1VeXWXITPZdWW7j+pLCOJwMyaA78BxgBrgH+a2XPu/lEm4mkq9t8/aOamKo80s+BgVrE7p0ULCvv2rbprr0MHGDeu6u3uuy/ccUfV9YceCi+/XK5ocg5hEihTSDsm7zedvJXTgwNj4tv7aafB2rW8tWABI3NzK7eYxowJDtQVE1Hv3kH9iBFBF1rF+kQLs3t3GDSo6m7FtWvhH/+oXP+LXwT1zzwDDz1EUvst6DZLJII776w8OFO3bkFLBoIk9uyz5esPOqisq/Ccc+CtCuefBg0KkjMEF1Qkuu4SrbJDD4Wbbw7qf/nL3d1++69ZEyTWgw4Kuksh6H4sKSm//Z49y1rYy5YF70XFRNa2bYoPu0zW/p80oEy1CIYDy919BYCZzQbOBJQI6mHq1PJ9n6Dr6vdUjd8Wk78Rt2kD++3Hzn33Ld9dltC7d9lBP5Vhw4KpKmefHUxVuemmYKpKeKHCG/Pmcezw4eUvgoDgIofEhRSJKfn81OWXB8ksuT5xkQUE3Xd77VX+/FaitQVQUBC0OJLX/+absvrp04OWSmkpByTKTj+9LBFccQV88UX513TeeUFLCCA3N+h2THbppcE1xO7B51OxtXT55Uyd+l9MvNQp/KasO65dsx1M3fvX8FQf+O53g4s9Jk0qt27OunVBQhs5MqifNatya23IkOA65oKCoOuxYn337tC+fZDgiosj08oyd2/8nZqdA5zs7t8L5y8Ajnb3q5OWmQhMBOjRo0fu7NmzK21n27ZtdEiclI2gTMT32mvdufPOgykubkaPHjv53vdW8K1vbahyeb2H5Y0bN4Ivv2xTqbxHjx3Mnr0w5Tp6D+uptJTCrVvpGN4QUtq+PQBt166lWVERVlxMs5ISrLiY4k6dKAyvaOv2+uuV6gtzcvg6Nxd27aLfo4/SrLgYKynZXf/VsGFsPOF5+vWWAAAI7klEQVQE5r3YhWn39WBtUQ96tfyCn3a9m++2f5Z1p5/O+jPOoNXGjeRefvnu9ZqVlNCsuJjlV13FmnPOod3KlQy/+OJKL+XjH/yAL049lY4ffUTuVVdVqv/wJz9h4wkn0GXRIgb/4AcAeLNm7GrRAm/Zkg+nTOHroUPZ6+236X///exq0YLFDz7IrsTNMtVI9Tkff/zxi9x9aE3rZioRnAuMrZAIhrv7NamWHzp0qL/99tuVyvPz8xmdrltoG0Cm4tuT8fn1HpZX8YoSCFpV1d2opvew/rIixvnzGT1qVNBqKimpfBFFUVFwBWDXrkGL4e9/L1+3c2dwRV6/fvDZZ/D00+VbU0VFwRWKhxwSdPklLqSYPbtWV8Gleg/NrFaJIFNdQ2uAPknzvYF1GYpFZLfEwV7XmEslZmVdZy1alO8Gq6hLFxg7tur6fv2CrqeqDB9edq6lEWQqEfwT6G9m/YC1wDjgPzMUS5MT5TuBs4HuVpa4yUgicPcSM7saeJng8tHH3P3DTMQiIhJ3GbuPwN1fBF7M1P5FRCSgQWNERGJOiUCkEaTzNyIaYnsSb0oEIg1s5kxYuBBefz0YHaPizbwimaZEINKAEvclJEYvT4x0qWQgUaJEINKAqhvpUiQqlAhEGpBGupRsoEQg0oCqGtGyPiNd6pyDpJsSgUgDmjo1GKsoWX1GhG3Icw66Eim+lAhEGlBeXjBgXWLwyL59qx/AriY65yANQT9VKdLA0jl2kc45SENQi0AkizTEOQcRJQKRLJLucw4NSeccsocSgUgWSfc5BxHQOQKRrJOXB488EjzXb09IOqhFICJZcW+CupoajhKBSMxpPCRRIhDJQvn56esW0r0JokQgEnO6N0GUCERiTvcmiBKBSMxlw70J2XAyO5spEYjEXEPcm5DOA7dOZjc8JQIRIS8PRoyAUaNg5cr6J4F0Hrh1MrvhKRGISFql+8Ctk9kNT4lARNIq3QduncxueEoEIpJW6T5wZ8PJ7GynRCAiaZXuA7cG2mt4GnRORNIqcYCeMCE4Ydy3b5AE6nPg1kB7DUuJQESA9B5gdeDOLuoaEhGJOSUCEZGYUyIQEYk5nSMQkaygcw0NRy0CEZGYUyIQEYk5dQ2JSINQV072UItARCTmlAhERGKuXonAzH5hZh+b2ftm9oyZdUmq+5GZLTezZWY2Nqn85LBsuZndXJ/9i4hI/dW3RfAqMMjdjwA+AX4EYGaHAuOAw4CTgQfMrLmZNQd+A3wbOBQ4L1xWREQypF6JwN1fcfeScHYh0Dt8fiYw2913uvtnwHJgeDgtd/cV7l4EzA6XFRGRDDF3T8+GzP4XeMrdZ5jZr4GF7j4jrHsUeClc9GR3/15YfgFwtLtfnWJ7E4GJAD169MidPXt2pX1u27aNDh06pCX+hhD1+CD6MUY9Poh+jFGPDxRjOqSK7/jjj1/k7kNrWrfGy0fN7DVg3xRVk9392XCZyUAJkPhVUkuxvJO6BZIyE7n7w8DDAEOHDvXRo0dXWiY/P59U5VER9fgg+jFGPT6IfoxRjw8UYzrUJ74aE4G7f6u6ejMbD5wGnOhlzYs1QJ+kxXoD68LnVZWLiEgG1PeqoZOBScAZ7p78c9XPAePMrLWZ9QP6A/8A/gn0N7N+ZtaK4ITyc/WJQURE6qde5wjMbDnQGtgcFi1098vDusnAJQRdRte7+0th+SnAvUBz4DF3r/EH7MxsI7AqRVVXYFOdX0DDi3p8EP0Yox4fRD/GqMcHijEdUsXX19271bRi2k4WZ4KZvV2bEyGZEvX4IPoxRj0+iH6MUY8PFGM61Cc+3VksIhJzSgQiIjGX7Yng4UwHUIOoxwfRjzHq8UH0Y4x6fKAY06HO8WX1OQIREam/bG8RiIhIPSkRiIjEXJNIBGb2AzNzM+ua6Vgqqm6o7kyK+nDgZtbHzOab2VIz+9DMrst0TKmEo+q+Y2bPZzqWVMysi5nNCf8Gl5rZyEzHVJGZ3RB+xkvMbJaZtclwPI+Z2QYzW5JUtreZvWpmn4aPe0Uwxjofa7I+EZhZH2AMsDrTsVQh5VDdmZQlw4GXADe6+yHACOCqCMYIcB2wNNNBVOM+4C/uPhA4kojFama9gGuBoe4+iOBG03GZjYrpBMPnJ7sZmOvu/YG54XwmTadyjHU+1mR9IgDuAX5IFYPXZVo1Q3VnUuSHA3f39e6+OHxeQHAA65XZqMozs97AqcDvMh1LKmbWCTgOeBTA3YvcfUtmo0qpBdDWzFoA7cjw+GPuvgD4qkLxmcDj4fPHgbMaNagKUsVYn2NNVicCMzsDWOvu72U6llq6hLLhuDOpF/B50vwaInaQTWZmOcBRwN8zG0kl9xJ8CdmV6UCqcACwEZgWdl/9zszaZzqoZO6+FriLoEW/Htjq7q9kNqqUerj7egi+pADdMxxPTfboWBP5RGBmr4V9hxWnM4HJwH9HPMbEMhWH6s6kqoYJjxwz6wD8kWC8qn9lOp4EMzsN2ODuizIdSzVaAEOAB939KGA7me/SKCfsaz8T6AfsB7Q3s/MzG1V2q8uxpsZhqDOtqmGwzexwgj+e98wMgmbQYjMb7u5fNGKIdR2qO5OqGyY8MsysJUESmOnuf8p0PBUcA5wRDqLYBuhkZjPcPUoHsTXAGndPtKTmELFEAHwL+MzdNwKY2Z+A/wfMyGhUlX1pZj3dfb2Z9QQ2ZDqgVOp6rIl8i6Aq7v6Bu3d39xx3zyH4ox/S2EmgJtUM1Z1JkR8O3ILs/iiw1N3vznQ8Fbn7j9y9d/i3Nw6YF7EkQPi/8LmZHRwWnQh8lMGQUlkNjDCzduFnfiIRO6Edeg4YHz4fDzybwVhSqs+xJmsTQRb5NdAReNXM3jWzhzIdUHhC6WrgZYJ/uqfd/cPMRlXJMcAFwAnh+/Zu+O1b9sw1wEwzex8YDPxPhuMpJ2ytzAEWAx8QHJMyOpSDmc0C3gIONrM1ZjYBuB0YY2afElyleHsEY6zzsUZDTIiIxJxaBCIiMadEICISc0oEIiIxp0QgIhJzSgQiIjGnRCAiEnNKBCIiMff/AWE/G+BCVJKMAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fcd8f06e860>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# test data\n",
    "x = [ -3.5, -2.5, -1.5, -0.5, 0.5, 1.5,  2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5]\n",
    "y = [ 219., 155., 89., 11., 54., 80., 48., 5., -47., -64., -101., -125., -124., 13., 239., 525. ]\n",
    "ey = [ 65., 41, 30., 26., 26., 24., 18., 4., 17., 33., 45., 49., 38., 6., 74., 160.]\n",
    " \n",
    "ey = np.array(ey)\n",
    "w=1./ey\n",
    "\n",
    "n=1 #polynom degree\n",
    "\n",
    "pars=sp.polyfit(x,y,n,w=w) # do the fit\n",
    "\n",
    "yr=sp.polyval(pars,x) # compute residuals based on fitted line\n",
    "chi2=sp.sum((yr-y)**2/ey**2) #compute Chi2\n",
    "\n",
    "print (pars)\n",
    "print ('Chi2, n, chi2/ndf:', chi2, n, chi2/(len(x)-n))\n",
    "\n",
    "plt.title('Polynom Regression Example')\n",
    "plt.errorbar(x,y,ey,fmt='ob')\n",
    "#plt.plot(x,yr,'r--')\n",
    "\n",
    "xr = np.linspace(x[0],x[-1],200)\n",
    "yr2 = sp.polyval(pars,xr) \n",
    "plt.plot(xr,yr2,'r--')\n",
    "\n",
    "plt.legend(['regression','values']);\n",
    "plt.grid(True);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [Root]",
   "language": "python",
   "name": "Python [Root]"
  },
  "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.5.2"
  },
  "toc": {
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": true,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}