{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Ejemplo de calibración " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Ejemplo.** Un sistema de medida de altura usando pulsos de luz. La tabla muestra los valores reales \n", "y los medidos (con error) cuando se incrementa la distancia y cuando se disminuye. Desarrolle la calibración y establezca el margen de error con un nivel de confianza del 99.7% para un valor medido de 4.32mm.\n", "\n", "| X ref(mm) | Y medido (Inc.) | Y medido (Dism.) |\n", "|-------------:|------------------:|--------------------:|\n", "| 0 | -1.12 | -0.69 |\n", "| 1 | 0.21 | 0.42 |\n", "| 2 | 1.18 | 1.65 |\n", "| 3 | 2.09 | 2.48 |\n", "| 4 | 3.33 | 3.62 |\n", "| 5 | 4.50 | 4.71 |\n", "| 6 | 5.26 | 5.87 |\n", "| 7 | 6.59 | 6.86 |\n", "| 8 | 7.73 | 7.92 |\n", "| 9 | 8.68 | 9.10 |\n", "| 10 | 9.88 | 10.20 |" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Solucion\n", "\n", "El ejercicio no lo indica, pero podríamos suponer que cada valor medido y (respuesta instrumental), viene de realizar varias mediciones para cada valor de entrada $X_{ref}$ y se calcula un promedio $\\bar{y}$ y una varianza $S$. Como tal cada valor se reporta con [un margen de error y nivel de confianza](http://www.stat.wmich.edu/s216/book/node79.html \"t-based Confidence Interval for the Mean\") $(1-\\alpha)100\\%$: \n", "\n", "$$\\bar{y}\\pm t \\left( \\tfrac{S}{\\sqrt{n}} \\right)$$\n", "\n", "donde $t$ es el valor crítico de la distribución $t_{n-1}$. Grafiquemos los datos y observemos su comportamiento." ] }, { "cell_type": "code", "execution_count": 83, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAESCAYAAAASQMmzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAF4JJREFUeJzt3X2UZHV95/F30T1sHOzVcS2fQWCPfIMEEXBkbDWADudg\nD6KebExYXQmQeHajhrCrJz5slPVkz+q6KCZudqM8CFFRw2YdkxkcBwlqzjgwDESMmK8b48hTgJ7Z\nwQUnCMPU/nGrLjXDdHf1Q9Wtuv1+ncOh6vadut/f6Zn+9O/e30Oj1WohSRLAIVUXIEkaHoaCJKlk\nKEiSSoaCJKlkKEiSSoaCJKk03u8LRMQVwDrggcw8vn3sY8BZwKPAj4DzMvOn/a5FkjS7QfQUrgTO\nPODY14HjMvME4IfA+wZQhyRpDn0Phcz8NrD7gGObM3Nf++1NwAv6XYckaW7D8EzhfGBj1UVIkioO\nhYj4APBoZn6hyjokSYW+P2ieSUT8BjAFvLaX8/fufbw1Pj7W15okqYYa8zm5klCIiDOB9wCnZuYj\nvfyZ3bv39LeotmZzgunphwZyrX6rU1ugXu2pU1ugXu2pU1ugaM98DGJI6jXAqcAzI+Iu4EMUo40O\nBTZHBMB3MvO3+12LJGl2fQ+FzDznIIev6Pd1JUnzNwyjjyRJQ8JQkCSVDAVJUslQkCSVDAVJUslQ\nkCSVDAVJUslQkCSVDAVJUslQkCSVDAVJUslQkCSVDAVJUslQkCSVDAVJUslQkCSVDAVJUslQkCSV\nDAVJUslQkCSVDAVJUslQkKS2VqtFq9WquoxKGQqSRBEIU1NrmZycXNbBMN7vC0TEFcA64IHMPL59\n7BnAl4AXAjuAN2fmg/2uRZIOphMI27dvA2Bqai0bN15Po9GouLLBG0RP4UrgzAOOvRfYnJnHAN9o\nv5ckVazvoZCZ3wZ2H3D4bOCq9uurgDf2uw5Jmkmj0WDjxus5+eTVrFmzZtn2EmAAt49m8OzMvL/9\n+n7g2RXVIUlAEQwbNmym2Zxg166fVV1OZaoKhVJmtiJizqc6q1atZHx8bBAl0WxODOQ6g1CntkC9\n2lOntsDot6fVajE5OQnAli1b7CkM2P0R8ZzMvC8ings8MNcf2L17zwDKKv5iT08/NJBr9Vud2gL1\nak+d2gKj354DHzSvXn1KbW4hzTesqxqS+lXg3Pbrc4GvVFSHJKnLIIakXgOcCjwzIu4CPgh8BPhy\nRFxAe0hqv+uQpJl0HjRPTa1lxYox1q/fVItewkL0PRQy85wZvrS239eWpF51gqHZnGDnzoerLqcy\nzmiWpLZGo7FsewgdhoIkqWQoSJJKhoIkqWQoSJJKhoIkqWQoSJJKhoKkoeeOaINjKEgaap11iaam\n1hoMA2AoSBpa3QvVbd++zWAYAENBklSqfD8FSZpJ90J1QG2Wsx5mhoKkodYJhs5r9ZehIGnoGQaD\n4zMFSUPPIamDYyhIGmoOSR0sQ0HS0HJI6uAZCpKkkg+aJQ0th6QOnqEgaag5JHWwDAVJQ88wGByf\nKUiSSoaCJKlU6e2jiLgIuABoAd8DzsvMn1dZkyQtZ5X1FCLi+cC7gJMz83hgDPj1quqRJFX/oHkc\nWBkRjwMrgXsqrkeSlrXKegqZeQ9wCXAncC/wYGZeX1U9khbGdYnqpVHVNzMiVgHXAm8Gfgr8GXBt\nZn7+YOfv3ft4a3x8bIAVSppLq9VicnISgC1btjh0dDjN65tS5e2jtcCPM3MXQET8OTAJHDQUdu/e\nM5Cims0JpqcfGsi1+q1ObYF6tacObelelwhg9epTajHjuA7fm27N5sS8zq8yFH4CrImIpwCPUITE\nzRXWI0nLXpXPFG6muH10K3B7+/Cnq6pH0vw0Gg02bNjMypUrOeyww9iwYfPI9xJU8eijzLwYuLjK\nGiQtTKvVYt26M9izp7i1u27dGbW4fbTcOaNZklSqep6CpBHVvaz1ihVjrF+/yV5CDRgKkhasEwzN\n5gQ7dz5cdTlaAt4+krQojUbDHkKNGAqSpJKhIEkqGQqSpNKMD5oj4lyKfQ5mu1nYysyrl7wqSVIl\nZht99HHgq7N8vQGcBRgKklQTs4XC1zLzvNn+cEQcdPE6SdJomjEUMvMtndcRcTyw6oCvf6v7HEnS\n6Jtz8lpEfBE4iSfvinZ6XyqStCQ6e6U4h0Dz0cuM5hOAYzPz8X4XI2lpdPY6AFykTvPSy5DUm4AX\n9bsQSUuje/Ob7du3MTW11u0y1bNeego3AH8bEf8I7G0fa2Xm0f0rS5JUhV5C4Q+A1wB39rkWSUug\ne/VS8PaR5qeXUHgA+OvM3NfvYiQtHW8ZaSF6CYXbge9ExGbgsfaxVmZ+uH9lSVqozjOFW2+9BYCp\nqbX2FtSzXkLhTp64deTfKkmqsTlDob2PsqQR4TMFLUYvk9d+F/gg8PSuw63MHOtbVZIWpRMMnddS\nr3q5fXQR8NLMdPSRNEIMAy1EL5PX7qAYgSRJqrleegqfBL4XEVvZf/La+Yu9eEQ8HbgMOI5i74bz\nM3PrYj9XkrQwvYTCHwF/yv6T15ZqAPQngY2Z+a8iYhw4bIk+V5K0AL2Ewj/1Y05CRDwNeHVmnguQ\nmXuBny71dSRJveslFK6PiEuA64BHOwcz81uLvPZRwHREXEmxEut24MLM3LPIz5WGWqvVcraxhlYv\noXASxe2ikw44vtj9FMbbn/nOzNwWEZcC76UY/vokq1atZHx8MKNgm82JgVxnEOrUFhj99rRaLSYn\nJwHYsmVLrUYIjfr3plud2jJfM4ZCRKzJzK2Zedpc5yzw2ncDd2fmtvb7aylC4aB27x5MB6LZnGB6\n+qGBXKvf6tQWGP32dC9pDbB69Sm1mVg26t+bbnVqC8w/4GbrKVweEa+b5esNipFDvzSvK7Zl5n0R\ncVdEHJOZPwTWAt9fyGdJkpbGbKFwGPDNPl//XcDnI+JQ4EfAeX2+nlSZ7uUnVqwYY/36TbXoJahe\nZgyFzDyy3xfPzO8Cq/t9HWlYdIKh2Zxg586Hqy5HepJeZjRLWkKNRsMegoaWoSBJKhkKkqRSL0tn\nPwv4FPDa9vk3AP82M+/vc22SpAHrpafwJ8DNFDOQXwh8B7i8n0VJdeaMZg2zXmY0H52Zb+p6/18j\n4m39Kkiqs84ENoekalj10lPYFxFHdN5ExAvpWgNJUm+6ZzRv3bqVqam19hg0dHrpKfw+sCUibm6/\nXwO8vX8lSZKqMmcoZOZfRsRJFJPMDqF4yOxObNI8OaNZo2C2BfE+1PW2RbHWEcCJEUE/9liQ6s4Z\nzRp2sz1T6ITAy4FfAR6neJZwFsX2mZIWwBnNGmazrX10MUBEbAFe0dn8JiI+Adw4iOIkSYPVy+ij\nZx7w/lBgVR9qkSRVrJfRR58BbomIDcAYxe2jS/talVSBzvBQb+1oOZuzp5CZHwPeBtxHsVvar2bm\nH/e7MGmQOnMInDug5W7GUIiI17f/fy7wYmAnsAt4qTOaVSfdk8q2b99mMGhZm+320cuAvwBOpxiS\neqCr+1KRJKkyjVH5jWh6+qGBFFqnTbvr1Bbob3v27dvH0Uc/D4B/+Id7OeSQ/q4q7/dmeNWpLQDN\n5sS8HpLNNnntx11vuyevAbQy8+h51iYNpVarxbp1Z7Bnzx4A1q07g40br/eBs5al2W4fnd7+/4eA\nHwGfpZjA9q8BA0GSami2yWs7ACLiJZl5XteXLomIW/tdmDQo3WsSAfYStKz1Mk+hERGvycwbACJi\nCnisv2VJg9UJhs5rabnqJRQuAK6OiOdSPFf4CfDWvlYlVcAwkHpbOvs24PiI+Bft97uWsoCIGANu\nAe7OzNcv5WdLkuZnznF3EXFkRGwGbgIOjYi/ioijlrCGC4E7OPhcCEnSAPUyGPtPgP8GPATcD3wB\nuGopLh4RLwCmgMvYf8irJKkCPa2SmpmbADJzX2Z+BnjaEl3/E8B7gH1L9HmSpEXo5UHznvZv9ABE\nxKuARxZ74Yg4C3ggM2+LiNPmOn/VqpWMj48t9rI9aTYnBnKdQahTW6Be7alTW6Be7alTW+arl1D4\n98AG4OiI+C7wDOBXl+Dak8DZ7SGuvwD884i4OjMPutje7t17luCSc6vTFPc6tQXq1Z46tQXq1Z46\ntQXmH3C9jD7aFhEvA46h2E/h7zLz0YWVt9/nvh94P0BEnAq8e6ZAkCQNxpyhEBG/CLydrt3WIqKV\nmecvcS2OPtKMWq2Wy1lLA9DL7aP/DVwDfJcnRggt6b/OzPwm8M2l/EzVR2e/gxUrxli/fpOTzKQ+\n6iUUdmfmh/teiXQQ3RvgAExNrXVtIqmPegmFz0bEfwa+AeztHMzMb/WtKklSJXoJhdOA1RSjhbqd\n/uRTpaXVvYKpt4+k/uslFF4GHJOZPuVTJRqNBhs2bKbZnGDXrp9VXY5Ua73MaP4e8JJ+FyLNpLMz\n2itf+UpHIEl91ktP4V8Ct0bEfUBnfoLbcWogfNAsDVYvofDGvlchSRoKvcxo3jGAOqSD8kGzNFi9\n9BSkSnWCodmcYOfOh6suR6q1GR80R8RhgyxEmk2j0bCHIA3AbKOPbo+IXx5YJZKkys0WCv8OuCIi\nLomIfzaogiRJ1ZkxFDLz68AJ7bc3R8QvR8QRnf8GU54kaZBmfdCcmT+LiN8HDgfWAw92ffmofhYm\nSRq8WUOhvWXmfwc2AUdkZn22I5IkPcmMoRARfwacDJyfmd8YXEmSpKrM1lO4Hzg+M12BTDPqrEXk\ncFGpHmYMhcx85yAL0ejprEsEuB6RVBO9rJIqPUn3QnXbt29jamqtK5hKNWAoSJJKrn2kBeleqA68\nfSTVhaGgBesEQ+e1pNFnKGhRDAOpXioLhYg4HLgaeBbQAj6dmX9YVT1aGIekSvVS5YPmx4CLMvM4\nYA3wjog4tsJ6NE+dEUiOPJLqo7JQyMz7MvNv2q8fBn4APK+qejQ/DkmV6mkohqRGxJHAicBNFZci\nSctao+rf7iLiqcCNwB9k5ldmOm/v3sdb4+NjA6tLc2u1WkxOTgKwZcsWnytIw2le/zArDYWIWAH8\nJXBdZl4627nT0w8NpNBmc4Lp6XosBjuItgzyQbPfm+FVp/bUqS0AzebEvP5xVjn6qAFcDtwxVyBo\neNk7kOqlynkKrwTeSrEX9G3tY+/LzK9VWFNttFotH/xKmrfKQiEz/5ohedBdN52RQStWjLF+/SZ/\nm5fUM38o10z3UNGtW7c6VFTSvBgKkqSSax/VTPfqpd4+kjRfhkINdYKh2Zxg586Hqy5H0gjx9lFN\nNRoNewiS5s1QqCmHpEpaCEOhhjojkCYnJw0GSfNiKNSMQ1IlLYahIEkqOfqoZhySKmkxDIUackiq\npIXy9lFNOSRV0kIYCpKkkqEgSSoZCpKkkqEgSSoZCpKkkqEwYK5JJGmYGQoD1FmCwqUnJA0rQ2FA\nutck2r59m8EgaSgZCpKkkstcDEj3mkQAGzde74xjSUPHUBgwbxlJGmaVhkJEnAlcCowBl2XmR6us\np586zxRuvfUWAKam1tpbkDR0KnumEBFjwKeAM4EXA+dExLFV1SNJqran8HLg7zNzB0BEfBF4A/CD\nCmvqG58pSBoFVYbC84G7ut7fDZxSUS0D0QmGzmtJGjZVhsK8nriuWrWS8fGxftWyn2ZzYiDXGYQ6\ntQXq1Z46tQXq1Z46tWW+qgyFe4DDu94fTtFbOKjdu/f0vSAo/jJMTz80kGv1W53aAvVqT53aAvVq\nT53aAvMPuCpD4RbgRRFxJHAv8GvAORXWI0nLXmWjjzJzL/BOYBNwB/ClzKzlQ2ZJGhWVzlPIzOuA\n66qsQZL0BNc+kiSVDAVJUslQkCSVDIUu7oomabkzFNo6C9ZNTk4aDJKWLUOB/XdF27p1q7uiSVq2\nDAVJUslNdth/BdMVK8ZYv36TC9ZJWpYMhbZOMDSbE+zc+XDV5UhSJbx91KXRaNhDkLSsGQpdHJIq\nabkzFNockipJhgLgkFRJ6jAUJEklRx/hkFRJ6jAU2hySKknePtqPQ1IlLXeGgiSpZChIkkqGgiSp\nZChIkkqGgiSpVMmQ1Ij4GHAW8CjwI+C8zPxpFbVIkp5QVU/h68BxmXkC8EPgfRXVIUnqUklPITM3\nd729CfiVKuqQJO1vGJ4pnA9srLoISVIfewoRsRl4zkG+9P7M/Iv2OR8AHs3ML/SrDklS7xpVLREd\nEb8B/Bbw2sx8pJIiJEn7qWr00ZnAe4BTDQRJGh6V9BQi4v8AhwL/t33oO5n52wMvRJK0n8puH0mS\nhs8wjD6SJA0JQ0GSVDIUJEklt+Nsa4+IuhQYAy7LzI9WXNKCRcThwNXAs4AW8OnM/MNqq1qciBgD\nbgHuzszXV13PYkTE04HLgOMovj/nZ+bWaqtamIi4CLiAoh3fo1jH7OfVVtW7iLgCWAc8kJnHt489\nA/gS8EJgB/DmzHywsiJ7NENb5r3OnD0Fyh84nwLOBF4MnBMRx1Zb1aI8BlyUmccBa4B3jHh7AC4E\n7qD44TPqPglszMxjgZcAP6i4ngWJiOcD7wJObv8QGgN+vdqq5u1Kin/33d4LbM7MY4BvtN+PgoO1\nZd7rzBkKhZcDf5+ZOzLzMeCLwBsqrmnBMvO+zPyb9uuHKX7oPK/aqhYuIl4ATFH8dj3Sm2hHxNOA\nV2fmFQCZuXfEVwgeB1ZGxDiwErin4nrmJTO/Dew+4PDZwFXt11cBbxxoUQt0sLZk5ubM3Nd+exPw\ngrk+x1AoPB+4q+v93e1jIy8ijgROpPgLMao+QTHZcd9cJ46Ao4DpiLgyIm6NiM9ExMqqi1qIzLwH\nuAS4E7gXeDAzr6+2qiXx7My8v/36fuDZVRazhHpaZ85QKNThlsSTRMRTgWuBC9s9hpETEWdR3CO9\njRHvJbSNAycBf5yZJwE/Y3RuT+wnIlZR/FZ9JEVP9KkR8ZZKi1pimdmiBj8f5rPOnKFQuAc4vOv9\n4RS9hZEVESuA/wV8LjO/UnU9izAJnB0RPwauAV4TEVdXXNNi3E3xsHxb+/21FCExitYCP87MXZm5\nF/hziu/XqLs/Ip4DEBHPBR6ouJ5Faa8zNwX0FNiGQuEW4EURcWREHAr8GvDVimtasIhoAJcDd2Tm\npVXXsxiZ+f7MPDwzj6J4iHlDZr6t6roWKjPvA+6KiGPah9YC36+wpMX4CbAmIp7S/ju3lmIwwKj7\nKnBu+/W5wMj+UtW1ztwbel1nzmUu2iLidTwxJPXyzPwvFZe0YBHxKuBbwO080fV9X2Z+rbqqFi8i\nTgX+Q2aeXXUtixERJ1A8ND+UEd+ONiIupvglai9wK/Cb7cEaIyEirgFOBZ5J8fzgg8B64MvAEYzW\nkNQD2/IhitFG81pnzlCQJJW8fSRJKhkKkqSSoSBJKhkKkqSSoSBJKhkKkqSSoaBai4jTIuLeiGh2\nHXt3RFy7iM/8rYjYEREfPeD4joj4fkSsWUzNB3zm5yNiV0ScO/fZ0uIZCqq1zLwR+BzwGYD2D+y3\nUywOtlDnUEzS+r0DjreA1y3l3giZ+RaKGbZOKNJAuMmOloMPADdHxO8A7wT+TWb+vwNPiogdwFbg\npcCrgddR7ONwCLAdeAfwe8Bq4H9ExO9k5nUHu2BE3Egxw3ct8BSKfQcupNiv4xOZeWl7NvARFHsq\nPAv4j8BrgFOA72Zm994EdVgMUCPAnoJqr73swluAjwPXZOZMy4i3KDa/+UWKH9K/CbwiM08EpoF3\nZ+aHKdbKumCmQOj6rFZmvgT4U+CPgDdRhM0Hu847jmI/j7cCVwAfAX4JOCkijl9Ie6XFMBS0XLyK\n4gf7Ge2d9mbSCYzTgRcBN0XEbRRLREfXeb385t4JjTuBrZn5SGbeCTy965zOJih3Av+YmX+XmY9T\nrNy7qodrSEvKUFDtRcSLgYuBVwA/p7hNM5N/av//EODLmXliu6fwcopbQB293ON/tOv13hnOeayH\nc6SBMRRUaxHxCxSbsL87M3dQLIX8rog4ZY4/eiPwpohotpeF/p8UzwSkWjMUVHcfp3ho+wWA9u2b\n3wU+N9s2mJl5O/CfgBuAv20f/sgCazhw967WHMelyrh0trRE2rvDnZaZP1niz/0s8FeZedVc50qL\nZU9BWloblnryGnAW9iI0IPYUJEklewqSpJKhIEkqGQqSpJKhIEkqGQqSpJKhIEkq/X9ASPTpWak3\nFQAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline\n", "from scipy import stats\n", "import seaborn as sns\n", "\n", "# X referencia: se asume su incertidumbre muy baja\n", "X = np.arange(0,11)\n", "x = np.hstack([X,X])\n", "\n", "# y\n", "y_inc = np.array([-1.12, 0.21, 1.18, 2.09, 3.33, 4.5, 5.26, 6.59, 7.73, 8.68, 9.88])\n", "y_dis = np.array([-0.69, 0.42, 1.65, 2.48, 3.62, 4.71, 5.87, 6.86, 7.92, 9.10, 10.20])\n", "# y = np.mean([y_inc, y_dis], axis=0)\n", "y = np.hstack([y_inc,y_dis])\n", "\n", "plt.plot(X,y_inc,'.k')\n", "plt.plot(X,y_dis,'.k')\n", "plt.ylabel('Y medido [mm]')\n", "plt.xlabel('X ref [mm]')\n", "plt.axis([-1,12,-2,12])\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Se observa que no hay mucha dispersión de los datos y la tendencia claramente es lineal. \n", "Ahora procedemos a realizar la calibración mediante regresión lineal no ponderada. Para $Y$ tomaremos \n", "las dos columnas. No analizaremos el error por histéresis en este momento." ] }, { "cell_type": "code", "execution_count": 78, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The coefficients are b1 = 1.0839 and b0 = -0.8525\n" ] } ], "source": [ "# Numero de datos\n", "n = len(x)\n", "\n", "# polyfit computes the coefficients b1 and b0 of degree=1 via Least Squares\n", "p, V = np.polyfit(x,y,1,cov=True)\n", "\n", "b1 = p[0] # slope\n", "b0 = p[1] # intercept\n", "\n", "print 'The coefficients are b1 = %2.4f and b0 = %2.4f' % (b1, b0)" ] }, { "cell_type": "code", "execution_count": 79, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "the coefficients are b1 = 1.0839, and b0 = -0.8525\n", "The residual standard deviation: 0.2043\n", "The standard deviation of b1: 0.0138\n", "The standard deviation of b0: 0.0689\n" ] } ], "source": [ "# Manual estimation\n", "\n", "b1 = np.sum((x-x.mean())*y)/np.sum((x-x.mean())**2)\n", "\n", "b0 = y.mean() - b1*x.mean()\n", "\n", "print 'the coefficients are b1 = %2.4f, and b0 = %2.4f' % (b1,b0)\n", "\n", "\n", "S2yx = np.sum((y-(b1*x+b0))**2)/(n-2)\n", "print 'The residual standard deviation: %6.4f' % np.sqrt(S2yx)\n", "\n", "S2b1 = S2yx/np.sum((x-x.mean())**2)\n", "print 'The standard deviation of b1: %6.4f' % np.sqrt(S2b1)\n", "\n", "S2b0 = S2yx*(1/n+(x.mean()**2)/np.sum((x-x.mean())**2))\n", "\n", "print 'The standard deviation of b0: %6.4f' % np.sqrt(S2b0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Recall that the residual standard deviation $S_{y/x}$ represents the dispersion of the residue against the model. It should be normal." ] }, { "cell_type": "code", "execution_count": 80, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Con un 99.7% de confianza los valores de b1 y b0:\n", "b1 = 1.0839 ± 0.0465\n", "b0 = -0.8525 ± 0.2325\n" ] } ], "source": [ "# Para reportar los parametros b1 y b0 con su respectivo margen de error 99.7% confianza\n", "# n-2 grados de libertad porque se estimaron dos parametros \n", "\n", "t_alpha = stats.t.ppf(1-0.003/2,n-2)\n", "print 'Con un 99.7% de confianza los valores de b1 y b0:'\n", "print 'b1 = %2.4f ± %2.4f' % (b1, t_alpha*np.sqrt(S2b1))\n", "print 'b0 = %2.4f ± %2.4f' % (b0, t_alpha*np.sqrt(S2b0))" ] }, { "cell_type": "code", "execution_count": 82, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The estimated real value is: 4.7723\n", "Estimated x0 variance is: 8.3729e-06\n", "(4.7625092630510446, 4.7820493446235215)\n", "The x0 value with uncertainty 99.7 is: 4.7723 ± 0.0098\n", "The x0 range is: (4.7625, 4.7820)\n", "The bias error is: 0.4523, and the random error is: ±0.0098\n" ] } ], "source": [ "# Calculate real value with uncertainty based on measurement\n", "\n", "y0m = 4.32 # valor medido en mm, respuesta instrumental\n", "\n", "# usando regresion lineal para estimar valor verdadero\n", "x0 = (y0m-b0)/b1\n", "print 'The estimated real value is: %6.4f' % x0\n", "\n", "# estimar la varianza de x0 estimado\n", "\n", "S2x0 = S2yx/b1**2 * (1/2+1/n+(x0-x.mean())**2/(np.sum((x-x.mean())**2)))\n", "print 'Estimated x0 variance is: %2.4e' % S2x0\n", "\n", "\n", "# Assuming x0 approximately normal the limits of the (1-alpha)100% confidence interval\n", "# for the true value of x0 corresponding to the response average y0m are\n", "# x0m± = x0 ± t(1-alpha/2,n-2)Sx0\n", "# alpha is 0.05 for 95% confidence\n", "# n-2 degrees of freedom\n", "\n", "t_alpha = stats.t.ppf(1-0.003/2,n-2)\n", "\n", "print stats.t.interval(1-0.003,n-2,loc=x0,scale=np.sqrt(S2x0))\n", "\n", "\n", "print 'The x0 value with uncertainty 99.7 is: %4.4f ± %4.4f' % (x0, t_alpha*np.sqrt(S2x0))\n", "print 'The x0 range is: (%4.4f, %4.4f)' % (x0-t_alpha*np.sqrt(S2x0), x0+t_alpha*np.sqrt(S2x0))\n", "\n", "print 'The bias error is: %2.4f, and the random error is: ±%2.4f' % (x0-y0m,t_alpha*np.sqrt(S2x0))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Nótese que el error *bias* (sesgo) es relativamente alto (0.4523 mm) y esto se debe a que el modelo de la respuesta instrumental no pasa por (0,0). El intersecto es en -0.8525. Esto se puede corregir en la calibración y es normal que ocurra. Por otra parte, el error aleatorio con 99.7% de confianza es relativamente pequeño y tiene sentido, pues hay poca dispersión en los datos." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.10" } }, "nbformat": 4, "nbformat_minor": 0 }