oseiskar/rc_discharge.ipynb

## rc_discharge.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from numpy.random import randn\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulate RC discharge measurement"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--- Actual ---\n",
      "\n",
      "        Resistance: 6.27323\n",
      "        Capacitance: 0.0105955\n",
      "        Initial voltage: 6.7069\n",
      "        \n",
      "--- Measurement ---\n",
      "\n",
      "        Capacitance, C: 0.0116087\n",
      "        Initial voltage, V0: 6.86521\n",
      "        Times, t: [0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1  0.11 0.12 0.13 0.14\n",
      " 0.15]\n",
      "        Voltages: [5.64289449 4.19601466 3.89185733 3.67108154 3.52121243 2.29449239\n",
      " 2.43223162 2.44059144 1.66906039 1.32721103 1.31180082 1.07327804\n",
      " 1.07107477 0.86220034 0.64562205]\n",
      "        \n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAHpZJREFUeJzt3Xd0lPed7/H3d0YNgYQACaEKovcqmik2duyAe2wwttNcYu9u+iZO7Ozu2U3uSe7d2Cnr7N6bXfc4iY0NOG64Oy5gMCDRbIppAiGKJLooqvO7f2hMZMAwQho9Uz6vc3SmPTP6SEgffvo9v2cec84hIiLRw+d1ABERaR0Vt4hIlFFxi4hEGRW3iEiUUXGLiEQZFbeISJRRcYuIRBkVt4hIlFFxi4hEmYRwvGhmZqbr06dPOF5aRCQmlZaW7nfOZYWybViKu0+fPpSUlITjpUVEYpKZ7Qx1W02ViIhEGRW3iEiUUXGLiEQZFbeISJRRcYuIRBkVt4hIlFFxi4hEmYgp7tqGJh56fxtLt+33OoqISESLmOJO8BkPLy7jsSU7vI4iIhLRIqe4/T5uGJvHO59UUV1T53UcEZGIFTHFDTBnXAFNAcfzq3d7HUVEJGJFVHH379mFMYUZPFuyC+ec13FERCJSRBU3NI+6t1QdY23FEa+jiIhEpIgr7qtH5ZCS6GN+yS6vo4iIRKSIK+70lERmDc/hxbV7qG1o8jqOiEjEibjiBpgzLp+a2kZeX7/P6ygiIhEnIot7Ut8e5HfrxPySCq+jiIhEnIgsbp/PmD0unw+27Wf34ZNexxERiSgRWdwAN47NxzlYWKpRt4hISxFb3AXdU7moXw/ml+4iENCabhGRT0VscQPMKc5n18GTLC876HUUEZGIEdHFPXNYDmnJCcwv1ZpuEZFPRXRxd0ryc/WoXF79aB81tQ1exxERiQgRXdzQPF1ysqGJRev2eh1FRCQiRHxxjynIoF9WZ+ZrdYmICBAFxW1m3FRcQOnOQ2yrPuZ1HBERz0V8cQN8aWwefp+xQKNuEZHoKO6eaSlcMjCLhaUVNDYFvI4jIuKpqChuaN5JWVVTx+ItOpmwiMS3qCnuSwdn071zktZ0i0jcC6m4zWyHmX1kZmvMrCTcoc4mKcHH9aPzeHNDJQeP13sRQUQkIrRmxD3DOTfaOVcctjTnMac4n4YmxwtrdDJhEYlfUTNVAjAkJ53heel6n24RiWuhFrcD3jCzUjO7O5yBzuem4gI27D3K+j06mbCIxKdQi3uqc24sMAv4lplNP30DM7vbzErMrKS6urpdQ7Z07ahckvw+jbpFJG6FVNzOud3ByyrgL8CEs2zzkHOu2DlXnJWV1b4pW8hITeLyYdk8v2Y3dY06mbCIxJ/zFreZdTaztE+vA1cAH4c72LnMGZfP4RMNvL2xyssYIiKeCGXEnQ0sMbO1wApgkXPutfDGOrdpA7LolZ7C/BKt6RaR+JNwvg2cc9uBUR2QJWR+n3HjuDx+/+42Ko/Wkp2e4nUkEZEOE1XLAVuaPa6AgIOFq7STUkTiS9QWd1FmZ8b36caCkgqc08mERSR+RG1xA8wpLmD7/uOsKj/kdRQRkQ4T1cV91YgcUpP8WtMtInElqou7c3ICV47I4aW1ezhR3+h1HBGRDhHVxQ3Na7qP1zfx6kf7vI4iItIhor64JxR1p0+PVL1Pt4jEjagvbjNj9rh8Ptx+kPIDJ7yOIyISdlFf3AA3jM3HDBZo1C0icSAmijs3oxNT+2eycNVuAgGt6RaR2BYTxQ3N79O9+/BJlm474HUUEZGwipnivnxoNukpCdpJKSIxL2aKOyXRz3Wj83jt430cOdngdRwRkbCJmeKG5pMJ1zUGeGntHq+jiIiETUwV94i8rgzulcb8Uh0CLyKxK6aK+9M13Wt3HWZzZY3XcUREwiKmihvgS2PySPCZzo4jIjEr5oq7R5dkLh3ck7+s3k1DU8DrOCIi7S7mihua13TvP1bPu59Uex1FRKTdxWRxXzIoi8wuyZouEZGYFJPFneD3ccPYPP66qYr9x+q8jiMi0q5isrih+X26GwOO51fv9jqKiEi7itniHpCdxuiCDObrZMIiEmNitrih+UjKTypr+Gj3Ea+jiIi0m5gu7mtG5ZKc4ONZ7aQUkRgS08WdnpLIzOG9eHHNHmobmryOIyLSLmK6uAFuHl/I0dpGrv7PJfx1U6Xmu0Uk6sV8cU/u14OHvjqOpoDjjidKuPXh5XysOW8RiWIxX9wAVwzrxRv/OJ2fXTuMTypruPo/l/CDZ9aw+/BJr6OJiLSahWPqoLi42JWUlLT767aHo7UN/P7dbTy6pAyAO6cW8Q+X9CM9JdHjZCISz8ys1DlXHMq2cTHibik9JZF7Zw7mnXsu4eoROfz+3W1c8sC7/GHpDr0plYhEhZCL28z8ZrbazF4OZ6COkpfRid/MHc1L357KoOw0/u3F9Xzxt+/z+vp92oEpIhGtNSPu7wEbwxXEKyPyu/LUXRN57LZifD7j7/5Yyk3/s4zV5Ye8jiYiclYhFbeZ5QNXAY+EN443zIxLB2fz2vem8YsvDads/3G+9P+W8p2nV7Pr4Amv44mIfEaoI+7/AH4MfO4ksJndbWYlZlZSXR2d74Od4Pfx5Ym9efdHM/jupf15c8M+Lvv1e/xi0QaOnNCZ40UkMpy3uM3saqDKOVd6ru2ccw8554qdc8VZWVntFtALXZIT+MEVg3j3nhlcNzqXR5aUMf2Bd3hk8XbqGnUEpoh4K5QR9xTgWjPbAcwDLjWzP4U1VYTo1TWFB+aMYtF3pjEyvys/X7SRy3/zPovW7dUOTBHxTKvWcZvZJcA9zrmrz7VdJK/jbov3Nlfzf17ZyKZ9NYwpzOC+mYOZUNQdM/M6mohEudas404Id5hYcvHALKb2z2RhaQW/euMT5j70IUNy0rl1YiHXj84lTQfxiEgHiLsjJ9vLifpGnlu1m6eWl7Nh71FSk/xcOyqXWyYUMjK/q0bhItIqrRlxq7jbyDnH2oojPLV8Jy+t3cvJhiaG5TaPwq8bnUeXZP1RIyLnp+L2yNHaBl5YvZs/Ly9n074aOif5uXZ0Hl+eWMjwvK5exxORCKbi9phzjtW7DvPU8nJeXreH2oYAI/O7csuEQq4dlUtnjcJF5DQq7ghy5GQDf1lVwVMrytlceYwuyQlcNzqXWycWMixXo3ARaabijkDOOVaVH+LPy8tZtG4vdY0BRhVkcOuEAq4ZlUtqkkbhIvFMxR3hDp+ob16RsqKcrVXHSEtO4Poxedw6sZAhOelexxMRD6i4o4RzjpU7DvHU8p288vE+6hsDjCnM4IYxeVw8sCeFPVK9jigiHUTFHYUOHa9n4aoKnl5Rzrbq4wAUdk9l2oBMpg3IYnK/HnTtpAN8RGKVijuKOeco23+cxVv2s3hLNcu2HeB4fRN+nzG6IIOp/TOZPjCTUfkZJPjj7gRGIjFLxR1D6hsDrC4/1FzkW/ezruIwzkFacgIX9e/BtAFZTB+QpWkVkSin4o5hh47Xs3TbARZvqWbxlv2nzlSvaRWR6KbijhOaVhGJHSruOHW+aZV7Zw6mb1YXr2OKyFmouAX47LTKy+v2MqYwgz/eOdHrWCJyFno/bgGgW+ckrhqZw1UjcyjonsoDr3/Cpn1HGdxLB/mIRDNNfMaJL08spFOin0cXl3kdRUTaSMUdJzJSk5hTnM8La/ZQVVPrdRwRaQMVdxy5Y0oRDYEATy7d6XUUEWkDFXcc6ZPZmcuHZPOn5Ts5Wd/kdRwRuUAq7jhz1/S+HD7RwIJVFV5HEZELpOKOM8W9uzEqvyuPLSkjEGj/paAiEn4q7jhjZnxjWl/K9h/nrY2VXscRkQug4o5Ds4b3Ii+jE48s0dJAkWik4o5DCX4ft0/pw4qyg6yrOOx1HBFpJRV3nJo7voC05AQe1gE5IlFHxR2n0lISuXlCAa98tPfUW8OKSHRQccex26YUAfDEBxp1i0QTFXccy8voxJUjcpi3Yhc1tQ1exxGREKm449xd04qoqWvkmZW7vI4iIiFScce5kfkZTCjqzuMf7KCxKeB1HBEJwXmL28xSzGyFma01s/Vm9rOOCCYd5xtTi9h9+CSvfrzP6ygiEoJQRtx1wKXOuVHAaGCmmU0KbyzpSF8Ykk1RZmceWbydcJwRSUTa13mL2zU7FryZGPzQb3cM8fmMO6YWsbbiCCU7D3kdR0TOI6Q5bjPzm9kaoAp40zm3PLyxpKPNHptPRmoiD7+/3esoInIeIRW3c67JOTcayAcmmNnw07cxs7vNrMTMSqqrq9s7p4RZpyQ/X5nYmzc3VlK2/7jXcUTkHFq1qsQ5dxh4B5h5lscecs4VO+eKs7Ky2iufdKCvXdSbRJ+Px3VAjkhEC2VVSZaZZQSvdwIuBzaFO5h0vJ5pKVw3Opf5JRUcPlHvdRwR+RyhjLhzgHfMbB2wkuY57pfDG0u8cue0Ik42NPHn5eVeRxGRz5Fwvg2cc+uAMR2QRSLA4F7pTBuQyRNLd/CNaUUkJ/i9jiQip9GRk3KGu6b1pbqmjpfW7vU6ioichYpbzjBtQCaDstN0QI5IhFJxyxnMjDunFbFpXw1Ltu73Oo6InEbFLWd13ehcMrsk84jOkCMScVTcclbJCX6+Prk3722uZnNljddxRKQFFbd8rq9M6k1Koo9HFusweJFIouKWz9WtcxKzx+Xz/Oo9VNfUeR1HRIJU3HJOd0wpoiEQ4I/LdngdRUSCVNxyTn2zunDZ4Gz++OFOTtY3eR1HRFBxSwjumlbEoRMNLFxV4XUUEUHFLSGYUNSdkfldeWxJGYGADsgR8ZqKW87LzLhzahHb9x/nr5uqvI4jEvdU3BKSK0fkkNs1hYe1NFDEcypuCUmi38ftU4pYXnaQjyqOeB1HJK6puCVkcycU0CU5gUeWaNQt4iUVt4QsPSWRueMLeHndXvYcPul1HJG4dd4TKYi0dPuUPjz+QRl/WLqDn1w5pN1fv/JoLR9uP8C2qmM0BByBgKMx4GgKOBoDgebLpk9vf/b+z97X8jJw6jnXjMrlu5cNaPfcIh1JxS2tkt8tlVkjcnhqRTnfuWwAXZLb9iNUdbSWZdsP8OH2gyzffoDtwTPMm0Giz4ffZyT4DL/f8Jt95nZC8PFT9/tbPO4zkhL9+FrcPnCsjt+8uZmB2WnMHN6rPb4dIp5QcUur3TWtL4vW7eWZlbu4c2pRq55bdbSWD8sO8uH2A3y4/QDbq5uLOi05gQlF3bllQiGT+vZgaG46fp+1a+76xgA3/n4p9z23jjGFGWSnp7Tr64t0FAvHGU6Ki4tdSUlJu7+uRI45/72UPYdree9Hl5Dg//xdJVU1tSzf/rei3hYs6i7Bop7UtzuT+2aGpajPZlv1Ma7+3RLG9e7Gk3dMwNcBn1MkFGZW6pwrDmVbjbjlgtw5tS9//6dSXl9fyVUjc07dX11Td6qkTy/q8X26MXd8QfOIOif9nIUfLv2yuvCv1wzlJ899xKNLyrhret8OzyDSVipuuSCXD82md49UHnp/Gw4XLOqDbK06BvytqG8qbi7qYbneFPXZ3Dy+gHc/qeL+1zdxUf8eDMvt6nUkkVbRVIlcsCeX7eBfX1gPQOckP+OLujOpbw8m9e3B8Agq6rM5dLyemQ++T1pKIi99eyqdkvxeR5I4p6kS6RC3TCgkNSmBflmdGZHXNaKL+nTdOifx6zmj+cqjy/nFKxv4+fUjvI4kErLo+U2TiJPo9zF7XD5jCrtFVWl/auqATO6e3pc/fVjOWxsqvY4jErLo+20TaUc/vGIgQ3PS+fHCdVQdrfU6jkhIVNwS15IT/PzultGcqG/kh/PX6v3GJSqouCXu9e+Zxr9cNZTFW/bzxNIdXscROS8Vtwjw5YmFfGFINv/+6iY27j3qdRyRc1Jxi9B8lp9f3jiCrqmJfG/eamobdGJkiVwqbpGgHl2S+dWcUWyuPMa/v7rJ6zgin+u8xW1mBWb2jpltMLP1Zva9jggm4oWLB2Zxx5Qinli6g3d0fk2JUKGMuBuBHzrnhgKTgG+Z2dDwxhLxzo9nDmJwrzR+tGAt1TV1XscROcN5i9s5t9c5typ4vQbYCOSFO5iIV1IS/fzuljHU1Dby4wVrCcfbQoi0RavmuM2sDzAGWH6Wx+42sxIzK6murm6fdCIeGZidxj9dOYR3PqnmyWU7vY4j8hkhF7eZdQEWAt93zp2xXso595Bzrtg5V5yVldWeGUU88bXJvZkxKItfvLKRzZU1XscROSWk4jazRJpL+8/OuefCG0kkMpgZD8wZRXpKAt99WksEJXKEsqrEgEeBjc6534Q/kkjkyOySzANzRrFpXw33v/aJ13FEgNBG3FOArwKXmtma4MeVYc4lEjFmDOrJbRf14bEPynhvs/bfiPdCWVWyxDlnzrmRzrnRwY9XOiKcSKS4b9ZgBmWncc/8tRw4piWC4i0dOSkSgpREPw/eMpojJxu4d+E6LREUT6m4RUI0uFc6980czFsbq/jz8nKv40gcU3GLtMJtF/Vh+sAsfr5oA1urtERQvKHiFmkFn8/41ZyRpCYl8N2n11DXqCWC0vFU3CKt1DMthQdmj2TD3qP86nUtEZSOp+IWuQCXDcnmq5N68/DiMpZs2e91HIkzKm6RC/RPVw6hf88u/ODZNezXEkHpQCpukQvUKcnPgzeP5vCJBmY9uJgX1+7RMkHpECpukTYYltuV5755Eb3SU/ju06v52mMr2HnguNexJMapuEXaaHheV57/1hR+es1QVpcf5orfvs//fWcr9Y0Br6NJjFJxi7QDv8+4bUoRb/3gYi4d3JMHXv+Eq363mBVlB72OJjFIxS3Sjnp1TeH3XxnHY7cVc6K+iZv+Zxn3LljHoeP1XkeTGKLiFgmDSwdn8+YPpvN3F/dlwaoKLvvNeywsrdDOS2kXKm6RMElNSuAns4bw8nem0rtHKj+cv5ZbH17OtupjXkeTKKfiFgmzITnpLPz7i/jFl4bz8Z4jzPqPxfz2zc06o45cMBW3SAfw+YwvT+zN2z+8mJnDe/Hg21uY9eBilm7VUZfSeipukQ7UMy2F390yhifvmEDAOW59ZDn/+IyOvJTWUXGLeGD6wCxe//50vj2jPy+v28Nlv36PeSvKCQS081LOT8Ut4pGURD/3fHEQr35vGoOy07jvuY+Y+9AyNlfqfb7l3FTcIh7r3zONeXdP4v4bR7Kl6hhXPriYX762iZP12nkpZ2fhWFdaXFzsSkpK2v11RWLdgWN1/O9XNrFwVQU905Ip6J5Kkt9HUoKP5ITmy1PX/T6SE/2nHk86dZ/vjOckJ/hPPZ6dnkKvrilef6lyGjMrdc4Vh7JtQrjDiEjoenRJ5tc3jWL2uHyeWFrGsbpG6hsDnDjRSF1jgPqmAHUNzZf1jc0fdY1NtGZq3AwuG5zNHVP6MLlfD8wsfF+QhIWKWyQCTe7Xg8n9eoS8fWPT38q8rvGzl81l33Tq8dXlh3lqRTlvbaxkcK80bruoD9ePySMl0R/Gr0jak6ZKROJQbUMTL67Zw2MflLFpXw3dUhO5ZUIhX53cm5yunbyOF5daM1Wi4haJY845lm0/wOMf7OCtjZX4zZg1Iofbp/RhbGE3r+PFFc1xi0hIzIyL+mVyUb9Myg+c4A/LdvDsyl28tHYPowoyuGNKH2YNzyEpQQvQIolG3CLyGcfqGllYWsETS3dQtv84PdOS+eqk3tw6sZAeXZK9jhezNFUiIm0WCDje21zNYx+UsXjLfpISfFw3KpfbpxQxNDfd63gxR1MlItJmPp8xY3BPZgzuyZbKGp5YuoPnVu1mfmkFk/p25/YpRXxhSDZ+n5YTdjSNuEUkZIdP1DNv5S6eXLqDPUdqKejeia9P7sNN4wtIT0n0Ol5Ua9epEjN7DLgaqHLODQ/lRVXcIrGtsSnAGxsqefyDMlbuOERqkp/heV0xmg/wMQyfr/ny0+N7zKzF45+9TXC7ls83az6X5+R+PbhmVG7M/8fQ3sU9HTgGPKniFpHTfbz7CE8u28GOAyfAgcPhHDialxs2Xzbf5jO3g9u12JbTHjtR38TuwydJTvAxc3gv5owr4KJ+PfDF4PRMu85xO+feN7M+bQ0lIrFpeF5X7p89Kiyv7Zzjo91HmF9SwQtrdvPCmj3kZXTixnH5zBmXT0H31LB83kgX0hx3sLhf1ohbRLxS29DEmxsqebZkF0u27sc5mNy3B3OK85k1PIdOSdF9yH67LwcMpbjN7G7gboDCwsJxO3fuDCmsiEhr7Tl8kudWVTC/tIKdB07QJTmBq0fmMKe4gLGFGVH5xlmeFHdLGnGLSEdwzrGi7CDzSytYtG4vJxua6JvVmTnjCrhhbB7Z6dHz9rUqbhGJO8fqGnll3V7ml+5i5Y5D+AwuHpjFTcUFXDYkO+IP22/vVSVPA5cAmUAl8G/OuUfP9RwVt4h4aXv1MRaUVvDcqt3sO1pLt9RErhudx5zifIbldvU63lnpkHcREaAp4Fi8pZr5pRW8ub6S+qYAw3LTmTMun+vH5JGRmuR1xFNU3CIipzl0vJ4X1+5hfukuPt59lKQEH1cO78XNEwqZWNTd8x2aKm4RkXNYv+cI81bs4vnVu6mpa6RvZmfmji/gxnH5ZHr0DogqbhGREJysb2LRR3uZt6Kckp2HSPQblw/NZu74Qqb1z+zQIzRV3CIirbSlsoZ5K3fx3KoKDp1oIC+jE3PHF3BTcQG9uoZ/WaGKW0TkAtU1NvHG+krmrSzng60H8BnMGNSTmycUMmNQFgn+8CwrVHGLiLSDnQeO88zKXcwvraC6po7s9GTmjCtg7viCdn+fFBW3iEg7amgK8NdNVcxbUc57m6sJOJjaP5ObJxRwxdBe7XJwj4pbRCRM9hw+ybMlu3h25S72HKmle+ckbhybx9zxhfTv2eWCX1fFLSISZk0Bx/tbqpm3opy3N1bRGHBMLOrOk3dOIDmh9e9UqHNOioiEmd9nzBjUkxmDelJVU8vC0t3sPHD8gkq7tVTcIiJt1DMthX+4pF+Hfb7IfrssERE5g4pbRCTKqLhFRKKMiltEJMqouEVEooyKW0Qkyqi4RUSijIpbRCTKhOWQdzOrBna2+wu3TSaw3+sQIVLW8ImmvNGUFaIrbyRm7e2cywplw7AUdyQys5JQ3wfAa8oaPtGUN5qyQnTljaasZ6OpEhGRKKPiFhGJMvFU3A95HaAVlDV8oilvNGWF6MobTVnPEDdz3CIisSKeRtwiIjEh6ovbzGaa2SdmttXM7jvL48lm9kzw8eVm1id4/+VmVmpmHwUvL43kvC0eLzSzY2Z2TyRnNbORZrbMzNYHv8cpkZrXzBLN7A/BnBvN7CcRkHW6ma0ys0Yzm33aY183sy3Bj69HalYzG93iZ2Cdmc0Nd9a25G3xeLqZVZjZf3VE3gvinIvaD8APbAP6AknAWmDoadt8E/jv4PWbgWeC18cAucHrw4HdkZy3xeMLgPnAPZGaleYTdKwDRgVv9wD8EZz3VmBe8HoqsAPo43HWPsBI4Elgdov7uwPbg5fdgte7RWjWgcCA4PVcYC+QEQE/B2fN2+LxB4GngP8KZ9a2fET7iHsCsNU5t905Vw/MA647bZvrgD8Ery8ALjMzc86tds7tCd6/HuhkZsmRmhfAzK4HyoJ5w60tWa8A1jnn1gI45w4455oiOK8DOptZAtAJqAeOepnVObfDObcOCJz23C8CbzrnDjrnDgFvAjMjMatzbrNzbkvw+h6gCgjpABMv8gKY2TggG3gjzDnbJNqLOw/Y1eJ2RfC+s27jnGsEjtA8AmzpRmCVc64uTDnPyBIUcl4z6wLcC/wszBnPyBHUmu/tQMCZ2evBP0l/HOF5FwDHaR4RlgO/cs4d9DhrOJ57Idrl85nZBJpHwNvaKdfnueC8ZuYDfg2EfRqyreL+nJNmNgz4Jc2jxEj2U+C3zrljwQF4JEsApgLjgRPA28EzWL/tbazPNQFoovnP+W7AYjN7yzm33dtYscHMcoA/Al93zp0xyo0g3wRecc5VRPrvWLQX926goMXt/OB9Z9umIvincFfgAICZ5QN/Ab7mnAv3SKBllk+1Ju9EYLaZ3Q9kAAEzq3XOhWsHSluyVgDvO+f2A5jZK8BYIJzF3Za8twKvOecagCoz+wAopnn+2Kus53ruJac99912SfX5n+9Cs2Jm6cAi4J+dcx+2c7azaUveycA0M/sm0AVIMrNjzrkzdnB6zutJ9rZ80Pwfz3agiL/tiBh22jbf4rM7pJ4NXs8Ibn9DNOQ9bZufEv6dk2353nYDVtG8oy8BeAu4KoLz3gs8HrzeGdgAjPQya4ttn+DMnZNlwe9xt+D17hGaNYnm/6y/H85/+/bKe9pjtxHBOyc9D9AO/1BXAptpnjv75+B9/wu4Nng9heZVGFuBFUDf4P3/QvO85poWHz0jNe9prxH24m5rVuArNO9E/Ri4P8J/FroE718fLO0fRUDW8TT/5XKc5r8K1rd47h3Br2ErcHukZg3+DDSc9js2OlLznvYatxHBxa0jJ0VEoky0ryoREYk7Km4RkSij4hYRiTIqbhGRKKPiFhGJMipuEZEoo+IWEYkyKm4RkSjz/wE5cYjEvUOPqQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f1cc4c10090>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def log_normal(mean, log_err):\n",
    "    return np.exp(np.log(mean) + randn()*log_err)\n",
    "\n",
    "class Measurement:\n",
    "    def __str__(self):\n",
    "        return \"\"\"\n",
    "        Capacitance, C: %g\n",
    "        Initial voltage, V0: %g\n",
    "        Times, t: %s\n",
    "        Voltages: %s\n",
    "        \"\"\" % (self.C, self.V0, str(self.t), str(self.V))\n",
    "\n",
    "class Simulation:\n",
    "    def __init__(self):\n",
    "        \"Randomly select true values\"\n",
    "        self.R = log_normal(5.0, 0.5)              # resistance \n",
    "        self.C = log_normal(0.01, 0.1)             # capacitance\n",
    "        self.V0 = log_normal(7.0, 0.1)             # begin voltage\n",
    "    \n",
    "    def simulate_measurement(self, n_samples = 15):\n",
    "        \"Simulate measurement of V(t) = V0 * exp(-t / (R*C))\"\n",
    "        sample_rate = 1e-2\n",
    "        # multiplicative independent errors in voltage measurement\n",
    "        log_v_error = 0.1\n",
    "        # additive error in capacitance measurement\n",
    "        c_error = 0.001\n",
    "        # random clock drift\n",
    "        clock_drift = randn() * 1e-6\n",
    "        \n",
    "        m = Measurement()\n",
    "        m.V0 = self.V0 + randn() * (1 + log_v_error * randn())\n",
    "        m.C = self.C + randn() * c_error\n",
    "        \n",
    "        # percieved measurement time, t\n",
    "        m.t = np.arange(1, n_samples+1) * sample_rate\n",
    "        \n",
    "        # actual measurement time with clock drift\n",
    "        t_true = m.t * (1 + clock_drift)\n",
    "        \n",
    "        # measured voltages V(t)\n",
    "        m.V = self.V0 * np.exp(-t_true / (self.R * self.C)) * (1 + randn(n_samples) * log_v_error)\n",
    "        \n",
    "        return m\n",
    "    \n",
    "    def __str__(self):\n",
    "        return \"\"\"\n",
    "        Resistance: %g\n",
    "        Capacitance: %g\n",
    "        Initial voltage: %g\n",
    "        \"\"\" % (self.R, self.C, self.V0)\n",
    "\n",
    "actual = Simulation()\n",
    "meas = actual.simulate_measurement()\n",
    "\n",
    "print('--- Actual ---')\n",
    "print(actual)\n",
    "print('--- Measurement ---')\n",
    "print(meas)\n",
    "plt.plot(meas.t, meas.V);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reconstruction of R assuming all errors are independent\n",
    "\n",
    "Simple, but not really accurate\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xd8VGW+x/HPk95DAgESAgmggEgnQADXugXRK6ioq6AgSmR1i9tcy3VX3d1717LVXZUOKgKKvYtldaUn0nuH0CGQSvpz/5iRC4gQksycmZzv+/WaV05mzsz5MUzme57zPOc5xlqLiIi4V4jTBYiIiLMUBCIiLqcgEBFxOQWBiIjLKQhERFxOQSAi4nIKAhERl1MQiIi4nIJARMTlwpwuoC5atGhhMzMznS5DRCSo5OXlHbLWppxtvaAIgszMTHJzc50uQ0QkqBhjdtRlPR0aEhFxOQWBiIjLKQhERFxOQSAi4nIKAhERl1MQiIi4nM+CwBgz1RhzwBiz+oT7ko0x84wxm7w/k3y1fRERqRtftgimA0NOue9+4BNr7fnAJ97fRUTEQT4LAmvtF0DBKXcPA2Z4l2cAw321ffGtmyYs5KYJC50uQ0Qagb/7CFpZa/d6l/cBrfy8fREROYVjncXWWgvYb3vcGJNjjMk1xuQePHjQj5WJiLiLv4NgvzEmFcD788C3rWitnWitzbLWZqWknHXOJBERqSd/B8FbwGjv8mjgTT9vX0RETuHL4aOzgIVAZ2NMvjHmDuBPwPeMMZuA73p/FxERB/lsGmpr7c3f8tAVvtqmiIicO51ZLCLicgoCERGXUxCIiLicgkBExOUUBCIiLqcgEBFxOQWBiIjLKQhERFxOQSAi4nIKAhERl1MQiIi4nIJARMTlFAQiIi6nIBARcTkFgYiIyykIROrppgkLuWnCQqfLEGkwBYGIiMspCEREXE5BICLicgqCJk7HsUXkbBQEIiIupyAQEXE5BYGIiMspCEREXE5BICLicgoCERGXUxCIiLicgkBExOUUBCIiLqcgEBFxOQWBiIjLKQhERFxOQSAi4nKOBIEx5ufGmDXGmNXGmFnGmCgn6hAREQeCwBjTBvgpkGWt7QaEAj/0dx0iIuLh1KGhMCDaGBMGxAB7HKpDRMT1/B4E1trdwFPATmAvUGit/ejU9YwxOcaYXGNM7sGDB/1dpoiIazhxaCgJGAa0B9KAWGPMqFPXs9ZOtNZmWWuzUlJS/F2miIhrOHFo6LvANmvtQWttFfAaMMiBOkREBGeCYCeQbYyJMcYY4ApgnQN1iEvous0iZ+ZEH8FiYC7wFbDKW8NEf9chIiIeYU5s1Fr7O+B3TmxbREROpjOLRURcTkEQQHQsW0ScoCAQEXE5BYGIiMspCEREXE5BICLicgoCERGXUxCIiGtoZN7pKQhERFxOQSAi4nIKAhERl1MQiIi4nIJAROpNna9Ng4JARMTlFAQiIg0U7C0jBYFIgAn2LxUJPgoCERGXUxCIiLicgkBExOUUBCIScNRP4l9NOgj0YRIRObswpwsQCVZHyyoprqimuqaWsNAmvU8lTZw+vSLnyFrL059sYsP+EvYcLefhN9dgrXW6LJF6U4tA5ByUVVbz61dW8u6qvTSPjSAiLIRZS3aSlhjFT6443+nyROpFQSBSR7sKyhj3fC4b9xfz4NAufLx2PwDZHZrz53kbaZ0YxQ1ZbR2uUuTcKQikXqpqap0uwa8WbjnM3TPzqK61TB3Tj0s7t+STdQcAePz6HhwsruCB11bRMiGKSzqlOFzt6X09cGLOXQMdrkTqwp//X+ojkHO2bOcRVuYXsnZPEeVVNU6X41PWWp5fuJ1RUxaTHBvBm/cM5tLOLU9aJyIshGdH9eH8VvH86MU8Vu8udKZYkXpSENSDm4elfrbhALdMWowxUF5dy8Qvtjpdks9UVtfy4Our+O2ba7i0Uwpv3DOYDilxp103Piqc6bf3IykmgjHTlrKroMzP1YrUn4JA6uzVvHzGzcilfYtYuqUlkhwbwb8+29wkv/QOFldwy6RFzFqyi3su68jE27KIjwo/43NaJUQxY2w/qmpqGT1tCUdKK/1UrUjDKAjkrKy1TPh8C798ZQX92ycz565sIsJCaJccQ2iI4dG31zTatgKhtbUy/yjX/PNLVu8p5Ombe/PrH3QhNMTU6bnntYxn8ugs8o8c444ZS5v8oTNpGhQEcka1tZY/vLuO/31/PVf1SGXa7f2O7xlHhoVw73fP5+N1B46PoAl2byzbzQ3PLSTEGF790SD+q2faOb9Gv8xk/n5TL5btOspPZy2jplbnGEhgUxDIt6qsruXnLy9nypfbGDMok6d/2JvIsNCT1rl9cHvObxnHI2+v4Vhl8O791tRa/ve9ddw7Zzk92zbjrR8P5sK0xHq/3pXdU/nt1V35aO1+HnlLJ5xJYHMkCIwxzYwxc40x640x64wxGs8WYEoqqrljxlLeXL6H+4Z05nf/1ZWQ0xweCQ8N4bFh3cg/coxn/73ZgUobrvBYFWOnL2XCF1sZld2OmXcOoHlcZINf9/bB7cm5uAMvLNrBs59vaYRKRXzDqfMI/g58YK0dYYyJAGIcqkNO41BJBbdPW8ravUU8OaLHWU+SGtixOcN7pfHc51u5tk867VvE+qnShtt8oIRxz+eyq6CMP17bjZEDMhr19e8f0oW9heU88cEGWidEcV2f9EZ9fZHG4PcWgTEmEbgYmAJgra201h71dx1yejsPlzHi2QVsOlDMpNv61vlM2QeHXkBkWAi/C6LDIJ+s28+1/5pP0bEqXhqX3eghABASYnjqhh4M7NCc++au5MtNhxp9G1I3FdU1HC6p4GBxBev3FVHtspMiz8SJFkF74CAwzRjTE8gDfmatLXWgFjnB6t2FjJm2lOraWmbemU3fjKQ6P7dlQhQ//14nHntnLR+u2ceQbqk+rLRhrLU88+8tPPXRBi5MS2DCrVm0aRbts+1FhoXy3K19ufG5hYx/MY85d2U3qP9Bzs2hkgpmLtrJC4t2cKikAoAhf/sP0eGhXJiWQI/0ZvRIT6RHeiKZzWNPewi0qXMiCMKAPsBPrLWLjTF/B+4HHj5xJWNMDpAD0K5dO78X6TYLNh8i54U8EqLCmJ0zkPNaxp/za9w2MIOXc3fx2NtrubhTCjERgTeDybHKGn49dwXvrNzLNT3TePz6HkRHhJ79iQ2UGB3O9LH9uO6ZBdw+bSmv3T2I9CQdEfWlNXsKmTZ/O28t30NlTS2XdU5hb2E5EaEhjL2oPSvyj7Iqv5CXluxg6nxP6yA+KozubRLpnp5Iz/RmdG+TSHpSNMY07XBw4i81H8i31i72/j4XTxCcxFo7EZgIkJWVFRzHGoLUOyv38PM5y2nfIpYZY/uTmli/veOw0BD+MLwbI55byNOfbuY3Q7o0cqUNs/voMXKez2Xt3iJ+M6QL4y/p4Nc/8NTEaKbf3p8Rzy1gzLSlzB0/kGYxEX7bvhvU1Fo+WbefqfO3sWhrAdHhodzUry1jBmfSMSXu+Dkqw3u3YXjvNgBU19Sy6UAJq/ILPeGwu5CpX26jqsbztZMcG+FpMbRJPN56aJkQ5di/0Rf8HgTW2n3GmF3GmM7W2g3AFcBaf9chHtPnb+PRd9aSlZHE5Nv6kRhz5rNnzyYrM5kRfdOZ/J+tXN8nnfNann5KBn8rOlbFNU9/SWV1LVNH9+OyLi3P/iQf6Nw6nom3ZjF66hJyns/j+Tv6ExXu+xZJU1dcXsXLufnMWLCdnQVltGkWzYNDu3BTVruzfqbDQkO4IDWBC1ITuLGfp0+sorqGDfuKWZFfyKr8o6zML+SLjQf5+pSQ1glR3lZDIt3Tm1FVU0t4EF+cyKm2+0+Amd4RQ1uB2x2qw7WstTz10Qb+9dkWvt+1Ff+4uXejfSHdf2UXPlqzj9++uZqZdw5wvFldUFrJ5gMltG8Ry6TRWXT8lvmC/GVgx+Y8dWNPfjprGb94eTn/vLmPK49LN4Ydh0uZvmA7r+TmU1JRTb/MJB64sgvf69qqQVeNiwwL9e79NwM8gwjKKqtZu6fopHCYd8KJlLGRoUz8YgtXdkulbXJwHfZzJAistcuBLCe2LZ6m8IOvr+Ll3Hxu7t+O3w+7sFEvtdgiLpJf/6AzD7+5hndW7q3X2bmN5bP1B9h8oITYyDBev2cwidENa/E0lmt6pnGgqJw/vLuO3yes5bdXd3U8MOvDWkt5VS1F5VUknGUupsbc5sKth5n65XY+Wb+fsBDD1T3SuH1wpveL2zdiIsLIykwmKzP5+H1F5VWszi/kvldXcqS0kv95bz3/8956eqQnMrR7KkO7pdKueeCHQuD15olPHaus4ccvfcUn6w/wsyvO597vnu+TL6BbBmQwJ3cXf3h3LZd1aUlcpP8/agu2HGL8i3lER4TSuVVcwITA1+78Tgf2HC1n6vxtpCVGM+7iDk6XVGcV1TW8sWw3K3cXUl5VS49HPqJFXATtW8R6b3HHlzOaxzRKa7O8qoa3Vuxh6pfbWL+vmOTYCH582XmMys6glUPH7BOiwhl0XgvaNIumTbNonhzRk/dW7+X9VXv50/vr+dP76+nWJoGh3VO5qnsqGc0D8xwbBYGLHCmt5I4ZS1m26yh/GN6NUdmNP27+a6Ehht8P68Z1zy7g7x9v5KGruvpsW6eTt6OAO2fkktE8hrjIsIC9uPx/X3UB+4vK+eN762iZEMmwXm2cLumMisqrmLloJ9Pmb+NAcQUxEaFkNo/h5v7t2HaolK2HSvlsw0Fezs0//hxjIC0xmg4psScEhefWpln0Wf9vDhSV8+KiHcxcvJPDpZV0aR3PE9f34JpeaQHXv9KueQzjL+nI+Es6squgjPdX7+XdVft44oMNPPHBBrqmJnBVj1SGdk8NqBMvFQQusefoMW6buoSdBWU8O7KPX8b5926XxA/7tWXq/O2M6NuWzq3PfUhqfazeXciYqUtpGR/Ji3cM4Cezlvllu/UREmL48409OVhSwa9eWUFKfMOntvCFfYXlTJu/jZmLd1JSUc13zm/BX27sxT8+2Ygxhrsu6XjS+sXlVWw/VMa2w6VsO1jKtkMlbDtUyutf7aa4ovr4euGhhnbJMbRvEXdSUFRW11JVU8sv5izn7ZV7qK61XNGlFWMHZzKwY/OgOIzWNjmGnIs7knNxR/KPlPHB6n28u2ovT364gSc/3ECX1vFc1T2VoT1SHe+3UhC4QFllNdc9s4DSimqeH9uf7A7N/bbt+37QhfdX7+PhN1czJyfb53/AG/YVc+uUxSREhzNzXHZQDPOLCg9l0q1ZjHhuAXc9n0dmi5iAOQdj0/5iJn6xlTeW76am1nJVjzTuurgD3dp4Toh7+tNNp31efFQ43dM94/FPZK3lcGkl2w55A+J4UJTyxaaDVFaffLbvtkOljByQwZhBmWQG0B70uUpPiuHO73TwHg48xvur9/Heqr38ed5G/jxvI51bxXsOH/VoXa9zeBoqMD5t4jPF5VVs2F9C89gIXh4/kAtSE/y6/aTYCH4zpAsPvLaKN5bv5trevptrZ9uhUkZOXkx4aAgvjRvg07OFG1tiTDjTx/bnumfms35fMRf6+f/pVEu3FzDh8y18vO4AUeEh3NK/HXd+p0ODR8MYY2gRF0mLuEj6ndDpCp4pz/cUHmPboVL++43VGOCtn1zkt05of0lrFs0dF7Xnjovas7fwGO+v2sf7q/fyt0828tePN9KpVRxXdkulrLLabzsECoImbN7a/azbV0xkWAiv/miQY0Pabspqy+ylu/jju+u5vEsrn3Ta5h8pY+SkRdRay6xx2QHbKXcmbZpFM21Mf65++j8szy/k+3/9nL4ZyfTNSCIrI4mM5jE+bVHV1lo+XrefCV9sJW/HEZJiwrn3u+dz28BMkmN9f+JbSIghPSmG9KQYWntbck0tBE6VmhjN2IvaM/ai9uwvKuf9VXt5b/U+/vHpJqyFqPAQNu4vplMr37YSFARN1JylO3ngtVXERITSuVW8o+OaQ0IMfxjWjWv+9SV/nbeRR665sFFff39RObdMWkxJRTWzcrI538d/NL7UNS2BbmmJFJRVktYsmndX7mHWkp0AtIiLoE+7JLIyk+ibkUy3NgnfuD5EfXw9AmjCF1vZerCU9KRoHht2ITf0beuX6TfEo1VCFGMGt2fM4PYcKCrnhgkLOVJa6ZeWrYKgifl6QrUnP9zAxZ1SKKuorvNlFn2pe3oiowZk8PzC7dyQld5ok64dLqlg5OTFHC6p4MU7BzSJydyiI0JpE+GZjqK21rL5YAm524+Qu6OAvB1H+Mh7ElNEWAg92iTSNzOJrIxk+rRrdk7XUTh1BNCFaQn84+beDO3WOmBHWblFy4QoWntvsX4Yel2vLRhjQoCbrbUzG7keaYDaWstj76xl+oLtDO+VxhMjenLrlMVnf6Kf/Or7nXlv1V4efmM1c8cPavDZtIVlVYyasoT8I2VMv70/vdvVfbbUYBESYujUKp5OreK5ZYBn8sWDxRXk7ThC3o4CcnccYeqX25jw+VYAOrSIpW9GkudwUmYSHVrEfeN9/rYRQIPPC47RONL4zhgExpgE4B6gDfAWMA/4MfBLYAWgIAgQldW1/PKVFby9Yg93XNSeh4ZeEHDTFiTGhHP/lV349dyVzP0qnxvreK2D0ympqGb0tCVsOVDCpNFZfh0J5bSU+EiGdGvNkG6tAc+JVqt2F5K73RMOH6/bzyt5nnH8zWLC6dPOEwxHyyopKK3kO098Sk2t5eoeaeScMAJI3OtsLYIXgCPAQuBO4EHAAMO900RIACipqGb8C3l8ufkQ91/Zhbsu9u+smufi+j7pzFm6iz+9v57vd21Vr9k3j1XWMHb6UlbtLuTZkX24pFOKDyoNHlHhofTLTPaOwumItZath0rJO+Fw0qfrDwAQYuDW7IxGGQEkTcfZgqCDtbY7gDFmMrAXaGetLfd5ZVInh0oqGDt9KWv21O2ykk4LCTH8fng3rn76S578cAN/vLb7OT2/orqGnBdyWbq9gL/d1IvvX9jaR5UGL2MMHVPi6JgSd3w2zYLSSm6ZtIiYiFAeHdbN4Qol0JytR6jq6wVrbQ2e6wgoBALEroIybnhuIRv3FzPx1rpfVtJpF6QmcNvADF5aspOV+XW/SmlVTS0/fmkZ/9l0iMev6xHw0zEEkuTYCBKjw4N6qmTxnbN9KnoaY4q8t2Kgx9fLxpgifxQop7d2TxHXPbuAgtJKZt45gCsuaOV0Sefk59/rRIu4SB5+YzU1tWe/7lBNreUXL69g3tr9PHrNhcf3dEWk4c4YBNbaUGttgvcWb60NO2HZ2VMfXWzR1sPcNGEhocbwyviB9M1IPvuTAkxCVDgPDb2AFfmFzF6684zr1tZaHnhtJW+v2MNvhnRh9KBM/xQp4hJqJwaZD1bv47apS2iZEMmrdw/y+RmHvjSsVxrZHZJ54oMNFJRWnnYday2Pvr2Gl3Pz+enl5/GjSzuedj0RqT8FQRB5afFO7p6ZR9fUBOaOHxRUc+mcjjGGx4Z1o7SimsffX/+Nx621PP7BBmYs3MGdF7Xn59/r5ECVIk2fgiAIWGv5xyebePD1VVzcKYWXxg0gyQ9zv/hDp1bxjL2oPXNyd/HVziMnPfbPTzfz3OdbGDmgHQ9ddUHADokVCXYKggBXU2v53Vtr+Mu8jVzXuw2TbssKmCmKG8vPrjif1glRPPzGaqz1dBxP/s9W/jxvI9f1acPvh3VTCIj4kIIggFVU1/DTWct4fuEOci7uwFM39GySw/9iI8N4+OqurNlTxP7iCvZ7r+V7VfdUnri+R8CdIS3S1DStXcsmpLi8irteyGPBlsM8OLQLORc37U7Sod1bc9F5LViw5RC1Fi7v0pK/3tRLk5+J+IH+ygLQweIKfjhxEUu2FfCXG3s2+RAAT8fxo8MuxAIJUWE8M7IPEWH6eIr4g1oEAaa8qoYRzy3gQFEFk0ZncVnnlk6X5DcdU+Lold6M8FATcBclF2nKFAQBpLSimg37i4mNDGPmuAH0aYLTKp+NWgHiS3PuGuh0CQFJQRAgcrcXsHZvEWEhIcwdP9CRC1iLSP0Ee8AoCAJA3o4CRk9dQkRoCF1SExQCIuJXCgKH5e04wuipS2mZEEWz6HAdGhERv2vS3zoHiyvYVxi4s2Z/tfMIo6cuoUVcBLPGZSsERMQRTbZFYK31XJqvrIpH3lrDw1d3DYiLuH9t2c4jjJ6yhOZxEczKyaZ1YpTTJYlIAPFnv0OT3QU1xnBeyzhaJ0QxfcF2xj2fS0lFtdNlAbB811Fum7KEpFhPSyA1MbgnjxOR4NZkWwTgCYOM5jHcc/l5PPLWGm54biFTRmeR5uCsnSt2HeXWKYtJio1gdk62o7VIwwT7SJHG4Kv3QO+tfzXZFsGJbs3OYOqYfuwqKGP4v+azKr/QkTpW5h9l1JTFNIsJZ5ZCQEQChCuCAOCSTim8+qNBhIeGcOOEhXy4Zp9ft78qv5BRkxeTGB3OrHHZQX8tARFpOhwLAmNMqDFmmTHmHX9ts3PreN64ZzCdW8cz/sU8Jn2x9fi0x760enchIycvIiE6nNk52aQnxfh8myIideVki+BnwDp/bzQlPpLZOdkM7ZbKH99bx0NvrKaqptZn2/OEwGLiozwtAYWAiAQaR4LAGJMOXAVMdmL7UeGhPH1zb+6+tCMvLd7J2OlLKSqvavTtfB0CcZFhzM7Jpm2yQkBEAo9TLYK/AfcB37orbozJMcbkGmNyDx482OgFhIQY7hvShSdG9GDhlsNc/8wCdhWUNdrrr9lTyKgpi4mNCGXWOIWAiAQuvw8fNcZcDRyw1uYZYy79tvWstROBiQBZWVk+O5B/Y1Zb0pOiGf9CHtc+M5+Jt2U1eNbPtXuKGDl5MTHhoczOGUi75goBqTsNyRR/c6JFMBi4xhizHZgNXG6MedGBOo4b1LEFr98zmNjIMG6euIh3Vu6p92ut21vEyMmLiA4PZVZOtkJARAKe34PAWvuAtTbdWpsJ/BD41Fo7yt91nKpjShyv3z2Y7m0S+fFLy/jXZ5vPeUTR+n2elkBkmOdwUEbzWB9VKyLSeFxzHkFdJMdGMHPcAIb3SuPJDzfwq1dWUlldtxFFG/YVc8ukxYSHGmbnZJPZQiEgIsHB0SkmrLX/Bv7tZA2nigwL5a839SKzRSx/+3gT+UfKmHBrX5rFRHzrczwhsMgbAgMVAiISVJr0XEP1ZYzh3u92IrN5LPfNXcm1zyxg6ph+tD/NF/zG/Z4QCA0xzBqXfdp1xFnqJBU5Mx0aOoPhvdswc9wAjpZVcu0z81myreCkxzedGAI52XRIiXOoUhGR+lMQnEW/zGTeuGcwybERjJy8iNe+ygfgWGUNN09ajDGeEOioEBCRIKUgqIOM5rG8/qPBZGUk84uXV7D9cCnr9hUBMGucQkBEgpuCoI4SY8KZMbY/N2als7+oAoDZOQM4r6VCQESCmzqLz0FEWAiPX9+Dr3YeJTYilPNaxjtdkohIgykIzpExhuax3z6UVEQk2OjQkIiIyykIRERcTkEgIuJyCgIREZdTEIiIuJyCQETE5TR8VAKKJogT8T+1CEREXE5BICLicgoCERGXUxCIiLicgkBExOU0akjqRaN7RJoOtQhERFxOQSAi4nIKAhERl1MQiIi4nIJARMTlFAQiIi6nIBARcTkFgYiIyykIRERcTkEgIuJyCgIREZdTEIiIuJzfg8AY09YY85kxZq0xZo0x5mf+rkFERP6fE7OPVgO/tNZ+ZYyJB/KMMfOstWsdqKXJ0yyhInI2fm8RWGv3Wmu/8i4XA+uANv6uQ0REPBztIzDGZAK9gcVO1iEi4maOXZjGGBMHvArca60tOs3jOUAOQLt27eq1DR0WERE5O0daBMaYcDwhMNNa+9rp1rHWTrTWZllrs1JSUvxboIiIizgxasgAU4B11tq/+Hv7IiJyMidaBIOBW4HLjTHLvbehDtQhIiI40Edgrf0SMP7eroiInJ7OLBYRcTkFgYiIyykIRERcTkEgIuJyCgIREZdTEIiIuJyCQETE5Ryba0i+SXMjiYgT1CIQEXE5BYGIiMspCEREXE5BICLicgoCERGXUxCIiLicgkBExOUUBCIiLqcgEBFxOQWBiIjLKQhERFxOQSAi4nIKAhERl1MQiIi4nIJARMTlFAQiIi6nIBARcTldoawedCUxEWlK1CIQEXE5BYGIiMspCEREXE5BICLicgoCERGXUxCIiLicgkBExOUUBCIiLqcgEBFxOWOtdbqGszLGHAR2OF3HKVoAh5wuoo6CqVYIrnqDqVYIrnqDqVYIzHozrLUpZ1spKIIgEBljcq21WU7XURfBVCsEV73BVCsEV73BVCsEX70n0qEhERGXUxCIiLicgqD+JjpdwDkIplohuOoNplohuOoNploh+Oo9Tn0EIiIupxaBiIjLKQhOYYwZYozZYIzZbIy5/zSPRxpj5ngfX2yMyfTe/z1jTJ4xZpX35+WBXO8Jj7czxpQYY34VyLUaY3oYYxYaY9Z43+OoQK3XGBNujJnhrXOdMeaBAKj1YmPMV8aYamPMiFMeG22M2eS9jfZ1rQ2p1xjT64TPwUpjzE2BWusJjycYY/KNMf/0da31Zq3VzXsDQoEtQAcgAlgBdD1lnbuB57zLPwTmeJd7A2ne5W7A7kCu94TH5wKvAL8K1FrxXElvJdDT+3tzIDSA670FmO1djgG2A5kO15oJ9ACeB0accH8ysNX7M8m7nBQA7+231dsJON+7nAbsBZoFYq0nPP534CXgn758XxtyU4vgZP2BzdbardbaSmA2MOyUdYYBM7zLc4ErjDHGWrvMWrvHe/8aINoYExmo9QIYY4YD27z1+lpDav0+sNJauwLAWnvYWlsTwPVaINYYEwZEA5VAkZO1Wmu3W2tXArWnPPcHwDxrbYG19ggwDxjiw1obVK+1dqO1dpN3eQ9wADjrCVNO1ApgjOkLtAI+8mGNDaYgOFkbYNcJv+d77zvtOtbaaqAQzx7qia4HvrLWVviozm/U4lXneo0xccBvgEd9XOM36vA6l/e2E2CNMR96m+D3BXi9c4FSPHurO4EaOJjLAAACs0lEQVSnrLUFDtfqi+fWV6Ns0xjTH89e+pZGqut06l2rMSYE+DPg88OuDaWL1zcyY8yFwON49mID2SPAX621Jd4GQiALAy4C+gFlwCfGmDxr7SfOlvWt+gM1eA5dJAH/McZ8bK3d6mxZTYcxJhV4ARhtrf3GnniAuBt4z1qbH+h/YwqCk+0G2p7we7r3vtOtk+9t+icChwGMMenA68Bt1lpf7qWcWsvXzqXeAcAIY8wTQDOg1hhTbq31VYdWQ2rNB76w1h4CMMa8B/QBfBkEDan3FuADa20VcMAYMx/IwnP83alaz/TcS0957r8bpaozb7O+9WKMSQDeBR6y1i5q5NpO1ZBaBwLfMcbcDcQBEcaYEmvtNzqcHed0J0Ug3fAE41agPf/fMXThKevcw8kdhC97l5t5178uGOo9ZZ1H8H1ncUPe2yTgKzwdr2HAx8BVAVzvb4Bp3uVYYC3Qw8laT1h3Ot/sLN7mfY+TvMvJTr+3Z6g3As8OwL2+rLExaj3lsTEEcGex4wUE2g0YCmzEc9zxIe99jwHXeJej8Iyy2QwsATp47/9vPMeFl59waxmo9Z7yGj4PgobWCozC06m9GngiwD8Lcd7713hD4NcBUGs/PC2rUjytljUnPHes99+wGbg9QN7b09br/RxUnfJ31isQaz3lNcYQwEGgM4tFRFxOo4ZERFxOQSAi4nIKAhERl1MQiIi4nIJARMTlFAQi9WSMaeY9WUgkqCkIROqvGZ5pBESCmoJApP7+BHQ0xiw3xjzpdDEi9aUTykTqyXshmnestd0cLkWkQdQiEBFxOQWBiIjLKQhE6q8YiHe6CJGGUhCI1JO19jAw3xizWp3FEszUWSwi4nJqEYiIuJyCQETE5RQEIiIupyAQEXE5BYGIiMspCEREXE5BICLicgoCERGX+z/RXIdVP3ww5AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f1cc2b1add0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "reconstruction: R = 5.49677 +- 0.42 (actual: 6.3)\n"
     ]
    }
   ],
   "source": [
    "class Reconstruction:\n",
    "    # assumed error levels in the measurement\n",
    "    c_error = 0.001 # standard deviation\n",
    "    clock_drift = 1e-6 # clock_drift standard deviation\n",
    "    \n",
    "    # assume multiplicative error, i.e., v_meas = v_true * (1 + log_v_error)\n",
    "    # instead of additive: v_meas = v_true + v_error\n",
    "    # avoids negative logs\n",
    "    log_v_error = 0.1\n",
    "    \n",
    "class SimpleReconstruction(Reconstruction):\n",
    "    \n",
    "    def reconstruct_intervals(self, meas):\n",
    "        \"Compute R and confidence intervals for each t separately using interval arithmetic\"\n",
    "        log_voltage = np.log(meas.V0/meas.V)\n",
    "        r_mid = meas.t / (meas.C * log_voltage)\n",
    "\n",
    "        r_min = (meas.t * (1-self.clock_drift)) / ((meas.C + self.c_error) * (log_voltage + 2*self.log_v_error))\n",
    "        r_max = (meas.t * (1+self.clock_drift)) / ((meas.C - self.c_error) * (log_voltage - 2*self.log_v_error))\n",
    "\n",
    "        r_err = np.maximum(r_mid - r_min, r_max - r_mid)\n",
    "\n",
    "        return r_mid, r_err\n",
    "\n",
    "    def compute(self, meas):\n",
    "        \"Reconstruct final R and confidence intervals assuming the errors for each t are independent (not true)\"\n",
    "        r_recs, r_rec_errs = self.reconstruct_intervals(meas)\n",
    "\n",
    "        # inverse-variance weighting\n",
    "        r_rec_var = 1.0 / np.sum(1.0/r_rec_errs**2)\n",
    "\n",
    "        r_rec = np.sum(r_recs / r_rec_errs**2) * r_rec_var\n",
    "        r_rec_err = np.sqrt(r_rec_var) # this assumes the errors are independent (not really true)\n",
    "        \n",
    "        return r_rec, r_rec_err\n",
    "    \n",
    "    \n",
    "rec = SimpleReconstruction()\n",
    "\n",
    "r_recs, r_rec_errs = rec.reconstruct_intervals(meas)\n",
    "plt.errorbar(meas.t, r_recs, r_rec_errs)\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('R')\n",
    "plt.show()\n",
    "\n",
    "r_rec, r_rec_err = rec.compute(meas)\n",
    "print(\"reconstruction: R = %g +- %.2g (actual: %.2g)\" % (r_rec, r_rec_err, actual.R))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reconstruction of R using linear regression\n",
    "\n",
    "The dependency between the errors is handled more correctly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "reconstruction: R = 5.88534 +- 0.8 (actual: 6.3)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEKCAYAAADuEgmxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xd4VNXaxuHfOymkAEkIQekdRJBmqIoFLCgeQVCwg4CIgFiP5ajHfuz6HQQERUGwgWJBxYYISifSiyigIEUJvQRS1/fHDJwAAQaTzEyS576uubJnZs2eZwLJm73XXmuZcw4RERF/eIIdQEREig4VDRER8ZuKhoiI+E1FQ0RE/KaiISIiflPREBERv6loiIiI31Q0RETEbyoaIiLit/BgByho5cuXdzVq1Ah2DBGRIuWnn37a6pxLOlG7Ylc0atSoQUpKSrBjiIgUKWa2zp92Oj0lIiJ+U9EQERG/qWiIiIjfVDRERMRvKhoiIuI3FQ0REfGbioaIiPgtaEXDzKqa2fdmtsLMlpvZ7Xm0MTMbYmarzWyJmTUPRlYREfEK5uC+LOBu59wCMysD/GRm3zrnVuRqcwlQ13drBbzq+1ooeoycDcD4W9oU1luIiBRpQTvScM5tds4t8G3vAVYClY9o1hkY67zmAPFmVjHAUUVExCck+jTMrAbQDJh7xFOVgT9y3d/A0YUFM+tnZilmlpKamlpYMUVESrygFw0zKw1MBO5wzu3+O/twzr3mnEt2ziUnJZ1wvi0REfmbglo0zCwCb8F4xzn3UR5NNgJVc92v4ntMRESCIJhXTxnwBrDSOffSMZpNAm70XUXVGtjlnNscsJAiInKYYF49dRZwA7DUzBb5HvsXUA3AOTcCmAxcCqwG0oCbCjOQcw5vLRMRkbwErWg452YAx/0N7ZxzwMBA5Nm1P5OlG3dTMS4qEG8nIlIkBb0jPFSkZ2UTHmas3bqPez5YTFpGVrAjiYiEHBUNnwplomhwahkqx0cxccEGOg+dyS9/7Ql2LBGRkKKikYuZUSUhhnG9W7EjLZPLh85g/Pz1eM+SiYiIikYuES4DgLPrlmfy7WfTvFoC901cyp3jF7E3XaerRERUNA7av4MhW3px3e5RsH8nFcpEMa5PK+66sB6TFm/i8ldmsGLT3xp7KCJSbKhoHJSTzZJSzbls30QY0gzmjiTMZTG4Q13e6duavelZdBk+k3fmrtPpKhEpsVQ0Dootz6vx9/BA+Vfg1Ebw5b0wrBWs/Jw2tcox+fZ2tK6VyIMfL2PQewvZcyAz2IlFRAJOReMIv0fUgRsnwbUTwBMO46+D0ZdSftcyxvRqwb0d6/PVsj+57JUZLN2w67j76jFy9qHp1kVEigMVjVzG39LGu5aGGdS7GG6dBZ1egq2/wOvt8Xx8MwOaRvJ+v9ZkZOXQ7dVZjJn5m05XiUiJoaJxPGHh0KIPDF4I7e6GlZ/BK8m0+PX/mHxLE9rVLc+jn63g1rcXsGu/TleJSPGnouGPqLLQ4d9w20/Q8AqY+V8SRrVk1GkLefCSekxZ+RedhvzIoj92BjupiEihUtE4GXFVoOtI6DcdKpyOffVPbl5yNRM6OpyDK1+dxagf1+p0lYgUWyoaf0elptDzM7jmfcBoPvVaJpcfSvsakTz5xUpuHpvCzrSMYKcUESlwKhp/lxnUvwQGzIZOLxK3fTEjN3bjkWpLmL4qlUv/+6MuyxWRYkdFI7/CIqBFXxi8AGt3Fzdt/z8mRj5KeMZOVmzew1+7DwQ7oYhIgVHRKChRcXDBIzAohcaNm/J5zkDODVvC79vS+HTB+mCnExEpECoaBS2+KnR9jbK3fMmdMV/RylZwz4SFzJz6BaiDXESKOBWNwlKpGU8nPk2jChHUDN9G/2/SWDmyJ2xaGOxkIiJ/m4pGYTJjWUwrxtx5FbFRpei17mI2juwGH/WDXRuCnU5E5KSpaARApcQyjOl/HmmRifQMf45dy76BV86EKY/BAU23LiJFh4pGgJx2allG3pjM+gOx3FxuNAfqd4EZL3mnYZ8/CrK1yJOIhD4VjQBqW7s8L3RvwrwNadyV0Y+cvlMhqT58cTe82gZWfaXOchEJaSoaAXZ5k0o81KkBk5f+yRMLo3A9P4er3wWXA+/1gLf+AZsXBzumiEiewoMdoDgbf0ubPB/v264Wm3Ye4M2Zv1EpLpqbz+kEdS+Cn8bAtKdh5LnQ5Gpo/5B3visRkRChI40geahTAzo1rshTk1fy6aKN3pHlLW/2TsN+1u2w7CNvZ/l3T0D6nmDHFREBVDSCxuMxXryqCS1rluOeDxYza81W7xNRcXDhY3BbCjT4B/z4grezPOVNdZaLSNCpaARRVEQYr9+QTM3ysdwy9id+/jPX5bfx1aDbKLh5KiTWhc/vhFfbwi9fq7NcRIJGRSPI4mIiGHNTS2JLhdPrzfls2rn/8AaVz4SbJkOPdyAnC97tDmM7w+YlwQksIiWaikYIqBQfzZjeLdiXnkWv0fPYlXbElOpm0OAyGDgXLnke/lwKI8+BTwbA7k3BCS0iJVJQi4aZvWlmW8xs2TGeP8/MdpnZIt/t34HOGCjewX9n8tvWfdw8LoUDmdlHNwqLgFb9fJ3lg2HpBzCkOUx9Up3lIhIQwT7SGAN0PEGbH51zTX23xwOQKWja1i7Pi92bMu+37dz9wWJycg7vu+gxcjY9Rs6G6Hi48HEYlAKndYIfnvcWj5TR6iwXkUIV1KLhnPsB2B7MDKHm8iaVePDSBnyxZDNPTV55/MYJ1eHKN6DvVEisDZ/fASPOhl+/VWe5iBSKYB9p+KONmS02sy/NrGGwwwRC33Y1uemsGrwx4zdG/bj2xC+ocibc9CX0eBuy0+GdK2FcF2/fh4hIAQr1orEAqO6cawK8AnySVyMz62dmKWaWkpqaGtCAhcHMeLjT6Vx6xqk8+cVKJi32o7PbzDuuY8Bc6PisdyqSEe3gk4HqLBeRAhPSRcM5t9s5t9e3PRmIMLPyebR7zTmX7JxLTkpKCnjOwuDxGC91b0rLGuW4Z8JiZq/Z5t8LwyOhdX8YvAjaDoKlE7wjy7//D6TvLdzQIlLshXTRMLNTzcx82y3x5vXzt2fRFxURxus3JlM9MYZ+41JIyziJTu7oeLjoSRg0H+pfAtOf9Y4s/2kM5ORxZZaIiB+Cfcnte8BsoL6ZbTCzPmbW38z6+5pcCSwzs8XAEOBq50pWD29cTARjerckJjKMn//cQ3rWSf7CT6gBV74Jfb+DcrXgs9t9neVTgFxXZImI+MGK2+/g5ORkl5KSEuwYBW7l5t1cNmQGYWFGv3a1OP+0JJpWTSDMY/7vxDlYOQm+fQR2/Aa12/PPXVeyPqLWMWfkFZGSwcx+cs4ln7CdikbRcel/f+CP7ftJy8wmO8cRFx1Bu7rlOa9+Bc6tl0RSmVL+7SgrA1LegOnPkrN/J9OiL6T9rUOgbMXC/QAiErL8LRpaT6MIKRMVwemVInjtxmRm/LqVaau2MO2XVD5fshmARpXLcn79CpxX/wRHIeGR0PpWaHI1Xwy9m477JsErzaHtYGh7G5QqHcBPJSJFiYpGERQXHUGnxhXp1LgizjmWb9rN9F9SmbZqC8OnreGVqav9OwqJTuDtsjfzTcxlvFJhEkx/xttR3v5BaHodeMIC/tlEJLSpaBRxZkajynE0qhzHwPPrsCstkxmrt/L9qi1Mz3UUckblOM6rn5TnUciW8Ipw1RhoPQC+fhAm3QZzXoWLnoA6F5xUnoOd6uojESmeVDSKmbiY/x2F5OQ4Vmz+31HIsO9XH3UUkpmdQ0SY7yK6qi2hzzew4lOY8gi83Q1qd/AWj1NKxGB8ETkBFY1izOM5+ijkx9WpTFuVethRSFSEh7snLKZ59XiaV0ugXoPOhNW/BOaPgunPeS/RbXodnP+gOstFSjgVjRIkLiaCyxpX4rLGlQ4dhdz69k/sSc9i2qotTFywAYDYyDCaVI2nebULaN6xI802jCVhwTBYNtG7fnnb2yAyNiCZdbpLJLSoaJRQB49CKsVHA/B+v9as357GwvU7WbB+BwvW7+DV6WvIznFAS2omnE0z+4XmU76j+ZzLqH9Bb8KaX6vOcpESRkWjCCnMv7bNjOqJsVRPjKVLs8oApGVksWTDrkOF5If19fgoqwbshJgPD9DksxE0P602zZs2o1m1BMrFRhZaPhEJDSoackwxkeG0rpVI61qJADjn+GP7fhas287CxQtYsCacEYszyF7sHUxZIzGGvelZxESG8/2qLdRJKk3l+Gg8JzNqXURCmoqG+M3MqJYYQ7XEGLo0rwJZ6eyfNYolP0xk4YFKLHDtmZZWga17M7hp9HzA28leq3xp6lQoTe0k39cKsdRIjCUqQqe2RIoaFQ35+8JLEX3OQFolX0OrH56HefdyoFQY70VfS6Mr7mH1jmzWbNnL6tS9LPxjB58t2XRoQUGPQdVyMdRJKk3tCqUP+xoXExHczyUix6SiIfkXUw46Pg0t+rLo9cHclPYmTPqSFu0fgkuvOdRZvj8jm7Vb97ImdR+rt+xlTepe1mzZy4+rt5KRlXNod+VLl6J2Uix1KpTmz10HiItWEREJFZqwUApUj5GzqZ+xnMej3oONKXBKI+/gwNrtj/ma7BzHhh1phwqJ96u3sOzan4nHYNb9HTg1LiqAn0SkZNGEhRI0qyIbQt8psPwjmPIojLsC6lzoLR4VGhzVPszzvyu3OjQ45dDjzjkuHzqDZZt288QXKxh2bfMAfgoRyUtIr9wnRZgZNOoGg1K8Kwj+MQ9ebetdBGrPX37uwoiJDKdyXDRfLNnM9F+K/vrvIkWdioYUrvBS3hHkty+ClrfAwre9y85Ofw4y0vzaRcX4KGqVj+Xfny7jQKaWqhUJJhUNCYyYcnDJMzBwHtRpD98/5V3DY+E7J1yz3GPGk10asW5bGsOnrQlQYBHJi4qGFKjxt7Q5/sj1xNrQ42246SsoWwk+HQCvnQtrpx13v23rlKdz00qMmLaGtal7Cza0iPhNRUOCo3ob6DMFur0B+3fB2M7wzlWwZeUxX/JgpwaUivDw8KfLKG5X/YkUFSoaEjweD5xxJQyaDxc+Duvn+jrL74C9W45qXqFMFPdeXJ+Zq7cxafGmIAQWEY3TkNCxbxtMfxZS3oDwKDj7Dmg9ECJjDjXJznFcMXwmm3Ye4Lu7z9XAP5EC4u84DR1pSOiITYRLn4MBc6HWeTD1SXjlTFj0LuR4R4yHeYynupzB9n3pvPjNqqDGFSmJVDQk9JSvA1e/A70mQ5lT4JNb4bVzYO10AM6oEseNbWowbs46lmzYGeSwIiWLioaErhpnQd+p0HUU7N8JYy+Hd3tA6iruuqge5UuX4sGPl/kWihKRQFDRkNDm8UDjq7wjyy94DNbNguFtKDvlXh6+oDJLN+7i7Tnrgp1SpMRQ0ZCiIcLXMT54IbToAwvG8o8pF9Ku/F5e+Ppntuw+EOyEIiWCioYULbHl4dLnYcAcrNY5PL7736SnH+CJt7861FkuIoVHRUOKpvJ14Zp3qdn7DW6Nn8dn6yP4cUgv+O2HYCcTKdZUNKRoq3E2t975GDVK5/Dv1PYcGNMV3r0aUn8JdjKRYimoRcPM3jSzLWa27BjPm5kNMbPVZrbEzLSgghwlqlQET/RozW/ZSYyo/iKsmwnDW8MXd8NeTacuUpCCfaQxBuh4nOcvAer6bv2AVwOQSYqgdnWTuKxxRYavrcDv18+G5N6QMto7DfuPL0Hm/mBHFCkWglo0nHM/ANuP06QzMNZ5zQHizaxiYNJJUfPwZacTGebh4W824i59HgbOhZrnwHePwSvJsHi8OstF8inYRxonUhn4I9f9Db7HDmNm/cwsxcxSUlN1OqKkOqVsFPdcVI8ff93K50s2H+osp+fn3quuPu4Hr58Pv88IdlSRIivUi4ZfnHOvOeeSnXPJSUlJwY4jQXRDmxo0qlyWJz5fwZ4Dmd4Ha7aDm7+Hrq/Dvq0wphO8dw1s/bVA3rPHyNn0GDm7QPZVGPsTKUihXjQ2AlVz3a/ie0wkTwcnNEzdm86L3+S6gsrjgcbd4bYU6PAI/PYjDGsFX9zjLSQi4pdQLxqTgBt9V1G1BnY55zYHO5SEtiZV47m+VXXGzv6dZRt3Hf5kRDS0u8s7sjz5Jkh509tZPuNlyNSocpETCfYlt+8Bs4H6ZrbBzPqYWX8z6+9rMhlYC6wGXgcGBCmqFDH3XFyfcrGlePDjpXlPaFg6CTq9CAPmQI2zYcqjMDQZlnygznKR4wj21VPXOOcqOucinHNVnHNvOOdGOOdG+J53zrmBzrnazrkznHNaXUn8EhcdwUOdGrB4wy7enbf+2A2T6sE170HPzyCmHHzUF0a1h99nBi6sSBFyUkXDzGLNLKywwogUpM5NK9G2diLPffUzW/Z4Tz0ds5O55jlw8zS4YqR3qdkxl8L718HW1YENLRLijls0zMxjZtea2RdmtgX4GdhsZivM7HkzqxOYmCInz8x4oksj0jNz+M8XK0/8Ao8HmlwNt/0E7R+GtdNgeCuY/E/vUrSFJCfHkZmdw4HMbPalZ5GVk0NxW4ZZio/wEzz/PTAFeABY5pzLATCzcsD5wLNm9rFz7u3CjSny99ROKs0t59bilamr6Z5c9cQvAG9n+Tn3QPMbYdozMP8NWPw+tLsbWvWHiCg279rPgnU7WbB+Bys27SbbOS5++QeynSM754ibc+TkOLJyvF+z3eHbedWHslHhZOc4wjxWsN8QkXyy4/1FY2YRzrnM4+7AjzaBlJyc7FJS1PUh/3MgM5uLXv6B8DAjMTYSjxnjb2nj/+s3r2T558NYuG4bC8Ias8BzBn/u9/4yLxXuISLMQ0SY0apmImEeO3TzmBHuMTweI8wD4R4PHvNuezze58LMDm17fPdHz/yNP3enc2/H+gw4TwfzEhhm9pNzLvlE7Y57pHGwGJjZOOfcDUe8wTjn3A2hVDBE8hIVEcbjnRvSa/R8MjJzqJwQfcy2zjk27TrAgnU7WLj+f0cSGdmXAFCFnbTMnEWzxAM0P68LDZqdzQ1vzAVgxA1nFkje71b+RUa24+Vvf+Hcekk0rBRXIPsVKQgnOj11UMPcd3yd4QXzEyISAOfVr8ClZ5zKl8v+JLF05KHHD2Rms2zjLhas38GCdTtZ+McO/tqdDkBUhIfGleO56ewaNKuaQPPq8VSIjYQl4+G7x+GLN2DNZVTM6sLm8CoFltXMqJEYw4Yd+7lz/CImDTqbqAhdfyKh4bhFw8weAP4FRJvZ7oMPAxnAa4WcTaRA/fuyhny17E/Wpu7j0UnLWbh+Bys27yYz23uKtmq5aFrXSqR5tQSaV0vgtIpliAjL41qRptfA6Z1hzjCY8X+8kPEl38RcBvvqQWxigWSNCPPw7JWNuWn0fF78ZhUPdjq9QPYrkl8nOj31NPC0mT3tnHsgQJlECsWpcVFUTYhh3fY0xs//g8ZV4uhzdi2aV4unWbUEksqU8n9nkTFwzj+heU++H347HdM+gyHTvB3oLft51zTPp/PrV+D61tUYNeM32p92Cm1qF0xBEsmPEx1p1HDO/X6sgmFmBlR2zm0olHQiBeyUsqWIj4ngw1vb5n0UcbJKV2BU3GC+jOnMSwkT4duHYf7r3vmtGnUDy9/VT/+6tAEzV2/jng8W8+Ud7SgbFZH/zCL5cKKfmufNbKKZ3WhmDc2sgplVM7P2ZvYEMBNoEICcIgXCzIiKCCuYgpHLxojqcN0HcMMnUCoOJvaBUR1gXf5mq42JDOfF7k3YvGs/j01aUUBpRf6+4/7kOOeuAh4G6gPDgB+BT4G+wCqgvXPu28IOKVJk1D4fbpkOXV6F3ZthdEcYfz1sW+P3Lsbf0uawS4KbV0tg0Pl1mLhgA18t03ydElwnOj0V4ZxbATwYoDwiRc5RYz48YdD0Wji9C8we5p1Bd9WX0OJmOPde7xxXJ+m2DnX5flUqD3y0lObVE6hQJv99JiJ/x4mO0Tea2Sgz6+DrvxARf0XGwLn/9E7D3ux6mDcShjSFWa9AVvpJ7SoizMPLPZqQlpHN/ROXapoRCZoTFY0GwHzgIeAPM/uvb10LEfFXmVPgH/+FW2dBlZbwzUMwtAUsm0iec4gcQ50KZbiv42lM/XkL78//48QvECkEJ+rT2OacG+mcOx9oiXdti5fNbI2ZPRWQhCLFRYUGcP2HcMPHUKoMfNgbRl0A6+f4vYtebWtwVp1Envh8Beu27SvEsCJ5O+7cU0c1NisNdAXuAio6504prGB/l+aekiIhJxsWvwdTn4Q9m6HB5XDBo5BY+4Qv3bRzPxf/3w/UO6UME25po0kNpUD4O/fUCa87NLMoM7vKzD7Cu4Jee+B+oFL+Y4qUUJ4wbz/HbT/Bef+C1d951yz/6gFI237cl1aKj+aJzo34ad0ORkz3/6oskYJwovU03gXWA92Bd4AazrlezrmvnHPZgQgoUqxFxsJ598HgBd4rruaO8KuzvHPTSnRqXJGXv/3l6HXQRQrRiY40vgJqO+eucs5NdM4dCEQokRKnzKlw+RDoPwMqJ+fqLP8oz85yM+OpLo0oFxvJneMXcSBTf8NJYJyoI3ysc25PoMKIlHinNIQbPoLrP4LI0vDhTfDGhbB+7lFN42Miee7Kxvy6ZS8vfL0qCGGlJCrYuRREpGDU6QD9f4TLh8LOP+DNi2BCT9i+9rBm59WvwA2tqzNqxm/MWrM1SGGlJFHREAlVnjBofoOvs/wB+PUbGNoSvvrXYZ3lD1x6GjXLx3LPhMXsPqA10aRw+V00zKytmV3rm7zwRjO7sTCDiYhPqdJw3v1w2wJocjXMGQ5DmnmnKMlKJyYynJe6N+GvPek8Oml5sNNKMedX0TCzccALwNlAC9/thNfzikgBKlsROg/1dZY3h6//BcNawvJPaFY1noHn1eajBRv5cqkmNZTC49fgPjNbCZzuisCENxrcJyXG6inwzcOwZQVUaUnmBU/S9bNMNuxI4+s7zqFCWe+khj1GeqdnP2piRZFcCmxwn88y4NT8RRKRAlXnAu9Rxz+GwM51RIy5iJfjx5OWkcV9E5doUkMpFP4WjfLACjP72swmHbwVZjAR8YMnDM7s6e3vOPd+6qz/kPttHN+vSuW9mT8HO50UQ8ddTyOXRwszhIjkU6nScP4DcGZPek59iu/mLeWJz9Npu/9HwtxpZJuWiZWC4VfRcM5NL+wgIlIAylbC02UYzzdcyMWjf+OuqfsYEt2f8WVvAtc632uWi5xo7qk9ZrY7j9seM9ud3zc3s45mtsrMVpvZ/Xk838vMUs1ske/WN7/vKVISVKzbjCeuaskCV4+3sjpw984n4c2OsEEXiUj+HPdIwzlXprDe2MzC8K47fiGwAZhvZpN8y8vmNt45N6iwcogUV5c3rcw3K7cwaklHsspV5pHtb8CoDtCwK1zwCCTUCHZEKYKCOSK8JbDaObfWOZcBvA90DmIekWLl4KSG4WHGe3uakNp7Dpxzr3e98qEtvJMi7t8R7JhSxASzaFQGcq9ZucH32JG6mdkSM/vQzKoGJppI8RAfE0ntpNJkZOXQ9Y3FrG50u3ca9jOugllDvSPL54yArIxgR5UiItTnnvoM7xoejYFvgbfyamRm/cwsxcxSUlNTAxpQJNTFRUfQoGJZ9mdk0+3VWczdWgq6DIdbfoBTG8NX98HwVrDys5Nas1xKpmAWjY1A7iOHKr7HDvGtUX5wJZpRwJl57cg595pzLtk5l5yUlFQoYUWKstKlwvl4wFmULx3JDW/M49NFG6FiY7jxU7j2AwiLhPHXw+hLYMNPwY4rISyYRWM+UNfMappZJHA1cNiAQTOrmOvu5cDKAOYTKVaqlovho1vPolm1eG5/fxHDvl+NA6h3EfSfCZe9DNtWw6j28GFv2LEu2JElBAWtaDjnsoBBwNd4i8EE59xyM3vczC73NRtsZsvNbDEwGOgVnLQixUNcTARj+7Skc9NKPP/1Kv718VKysnMgLBySe8PghdDuHvj5Cxia7J3bav/OYMeWEOLXhIVFiSYsFDlcXhMWOud48ZtfGPr9as6tl8Sw65pTulSuK/B3bYSpT8Di9yE6wTs1e3JvCNPI8uKqoCcsFJFixMy45+L6PNP1DGas3kr3EbP5c9eB/zWIqwxXjIB+07xL0H55LwxrBSs/V2d5CaeiIVKCXd2yGm/2asG6bfu4YvhMfv7ziIkeKjWFnp/BtRPAEw7jr4MxnWCjOstLKp2eEhFWbNpN7zHz2ZuexavXN6dd3TyuQszOggVvwff/gbSt3rEeHf4N8dUCH1gKnE5PiYjfTq9Ulo8HtqVKQjQ3jZ7PhJQ/jm4UFg4t+ng7y8++yzuu45Vk+PYROLAr8KElKFQ0RASAinHRfNC/DW1qJ3Lvh0t46ZtVeS/kFFWWHmsuYkC516HhFTDz/7wjy+e+BtmZgQ8uAaWiISKHlImK4M1eLeieXIUhU1dz94TFZGTl5Nl2W1gSdB0J/aZDhdPhy3/C8Nbey3VP8rR3j5GzD13lJaFNRUNEDhMR5uHZbo25+8J6fLRwIz3fnMeu/cc5gjjYWX7N+4DB+9fCmMtg44KAZZbAUdEQkaOYGbd1qMvLPZqQsm47V746iw070o73Aqh/CQyYDZ1ehNSf4fXzYeLNsDOP/hEpslQ0ROSYrmhWhbG9W/HX7gNcMXwWSzecoMM7LAJa9M3VWT4JXjkTpjwKB/K9bpuEABUNETmuNrUTmXhrWyLDPHQfOZvvVv514hdFlfUu9DQoBRp2gRkvezvL572uzvIiTkVDRE6o7ill+HhgW+pUKM3NY1P4c/cBsrJzyMzOu5P8kPiq0PU178jypNNg8j0wvA38PFkjy4uo4y73KiJyUIUyUYy/pTWD31vIlJVbWLctjboPfklkmIfoyDBiI8OIjgwjJjKcmMgw3+3gdgTRp75AbMxaotd+RczbrxOb9C3Rza8i5tR67DmQRUxkWLA/ovhBRUNE/BYTGc7IG5K5+OXpZGQ7rjqzCmmZ2aSlZ5GWke27ebe37s0gLSON/RnZ7MvIZn9GNhnZAB29O9sMfLEDmAtAqXAPyzftomGluCB9OvGHioaInJSPy/hmAAAQ/0lEQVQwj5FYuhQAt3Woe1KvzczOIc1XQPbt2cH+eeNIW/IpG7ITeCyrF1cMm8njnRvRo0VVzKww4ks+qU9DRAImIsxDXHQEp8ZFUbtKRRp1vZeWd7xHUmw430XeSUtbyf0fLeWeCYvY7z0skRCjIw0RCa74agxNuI9asb/wVuwEhqxZzpCFXVn22yaG9z6H2hXKBDuh5KIjDREJCWsj6xF20+fceUM33kocS+rO3Vz+f1P4bLqmFwklKhoiEjrM4LROnHPXO3zRMY3TPBu47cvtPPLyUNK3rQ92OkFFQ0RCUVgEFc/rw/sPXE/fapt566+adH/xUzZ8/iyk7wl2uhJNfRoictJyrzdemCJiE3hoQF+S5y7ln5My6TQjnZcXX0f7i7tAsxu9a3xIQGnlPhEpEtZt28eto2ewYmsWt4Z9yt0VlxB+0eNQ90Lvaa1cDk6zXlDFraD3F4q0cp+IFCvVE2P56PYLuKZFVV7N7sx1W65nyzs3w9jOsHlJsOOVGCoaIlJkREWE8XS3xrzUvQlLXG0u5RVm/3EARp4DnwyA3ZuCHbHYU9EQkSKna/MqfDrobMqWLct1+25nWKWnyVnyIQxpDlOfIirnOGt/SL6oaIhIkVTvlDJMGnQ2nRpX4vm11ehz6gfsqNMFfniO/6b2oUPaZMjOCnbMYkdFQ0SKrNKlwhlydVOe6NyQmevSuOz37iz6x1dsDqtEv11DYMTZ8Ou3moa9AKloiEiRZmbc0KYGH97aBjO46qMd9I94khfiHoLsdHjnShjXBf5cGuyoxYIuchaRYqFxlXi+uK0dd01YxHc/b2F05Bksqf4qnuwNeH5djefnj7G4eXiS6uOJjMZjhhl4zPD4vlqubY+HQ/d/37qPMI+xMy2D+JjIYH/UoNI4DREpVnJyHOe98D3b92VSp0JpnHPkZGeRs+cvctK24/CQE5NITnQiORjOQY5z3lsO3vaHHvPe37U/k6wcR43EGEb1TKZOMZxE0d9xGjrSEJFixeMxKsZFUzEu+ujBeDt+hymPwfKPIKcCtH8Qmt0AnuOvGthj5Gz2HMhky54Mrhg2iyHXNOP80yoU3ocIYUHt0zCzjma2ysxWm9n9eTxfyszG+56fa2Y1Ap9SRIqNhBpw1WjoMwXK1YTPbvd1lk85YWd5magIJg06i2qJMfR+az6v/bCG4namxh9BKxpmFgYMAy4BTgeuMbPTj2jWB9jhnKsDvAw8G9iUIlIsVW0Bvb+G7mMhcz+80w3GXXHCzvJK8dF80L8NlzaqyH8m/8zdHyzmQGbJWiwqmEcaLYHVzrm1zrkM4H2g8xFtOgNv+bY/BDqY1oAUkYJgBqd3hoHz4OKnYdNCGNEOPh0Iuzcf82UxkeEMvbYZd15Qj48WbOSa1+ewZc+BAAYPrmAWjcrAH7nub/A9lmcb51wWsAtIPHJHZtbPzFLMLCU1NbWQ4opIsRQeCW0GwO2LoM1AWDweXmkO3/8H0vfm+RIz4/YL6vLqdc35efMeOg+dybKNuwIcPDiKxTgN59xrzrlk51xyUlJSsOOISFEUnQAXPwWD5kO9i2H6s97i8dNbmMv7FNQlZ1Tkw1vb4DHjyhGz+HxJ8Z/7KphFYyNQNdf9Kr7H8mxjZuFAHLAtIOlEpGQqVxOuGgN9voX46vDZYJ7bOpAmB/K+lL9hpTg+HXQWjSrFMejdhbz0zSpycopvB3nQxmn4isAvQAe8xWE+cK1zbnmuNgOBM5xz/c3saqCrc6778farcRoiUmCcgxWfwpRHvJfr1m4PFz0JpzQ8qml6VjYPf7KMCSkbuLjhKbzUvSmxpYrOqIaQX0/D10cxCPgaWAlMcM4tN7PHzexyX7M3gEQzWw3cBRx1Wa6ISKExg4ZdfJ3l/4GNC7yX6H466KjO8lLhYTzbrTEPX3Y63674i26vzmLDjmPPtttj5OxDizsVJRoRLiLir7Tt8MMLMO81CIuAtoPhrMEQGXtYs+m/pDLo3QVEhnkYccOZtKhR7qhdhdpqgCF/pCEiUuTElIOO/4FB87zLzE5/xruGx4KxkPO/zvJz6yXxycCziIuO4NrX5zB+/voghi5YKhoiIierXC3vwMDe30B8VZh0m3eMx+rvDjWpnVSajwecRetaidw3cSmPfbacrOycIIYuGCoaIiJ/V7VW3qusrhwNGXvh7a7wdjf4awUAcTERjO7Vgt5n1WT0zN+5acx8dqVlBjl0/qhoiIjkhxk06uod33HRU7BhPow4y3v0sedPwsM8/Psfp/Nct8bMWbuNK4bPZE1q3oMGiwIVDRGRghBeCtoOgsGLoFV/WPSet79j2rOQsY/uLary7s2t2bU/ky7DZrIzLSPYif8WFQ0RkYIUUw46Pg0D50KdDjDtP/DKmbDwbVpU8w4ErJIQw6q/9pK6Jz3YaU+aioaISGFIrA09xnln0y1b2TsR4shzqbJ9Lh/2b0NcdDhrt+7j6+V/BjvpSVHREBEpTNVaQ98pcOWbkL4LxnUh9oOrOS9hO7GlwrjtvYXMXVt0ZkdS0RARKWxm0KgbDJwPFz4Bf8zj5e0DeCr6XarGRdJ3bAorN+8Odkq/qGiIiARKRJR3BPnti/gq5h9cduBzxqbfQazbT8835/LH9mNPOxIqVDRERAItphxvxd3K3UmvUbluE8a6f5G+dyc3Dp/C1t2hXThUNEREguTP8MrQ423q9RnFmxU/YfPeLG564V32/vx9sKMdk4qGiEiwVW/LmYPGMfw8WJFRgf5j55I+rgds+TnYyY6ioiEiEgo8Htp37MqzXc9gRs4Z3P3L6eQMbwuf3QF7twQ73SFFZ4UQEZFi5FhTol/Zsibb9ufw9JeQmPQMjy54AFv6AZx9B7QeCJExAU56OB1piIiEmFvOrc3N7Wry1qbKDG0+GWqdB1OfhKHJ3ulJcoI3W66KhohICHrgkgZ0bVaZF2du5d2aT0OvyVC6AnzSH147B9ZOD0ouFQ0RkRDk8RjPXtmY8+sn8dAnS/lqb23oOxW6joL9O2Hs5fBuD0hdFdhcAX03ERHxW0SYh2HXNadJ1XgGv7+QOb/vgMZXwaAUuOAxWDcLhreBz++EvakByaSiISISwmIiw3mzZwuqlYvh5rdSWLFpt3dk+dl3wOCF0KKPd7nZIc2865c7V6h5VDREREJcQmwkY3u3pHRUOD1Hz2P9Nt+o8djycOnzMGAO860R82Z+553nqhCpaIiIFAGV4qMZ27slmdk53PjmXLbuzbUWR/m6vFDuEf6bcH+h51DREBEpIuqeUoY3erbgz90H6DV6HnsOHL7eeJZFFnoGFQ0RkSLkzOoJvHrdmazcvIdbxv1EelZ2QN9fRUNEpIg5/7QKPH9lY2at2cZd4xeTnVO4nd+5aRoREZEiqGvzKmzbm8FTk1dSLjYS5xxWyJ3goKIhIlJk3XxOLbbuTWfkD2upHB9NlYToQn9PFQ0RkSLs/ktOY+veDCYu2EBEWOEfaahPQ0SkCDMznul2BvHREWzbl1Ho/RtBOdIws3LAeKAG8DvQ3Tm3I4922cBS3931zrnLA5VRRKSoiAjzUKdCaQwI8xTPwX33A9855+oC3/nu52W/c66p76aCISJyDGEew1PIBQOCVzQ6A2/5tt8CugQph4iInIRgFY1TnHObfdt/Aqcco12UmaWY2RwzU2EREQmyQuvTMLMpwKl5PPVg7jvOOWdmx+q5qe6c22hmtYCpZrbUObcmj/fqB/QDqFatWj6Ti4jIsRRa0XDOXXCs58zsLzOr6JzbbGYVgTxXTXfObfR9XWtm04BmwFFFwzn3GvAaQHJycuCGRoqIlDDBOj01Cejp2+4JfHpkAzNLMLNSvu3ywFnAioAlFBGRowRrcN8zwAQz6wOsA7oDmFky0N851xdoAIw0sxy8xe0Z55yKhohIHsbf0iYg7xOUouGc2wZ0yOPxFKCvb3sWcEaAo4mIyHFoRLiIiPhNRUNERPymoiEiIn5T0RAREb+paIiIiN9UNERExG8qGiIi4jcVDRER8Zs5V7ymajKzVLyjzP+u8sDWAopTGEI9H4R+xlDPB8pYEEI9H4RWxurOuaQTNSp2RSO/zCzFOZcc7BzHEur5IPQzhno+UMaCEOr5oGhkPJJOT4mIiN9UNERExG8qGkd7LdgBTiDU80HoZwz1fKCMBSHU80HRyHgY9WmIiIjfdKQhIiJ+KzFFw8w6mtkqM1ttZvfn8XwpMxvve36umdXI9dwDvsdXmdnFoZbRzC40s5/MbKnva/tQypfr+WpmttfM7imMfPnNaGaNzWy2mS33fS+jQimjmUWY2Vu+bCvN7IEg5TvHzBaYWZaZXXnEcz3N7FffreeRrw12RjNrmuvfeImZ9QilfLmeL2tmG8xsaGHkyxfnXLG/AWF41xavBUQCi4HTj2gzABjh274aGO/bPt3XvhRQ07efsBDL2Ayo5NtuBGwMpXy5nv8Q+AC4JwT/ncOBJUAT3/3EEPx3vhZ437cdA/wO1AhCvhpAY2AscGWux8sBa31fE3zbCUH6Hh4rYz2grm+7ErAZiA+VfLme/y/wLjC0oL9/+b2VlCONlsBq59xa51wG8D7Q+Yg2nYG3fNsfAh3MzHyPv++cS3fO/Qas9u0vZDI65xY65zb5Hl8ORB9cXz0U8gGYWRfgN1++wpKfjBcBS5xzi8G7uqRzLjvEMjog1szCgWggA9gd6HzOud+dc0uAnCNeezHwrXNuu3NuB/At0LGA8+Uro3PuF+fcr77tTcAW4IQD2gKVD8DMzgROAb4p4FwFoqQUjcrAH7nub/A9lmcb51wWsAvvX5v+vDbYGXPrBixwzqWHSj4zKw3cBzxWwJkKLCPev0CdmX3tO21wbwhm/BDYh/ev4/XAC8657UHIVxivPRkF8j5m1hLvkcCaAsp10N/OZ2Ye4EWg0E7h5ldQ1giXwmFmDYFn8f7VHEoeBV52zu31HXiEonDgbKAFkAZ8Z2Y/Oee+C26sw7QEsvGeVkkAfjSzKc65tcGNVfSYWUVgHNDTOXfUX/tBNACY7JzbEKo/KyXlSGMjUDXX/Sq+x/Js4zv8jwO2+fnaYGfEzKoAHwM3OucK+i+n/OZrBTxnZr8DdwD/MrNBIZZxA/CDc26rcy4NmAw0D7GM1wJfOecynXNbgJlAQU9BkZ//76H0s3JMZlYW+AJ40Dk3p4CzQf7ytQEG+X5WXgBuNLNnCjZePgW7UyUQN7x/Ra7F25F9sGOq4RFtBnJ45+ME33ZDDu8IX0vhdJDmJ2O8r33XUPweHtHmUQqvIzw/38MEYAHeDuZwYArQKcQy3geM9m3HAiuAxoHOl6vtGI7uCP/N971M8G2XC8b38DgZI4HvgDsK4/9gfvMd8VwvQrAjPOgBAvZB4VLgF7znLx/0PfY4cLlvOwrvlT2rgXlArVyvfdD3ulXAJaGWEXgI77nuRbluFUIl3xH7eJRCKhoF8O98Pd6O+mXAc6GWESjte3w53oLxzyDla4H3yGwf3iOg5ble29uXezVwUxC/h3lm9P0bZx7xs9I0VPIdsY9ehGDR0IhwERHxW0np0xARkQKgoiEiIn5T0RAREb+paIiIiN9UNERExG8qGiIBYGbxZjYg2DlE8ktFQyQw4vFOESFSpKloiATGM0BtM1tkZs8HO4zI36XBfSIB4FtI6XPnXKMgRxHJFx1piIiI31Q0RETEbyoaIoGxBygT7BAi+aWiIRIAzrltwEwzW6aOcCnK1BEuIiJ+05GGiIj4TUVDRET8pqIhIiJ+U9EQERG/qWiIiIjfVDRERMRvKhoiIuI3FQ0REfHb/wPTAConlC38RAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f1cc2ae9f50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "class LinearRegressionReconstruction(Reconstruction):\n",
    "    \n",
    "    def fit_linear_regression(self, meas):\n",
    "        # fit: y = ln V(t) = a t + b\n",
    "        # where a = -(1+true_clock_drift)/(RC), b = ln V(0)\n",
    "        \n",
    "        # use (t=0, V=V0) similarly to all other voltage measurements\n",
    "        x = np.array([0] + list(meas.t))\n",
    "        y = np.log([meas.V0] + list(meas.V))\n",
    "        \n",
    "        # solve least squares problem\n",
    "        M = np.ones((len(x), 2))\n",
    "        M[:,0] = x\n",
    "        inv_MTM = np.linalg.inv(np.dot(M.T, M))\n",
    "        \n",
    "        # linear regression: coeffs = [a, b]\n",
    "        coeffs = np.dot(inv_MTM, np.dot(M.T, y))\n",
    "        # cov = covariance matrix of [a, b]\n",
    "        # here the error level in the measurements is assumed to be known\n",
    "        # and not estimated from the data \n",
    "        cov = inv_MTM * self.log_v_error**2\n",
    "        \n",
    "        return coeffs, cov, x, y\n",
    "    \n",
    "    def compute(self, meas):\n",
    "        coeffs, cov = self.fit_linear_regression(meas)[:2]\n",
    "        \n",
    "        # R = -(1+true_clock_drift) / (a * C)\n",
    "        inv_rc = -coeffs[0]\n",
    "        # coeffs[1] would be a better approximation for ln V(0)\n",
    "        \n",
    "        # confidence intervals for a from the least squares fit\n",
    "        inv_rc_err = np.sqrt(cov[0,0])\n",
    "        \n",
    "        # final confidence intervals using interval arithmetic with known confidence intervals\n",
    "        # for C error and clock drift\n",
    "        r_min = 1.0 / (inv_rc + inv_rc_err) / (meas.C + self.c_error) * (1 - self.clock_drift)\n",
    "        r_max = 1.0 / (inv_rc - inv_rc_err) / (meas.C - self.c_error) * (1 + self.clock_drift)\n",
    "        r_mid = 1.0 / inv_rc / meas.C\n",
    "        \n",
    "        r_err = max(r_max - r_mid, r_mid - r_min)\n",
    "        return r_mid, r_err\n",
    "\n",
    "rec = LinearRegressionReconstruction()\n",
    "\n",
    "coeffs, cov, x, y = rec.fit_linear_regression(meas)\n",
    "plt.errorbar(x, y, rec.log_v_error)\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('ln V(t)')\n",
    "plt.plot(x, coeffs[0]*x + coeffs[1])\n",
    "\n",
    "r_rec, r_rec_err = rec.compute(meas)\n",
    "print(\"reconstruction: R = %g +- %.2g (actual: %.2g)\" % (r_rec, r_rec_err, actual.R))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare RMSE of the two methods and true vs. estimated errors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEWCAYAAABv+EDhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAF39JREFUeJzt3Xu0Z2V93/H3RxTS4AUMIyGADqWjKbYGzBRovCyyTLjaoEmlEKNgbEcTsJriMqhtQCldo9ZrYoiAs4AWQYxaJ0AExBhv4TIgAgMaRxwKU4RRFFBTGvDbP/Yz8PNwzpzrnDPnPO/XWr81ez/72Xs/z++cOZ/fvvyenapCktSfJyx0AyRJC8MAkKROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgAUlemeSKbbTtc5P8122xbWk2DAB1JckLk3w1yf1J7kvylST/qqouqKpDF7p90nx64kI3QJovSZ4KXAL8AXAxsCPwIuChhWyXtFA8AlBPng1QVRdW1SNV9Q9VdUVV3ZTkhCRf3lIxSSX5wyTfSvJgktOT7NuOHh5IcnGSHVvdQ5LcleRtSb6XZGOSV07UiCQvTXJjkh+27T1v23ddejwDQD35e+CRJOclOSLJrpPUPwz4VeBg4C3AWcDvAXsD/wI4bqTuLwK7AXsCxwNnJXnO2A0mOQBYA7wO+AXgI8DaJDvNpmPSTBgA6kZVPQC8ECjgbGBzkrVJdp9glXdX1QNVtR64Bbiiqm6vqvuBvwYOGFP/v1TVQ1X1t8ClwDHjbHMV8JGquqYdhZzHcArq4Nn3UJoeA0BdqarbquqEqtqL4VP8LwEfmKD6PSPT/zDO/JNH5n9QVT8emb+jbXusZwEnt9M/P0zyQ4YjivHqStuUAaBuVdU3gHMZgmC2dk2y88j8M4H/M069O4EzqmqXkdfPV9WFc9AGaVoMAHUjyS8nOTnJXm1+b4bz+FfP0S7ekWTHJC8CXgp8Ypw6ZwOvT3JQBjsnOSrJU+aoDdKUGQDqyYPAQcA1SX7M8If/FuDkOdj2d4EfMHzqvwB4fTvC+BlVtQ74D8CftfobgBPmYP/StMUHwkizk+QQ4H+26wrSouERgCR1ygCQpE55CkiSOuURgCR1arseDG633Xar5cuXL3QzJGlRuf76679XVcsmq7ddB8Dy5ctZt27dQjdDkhaVJHdMpZ6ngCSpUwaAJHXKAJCkThkAktQpA0CSOmUASFKnDABJ6pQBIEmdmjQAkuyd5G+S3JpkfZI3tvLTkmxKcmN7HTmyzluTbEjyzSSHjZQf3so2JDll23RJkjQVU/km8MPAyVV1Q3tq0fVJrmzL3l9V/320cpL9gGOB5zI85/RzSZ7dFn8Y+E3gLuC6JGur6ta56IgkzZXlp1z66PTG1UctYEu2rUkDoKruBu5u0w8muQ3YcyurHA1cVFUPAd9JsgE4sC3bUFW3AyS5qNU1ACRpAUzrGkCS5cABwDWt6KQkNyVZk2TXVrYnw4Ovt7irlU1UPnYfq5KsS7Ju8+bN02meJGkaphwASZ4MfBJ4U1U9AJwJ7Avsz3CE8N65aFBVnVVVK6tq5bJlkw5mJ0maoSmNBprkSQx//C+oqk8BVNU9I8vPBi5ps5uAvUdW36uVsZVySdI8mzQAkgT4KHBbVb1vpHyPdn0A4OXALW16LfCxJO9juAi8ArgWCLAiyT4Mf/iPBX53rjoiSbMxeuG3F1M5AngB8Crg5iQ3trK3Accl2R8oYCPwOoCqWp/kYoaLuw8DJ1bVIwBJTgIuB3YA1lTV+jnsiyRpGqZyF9CXGT69j3XZVtY5AzhjnPLLtraeJGn++E1gSeqUASBJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ0yACSpUwaAJHXKAJCkThkAktQpA0CSOmUASFKnDABJ6pQBIEmdMgAkqVNTeSi8JHVr+SmXPjq9cfVRC9iSuecRgCR1ygCQpE55CkhSl0ZP7fTKIwBJ6pQBIEmdMgAkqVMGgCR1ygCQpE4ZAJLUKQNAkjplAEhSpwwASerUpAGQZO8kf5Pk1iTrk7yxlT89yZVJvtX+3bWVJ8mHkmxIclOS549s6/hW/1tJjt923ZIkTWYqRwAPAydX1X7AwcCJSfYDTgGuqqoVwFVtHuAIYEV7rQLOhCEwgFOBg4ADgVO3hIYkaf5NGgBVdXdV3dCmHwRuA/YEjgbOa9XOA17Wpo8Gzq/B1cAuSfYADgOurKr7quoHwJXA4XPaG0nSlE3rGkCS5cABwDXA7lV1d1v0XWD3Nr0ncOfIane1sonKx+5jVZJ1SdZt3rx5Os2TJE3DlAMgyZOBTwJvqqoHRpdVVQE1Fw2qqrOqamVVrVy2bNlcbFKSNI4pBUCSJzH88b+gqj7Viu9pp3Zo/97byjcBe4+svlcrm6hckrQApnIXUICPArdV1ftGFq0FttzJczzwmZHyV7e7gQ4G7m+nii4HDk2ya7v4e2grkyQtgKk8EOYFwKuAm5Pc2MreBqwGLk7yWuAO4Ji27DLgSGAD8BPgNQBVdV+S04HrWr13VtV9c9ILSdK0TRoAVfVlIBMsfsk49Qs4cYJtrQHWTKeBkqRtw28CS1KnDABJ6pQBIEmdMgAkqVMGgCR1ygCQpE4ZAJLUKQNAkjplAEhSpwwASeqUASBJnZrKYHCSJGD5KZc+Or1x9VEL2JK54RGAJHXKAJCkTnkKSFI3Rk/hyCMASeqWASBJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ0yACSpUwaAJHXKAJCkThkAktQpA0CSOmUASFKnDABJ6tSkAZBkTZJ7k9wyUnZakk1JbmyvI0eWvTXJhiTfTHLYSPnhrWxDklPmviuSpOmYyhHAucDh45S/v6r2b6/LAJLsBxwLPLet8+dJdkiyA/Bh4AhgP+C4VleStEAmfSZwVX0xyfIpbu9o4KKqegj4TpINwIFt2Yaquh0gyUWt7q3TbrEkaU7M5hrASUluaqeIdm1lewJ3jtS5q5VNVP44SVYlWZdk3ebNm2fRPEnS1sw0AM4E9gX2B+4G3jtXDaqqs6pqZVWtXLZs2VxtVpI0xqSngMZTVfdsmU5yNnBJm90E7D1Sda9WxlbKJUkLYEZHAEn2GJl9ObDlDqG1wLFJdkqyD7ACuBa4DliRZJ8kOzJcKF4782ZLkmZr0iOAJBcChwC7JbkLOBU4JMn+QAEbgdcBVNX6JBczXNx9GDixqh5p2zkJuBzYAVhTVevnvDeSpCmbyl1Ax41T/NGt1D8DOGOc8suAy6bVOknSNuM3gSWpUwaAJHXKAJCkThkAktQpA0CSOjWjL4JJ0mKx/JRLF7oJ2y2PACSpUwaAJHXKAJCkThkAktQpA0CSOmUASFKnDABJ6pQBIEmdMgAkqVMGgCR1ygCQpE4ZAJLUKQNAkjplAEhSpwwASeqUASBJnfKBMJI0A6MPmtm4+qgFbMnMeQQgSZ0yACSpUwaAJHXKAJCkThkAktQpA0CSOmUASFKnDABJ6pQBIEmdmjQAkqxJcm+SW0bKnp7kyiTfav/u2sqT5ENJNiS5KcnzR9Y5vtX/VpLjt013JElTNZUjgHOBw8eUnQJcVVUrgKvaPMARwIr2WgWcCUNgAKcCBwEHAqduCQ1J0sKYdCygqvpikuVjio8GDmnT5wFfAP64lZ9fVQVcnWSXJHu0uldW1X0ASa5kCJULZ90DSRpjdJweTWym1wB2r6q72/R3gd3b9J7AnSP17mplE5U/TpJVSdYlWbd58+YZNk+SNJlZXwRun/ZrDtqyZXtnVdXKqlq5bNmyudqsJGmMmQbAPe3UDu3fe1v5JmDvkXp7tbKJyiVJC2SmAbAW2HInz/HAZ0bKX93uBjoYuL+dKrocODTJru3i76GtTJK0QCa9CJzkQoaLuLsluYvhbp7VwMVJXgvcARzTql8GHAlsAH4CvAagqu5LcjpwXav3zi0XhCVJC2MqdwEdN8Gil4xTt4ATJ9jOGmDNtFonSdpm/CawJHXKAJCkThkAktQpA0CSOmUASFKnDABJ6pQBIEmdMgAkqVMGgCR1ygCQpE4ZAJLUKQNAkjplAEhSpwwASerUpMNBS9Ji4IPgp88jAEnqlAEgSZ0yACSpUwaAJHXKAJCkThkAktQpA0CSOmUASFKn/CKYJM3S6JfQNq4+agFbMj0eAUhSpwwASeqUASBJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI6NasvgiXZCDwIPAI8XFUrkzwd+DiwHNgIHFNVP0gS4IPAkcBPgBOq6obZ7F9S33wK2OzMxRHAr1fV/lW1ss2fAlxVVSuAq9o8wBHAivZaBZw5B/uWJM3QtjgFdDRwXps+D3jZSPn5Nbga2CXJHttg/5KkKZhtABRwRZLrk6xqZbtX1d1t+rvA7m16T+DOkXXvamU/I8mqJOuSrNu8efMsmydJmshsB4N7YVVtSvIM4Mok3xhdWFWVpKazwao6CzgLYOXKldNaV5I0dbM6AqiqTe3fe4FPAwcC92w5tdP+vbdV3wTsPbL6Xq1MkrQAZhwASXZO8pQt08ChwC3AWuD4Vu144DNtei3w6gwOBu4fOVUkSZpnszkFtDvw6eHuTp4IfKyqPpvkOuDiJK8F7gCOafUvY7gFdAPDbaCvmcW+JUmzNOMAqKrbgV8Zp/z7wEvGKS/gxJnuT5I0t/wmsCR1ygCQpE4ZAJLUKQNAkjo12y+CSdK8cgC4ueMRgCR1ygCQpE4ZAJLUKa8BSNIcGr1GsXH1UQvYksl5BCBJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI65W2gkrZ7Dv+wbXgEIEmdMgAkqVMGgCR1ygCQpE4ZAJLUKe8CkrRd8s6fbc8jAEnqlAEgSZ3yFJCk7cZSO+2zvT8bwCMASeqUASBJnTIAJKlTBoAkdcqLwJIW1FK78LuYeAQgSZ3yCECS5sH2eEvovAdAksOBDwI7AOdU1er5boOk+eHpne3bvJ4CSrID8GHgCGA/4Lgk+81nGyRJg/k+AjgQ2FBVtwMkuQg4Grh1ntshaSv85L5tbS+ng+Y7APYE7hyZvws4aLRCklXAqjb7oyTfnMF+dgO+N6MWbv+Wat+War9g6fZtqfYL5rFvedc22eyzplJpu7sIXFVnAWfNZhtJ1lXVyjlq0nZlqfZtqfYLlm7flmq/YGn3bdR83wa6Cdh7ZH6vViZJmmfzHQDXASuS7JNkR+BYYO08t0GSxDyfAqqqh5OcBFzOcBvomqpavw12NatTSNu5pdq3pdovWLp9W6r9gqXdt0elqha6DZKkBeBQEJLUKQNAkjq1pAIgyWlJNiW5sb2OHFn21iQbknwzyWEL2c6ZSnJykkqyW5tPkg+1ft2U5PkL3cbpSnJ6a/uNSa5I8kutfFH3Lcl7knyjtf3TSXYZWbaofxeTvCLJ+iQ/TbJyzLLF3rfDW9s3JDlloduzzVXVknkBpwFvHqd8P+DrwE7APsC3gR0Wur3T7NveDBfP7wB2a2VHAn8NBDgYuGah2zmDfj11ZPo/An+xFPoGHAo8sU2/C3hXm14Kv4v/HHgO8AVg5Uj5ou4bw40p3wb+KbBj68t+C92ubflaUkcAW3E0cFFVPVRV3wE2MAxLsZi8H3gLMHrV/mjg/BpcDeySZI8Fad0MVdUDI7M781j/FnXfquqKqnq4zV7N8J0XWAK/i1V1W1WN9w39xd63R4eqqar/B2wZqmbJWooBcFI77F6TZNdWNt4QFHvOf9NmJsnRwKaq+vqYRYu6X1skOSPJncArgT9pxUuib83vMxzNwNLq11iLvW+Lvf3Ttt0NBTGZJJ8DfnGcRW8HzgROZ/gUeTrwXob/fNu9Sfr1NoZTCovS1vpWVZ+pqrcDb0/yVuAk4NR5beAMTdavVuftwMPABfPZttmaSt+0+C26AKiq35hKvSRnA5e02e1+CIqJ+pXkXzKcT/16EhjafkOSA1kE/YKp/8wY/khexhAA233fJutXkhOAlwIvqXaSmUXQL5jWz2zUoujbViz29k/bkjoFNOYc8cuBW9r0WuDYJDsl2QdYAVw73+2biaq6uaqeUVXLq2o5w2Hp86vquwz9enW7Y+Zg4P6qunsh2ztdSVaMzB4NfKNNL+q+tQcfvQX4rar6yciiRfu7OAWLvW/dDVWz6I4AJvHuJPsznALaCLwOoKrWJ7mY4bkDDwMnVtUjC9bKuXMZw90yG4CfAK9Z2ObMyOokzwF+ynCH0+tb+WLv258x3A1zZTtyu7qqXr8UfheTvBz4U2AZcGmSG6vqsMXet5q/oWq2Gw4FIUmdWlKngCRJU2cASFKnDABJ6pQBIEmdMgAkqVMGQKeS7JLkDxe6HXMpybHtm7ez2cbbxsx/dXatenQ7hyT5tRmst3HL6K8LpbX9kslrLuw2NX0GQL92AcYNgCTz8v2QsfuZ6n7bl8PG+909AvjsLJv1MwFQVdP+oz2BQ4C52taUJNlhzPxU39+l9v0gTcAA6NdqYN82Dv972ieyLyVZC9yaZHmSLd+kJsmbk5zWpvdN8tkk17d1fnnsxpPs3AbkuzbJ19qAdiQ5IcnaJJ8Hrhq731bnPyW5pb3e1MqWt3Haz2f4hvfeY/YXYH/ghjHlO7T+XdcGCXxdK98jyRdb/29J8qIkq4F/0souaPV+1P49JMnfJvlMktuTrE7yyta/m5Ps2+r9myTXtD5/LsnuSZYzfMHtj9q2X5RkWZJPtnZdl+QFbf1fyPBchPVJzmEYDvtxkhya5O+S3JDkE0me3Mo3JnlXkhuAVyT5QpIPJFkHvLG9j59v78VVSZ7Z1js3yV8kuQZ490S/NFv5uV6d5Lkj9b6QZOVE9bWdWOjxqH0tzAtYDtwyMn8I8GNgnwmWvxk4rU1fBaxo0wcBnx9n+/8N+L02vQvw9wzDPZ/AMJzF0yfY768CN7e6TwbWAwe09vwUOHiC/jyfYfjoseWrgP/cpncC1jGMrXQyw8BmMHzr8ylt+kdj1v/RSDt/COzRtrMJeEdb9kbgA216Vx77guW/B97bpk9j5FkVwMeAF7bpZwK3tekPAX/Spo9i+Fb7bmPatBvwRWDnNv/HI+tsBN4yUvcLwJ+PzP8VcHyb/n3gf7XpcxnGznrc+P2t75dM8nP9o5H3Yw/gm5PUf3Sbvhbu5aGeRl1bwzjuE2qfNH8N+MTwoRsY/iCOdSjwW0ne3OZ/juEPHcCVVXXfBPt9IfDpqvpx29+ngBcxjMlyRw3PBhjP4Tw25PLYdjwvyb9t809jGKPmOmBNkicx/BG8caI+j7iu2nhESb4NXNHKbwZ+vU3vBXw8w7hUOwITvZ+/Aew38h4+tb23LwZ+G6CqLk3yg3HWPZjh4StfaevvCPzdyPKPj6k/Ov+vt2wf+B/87Kf9T9TkQzdM9HO9mOH9OBU4BvjLSeprO2AAaNSPR6Yf5mdPEf5c+/cJwA+rav9JthXgd2rMg0OSHDRmP2P3O9X2jXUo8DsTtOMNVXX54xYkL2b4lH1ukvdV1fmT7P+hkemfjsz/lMf+L/0p8L6qWpvkEIZP/uN5AsPRzP8d06ZJmjBUYwjR4yZYvi3e39F9P+7nCpDk+0meB/w7HhvTaaLfg92n2CZtQ14D6NeDwFO2svwe4BntnPRODMMaU8MTvL6T5BXw6AXZXxln/cuBN7Rz8yQ5YIrt+hLwsiQ/n2RnhlFdv7S1FZI8jeHxi9+foB1/0D7pk+TZ7bz0s4B7qups4ByGU0gA/7il7gw9jceGED5+pHzs+30F8IaRPmwJ1C8Cv9vKjmA4pTTW1cALkvyzVm/nJM+eYvu+yjDKJQwP4NnqezuOrf1cP84wAurTquqmKdTXAjMAOtX+WH6lXQB9zzjL/xF4J8Nwvlfy2DDNMPzheG2SrzOcox/vwt7pwJOAm5Ksb/NTadcNDOejrwWuAc6pqq9NstpvAp+bYNk5DBeXb8hwUfsjDJ/WD2F4xsLXGD6xfrDVP6u1eaYPcDmN4fTY9cD3Rsr/Cnj5lovADM8/Xtkuxt7KY5+Y3wG8uL1nvw3877E7qKrNDNdSLkxyE8Ppn8ddiJ/AG4DXtPVexXD9Yjq29nP9S4ZwuXiK9bXAHA1Ui167W+acrVwfkDQOA0CSOuUpIEnqlAEgSZ0yACSpUwaAJHXKAJCkThkAktSp/w+I5Kt+tvusCgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f1cc4c17e10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RMSE: 1.63112\n",
      "confidence level: 0.4231\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEWCAYAAABollyxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAHDBJREFUeJzt3X2UHVWZ7/Hvj3cMDOGlzYQkEJZGXbnK2/SCCOIgCIsAGpjhTRkITLzRGeCiwIIo9ypex7viOAPCOINGYAgOA8QIJgIjYIBhVIJ0QgyQwBAgMQmBtIGEVxHkuX/UDimaPjl1uk/36d79+6x11qnatavqqdPJc/bZVbVLEYGZmeVri1YHYGZmfcuJ3swsc070ZmaZc6I3M8ucE72ZWeac6M3MMudEb/1G0iGSHm91HIOVpFMl3dnqOGzwka+jt2aTtBz4XET8vNWxmJlb9DYESNqqGXWavU+z/uJEb/1G0qGSVpXml0u6QNJiSRsk3SRpu9LyYyUtkrRe0q8k7V1aNk3Sk5JekrRE0vGlZWdI+qWkyyStAy7pJpZLJM2W9G+SXgTOkLRFabvrJM2StEtpndMlrUjL/k+K/5M92Z6k7VLdden4HpQ0ohT/U+nYnpZ0aqn8F6V4DkrrbUjvB5WW3SvpG+lzeEnSnZJ2691f0AYrJ3prtZOAo4C9gL2BMwAk7QdcA3we2BX4PjBX0rZpvSeBQ4CdgK8D/yZpZGm7BwJPASOAb9bY9yRgNjAcuB44BzgO+HNgd+AF4J9TPOOBfwFOBUam/Y7q6faAyWkbY9LxfQF4TdIw4ApgYkTsCBwELOoaePrCuC3V3RW4FLhN0q6lap8FzgTeC2wDXFDjc7DMOdFbq10REc9ExPPAT4F9U/lU4PsR8UBE/DEiZgKvAxMAIuJHab23IuIm4AnggNJ2n4mIf4qINyPitRr7vj8ifpK28RpFsr04IlZFxOsUvwROSN0wJwA/jYhfRMQfgK8CXU9wNbK9NygS9PvT8S2IiBfTdt4CPixp+4hYExGPdhP7McATEfHDdIw3AI8BnyrV+deI+O8Uy6zSZ2tDjBO9tdqzpelXgR3S9J7A+albY72k9RSt393h7W6URaVlHwbKXRMrK+y7a509gVtK21wK/JHiV8Hu5foR8Sqwrhfb+yFwB3CjpGck/b2krSPiFeBkii+JNZJuk/ShbmLfHVjRpWwF7/yVUeuztSHGid4GqpXANyNieOn1noi4QdKewA+As4FdI2I48Aig0vpVLifrWmclRZdJeZ/bRcRqYA0wemNFSdtTtMh7tL2IeCMivh4R4ym6Z44FTgeIiDsi4giKLqLH0rF29QzFF0nZHsDqCsdtQ4wTvfWVrdMJx42vRq9C+QHwBUkHqjBM0jGSdgSGUSTVTgBJZ1K06Hvre8A30xcJktokTUrLZgOfSidAt6HohlH3m6m/PUmfkPQRSVsCL1J05bwlaYSkSamv/nXgZYqunK5uBz4g6bOStpJ0MjAeuLXnh2+5cqK3vnI78FrpdUkjK0dEB/A/ge9SnMRcRjpRGxFLgH8E7geeAz4C/LIJMV8OzAXulPQSMJ/ipC6pn/wc4EaK1v3LwFqKZNzw9oA/pfjyeJGiS+c/KbpztgDOo2ixP09xIvdvum44ItZR/Ao4n6IL6ULg2Ij4Xc8O3XLmG6bMekDSDsB6YFxEPN3qeMw2xy16s4okfUrSe1K3yj8ADwPLWxuVWX1O9GbVTaLoUnkGGAecEv5JbIOAu27MzDJXqUUv6UuSHpX0iKQb0lUUe0l6QNIyFbeub5Pqbpvml6XlY/vyAMzMbPPqtugljQJ+AYyPiNckzaK4ouJo4OaIuFHS94DfRMSVkv4W2DsiviDpFOD4iDh5c/vYbbfdYuzYsc04HjOzIWPBggW/i4i2evWqXtu8FbC9pDeA91BcXnYYxVgaADMpLp+7kqIf85JUPhv4riRtri9z7NixdHR0VAzFzMwAJHW9O7pbdbtu0l2B/wD8liLBbwAWAOsj4s1UbRWbbr0eRboVPC3fwLvvIDQzs35SN9FL2pmilb4XxfgawyhGG+wVSVMldUjq6Ozs7O3mzMyshionYz8JPB0RnRHxBnAzcDAwvHRb+2g2jbGxmmLwqY0PX9iJdw/+RETMiIj2iGhva6vbxWRmZj1UJdH/FpiQbhQRcDiwBLiHYuhWKMbWnpOm56Z50vK7fa2xmVnrVOmjf4DipOpCijsBtwBmABcB50laRtEHf3Va5Wpg11R+HjCtD+I2M7OKBsQNU+3t7eGrbszMGiNpQUS016vnIRDMzDLnRG9mljknejOzzDX61B8zG0LGTrvt7enl049pYSTWG070ZplzsjZ33ZiZZc6J3swsc070ZmaZc6I3M8ucE72ZWeZ81Y2Z+cqczDnRmw1R5eRueXOiN7OGdf2S8K+Agc2J3mwIcSt+aPLJWDOzzLlFb5YJn1C1WtyiNzPLnBO9mVnm6iZ6SR+UtKj0elHSFyXtIukuSU+k951TfUm6QtIySYsl7d/3h2FmZrXU7aOPiMeBfQEkbQmsBm6heOj3vIiYLmlamr8ImAiMS68DgSvTu5kNYr5iZ/BqtOvmcODJiFgBTAJmpvKZwHFpehJwXRTmA8MljWxKtGZm1rBGE/0pwA1pekRErEnTzwIj0vQoYGVpnVWpzMzMWqDy5ZWStgE+DXy567KICEnRyI4lTQWmAuyxxx6NrGpmfchdNPlp5Dr6icDCiHguzT8naWRErEldM2tT+WpgTGm90ansHSJiBjADoL29vaEvCTMbWHwN/8DWSNfNZ9jUbQMwF5icpicDc0rlp6erbyYAG0pdPGZm1s8qteglDQOOAD5fKp4OzJI0BVgBnJTKbweOBpYBrwJnNi1aMzNrWKVEHxGvALt2KVtHcRVO17oBnNWU6MzMrNd8Z6yZWeY8qJlZhnzljJW5RW9mljknejOzzDnRm5llzonezCxzTvRmZplzojczy5wvrzQbxHwZpVXhFr2ZWeac6M3MMudEb2aWOSd6M7PMOdGbmWXOid7MLHNO9GZmmfN19GbWVH5+7MDjFr2ZWeYqJXpJwyXNlvSYpKWSPippF0l3SXoive+c6krSFZKWSVosaf++PQQzM9ucqi36y4GfRcSHgH2ApcA0YF5EjAPmpXmAicC49JoKXNnUiM3MrCF1E72knYCPA1cDRMQfImI9MAmYmarNBI5L05OA66IwHxguaWTTIzczs0qqnIzdC+gE/lXSPsAC4FxgRESsSXWeBUak6VHAytL6q1LZmlIZkqZStPjZY489ehq/2ZDjgcysUVW6brYC9geujIj9gFfY1E0DQEQEEI3sOCJmRER7RLS3tbU1sqqZmTWgSqJfBayKiAfS/GyKxP/cxi6Z9L42LV8NjCmtPzqVmZlZC9RN9BHxLLBS0gdT0eHAEmAuMDmVTQbmpOm5wOnp6psJwIZSF4+ZmfWzqjdMnQNcL2kb4CngTIoviVmSpgArgJNS3duBo4FlwKuprpmZtUilRB8Ri4D2bhYd3k3dAM7qZVxmZtYkHgLBzPqMh0MYGDwEgplZ5pzozcwy50RvZpY599GbDQK+G9Z6wy16M7PMOdGbmWXOid7MLHNO9GZmmXOiNzPLnBO9mVnmnOjNzDLn6+jNrF943JvWcaI3G6B8k5Q1i7tuzMwy50RvZpY5J3ozs8w50ZuZZa5Sope0XNLDkhZJ6khlu0i6S9IT6X3nVC5JV0haJmmxpP378gDMzGzzGrnq5hMR8bvS/DRgXkRMlzQtzV8ETATGpdeBwJXp3cwM8KWW/a03XTeTgJlpeiZwXKn8uijMB4ZLGtmL/ZiZWS9UTfQB3ClpgaSpqWxERKxJ088CI9L0KGBlad1VqewdJE2V1CGpo7Ozswehm5lZFVW7bj4WEaslvRe4S9Jj5YUREZKikR1HxAxgBkB7e3tD65qZWXWVWvQRsTq9rwVuAQ4AntvYJZPe16bqq4ExpdVHpzIzM2uBuole0jBJO26cBo4EHgHmApNTtcnAnDQ9Fzg9XX0zAdhQ6uIxM7N+VqXrZgRwi6SN9f89In4m6UFglqQpwArgpFT/duBoYBnwKnBm06M2M7PK6ib6iHgK2Keb8nXA4d2UB3BWU6IzM7Ne852xZmaZc6I3M8ucE72ZWeac6M3MMudEb2aWOT9K0GwA8eMDrS840ZtZS3kky77nrhszs8w50ZuZZc6J3swsc070ZmaZc6I3M8ucE72ZWeZ8eaVZi/naeetrTvRmLeDkbv3Jid6snzi5W6u4j97MLHNO9GZmmauc6CVtKekhSbem+b0kPSBpmaSbJG2TyrdN88vS8rF9E7qZmVXRSIv+XGBpaf5bwGUR8X7gBWBKKp8CvJDKL0v1zMzqGjvttrdf1jyVEr2k0cAxwFVpXsBhwOxUZSZwXJqelOZJyw9P9c3MrAWqtui/A1wIvJXmdwXWR8SbaX4VMCpNjwJWAqTlG1L9d5A0VVKHpI7Ozs4ehm9mZvXUTfSSjgXWRsSCZu44ImZERHtEtLe1tTVz02ZmVlLlOvqDgU9LOhrYDvgT4HJguKStUqt9NLA61V8NjAFWSdoK2AlY1/TIzcyskrot+oj4ckSMjoixwCnA3RFxKnAPcEKqNhmYk6bnpnnS8rsjIpoatZmZVdab6+gvAs6TtIyiD/7qVH41sGsqPw+Y1rsQzcysNxoaAiEi7gXuTdNPAQd0U+f3wIlNiM3MzJrAd8aamWXOid7MLHNO9GZmmXOiNzPLnMejN2uy8jgty6cf08JIBjd/js3jFr2ZWeac6M3MMudEb2aWOSd6M7PMOdGbmWXOV92Y9SE/KckGArfozcwy50RvZpY5J3ozs8w50ZuZZc6J3swsc070ZmaZc6I3M8tc3evoJW0H3Adsm+rPjoivSdoLuJHiebELgNMi4g+StgWuA/4MWAecHBHL+yh+swHB18vbQFalRf86cFhE7APsCxwlaQLwLeCyiHg/8AIwJdWfAryQyi9L9czMrEXqJvoovJxmt06vAA4DZqfymcBxaXpSmictP1ySmhaxmZk1pFIfvaQtJS0C1gJ3AU8C6yPizVRlFTAqTY8CVgKk5Rsoune6bnOqpA5JHZ2dnb07CjMzq6lSoo+IP0bEvsBo4ADgQ73dcUTMiIj2iGhva2vr7ebMzKyGhgY1i4j1ku4BPgoMl7RVarWPBlanaquBMcAqSVsBO1GclDUz6xE/VrB36rboJbVJGp6mtweOAJYC9wAnpGqTgTlpem6aJy2/OyKimUGbmVl1VVr0I4GZkrak+GKYFRG3SloC3Cjp74CHgKtT/auBH0paBjwPnNIHcZuZWUV1E31ELAb266b8KYr++q7lvwdObEp0ZmbWa74z1swsc070ZmaZc6I3M8ucE72ZWeac6M3MMudEb2aWOSd6M7PMOdGbmWXOid7MLHMNDWpmZpv4qVI2WLhFb2aWOSd6M7PMuevGbDO6ds94LHQbjJzozRrgfnkbjJzozWxQ8dOmGuc+ejOzzDnRm5llzonezCxzVR4OPkbSPZKWSHpU0rmpfBdJd0l6Ir3vnMol6QpJyyQtlrR/Xx+EmZnVVqVF/yZwfkSMByYAZ0kaD0wD5kXEOGBemgeYCIxLr6nAlU2P2szMKqub6CNiTUQsTNMvAUuBUcAkYGaqNhM4Lk1PAq6LwnxguKSRTY/czMwqaaiPXtJYYD/gAWBERKxJi54FRqTpUcDK0mqrUlnXbU2V1CGpo7Ozs8GwzcysqsqJXtIOwI+BL0bEi+VlERFANLLjiJgREe0R0d7W1tbIqmZm1oBKiV7S1hRJ/vqIuDkVP7exSya9r03lq4ExpdVHpzIzM2uBunfGShJwNbA0Ii4tLZoLTAamp/c5pfKzJd0IHAhsKHXxmA14HubAclNlCISDgdOAhyUtSmVfoUjwsyRNAVYAJ6VltwNHA8uAV4EzmxqxmZk1pG6ij4hfAKqx+PBu6gdwVi/jMjOzJvGdsWZmmXOiNzPLnBO9mVnmnOjNzDLnB4+Y2aDlh5BU4xa9mVnmnOjNzDLnrhszfDes5c0tejOzzLlFb2ZZ8InZ2pzobUhxMrChyF03ZmaZc4vehiyfgLWhwi16M7PMuUVvZtnxuZh3cqK37LmLxoY6d92YmWXOid7MLHNVHg5+DXAssDYiPpzKdgFuAsYCy4GTIuKF9CDxyymeGfsqcEZELOyb0M1qc3eN2SZVWvTXAkd1KZsGzIuIccC8NA8wERiXXlOBK5sTppmZ9VTdRB8R9wHPdymeBMxM0zOB40rl10VhPjBc0shmBWtmZo3raR/9iIhYk6afBUak6VHAylK9VansXSRNldQhqaOzs7OHYZiZWT29PhkbEQFED9abERHtEdHe1tbW2zDMzKyGnib65zZ2yaT3tal8NTCmVG90KjMzsxbpaaKfC0xO05OBOaXy01WYAGwodfGYmVkLVLm88gbgUGA3SauArwHTgVmSpgArgJNS9dspLq1cRnF55Zl9ELPZ23yru1l9KrrYW6u9vT06OjpaHYYNQr5e3hqRW2NA0oKIaK9Xz3fGmpllzoOa2aDgLhqznnOit0HH3TVmjXHXjZlZ5tyiN7MhY6h2AbpFb2aWOSd6M7PMuevGBiyfdDVrDid6G1Cc3M2az4neWs7J3axvuY/ezCxzbtFbS7gVb9Z/nOitTzmh20A1lK6pd9eNmVnmnOjNzDLnrhtriqH0M9jyk/u/Xyd661aVf/i1+t/dL282sDjRW0OcxC13Obbu+yTRSzoKuBzYErgqIqb3xX6seZzAzfLV9GfGStoS+G/gCGAV8CDwmYhYUmsdPzO2+Wq1SpzQzXpmILbuqz4zti9a9AcAyyLiqRTIjcAkoGaiHyoa/UnYrGTt5G7WXLX+Tw3ELwPomxb9CcBREfG5NH8acGBEnN2l3lRgapr9IPB4UwPpvd2A37U6iB5w3P3LcfefwRgz9G3ce0ZEW71KLTsZGxEzgBmt2n89kjqq/CQaaBx3/3Lc/WcwxgwDI+6+uGFqNTCmND86lZmZWQv0RaJ/EBgnaS9J2wCnAHP7YD9mZlZB07tuIuJNSWcDd1BcXnlNRDza7P30gwHbrVSH4+5fjrv/DMaYYQDE3fSTsWZmNrB4UDMzs8w50ZuZZc6JvgJJ50sKSbu1OpYqJH1D0mJJiyTdKWn3VsdUhaRvS3osxX6LpOGtjqkeSSdKelTSW5IG/KV/ko6S9LikZZKmtTqeKiRdI2mtpEdaHUsjJI2RdI+kJenfyLmtisWJvg5JY4Ajgd+2OpYGfDsi9o6IfYFbga+2OqCK7gI+HBF7Uwyj8eUWx1PFI8BfAPe1OpB60vAk/wxMBMYDn5E0vrVRVXItcFSrg+iBN4HzI2I8MAE4q1WftxN9fZcBFwKD5qx1RLxYmh3GIIk9Iu6MiDfT7HyKezAGtIhYGhED7a7uWt4eniQi/gBsHJ5kQIuI+4DnWx1HoyJiTUQsTNMvAUuBUa2IxcMUb4akScDqiPiNpFaH0xBJ3wROBzYAn2hxOD3x18BNrQ4iM6OAlaX5VcCBLYplSJE0FtgPeKAV+x/yiV7Sz4E/7WbRxcBXKLptBpzNxR0RcyLiYuBiSV8Gzga+1q8B1lAv7lTnYoqfvdf3Z2y1VInZrBZJOwA/Br7Y5dd2vxnyiT4iPtlduaSPAHsBG1vzo4GFkg6IiGf7McRu1Yq7G9cDtzNAEn29uCWdARwLHB4D5CaPBj7rgc7Dk/QzSVtTJPnrI+LmVsUx5BN9LRHxMPDejfOSlgPtETHgR8+TNC4inkizk4DHWhlPVemBNRcCfx4Rr7Y6ngy9PTwJRYI/Bfhsa0PKl4oW4tXA0oi4tJWx+GRsnqZLekTSYoqup5Zd1tWg7wI7AnelS0O/1+qA6pF0vKRVwEeB2yTd0eqYakknujcOT7IUmDUYhieRdANwP/BBSaskTWl1TBUdDJwGHJb+PS+SdHQrAvEQCGZmmXOL3swsc070ZmaZc6I3M8ucE72ZWeac6M3MMudEnzFJwyX9bavjaCZJp6Q7Z3uzja90mf9V76J6ezuHSjqoB+stb/XIqCn2Wwf6Nq1nnOjzNhzoNtFL6peb5brup+p+Veju3+dE4Ge9DOsdiT4iGk7ONRwKNGtblaQRKcvzVT9f3yw5hDjR52068L50o8a3UwvrvyTNBZZIGlse41vSBZIuSdPvk/QzSQvSOh/qunFJw9JY4b+W9FAaBA5JZ0iaK+luYF7X/aY656Wbuh6R9MVUNjaNlX4dxfC/Y7rsT8C+wMIu5Vum43swjWX/+VQ+UtJ96fgfkXSIpOnA9qns+lTv5fR+qKT/lDRH0lOSpks6NR3fw5Lel+p9StID6Zh/LmlEGrTqC8CX0rYPkdQm6ccprgclHZzW31XFcwIelXQV0O2IeZKOlHS/pIWSfqRizJSNvwC+JWkhcKKkeyV9R1IHcG76HO9On8U8SXuk9a6V9D1JDwB/X+sfzWb+rvMl/Y9SvXsltdeqbwNIRPiV6QsYCzxSmj8UeAXYq8byC4BL0vQ8YFyaPhC4u5vt/z/gr9L0cIox5IcBZ1CMjLhLjf3+GfBwqrsD8CjFyH5jgbeACTWOZ3/gum7KpwL/O01vC3RQjFN0PsXAY1A8qH7HNP1yl/VfLsW5HhiZtrMa+Hpadi7wnTS9M5tuNvwc8I9p+hLggtJ2/x34WJreg+JWeIArgK+m6WMohpHerUtMu1GMcT8szV9UWmc5cGGp7r3Av5TmfwpMTtN/DfwkTV9L8XyCLbv5DA8Fbq3zd/1S6fMYCTxep/7b2/SrtS//fBt6fh0RT2+uQmo5HgT8SJuGZ962m6pHAp+WdEGa344ioQHcFRHlMcTL+/0YcEtEvJL2dzNwCDAXWBER82uEdhTwHzXi2FvSCWl+J2Acxdgu16gYWOonEbGo1jGXPBgRa1JcTwJ3pvKH2TTc82jgJkkjgW2AWp/nJ4Hxpc/wT9Jn+3GKh5UQEbdJeqGbdSdQPBzkl2n9bSiGAdio6xDO5fmPbtw+8EPe2Xr/UUT8sUa8G9X6u86i+Dy+BpwEzK5T3wYIJ/qh55XS9Ju8s/tuu/S+BbA+iidUbY6Av4wuD96QdGCX/XTdb9X4ujoS+MsacZwTEe8aZ0bSxylazddKujQirquz/9dL02+V5t9i0/+XfwIujYi5kg6laMl3ZwuKXye/7xJTnRCKahRflp+psbwvPt/yvt/1dwWQtE7S3sDJFF1VNetLGlExJutj7qPP20sUg4TV8hzw3tRnvC3F8MBEMWb205JOhLdPjO7Tzfp3AOekvnMk7Vcxrv8CjpP0HknDgONTWU2SdgK2ioh1NeL4m9RyR9IHUr/xnsBzEfED4CqKrh+ANzbW7aGd2DS87+RSedfP+07gnNIxbPzivI80aqSkiRRdQV3NBw6W9P5Ub5ikD1SM71cUI1MCnEqdz7Ybm/u73kQxwuhOEbG4Qn0bAJzoM5aS4i/Tichvd7P8DeD/Ar+meF5reTjjU4Epkn5D0Yfe3Qm2bwBbA4slPZrmq8S1kKK/+NcUT9y5KiIeqrPaEcDPayy7iuIk70IVJ5e/T9H6PpTieQIPUbRAL0/1Z6SYe/pgk0sourUWAOVhq38KHL/xZCzwv4D2dFJ0CZtawF8HPp4+s7+gm+cRR0QnxbmOG1SMQno/8K4T4jWcA5yZ1juNxkcv3dzfdTbFl8isivVtAPDolTYopKtTrtpM/72Z1eBEb2aWOXfdmJllzonezCxzTvRmZplzojczy5wTvZlZ5pzozcwy9/8BQhjWklAK3GwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f1cc29d98d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RMSE: 0.805382\n",
      "confidence level: 0.84405\n"
     ]
    }
   ],
   "source": [
    "recs = [SimpleReconstruction(), LinearRegressionReconstruction()]\n",
    "n_samples = 20000\n",
    "\n",
    "rel_err_sample = np.zeros((n_samples, len(recs)))\n",
    "err_sample = 0*rel_err_sample\n",
    "\n",
    "for i in range(n_samples):\n",
    "    actual = Simulation()\n",
    "    meas = actual.simulate_measurement()\n",
    "    \n",
    "    for j in range(len(recs)):\n",
    "        r, r_err = recs[j].compute(meas)\n",
    "        err = r - actual.R\n",
    "        \n",
    "        err_sample[i, j] = err\n",
    "        rel_err_sample[i, j] = err / r_err\n",
    "\n",
    "rmse = np.sqrt(np.mean(err_sample**2, axis=0))\n",
    "\n",
    "for j in range(len(recs)):\n",
    "    plt.hist(rel_err_sample[:,j], bins=100)\n",
    "    plt.xlabel('true error / estimated error level')\n",
    "    plt.title(['Simple', 'Linear regression'][j])\n",
    "    plt.show()\n",
    "    print('RMSE: %g' % np.sqrt(np.mean(err_sample[:,j]**2)))\n",
    "    print('confidence level: %g' % np.mean(np.abs(rel_err_sample[:,j]) < 1))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}