Skip to content

Instantly share code, notes, and snippets.

@joezuntz
Last active May 18, 2018 14:54
Show Gist options
  • Save joezuntz/f08bf7f5fc69cc7e60bf7851efb90906 to your computer and use it in GitHub Desktop.
Save joezuntz/f08bf7f5fc69cc7e60bf7851efb90906 to your computer and use it in GitHub Desktop.
Correlating noise can improve constraints
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Populating the interactive namespace from numpy and matplotlib\n"
]
}
],
"source": [
"pylab inline"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"import scipy.optimize"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"# Mock data - just two points\n",
"x = np.array([1.0, 2.0])\n",
"y = np.array([1.0, 2.0])\n",
"\n",
"# Fitting function - just fit the gradient of a straight line\n",
"f = lambda x, p0: p0*x\n",
"\n",
"# Standard deviations of individual data points\n",
"sigma_x = 1.0\n",
"sigma_y = 1.0\n",
"\n",
"# We will vary the covariance between the data points\n",
"def make_cov(rho):\n",
" return np.array([[sigma_x**2, rho*sigma_x*sigma_y], [rho*sigma_x*sigma_y, sigma_y**2]])\n",
"\n",
"# Loop through possible correlation values\n",
"Rho = np.linspace(-0.99, 0.99, 101)\n",
"results = []\n",
"for rho in Rho:\n",
" C = make_cov(rho)\n",
" mu, c = scipy.optimize.curve_fit(f, x, y, p0=[0.5], sigma=C, absolute_sigma=True)\n",
" results.append(c[0,0]**0.5)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(-1, 1)"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEOCAYAAAB1g0unAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xl8FPX9x/HXJyc5IJBwQwj3EeSO\ngPcFiidataJF631UtLb1Z7W2trX2tFV7oKLWE2+qFhUPwANEUW4h4Qp3uEm4Qsj9+f2xE11pjgF2\ndnaTz/Px2Ed2Zmdn3pkcn535zny/oqoYY4wxbsT4HcAYY0z0sKJhjDHGNSsaxhhjXLOiYYwxxjUr\nGsYYY1yzomGMMcY1KxrGGGNcs6JhjDHGNSsaxhhjXIvzO8Dhat26tXbt2tXvGMYYE1UWLFiwS1Xb\nHO16oq5odO3alfnz5/sdwxhjooqIbAjFeuz0lDHGGNesaBhjjHHNioYxxhjXrGgYY4xxzYqGMcYY\n1zwtGiIyRkRWiki+iNxdy+tXi8hOEVnsPK73Mo8xxpij49kltyISC0wERgMFwDwRmaqqeYcs+qqq\nTvAqhzHGmNDx8j6N4UC+qq4FEJFXgLHAoUXDGNNEqCoHyqsoLC5jV3E5+0orKC6tZH9pJQcrqiiv\nrKaiqprKqmpEhBgRYmOgWXwsKYlxJCfE0iIpnoyUBDJSE8lISaBZfKzf31aT4mXR6ARsCpouAEbU\nstzFInIysAr4iapuOnQBEbkRuBGgS5cuHkQ1xoSCqrKzuIz1u0rYUHiAjUUlFOw+yNa9B9m6t5Rt\ne0spq6wO6TbbNk8kMz2ZLunJ9GybSt/2zenTvjmdWiYhIiHdlvH/jvC3gZdVtUxEbgKeA04/dCFV\nfQJ4AiAnJ0fDG9EYU5u9BytYvnUfy7fuY8XW/azesZ/8HcXsK638ZpkYgQ5pSXRIa8bAzi05q38z\nWqcmkJGSSHpqAmlJ8bRoFkdqYjxJCbEkxsUQHxtDbIygqlQrVFZXU1pRzcHyKg6UV7L3YAWFxeUU\nHShj+74yNhWVsGl3CV+uLeTNRZu/2XZ6SgJDu7RiWFYrhndrxaDOLYmLtWt/jpaXRWMzkBk03dmZ\n9w1VLQyafAr4i4d5jDFH6GB5FUsK9rBk0x6+3ryXrwv2sKno4Devp6ck0KttKucP6kjPtql0a51C\nVkYKnVomkRB3ZP+oRYRYgdiYWBLjYklLim/wPftLK1i1fT/Lt+5nyaY9LNiwmxnLtwPQolkcJ/Vu\nw2l92jK6XzvSkhten/lfourNB3cRiSNwyukMAsViHnCFquYGLdNBVbc6zy8Cfq6qI+tbb05Ojlrf\nU8Z4a+f+MuatL+KrdUUs2LCbvK37qKoO/K/o1DKJQZlpHNMpjX4dWpDdoQVtmydG7KmgogPlfLGm\nkE9W7uCTVTvZub+MhNgYTunThrGDOzKqX7sm0S4iIgtUNedo1+PZkYaqVorIBOADIBZ4WlVzReR+\nYL6qTgVuF5ELgEqgCLjaqzzGmLrtKQn8Y52zZhdfrClkzc4DACTFxzIoM42bT+nO0C6tGJzZkozU\nRJ/THp70lATOHdiBcwd2oLpaWbp5L1OXbOHtJVuYnredVsnxXHZsF8aP7ELnVsl+x414nh1peMWO\nNIw5epVV1SzatIdPV+5k1uqdLN28F1VISYjl2G7pjOiWwYju6QzolEZ8I20HqKpWvlhTyOS5G/gw\nbxsAZ2a357YzetK/Y5rP6UIvVEcaVjSMaSKKDpTzycodzFyxg1mrdrK/tJLYGGFIZktO7NWaE3u2\nZlBmy0ZbJOqzec9BXpy7gRfmbmB/aSVj+rfnjtG96Nu+hd/RQsaKhjGmQet3HWB63nY+zNvGgg27\nqVZo0zyR0/oEGoSP79naVQNzU7H3YAX//mwdz3y2juLySi7LyeSuMX1JT0nwO9pRs6JhjKnV6u37\nmbZ0G+8t28qKbfsByO7QgtHZ7RjVrx39O7YgJiYyG60jxZ6SciZ+nM8zc9aTkhjHnWf14YrhXYiN\n4v1mRcMY8438HcW88/UW3vl6K/k7ihGBY7PSOeuY9pyZ3Y7MdGvgPRKrt+/n11Nz+XxNIUO7tOTh\nywaTlZHid6wjYkXDmCZuy56DTF2yhamLt5C3dR8iMLxrOucN7MBZ/dvTtkUzvyM2CqrKfxdv4Vf/\nXUZ1tfLr8/tzaU7niL3EuC5WNIxpgvaXVvDe0m28saiAL9cVoQpDurTkgkEdOWdAB9pZofDM5j0H\n+dlri5m7tohzBrTnwUsGkZLod6ca7kX8fRrGmNCoqlbm5O/iPwsL+CB3G6UV1XRrncIdZ/TmwiEd\no/Z0SbTp1DKJl64fyROz1/KX91ewducBnrwqp8md+rOiYUyE2lRUwuvzNzFlQQFb9pbSolkcFw/t\nzMXDOjMks2XUnR5pDGJihJtP6UF2hxZMeGkhYyfO4dEfDGVk9wy/o4WNnZ4yJoKUV1YzY/l2Xv5q\nI7NX70IETurVhkuHdWZ0dtPo7iJarNt1gOufm8eGwhL+Pm4I5w7s4HeketnpKWMakU1FJbz45Uam\nLNjEruJyOrVM4iejenNJTmc6tUzyO56pRbfWKbx56wlc9+w8bnt5IaUVg7h4WGe/Y3nOioYxPqmq\nVj5ZuYMX5m7g01U7EeCMfu24YkQXTu7VJqrvCWgqWjSL57lrh3PD8/P52etLKK2s4gcjsvyO5Skr\nGsaE2Z6Scl6bv4kX5m5gU9FB2jZP5LbTe3H58Ew6pNlRRbRJTojj3z88lh+9uJB731yGKowf2XgL\nhxUNY8Jk9fb9PD1nPW8uKqC0oprh3dK5e0w/zuzfrkn299SYNIuP5fHxw7hl8gLu++8y2jRP5Kz+\n7f2O5QkrGsZ4qLpa+XT1Tp7+bB2zV+8iMS6Gi4Z04qrjupLdsfF0hmcgIS6Gf14xhCue/JLbX17E\nSzeMYFhWut+xQs6unjLGA6UVVby1aDNPfbaO/B3FtGuRyFXHdeXy4V0aRed3pm6FxWVc8vgX7C4p\nZ8rNx9OzbarfkQC7I9zvGMbUak9JOZPnbuDZz9ezq7ic7A4tuOHkbpw7oOMRD3tqos/GwhK+99gc\nkhPiePu2EyOiJ2G75NaYCLJlz0Gemr2OV+ZtpKS8ilP7tOHGk7pzXI8MuwmvCeqSkcykK4dx2aS5\n3DVlCY+PH9Zofg+saBhzFPJ3FPP4p2t4a9FmAC4Y1JEbTu5Ovw7WXtHUDctK5+6z+/LAu8t5es56\nrjuxm9+RQsKKhjFHIHfLXiZ+nM97y7aRGBfD+JFZXH9SNxtj2nzHdSd248t1Rfxx2nIGZ7ZkWFYr\nvyMdNWvTMOYwLNy4m399lM9HK3bQPDGOq47P4toTupGRmuh3NBOh9pZUcN6/ZlNVpbx3x8m+tW9Y\nm4YxYbRgQxGPzFjN7NW7aJUcz51n9ubK47pGRAOniWxpyfH86/KhXPToHP44bTl/unig35GOihUN\nY+oRXCwyUhK4++y+XDkyK6rGUTD+G5TZkhtO6s6kWWu5YFBHju/Z2u9IR8x+842pxeJNe3h4+io+\nXbWTjJQEfnFOX8aPzCI5wf5kzJG5Y1RvPsjdxt1vLOWDO04mKSE6eyy2vwBjguRt2cdD01cyY/kO\nWiXHc8/ZfbnyOCsW5uglJcTyx+8N5PIn5/LwjFX84px+fkc6IvaXYAywdmcxD01fxTtfb6V5szju\nPLM3V5/QjVQ7DWVC6LgeGVw+PJOnZq/l/IEdGdA5ze9Ih63BW1RF5AU384yJRlv3HuSeN75m9MOz\n+GjFDiac1pPP7jqdCaf3soJhPHHPOf1IT0nkd+/mEW1Xr4K7I43+wRMiEgsM8yaOMeGxt6SCRz/J\n59nP11OtylXHZXHraT1pbZfOGo+1aBbPj0f14ldvLWPm8h2Mym7nd6TDUmfREJF7gF8ASSKyr2Y2\nUA48EYZsxoRcaUUVz32+nokf57O/rJKLhnTiJ6N6k5luN+WZ8Bl3bCbPfLaOP72/glP7tCEuirrG\nr7NoqOofgT+KyB9V9Z4wZjIm5KqrlbcWb+ZvH65i856DnNanDT8/uy9921t3Hyb84mNjuGtMX26e\nvIDXFxRw+fAufkdyrcHTU6p6j4h0ArKCl1fVWV4GMyZUPl+zi9+/u5zcLfsY2DmNv146iON6ZPgd\nyzRxZ/Vvx7CsVjw0fRVjB3eMmiv0GkwpIn8CxgF5QJUzWwErGiairdlZzB/eXc7MFTvo1DKJv48b\nzPkDOxJjY2+bCCAi/OKcvlz82Bc8NXsdt5/Ry+9IrrgpbRcBfVS1zOswxoTCnpJyHpmxmslzN9As\nPpafj+nLNSd0pVl8dN5MZRqvYVnpjM5ux78/W8d1J3aLip4G3CRcC8QDVjRMRKusqubFLzfy0PRV\n7C+tYNzwLvx0dG+7IspEtFtO7cH0vO28Om8T10ZB9+luikYJsFhEZhJUOFT1ds9SGXOY5uTv4rdv\n57JqezHH98jgvvOzrZHbRIWhXVpxbNdW/PuzdVx5XBbxEX4llZt0U4HfAZ8DC4IeDRKRMSKyUkTy\nReTuepa7WERURI66217TtBTsLuHmFxbwg6e+5GBFFZOuHMaL14+wgmGiyo0n92DznoNMW7rV7ygN\ncnP11HMikgR0UdWVblfs3AQ4ERgNFADzRGSqquYdslxz4MfAl4eV3DRppRVVTPp0LY9+kk+MCHee\n2ZvrT+pu7RYmKp3Rty092qQw6dNAL7iRPDSsm25EzgcWA+8704NFZKqLdQ8H8lV1raqWA68AY2tZ\n7nfAn4FS16lNk/bRiu2c+fAsHp6xilHZ7Zj5s1OYcHovKxgmasXECDee3J28rfv4LH+X33Hq5eb0\n1G8IFIA9AKq6GOju4n2dgE1B0wXOvG+IyFAgU1XfdRPWNG2bikq4/rn5XPvsfBLiYnjp+hFMvGIo\nHVsm+R3NmKN24ZBOtGmeyBOz1vodpV5uGsIrVHXvIYdL1Ue7YRGJAR4Crnax7I3AjQBdukTPnZMm\nNMorq3ly9lr++dFqYkS45+y+XHNCNxLiIrvB0JjDkRgXy9XHd+XBD1aSv6OYnm1T/Y5UKzd/dbki\ncgUQKyK9ROSfBBrFG7IZyAya7uzMq9EcOAb4RETWAyOBqbU1hqvqE6qao6o5bdq0cbFp01h8ubaQ\nc/4xmwc/WMmpvdsy46encNMpPaxgmEbp0pzOxMYIry/Y1PDCPnHzl3cbgZ5uy4CXgX3AHS7eNw/o\nJSLdRCSBwF3l37SFqOpeVW2tql1VtSswF7hAVecf5vdgGqHdB8q5a8oSLntiLqUVVTx9dQ6PXznM\nTkWZRq1t82ac1qcNbyzcTGXVUZ/Q8YSbq6dKgHudh2uqWikiE4APgFjgaVXNFZH7gfmq6qYx3TQx\nqoGOBR94Zzl7D1Zwy6k9uP30XlE7NKYxh+vSnExmLN/BrNU7Ob1v5HWbXl/X6I+o6h0i8jaBvqa+\nQ1UvaGjlqjoNmHbIvPvqWPbUBtOaRm1TUQm/eHMps1fvYnBmSyZ/bwD9Otj9FqZpOb1vWzJSEnht\nXkF0FQ2gZnS+v4YjiGm6qqqVZ+as428friJG4LcX9Gf8yCxirWNB0wTFx8Zw4ZBOPP/FeooOlJOe\nkuB3pO+obzyNBc7XT8MXxzQ1q7bv5/+mfM2STXs4vW9bHrjwGGu3ME3e93My+fdn63hr0eaI64+q\nvtNTS6nltFQNVR3oSSLTJFRUVfPYJ2v450erad4snr+PGxzxd8IaEy592jdnYOc0Xpu/iWtO6BpR\nfxf1nZ46z/l6q/O15nTVeOopJsY0JHfLXu58/WuWb93H+YM68pvzs8mwnmiN+Y5LczL51VvLyN2y\nj2M6pfkd5xt1XnKrqhtUdQMwWlXvUtWlzuPnwJnhi2gai4qqah6ZsYqx/5rDruIyJl05jH9ePsQK\nhjG1uGBgR+JihHe+jqxODN3cpyEickLQxPEu32fMN1Zs28eFE+fwyIzVnDewA9N/cjJn9W/vdyxj\nIlZacjwju2cwPW+b31G+w003ItcBT4tIGiDAbuBaT1OZRqOqWpk0aw0PT19FWlI8j48fxphjrFgY\n48aZ/dtx339zWbOzmB5tIqNbETc39y0ABjlFA1Xd63kq0yis23WAn722mIUb93DOgPY8cOGAiLt8\n0JhINqpfoGhMz9tOj1OipGgAiMi5BLoSaVbTiq+q93uYy0QxVeWlrzbywDvLiY8VuzLKmCPUsWUS\nAzql8WHuNm4+pYffcQAXRUNEHgeSgdOAp4BLgK88zmWi1M79Zdz9n6+ZuWIHJ/ZszYOXDqRDmt13\nYcyRGp3djodnrGLH/lLaNm/mdxxXDdrHq+pVwG5V/S1wHNDb21gmGs1cvp0xj8xidv4ufn1+Ns9f\nO9wKhjFH6cz+7VCFmct3+B0FcFc0akbUKxGRjkAF0MG7SCbalFZU8au3lnHdc/Np26IZ79x2Itec\n0I0Y6wbEmKPWp11zMtOT+DA3Mq6ictOm8baItAQeBBYSuLHvSU9TmaixfOs+bn95Eat3FHP9id34\nvzF9SIyzHmmNCRUR4czs9rwwdwPFZZWkJrpqivZMvUcazuh6M1V1j6r+B8gC+tbVU61pOlSV579Y\nz9iJc9hzsILnrx3OL8/LtoJhjAdGZ7ejvLKaWat2+h2l/qKhqtXAxKDpMrvk1uw+UM6NLyzgvv/m\nckKPDN7/8Umc3NtGVDTGKzlZrWiVHM+MvO1+R3F1emqmiFwMvKGq1udUEzd/fRG3vbyIXcVl/Oq8\nbK6NsM7UjGmM4mJjOKFnaz5fU4iq+vo356Yh/CbgdaBMRPaJyH4R2edxLhNhqquVRz/J57In5pIQ\nF8Mbt5zAdSd2s4JhTJiM7J7Btn2lbCgs8TWHmzvCm4cjiIlchcVl/PS1JXy6aifnDezAH783gObN\n4v2OZUyTMrJ7OgBz1xbStXWKbznc3Nw3tJbZe4ENqloZ+kgmksxfX8SElxZRVFLO7y86hiuGd7Gj\nC2N80KNNKq1TE/hyXRHjhnfxLYebNo1HgaHAUmd6ALAMSBORW1T1Q6/CGf+oKv/+bB1/em8FnVol\n8cYtx0dUn/7GNDUiwohuGcxd62+7hps2jS3AEFUdpqrDgMHAWmA08Bcvwxl/7C+t4JbJC3ng3eWc\n0a8tb992ohUMYyLAyO7pbN1byqaig75lcHOk0VtVc2smVDVPRPqq6lo7TdH4rNq+n5tfWMCGohLu\nPacf159kjd3GRIqR3TOAQLtGl4xkXzK4OdLIFZHHROQU5/EokCciiQS6FDGNxNtLtnDhxDnsK63g\nxetHcMPJ3a1gGBNBerZNJSMlgblrC33L4OZI42rgR8AdzvQc4E4CBeM0b2KZcKqsquYvH6zkiVlr\nGZbViolXDKV9mv+9aRpjvktEGNE9nS/XFfnWruHmktuDwN+cx6GKQ57IhNXuA+Xc9vIiPsvfxfiR\nXbjvvP4kxNlovsZEqhHdMpi2dBsFuw+SmR7+U1T+9nxlfLV86z5ueH4+O/aV8eeLB3DZsf5dxmeM\ncSe4XcOPomEfKZuo95dt5eLHPqeiqppXbxppBcOYKNGrbSqtkuOZu7bIl+3bkUYTU12t/OOj1Twy\nYzWDM1vyxJXDaNvC2i+MiRYxMd/er+GHOouGiLxNYOyMWqnqBZ4kMp45WF7Fna8v4d2lW/ne0E78\n4aIBNIu3rsyNiTY5XVvxfu42X4aAre9I46/O1+8B7YHJzvTlgP/985rDsm1vKTc8P59lW/Zyz9l9\nudEupzUmavXvGLjZNnfLPtr2iZCioaqfAojI31Q1J+ilt0VkvufJTMgsLdjL9c/Po7i0kievzGFU\ndju/IxljjkJ2xxYA5G3Zx2l92oZ1224awlNEpHvNhIh0A/zrYtEclg9yt3HppM+Ji4lhyi3HW8Ew\nphFIS4onMz2J3C3hHxPPTUP4T4BPRGQtIASGfL3J01TmqNV0OPj7acsZ2LklT12VQ5vmiX7HMsaE\nSP8OaeRtCf/QRm5u7ntfRHoBfZ1ZK1S1zNtY5mhUVlVz/zt5PP/FBsb0b8/Dlw0mKcEavI1pTPp3\nbMH7udvYX1oR1vFtGjw9JSLJwP8BE1R1CdBFRM7zPJk5IiXlldw8eQHPf7GBG0/uzqM/GGoFw5hG\nqH+nQLvG8q37w7pdN20azwDlwHHO9GbgATcrF5ExIrJSRPJF5O5aXr9ZRJaKyGIR+UxEsl0nN/9j\nV3EZlz/5JTNX7OD+sf35xTn9iImxK6SMaYy+vYIqvO0abopGD1X9C06PtqpaQqBto14iEgtMBM4G\nsoHLaykKL6nqAFUdTGBsjocOJ7z51vpdB7j4sc9ZsXUfj48fxlXHdfU7kjHGQ22bJ9I6NYHcMLdr\nuGkILxeRJJwb/USkB+CmTWM4kK+qa533vQKMBfJqFlDV4O82hXpuJjR1W7Z5L1c/8xVV1crLN45k\naJdWfkcyxnhMRMjumBb2ouHmSOPXwPtApoi8CMwE7nLxvk7ApqDpAmfed4jIrSKyhsCRxu0u1muC\nfLZ6F5dN+oLEuFim3HK8FQxjmpD+HVuwevt+yiqrwrbNeouGBG4ZXkHgrvCrgZeBHFX9JFQBVHWi\nqvYAfg78so4cN4rIfBGZv3PnzlBtOuq9vWQL1zz7FZnpybzxo+Pp0SbV70jGmDDq37EFldXK6u3h\nG6Wi3qKhqgpMU9VCVX1XVd9R1V0u170ZyAya7uzMq8srwIV15HhCVXNUNadNmzYuN9+4TZ67gdtf\nWcSQzFa8etNxtLNOB41pcvxoDHdzemqhiBx7BOueB/QSkW4ikgCMA6YGL+Dc/1HjXGD1EWynSVFV\nJn6czy/fWsZpfdry/HXDSUsK3zXaxpjIkZWeTGpiXFjbNdw0hI8AfiAiG4ADBK6cUlUdWN+bVLVS\nRCYAHwCxwNOqmisi9wPzVXUqMEFERhG4Mms38MOj+F4aPVXlT++tYNKstYwd3JG/XjqI+FgbEsWY\npiomRujXoXnEFY2zjnTlqjoNmHbIvPuCnv/4SNfd1FRXK/dNXcbkuRu5cmQWv72gv92DYYyhf8c0\nXpu/iapqJTYM/xMa/JiqqhsItE2c7jwvcfM+EzpV1cpd//mayXM3ctMp3bl/rBUMY0xAdocWlJRX\nsaHwQFi256YbkV8TuLLpHmdWPN+OrWE8VlFVzY9fWcSUBQX8ZFRv7h7T18bBMMZ8o6ab9HB1J+Lm\n9NRFwBBgIYCqbhGR5p6mMgCUV1Yz4aWFfJi3nV+c05cbT+7hdyRjTITp2jowUsX6MB1puLojXFVV\nRGruCLexNMKgrLKKW19cxIzl2/nN+dlcfUI3vyMZYyJQamIcrVMT2FhYEpbtuWmbeE1EJgEtReQG\nYAbwpLexmrayyipumbyQGcu387ux/a1gGGPqlZWRwoaiCDnSUNW/ishoYB/QB7hPVad7nqyJqikY\nH63Ywe8vOoYfjMjyO5IxJsJlpSczd21hWLbVYNEQkZ8Cr1qh8F55ZTW3vrjICoYx5rB0yUjmzcWb\nKa2oolm8t+PnuDk91Rz4UERmi8gEEbFBpj1Qc5XUjOXbuX9sfysYxhjXumakoAoFu71v13Bzn8Zv\nVbU/cCvQAfhURGZ4nqwJqapWfvLqYt5bto37zsu2sTCMMYelS0YyABvC0Bh+ODfp7QC2AYVAW2/i\nND3V1crd//mad77eyj1n9+XaE63R2xhzeLLSA0VjfSQUDRH5kYh8QmAcjQzghob6nTLuqCr3v5PH\n6wsK+PEZvbjpFLsPwxhz+NJTEkhNjGNjGO7VcHOfRiZwh6ou9jpMU/PXD1fy7Ofruf7EbtwxqlfD\nbzDGmFqICFkZyWwoioAjDVW9B1CnEXyCiAzyPFUT8MSsNUz8eA2XD8/k3nP7WdcgxpijkpWRHJYb\n/NycnrodeJFAO0ZbYLKI3OZ1sMbs9fmb+MO0FZw7sAMPXDjACoYx5qh1SU9h0+4SqqrV0+24OT11\nPTBCVQ8AiMifgS+Af3oZrLGanredu99Yykm9WvPw9weHpStjY0zj1zUjmYoqZcueg2Q6DeNecHP1\nlADBo5ZXOfPMYfpqXRETXlrIMR1b8Nj4YSTEWQ/zxpjQqLnsdqPH7RpujjSeAb4UkTed6QuBf3sX\nqXFavX0/1z83j04tk3j66mNJTXSz640xxp2sjG97uz2hZ2vPtuOm76mHnEtuT3RmXaOqizxL1Aht\n31fK1c/MIyEulueuHU5GaqLfkYwxjUyHFs1IiIvxvDHcTd9TI4FcVV3oTLcQkRGq+qWnyRqJ4rJK\nrnlmHrtLynn1xuM8PddojGm6YmKEzFZJnt8V7uak+mNAcdB0sTPPNKCyqpofvbiQldv3M/EHQxnQ\nOc3vSMaYRiwrI8XzwZhcNYSr6jfXcKlqNe7aQpo0VeU3b+cya9VOfn/hMZzWx3peMcZ4q0t6MhuL\nSgj6lx1yborGWhG5XUTincePgbWeJWoknpmznslzN3LTyd0ZN7yL33GMMU1A14xkSsqr2FVc7tk2\n3BSNm4Hjgc1AATACuNGzRI3AzOXbeeDdPM7MbsfPx/T1O44xpomouYJqo4ej+Lm5emoHMM6zBI3M\nim37uP3lRWR3bMEj4wYTYzfvGWPCpOZejfW7ShiWle7JNtx0I9JbRGaKyDJneqCI/NKTNFFu94Fy\nbnh+PimJcTx11bEkJ1jTjzEmfDq3SiJG8LTjQjenp54E7gEqAFT1a+zI439UVlVz60sL2b63jElX\nDqN9WjO/IxljmpjEuFg6pCVR4HPRSFbVrw6ZV+lFmGj2+2nL+XxNIX/43gCGdGnldxxjTBPVunki\nO4vLPFu/m6KxS0R6AAogIpcAWz1LFIWmLCjgmTnrufaEblwyrLPfcYwxTVhGSgJFB7y7esrNSfdb\ngSeAviKyGVgH/MCzRFEmd8te7n1zKcf3yOAX59iVUsYYf6WnJJC3ZZ9n66+3aIhIDJCjqqNEJAWI\nUdX9nqWJMntLKrh58gJaJSfwj8uHEBdrvdYaY/yVkRo40lBVT8bqqfe/nHP3913O8wNWML5VXa38\n5LXFbNtbyqPjh9LaOiE0xkSAjJQEyquqKS7zpunZzUfjGSJyp4hkikh6zcOTNFHk0U/y+WjFDn51\nXjZDreHbGBMhMlICH2ALPbqEo/roAAAUIklEQVQr3E2bxmXO11uD5inQPfRxosPctYU8NH0VYwd3\n5MqRWX7HMcaYb6SnJgBQeKCcrq1TQr5+N3eEdwv5VqNY0YFyfvzKIrIyUvj9RTa+tzEmsrT+5kjD\nm8tureX2MKgqd76+hN0HKvjXFUNs9D1jTMSpOdLw6rJbT4uGiIwRkZUiki8id9fy+k9FJE9Evna6\nKonocz3//mwdH63Ywb3n9qN/RxsbwxgTeTJSvj095QXPioaIxAITgbOBbOByEck+ZLFFBC7pHQhM\nAf7iVZ6jtbRgL39+fwVn9W/HVcdFdG0zxjRhzeJjSUmI9bUhHBHpBGQFL6+qsxp423AgX1XXOut4\nBRgL5AWt4+Og5ecC493FDq+D5VXc8eoiMlIS+fPFA60dwxgT0dJTEyg64E2bhpsxwv9M4AqqPKDK\nma1AQ0WjE7ApaLpmLI66XAe8V0eGG3HG8OjSJfwDGv3pveWs2XmAydeNoGVyQti3b4wxhyMjJdGz\n01NujjQuBPqoqmc9YInIeCAHOKW211X1CQJdmZCTk+PdOIa1+HTVTp77YgPXnNCVE3u1DuemjTHm\niGSkJLB1b6kn63Y13CsQfwTr3gxkBk13duZ9h4iMAu4FLvCyMB2J3QfK+b/Xl9CrbaqNwGeMiRoZ\nqQkU+nV6CigBFovITOCbFKp6ewPvmwf0EpFuBIrFOOCK4AVEZAgwCRjjjBAYUe6bmsvuknKevvpY\nmsXH+h3HGGNcSU9J9Kz/KTdFY6rzOCyqWikiE4APgFjgaVXNFZH7gfmqOhV4EEgFXne+sY2qesHh\nbssLH+Ru4+0lW/jp6N4c08kurzXGRI+MlAQqqpR9pZWkJR3JiaK6ubkj/DkRSQB6O7NWqmqFm5Wr\n6jRg2iHz7gt6PuowsobN3pIKfvnWMvp1aMEtp/bwO44xxhyWjKAb/EJdNNyMEX4qsJrAPRePAqtE\n5OSQpogwv3s3j6ID5Tx4yUDirbtzY0yUSU+pKRqhb9dwc3rqb8CZqroSQER6Ay8Dw0KeJgJ8vHIH\nUxYUMOG0nnZayhgTlWqGatjlwQ1+bj5Gx9cUDABVXcWRXU0V8UrKK/nlm8vo2TaV287o6XccY4w5\nIt8eaYS+aLg50pgvIk8Bk53pHwDzQ54kAvx95mo27znIlJuPIzHOrpYyxkSnmqLhRU+3borGLQTG\n0qi5xHY2gbaNRmXltv38e/Y6vp/TmZyuTX6MKWNMFGsWH0tqYpwnd4W7uXqqDHjIeTRK1dXKL99a\nSmqzOO4+u5/fcYwx5qilpyR40mmhXRoETFlYwLz1u7nn7L7fHNYZY0w0y0hN8KRNo8kXjb0lFfzp\nvRUMy2rFpcMyG36DMcZEgYyUBE9OT9VbNEQkVkT+GvKtRpC/z1zNnpJyfjf2GGJirMtzY0zjkJGS\n6ElDeL1FQ1WrgBNDvtUIsXZnMc9/sZ7Ljs0ku2MLv+MYY0zIpDunp1RD2zG4m6unFonIVOB14EDN\nTFV9I6RJfPCHaStoFh/LT0f38TuKMcaEVEZKApXVyr6DlaQlh+7WOjdFoxlQCJweNE+BqC4an+fv\nYsby7dw1pg9tmif6HccYY0Kqpv+pwgNl4S0aqnpNyLYWIaqqlfvfyaNzqySuPaGb33GMMSbk0lMC\nH4YLD5TTvU3o1uumw8LOIvKmiOxwHv8Rkc6hixB+/1lQwIpt+7n77L42ToYxplHK+Oau8NBeQeXm\nkttnCIyn0dF5vO3Mi0pllVX8feZqBnVO49wBHfyOY4wxngjuHj2U3BSNNqr6jKpWOo9ngRAe7ITX\nq/M2sXnPQe48q0/IR7QyxphI4VX/U26KRqGIjHfu2YgVkfEEGsajzsHyKv75UT7Du6VzYs/Wfscx\nxhjPJMbF0tyD/qfcFI1rge8D24CtwCVAVDaOT567gZ37y/jZ6N52lGGMafQyUkN/V3i9V0+JSCzw\nvUgZt/toFJdV8tinazipV2tGdM/wO44xxnguPSUh5KP3ubkj/PKQbtEnz85ZR9GBcu48027kM8Y0\nDekpiSG/esrNzX1zRORfwKt8947whSFN4qGS8kqe+mwdo/q1ZVBmS7/jGGNMWLROTWBJwZ6QrtNN\n0RjsfL0/aJ7y3TvEI9rr8wvYU1LBLaf28DuKMcaETWZ6MhkpCSHtf6qhNo0Y4DFVfS1kWwyzyqpq\nnpy9lpysVgzLshH5jDFNx62n9eTW03qGdJ0NtWlUA3eFdIthNm3ZNgp2H+SmU+wowxhjjpabS25n\niMidIpIpIuk1D8+ThYCqMunTNfRok8IZfdv6HccYY6KemzaNy5yvtwbNU6B76OOE1pz8QnK37OPP\nFw+wAZaMMSYE3PRyG7XdwE6atYY2zRO5cEgnv6MYY0yjUOfpKRG5K+j5pYe89gcvQ4XC6u37mb16\nF1cf35XEOOvJ1hhjQqG+No1xQc/vOeS1MR5kCamXv9pEfKww7thMv6MYY0yjUV/RkDqe1zYdUUor\nqnhjUQFn9m9PRqqNymeMMaFSX9HQOp7XNh1RPsjdxp6SCq4Y3sXvKMYY06jU1xA+SET2ETiqSHKe\n40w38zzZUXjpy41kZSRznHVMaIwxIVVn0VDVqGw9XrOzmC/XFXHXmD52ma0xxoSYm5v7osorX20k\nLka4ZFhUD2NujDERydOiISJjRGSliOSLyN21vH6yiCwUkUoRueRot1dWWcWUBQWMzm5H2+YRfQbN\nGGOikmdFwxnAaSJwNpANXC4i2YcsthG4GngpFNuckbeD3SUVjLMGcGOM8YSbbkSO1HAgX1XXAojI\nK8BYIK9mAVVd77xWHYoNTlu6ldapCTb+tzHGeMTL01OdgE1B0wXOPE+UVlTx8codnNm/PbHWAG6M\nMZ6IioZwEblRROaLyPydO3fWusysVTspKa9iTP/2YU5njDFNh5dFYzMQ3IdHZ2feYVPVJ1Q1R1Vz\n2rRpU+sy7+duIy0pnuN62L0ZxhjjFS+Lxjygl4h0E5EEAn1ZTfViQ+WV1czI286ofu2Ij42Kgydj\njIlKnv2HVdVKYALwAbAceE1Vc0XkfhG5AEBEjhWRAuBSYJKI5B7Jtr5YW8i+0krGHGOnpowxxkte\nXj2Fqk4Dph0y776g5/MInLY6Ku8v20ZKQiwn9bKrpowxxktRfy6nqlqZnreN0/q2pVl8VPZ8Yowx\nUSPqi8b89UXsKi63U1PGGBMGUV803s/dRmJcDKf1aet3FGOMafSivmh8saaQ4d3SSUn0tHnGGGMM\nUV409pdWsHL7foZ2aeV3FGOMaRKiumgs2bQXVRiaZUXDGGPCIaqLxqKNuwEYnNnS5yTGGNM0RHXR\nWLhxN73appKWFO93FGOMaRKitmioKos27bH2DGOMCaOoLRprdx1gT0kFQ7Ps1JQxxoRL1BaNhRsC\n7Rl2pGGMMeETvUVj4x6aN4ujR5tUv6MYY0yTEbVFY9HG3QzObEmMjdJnjDFhE5VFw27qM8YYf0Rl\n0fi6wG7qM8YYP0Rl0ahpBLeb+owxJryis2jYTX3GGOOLqCwaizbtYUgXO8owxphwi7qiUa3K6X3b\n2vgZxhjjg6gbhCJGhIe+P9jvGMYY0yRF3ZGGMcYY/1jRMMYY45oVDWOMMa5Z0TDGGOOaFQ1jjDGu\nWdEwxhjjmhUNY4wxrlnRMMYY45qoqt8ZDouI7AdW+p3DhdbALr9DuGA5QycaMoLlDLVoydlHVZsf\n7Uqi7o5wYKWq5vgdoiEiMt9yhk405IyGjGA5Qy2acoZiPXZ6yhhjjGtWNIwxxrgWjUXjCb8DuGQ5\nQysackZDRrCcodakckZdQ7gxxhj/ROORhjHGGJ9EZNEQkUtFJFdEqkWkzqsSRGSMiKwUkXwRuTto\nfjcR+dKZ/6qIJHiUM11EpovIaudrq1qWOU1EFgc9SkXkQue1Z0VkXdBrngwU4ians1xVUJapQfM9\n358u9+VgEfnC+d34WkQuC3rN031Z1+9a0OuJzr7Jd/ZV16DX7nHmrxSRs0KZ6why/lRE8pz9N1NE\nsoJeq/Xn71POq0VkZ1Ce64Ne+6Hze7JaRH7oc86HgzKuEpE9Qa+FZX+KyNMiskNEltXxuojIP5zv\n4WsRGRr02uHvS1WNuAfQD+gDfALk1LFMLLAG6A4kAEuAbOe114BxzvPHgVs8yvkX4G7n+d3AnxtY\nPh0oApKd6WeBS8KwP13lBIrrmO/5/nSTEegN9HKedwS2Ai293pf1/a4FLfMj4HHn+TjgVed5trN8\nItDNWU+sjzlPC/r9u6UmZ30/f59yXg38q5b3pgNrna+tnOet/Mp5yPK3AU/7sD9PBoYCy+p4/Rzg\nPUCAkcCXR7MvI/JIQ1WXq2pDN/ANB/JVda2qlgOvAGNFRIDTgSnOcs8BF3oUdayzfrfbuQR4T1VL\nPMpTl8PN+Y0w7s8GM6rqKlVd7TzfAuwA2niQ5VC1/q4dskxw/inAGc6+Gwu8oqplqroOyHfW50tO\nVf046PdvLtDZoyz1cbM/63IWMF1Vi1R1NzAdGBMhOS8HXvYoS51UdRaBD6N1GQs8rwFzgZYi0oEj\n3JcRWTRc6gRsCpoucOZlAHtUtfKQ+V5op6pbnefbgHYNLD+O//2l+r1zyPiwiCSGPGGA25zNRGS+\niMytOYVG+PbnYe1LERlO4NPfmqDZXu3Lun7Xal3G2Vd7Cew7N+8NZ85g1xH4BFqjtp+/F9zmvNj5\neU4RkczDfG8ouN6Wc5qvG/BR0Oxw7c+G1PV9HNG+9O2OcBGZAbSv5aV7VfW/4c5Tl/pyBk+oqopI\nnZeiOZV9APBB0Ox7CPyDTCBwOdzPgft9zJmlqptFpDvwkYgsJfDPLyRCvC9fAH6oqtXO7JDty6ZA\nRMYDOcApQbP/5+evqmtqX4Pn3gZeVtUyEbmJwFHc6T5lcWMcMEVVq4LmRdL+DBnfioaqjjrKVWwG\nMoOmOzvzCgkcfsU5n/hq5h+R+nKKyHYR6aCqW51/ZDvqWdX3gTdVtSJo3TWfrMtE5BngTj9zqupm\n5+taEfkEGAL8hxDtz1BkFJEWwLsEPlzMDVp3yPZlLer6XattmQIRiQPSCPwuunlvOHMiIqMIFOpT\nVLWsZn4dP38v/sk1mFNVC4MmnyLQ5lXz3lMPee8nIU/47bbc/uzGAbcGzwjj/mxIXd/HEe3LaD49\nNQ/oJYErexII/NCmaqCF52MC7QcAPwS8OnKZ6qzfzXb+53yn88+xpt3gQqDWqx9CoMGcItKq5pSO\niLQGTgDywrg/3WRMAN4kcH52yiGvebkva/1dqyf/JcBHzr6bCoyTwNVV3YBewFchzHZYOUVkCDAJ\nuEBVdwTNr/Xn72PODkGTFwDLnecfAGc6eVsBZ/Ldo/ew5nSy9iXQkPxF0Lxw7s+GTAWucq6iGgns\ndT5kHdm+DEfr/uE+gIsInF8rA7YDHzjzOwLTgpY7B1hFoHrfGzS/O4E/zHzgdSDRo5wZwExgNTAD\nSHfm5wBPBS3XlUBVjznk/R8BSwn8g5sMpPqVEzjeybLE+XpdOPeny4zjgQpgcdBjcDj2ZW2/awRO\nf13gPG/m7Jt8Z191D3rvvc77VgJne/y301DOGc7fVM3+m9rQz9+nnH8Ecp08HwN9g957rbOf84Fr\n/MzpTP8G+NMh7wvb/iTwYXSr87dRQKCt6mbgZud1ASY638NSgq5IPZJ9aXeEG2OMcS2aT08ZY4wJ\nMysaxhhjXLOiYYwxxjUrGsYYY1yzomGMMcY1KxrGGGNcs6JhjDHGNSsaJuKIiIrI34Km7xSR3zTw\nns89DxYiIlLcwOstReRHh8yLmu/PNG5WNEwkKgO+53S/4IqqHu9hnno53TPE1DV9BFoSGJ/jG35+\nf8YEs6JhIlElgZ5qf3LoCxIYeW6Z87gjaH6x8zVFRN4VkSXOMpc588eLyFcSGEVtkojE1rLuq5yu\nuJeIyAv1bVNEukpgRLfnCXRdctIh05kut/mWiCyQwGiENzqz/wT0cN73YPD310Ce5SLypLOuD0Uk\nqa4dLCI/ct6/QURuq/tHYcwhvOy3xR72OJIHUAy0ANYT6C32TgL9+wwj0HdOCpBKoG+iITXvcb5e\nDDwZtK40AiNBvg3EO/MeBa46ZJv9CfQx1NqZrun7qtZtEuhPrBoY6Sx36HSd2yRoRLeg7SQRKDYZ\nzrqWHbpPXOSp5Nu+uF4Dxtexfy8m0PFjPNAB2AXE+f1zt0d0POxIw0QkVd0HPA/cHjT7RALdyx9Q\n1WLgDeCkQ966FBgtIn8WkZNUdS9wBoF/tvNEZLEz3f2Q950OvK6qu5zt14yEVt82N2hQ9+yHTLvZ\nJsDtIrKEwCh6mQR6wa1PfXnWqepi5/kCAoWkNrcDP1fVCg30dlqBnXUwLvk2noYxLjwCLASecfsG\nVV0lIkMJ9E76gIjMBHYDz6nqPSHOd6CeaWlomyJyKjAKOE5VSyQw5kKzo8hTFvS8isDRy6HbjAcG\nqeoqZ7oDUKiB4UyNaZB9ujARy/m0/xqBrp4BZgMXikiyiKQQ6EJ/dvB7RKQjUKKqk4EHgaEEuly/\nRETaOsukS2B4zmAfAZeKSEbNMm63WQc320wDdjsFoy8w0pm/H2hex3qPNE+NbKCFiHR3Guv/CPzj\nMN5vmjg70jCR7m/ABABVXSgiz/LtIEZPqeqiQ5YfADwoItUETrvcoqp5IvJL4EPnH2UFgVHWNtS8\nSVVzReT3wKciUgUsAq6ua5si0rW+0G62CbwP3CwiywmMtTHXeW+hiMwRkWXAe6r6f0HrPaI8QYYA\nLxIYgyEFeENVn3D5XmNsPA1jmhIReQT4QlVf9TuLiU52esqYpmUwgRH7jDkidqRhjDHGNTvSMMYY\n45oVDWOMMa5Z0TDGGOOaFQ1jjDGuWdEwxhjjmhUNY4wxrlnRMMYY45oVDWOMMa79PxPUqYdthuUe\nAAAAAElFTkSuQmCC\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x10c14d0b8>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plot(Rho, results)\n",
"xlabel(r\"Noise correlation $\\rho$\")\n",
"ylabel(r\"Error on recovered gradient\")\n",
"xlim(-1,1)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"# Another version - fixed rho, varying sigma_y\n",
"\n",
"# Mock data - just two points\n",
"x = np.array([1.0, 2.0])\n",
"y = np.array([1.0, 2.0])\n",
"\n",
"# Fitting function - just fit the gradient of a straight line\n",
"f = lambda x, p0: p0*x\n",
"\n",
"# Standard deviations of individual data points\n",
"sigma_x = 1.0\n",
"rho = 0.8\n",
"\n",
"# We will vary the covariance between the data points\n",
"def make_cov(sigma_y):\n",
" return np.array([[sigma_x**2, rho*sigma_x*sigma_y], [rho*sigma_x*sigma_y, sigma_y**2]])\n",
"\n",
"# Loop through possible correlation values\n",
"#Rho = np.linspace(-0.99, 0.99, 101)\n",
"Sigma_y = np.linspace(0.1, 10.0, 1000)\n",
"results = []\n",
"for sigma_y in Sigma_y:\n",
" C = make_cov(sigma_y)\n",
" mu, c = scipy.optimize.curve_fit(f, x, y, p0=[0.5], sigma=C, absolute_sigma=True)\n",
" results.append(c[0,0]**0.5)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(0, 10)"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEOCAYAAACTqoDjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3XeYXXW56PHvO73PZGrqZFImPRBg\nIAkkIQmhiQIqIipHD0dFjxSRo1iORz14vNd+vV6VIyBF8CCIIgQj1UAKkB7Sy6RPybRM7+W9f+w1\nYTLszKxJsmft8n6eZz97r7V/a613djL7nfWroqoYY4wx/UV5HYAxxpjgZAnCGGOMX5YgjDHG+GUJ\nwhhjjF+WIIwxxvhlCcIYY4xfliCMMcb4ZQnCGGOMX5YgjDHG+BXjdQBDlZ2drQUFBV6HYYwxIWXT\npk3VqpozlGNCLkEUFBSwceNGr8MwxpiQIiJHhnqMVTEZY4zxyxKEMcYYvyxBGGOM8csShDHGGL8s\nQRhjjPErYAlCRB4RkUoR2XGa90VEfikixSKyTUQuDFQsxhhjhi6QdxCPAdcM8P61QKHzuB14IICx\nGGOMGaKAjYNQ1VUiUjBAkRuA36tvzdN3RCRDREapanmgYopkZXWtrNpXxbbSeoormqhpbqepvYvY\n6CiS42IYnZHAuMwkZo5O48L8EUzKSSEqSrwO2xjjIS8Hyo0BjvXZLnH2vS9BiMjt+O4yyM/PH5bg\nwkFXdw8v7TzOo2sPs+lILQCpCTFMH5nGtJFppMTH0NnTQ2NbF6W1rWw8XMvv3/aNpRmRFMuSqblc\nOSOPRVNySI4PuTGVxpizFBK/9ar6IPAgQFFRkXocTkhYs7+a7y3fSXFlExOyk/na1VO5emYek3JS\nEPF/Z6CqHKxuZtORWt45UMM/9lbyly2lJMRGce2sUdxcNI65EzLtzsKYCOFlgigFxvXZHuvsM2eh\nrbOb/1y+i6fWHyU/M4kHPnUhV88c6epLXUSYlJPCpJwUbi4aR1d3DxsO17J8WxnLt5bx3JZS8jOT\nuHVePh+/OJ/0xNhh+ImMMV4RXxNAgE7ua4N4UVVn+XnvOuBO4APAXOCXqnrJYOcsKipSm4vJv/L6\nVm57dAN7jjfyxcsncc+yQhJio8/JuVs7unl553H+Z/1R1h86QVJcNDcXjeO2ywoYn5V8Tq5hjAkc\nEdmkqkVDOiZQCUJEngIWA9lABfBdIBZAVf9bfPUcv8LX06kFuE1VB/3mtwThX3FlI5/+3Xoa2rr4\n1ScvYPHU3IBda0dpPY+sOcTybWV09SjXzR7FXUsLmToyNWDXNMacnaBKEIFiCeL9jtQ0c9N/v40q\nPHbbxcwakz4s161oaOPRtYd54u3DNHd0c+2skdy5dDIzRw/P9Y0x7lmCiEBVje189IG3aGjr5E9f\nmE9h3vD/FV/b3MGjaw/x6NrDNLZ3sWx6HndfMZnzxmYMeyzGGP8sQUSYzu4ePvXQOraV1vHH2+cz\nZ5y3X8j1rZ08tvYwj6w9RH1rJ1dMy+UrV04ZtjsaY8zpnUmCsLmYQtgP/rab9YdP8KOPnud5cgBI\nT4zly8sKWfP1JXz1qilsPFLLB//fGj7/+43sLKv3OjxjzBBZgghRK/dU8thbh7ntsgJumDPG63BO\nkZoQy51LC1n99SV8ZdkU3jlYw3W/XMMXntjI7vIGr8MzxrhkVUwhqK6lg6v+zyoykmJZftcC4mPO\nTVfWQKlv7eSRNYd4ZM0hGtu7uHbWSL68rJBpI9O8Ds2YiGFVTBHiP5fv4kRzBz+/eU7QJwfwVT19\n5coprPn6Uu5eOpnV+6u55herueMPm9lX0eh1eMaY07AEEWLeOlDNc1tK+dLiSSHX+JueFMu9V01l\nzdeXcMeSSbyxt5Krf7GKu57aQnGlJQpjgo1VMYWQru4ervvlGpo7unjt3svP2Shpr5xo7uCh1Qd5\n/K3DtHZ2c/35o7n7ikIm5aR4HZoxYceqmMLc/6w/yt6KRr593fSQTw4AmclxfP2aaay+bwm3L5rI\nKzsruPLnb/KVp7dyqLrZ6/CMiXh2BxEiGts6WfjjlcwYlcYfPjf3tDOyhrLqpnYeXHWQ3799mI6u\nHm68YAx3Ly2kINvmejLmbNkdRBh7dO1h6lo6+dYHpodlcgDITonnWx+Yzqr7lnDbZRP427Zyrvj5\nm3ztT+9ytKbF6/CMiTiWIEJAfUsnD60+yFUz8kKuYfpM5KYm8B8fnMHq+5bw6fnjef7dMpb+7A2+\n/uw2jp2wRGHMcLEEEQJ+t+YgjW1d3LNsitehDKvctAS++6GZrL5vCbfOG89zW0pZ8tM3+OZftlFS\na4nCmECzBBHkGto6eXTtYa6dNZIZoyNzYFleWgLfu34mb963mE9cks+fN72XKA5UNXkdnjFhyxJE\nkHt6/TEa27v40uLJXofiuVHpiXz/xlm88bXF3Fw0jj9vKmXZz9/kc49vZN3BGkKtw4UxwS4k1qSO\nVJ3dPTyy9hDzJmYye2z4tz24NTojkR98eDb3LJvCE28f5ol3jvDa7grOH5vO5xZO5NpZI4mJtr99\njDlb9lsUxFZsL6e8vo3bF030OpSglJMaz71XTeWtb1zBf904i4a2Lu56aguX/+QNHl59kIa2Tq9D\nNCak2TiIIKWqfOhXa2jr7OGVexYRFRWeXVvPpZ4e5bXdFTy8+hDrD58gMTaaGy8Yw6fnj2f6qMhs\nvzGm15mMg7AqpiC19VgdO0ob+K8bZ1lycCkqSrhq5kiumjmS7SX1PPHOYf6yuYSn1h/l4oIR/NP8\nAq6ZOZK4GLtxNsYNSxBB6n/WHSU5zvcXsBm62WPT+fFN5/OtD0znTxtLeHLdEe5+agvZKfF84pJx\n3HJJPmMyEr0O05igZgkiCNW3drJ8WxkfvmAsKfH2T3Q2MpLi+PyiiXx2wQTe3F/FE28f4Vcri/nV\nymIWTM7m5qJxXDUzLySmTTdmuNm3TxB6bnMJbZ09fGpuvtehhI2oKGHJ1FyWTM3l2IkWnt1UwrOb\nSrjrqS1kJMVy45wxfKxoLDNHW28xY3pZI3WQUVWu/sUqEmOjef7OBV6HE9Z6epS1B6p5ZmMJL+84\nTkd3D7PGpHHThWP54PmjyU6J9zpEY84Za6QOA5uP1rGvookffXS216GEvagoYWFhDgsLc6hr6eD5\nrWU8veEY31u+i+//bTcLJmdzw5zRXDVzpFX1mYhk/+uDzF82l5AQG8V15432OpSIkpEUx2cuLeAz\nlxaw93gjz28t5fmtZdz7zLvEx2xn2Yw8bjh/NIun5lovKBMxBk0QIvKEqv7TYPvM2Wvv6ubFbeVc\nbX+xemrqyFTuu2YaX7t6KpuP1vL81jJe3FbO37aVk54Yy5Uz8rh21kgWFGZb47YJa26+hWb23RCR\naOCiwIQT2VbuqaK+tZMPW9fWoCAiXDQ+k4vGZ/IfH5zBmuJqXthaxss7j/PsphJS4mNYMi2Xa2eN\nZPHUHJLiLKmb8HLa/9Ei8k3gW0CiiDT07gY6gAeHIbaI89yWErJT4lkwOdvrUEw/sdFRJ3tBdXT1\n8NaBal7acZxXdlWw/N0y4mOiuHxKDtfOHsmSqblkJMV5HbIxZ23QXkwi8r9V9ZvDFM+gwrUXU11L\nBxf/4DU+Pb+A//jgDK/DMS51dfew4XAtL+0o56Wdx6loaCdK4KLxI1gyLZel03KZmpcatqsAmtBx\nJr2YXHVzFZExwHj63HGo6qohR3gOhGuCePKdI3z7rzt48a4FEbFqXDjq6VHeLalj5Z5KXt9Tyc4y\n3433mIxEFk/NYem0XC6dlE1inLVbmOEXkG6uIvJD4BZgF9Dt7FbAkwQRrl7cVsbk3BRmRuiiQOEg\nKkq4IH8EF+SP4N6rplLR0MbKPZX8Y08lz20p5Q/rjhIfE8XciVksnJzNZZOzmTYy1ebaMkHLTava\nh4Gpqto+1JOLyDXA/wWigYdV9Yf93s8HHgcynDLfUNUVQ71OqKtqbGf9oRPcubTQqiLCSF5aArdc\nks8tl+TT3tXNuoMn+MeeSlbvr+IHK3YDkJ0Sx6WTsllQmM2CydmMtvmhTBBxkyAOArHAkBKE09vp\n18CVQAmwQUReUNVdfYp9G3hGVR8QkRnACqBgKNcJBy/tPE6PwnWzR3kdigmQ+JhoFk3JYdGUHADK\n61tZW1zDmv1VrCmu4YV3ywCYmJ3MgsJs5k7I4uIJI8hNTfAybBPh3CSIFmCriLxOnyShqncPctwl\nQLGqHgQQkT8CN+Crqjp5GqC3TiUdKHMZd1j5+/ZyJuYkMyUvxetQzDAZlZ7ITReN5aaLxqKq7Kto\nYvX+KtYWV/PsphJ+//YRwJcwLpmQefIxdkSSx5GbSOImQbzgPIZqDHCsz3YJMLdfme8Br4jIXUAy\nsOwMrhPSqpvaeedgDXcsmWzVSxFKRJg6MpWpI1P53MKJdHb3sKO0nvWHTrD+0AlWbC/njxt8v0pj\nMhJPJoui8SOYlJNibRgmYAZNEKr6uIgkAvmquvccX/8TwGOq+jMRmQ88ISKzVLWnbyERuR24HSA/\nP7xmOH3ZqV76gFUvGUdsdNTJxu4vXD6Jnh5lb0XjyYSxen8Vz20pBSA1PoY5+RlcMC6DC/JHMGdc\nBiOSbQyGOTfc9GL6EPBTIA6YICJzgPtV9fpBDi0FxvXZHuvs6+uzwDUAqvq2iCQA2UBl30Kq+iDO\n4LyioqLQmn52EH/ffpwJ2clMG5nqdSgmSEVFCdNHpTF9VBqfubQAVeVgdTObj9Sy9VgdW47W8auV\nxfQ4vxkFWUlOgsnggnEjmDoy1eaPMmfETRXT9/C1J7wBoKpbRWSii+M2AIUiMgFfYrgF+GS/MkeB\nK4DHRGQ6kABUuYo8DJxo7uDtgzV88fKJVr1kXBMRJuWkMCknhY8V+f4Ga27vYntpPVuO1rHlaC1r\niqtP3mXERUcxbVQqM0enM2tMGrPHpDMlL5WEWBuPYQbmJkF0qmp9vy+wntMV7qWqXSJyJ/Ayvi6s\nj6jqThG5H9ioqi8A/wY8JCJfwddg/c8aagtUnIXXd1fQ3aNcO8uql8zZSY6PYd7ELOZNzAJ864qU\n1rWy5WgdO0rr2VFWz9+2lfHU+qMAxEQJhXmpzB6Txqwx6cwcnc6MUWk2iM+cwk2C2CkinwSiRaQQ\nuBt4y83JnTENK/rt+06f17uAy9yHG15e213BqPQEGxxnzjkRYeyIJMaOSOJD5/umjldVSmpb2V5a\n7ySNBl7bXckzG0sAiBIoyEpm6shUpo1Mc55Tyc9MsobwCOUmQdwF/Du+Lq5P4bsj+H4gg4oEbZ3d\nrNpXzUcvGmPVS2ZYiAjjMpMYl5l0slOEqlJe38b20np2ljWw93gDu8obeGnncXrv5RNjo5mSl3Iy\ncUxzelxl2Yp7Yc9NL6YWfAni3wMfTuR4+0ANrZ3dLJue53UoJoKJCKMzEhmdkcjVM0ee3N/S0cW+\niib2Hm9gz/FG9h5vPOVuAyAzOY7JOSlMyk32tYnkpjA5J4UxGYl2xxEmBpru+xeqeo+ILMfXPnAK\nF72YzABe3V1Bclw08ydleR2KMe+TFBfDnHEZzBmXcXKfqlLV1M5eJ2EUVzZxoKqJl3Ycp7al82S5\n+JgoJuakMDk3hUk5yc5zChOyk61hPMQMdAfxhPP80+EIJJKoKq/vrmDRlBxbkcyEDBEhNzWB3NQE\nFhbmnPLeieYODlQ1+ZJGZRPFVU1sPVbLi9vK6NvtZFR6AuOzkijISmZ8VjIFWUmMz0pmfFYSybaK\nYtA57b+Iqm5ynt8cvnAiw47SBioa2q16yYSNzOQ4MpMzubgg85T9rR3dHKpupriqicPVzRyuaeZI\nTQuv7a6guqnjlLI5qfEnE0ZBVhIF2cmMz0xm7IhEMpJira3OAwNVMW3HT9VSL1U9LyARRYBXd1cQ\nJbBkWq7XoRgTUIlx0cwYncYMPz31Gts6OVLTwpGaFidxNHO4uoVV+6p4tvHUuUGT46KdXlmJziPp\nlGdLIIEx0D3dB53nO5zn3iqnWxkgcZjBvbargqLxmWTalAgmgqUmxDJrTLrfBbJaOrpOJo/SulZK\nalsoqW2lpLaV9YdO0NjedUp5fwlkdEYiI9MTGJWeQG5qPDHRNpp8qAaqYjoCICJXquoFfd76uohs\nBr4R6ODCUWldK7vKG/jWB6Z5HYoxQSspLubk9CL+1Ld2npI0BksgUQK5qQknE8ao9ERGpb+3PTI9\ngby0BGItiZzCTauQiMhlqrrW2bgUsE/xDK3a55tJZMlUq14y5kylJ8aSnugbAe5PfUsn5Q2tlNe3\nUV7XxvF63+vjDW3sq2jkzX1VtHR0n3KMCOSkxJ+SMHJT48lNTSAnLf7k66zkuIjpxusmQXwWeERE\n0gEBaoF/CWhUYWzN/mpGpiUwOdfWfjAmUNKTYklPimXaSP93IKpKY3sXx+vbKKtr5Xh9my+B1LdR\nVt/Kwapm3j5QQ0Nb1/uOjY4SslPiyHEShi9xxJOT9t7r3LQEclLiQ36SRDcD5TYB5zsJAlWtD3hU\nYaq7R1lTXM1VM/KsQc0YD4kIaQmxpCXEMiXv9DMpt3V2U9XYTmVjO1WNbVQ2tlPZ0E6l8/p4fRvb\nSuqpaW7H3yxyGUmxZCXHkZUS7zzHkZUcT3aKb19mcpzvdXI86YmxQXdn4qrjsYhcB8wEEnq/2FT1\n/gDGFZa2l9ZT39rJwik5gxc2xnguITb65PQkA+nq7qGmuYPKhnaqmtqcJNJOdVM7NU0dVDe1s7+y\niXcOtlPX2uk3mURHCZnJcackkqyUOLKdRJKZHMeIpDgyk2PJSIojIzE24A3vbtaD+G8gCVgCPAzc\nBKwPaFRhas1+X/vDZTZ62piwEhMdRV6ar93Ct3ry6XV191Db0klN83vJ40RzBzVNHdQ0t1Pd1EFN\nUzvv1tZxoqnjfQ3ufaUlxDAiOY6MpDgyk2IZkRTHiOQ4RiT5kkhmchwZSbFn3GPSzR3Epap6nohs\nU9X/FJGfAX8/o6tFuFX7q5k1Js0mOTMmgsVER5GTGk9OqrvvgbbObmqaO6ht7qC2pYPals6Tr+ta\nOjnhvK5qamdfRRN1LR0092uAP+NY3cTnPLeIyGigBrAFDIaoqb2LzUdq+fwiN2stGWOMT0JsNGMy\nEhmTkej6mPaubupaOqlt6eBEsy+RXPejoV/bTYJYLiIZwE+AzfgGyT009EtFtnUHa+jqURZOzvY6\nFGNMmIuPiSYvLdqp8jpzAyYIEYkCXlfVOuDPIvIikGA9mYZu9f5qEmKjuKhghNehGGOMKwM2gatq\nD/DrPtvtlhzOzKr9VcybmGWztxpjQoabPlKvi8hHxTrun7HSOt/AmwVWvWSMCSFuEsQXgD8B7SLS\nICKNItIQ4LjCSm/31kU2/sEYE0LcjKQ+/TBD48qq/dXkpcVTaNNrGGNCiJuBchf62V0PHFHV04/g\nMIBveo21xdUsm27TaxhjQoubbq6/AS4Etjvbs4EdQLqI/KuqvhKo4MLBzrJ66lo6WVho7Q/GmNDi\npg2iDLhAVS9S1YuAOcBB4Ergx4EMLhys3l8NwGXWQG2MCTFuEsQUVd3Zu6Gqu4BpqnowcGGFj1X7\nqpg5Oo1sm17DGBNi3CSInSLygIhc7jx+A+wSkXigM8DxhbTm9i42H61lgVUvGWNCkJsE8c9AMXCP\n8zjo7OvEN8OrOY11h2ro7FYWFVr3VmNM6HHTzbUV+Jnz6K/pnEcURlbtc6bXGG/TaxhjQk9or4cX\n5Fbvr2LuhCwSYm16DWNM6LEEESBlda0cqGq27q3GmJBlCSJA1jjdWxda+4MxJkSdtg1CRJbjW/vB\nL1W9PiARhYlV+6vITY1nSp5Nr2GMCU0D3UH8FF/D9CGgFd8iQQ/ha5g+4ObkInKNiOwVkWIR+cZp\nytwsIrtEZKeI/M/Qwg9OPc70GgsKs216DWNMyDrtHYSqvgkgIj9T1aI+by0XkY2DnVhEovGtJXEl\nUAJsEJEXnIF2vWUKgW8Cl6lqrYjknuHPEVR2ljVQ29Jp3VuNMSHNTRtEsoicXEhZRCYAyS6OuwQo\nVtWDqtoB/BG4oV+ZzwO/VtVaAFWtdBd2cFvlTO9t02sYY0KZm8n6vgK8ISIHAQHG41sjYjBjgGN9\ntkuAuf3KTAEQkbVANPA9VX2p/4lE5HbgdoD8/HwXl/bW6v1VTB+VRk6qTa9hjAldbgbKveRUBU1z\ndu1R1fZzeP1CYDEwFlglIrOdNbD7xvAg8CBAUVHRaRvOg0FLRxebjtTyL5dN8DoUY4w5K4NWMYlI\nEvA14E5VfRfIF5EPujh3KTCuz/ZYZ19fJcALqtqpqoeAffgSRshad/AEnd1q3VuNMSHPTRvEo0AH\nMN/ZLgX+y8VxG4BCEZkgInHALcAL/cr8Fd/dAyKSja/KKaRniV21v4r4mCiKCmx6DWNMaHOTICap\n6o9xZm5V1RZ8bREDclabuxN4GdgNPKOqO0XkfhHpHUPxMlAjIruAlcDXVLXmDH6OoLFmfzWXTMi0\n6TWMMSHPTSN1h4gk4gyaE5FJgKs2CFVdAazot+87fV4rcK/zCHnl9a3sr2zi5qJxgxc2xpgg5yZB\nfBd4CRgnIn8ALsM33bfpp3f1uIVTrHurMSb0DZggxDcMeA/wEWAevqqlL6tq9TDEFnJW768mJzWe\nqXmpXodijDFnbcAEoaoqIitUdTbwt2GKKST1Tq+xeEqOTa9hjAkLbhqpN4vIxQGPJMTtKm/gRHOH\nVS8ZY8KGmzaIucCnROQI0IyvmklV9byARhZibHoNY0y4cZMgrg54FGFgzf5qpo1MJTc1wetQjDHm\nnBi0iklVj+AbEb3Ued3i5rhI0tLRxcbDtSyaYqOnjTHhw81UG98Fvo5vWm6AWODJQAYVatYdOkFH\nd48tL2qMCStu7gQ+DFyPr/0BVS0DrB9nH2v2VxMXE8XFBZleh2KMMeeMmwTR4Yx47h1J7WYtiIiy\nen8Vc216DWNMmHGTIJ4Rkd8CGSLyeeA1fEuPGuB4fRv7KpqseskYE3bcrAfxUxG5EmgApgLfUdVX\nAx5ZiFjtdG9dMNkaqI0x4WXQBCEi9wJPW1Lwb01xNdkp8Uwbac0yxpjw4qaKKRV4RURWi8idIpIX\n6KBCRU+PsmZ/NQsLs4mKsuk1jDHhxc04iP9U1ZnAHcAo4E0ReS3gkYWAXeUN1DR3WPuDMSYsDWXA\nWyVwHKgBcgMTTmhZU+yb1HaBTa9hjAlDbgbKfUlE3gBeB7KAz9s8TD6r91f5ptdIs+k1jDHhx81c\nTOOAe1R1a6CDCSWtHd1sOFTLZy4d73UoxhgTEG66uX5TRM4XkTudXatV9d0AxxX01h2qoaO7hwWF\n1r3VGBOe3FQx3Q38AV+7Qy7wpIjcFejAgt2qfdXEx0Qxd4JNr2GMCU9uqpg+B8xV1WYAEfkR8Dbw\n/wIZWLBbtb+KS2x6DWNMGHPTi0mA7j7b3c6+iFVW10pxZROLrHrJGBPG3NxBPAqsE5HnnO0bgd8F\nLqTg1zu9hq3/YIwJZ24aqX/udHNd4Oy6TVW3BDSqILdqXzV5afFMyUvxOhRjjAkYN3MxzQN2qupm\nZztNROaq6rqARxeEunuUNcXVXDkjD5GIrmkzxoQ5N20QDwBNfbabnH0RaVtJHfWtnVa9ZIwJe64a\nqZ0FgwBQ1R7ctV2EpVX7qhGx6TWMMeHPTYI4KCJ3i0is8/gycDDQgQWrVfurmD0mnczkOK9DMcaY\ngHKTIL4IXAqUAiXAXOD2QAYVrOpbO9l6rM66txpjIoKbXkyVwC3DEEvQe/tANd09au0PxpiI4Gaq\njSki8rqI7HC2zxORbwc+tODz5r5qUuJjuCA/w+tQjDEm4NxUMT0EfBPoBFDVbUTgHYWqsmpfFfMn\nZREbPZRlNIwxJjS5+aZLUtX1/fZ1uTm5iFwjIntFpFhEvjFAuY+KiIpIkZvzeuFQdTOlda1WvWSM\niRhuEkS1iEwCFEBEbgLKBztIRKKBXwPXAjOAT4jIDD/lUoEvA0E98G7VPmd6DVte1BgTIdwkiDuA\n3wLTRKQUuAdfz6bBXAIUq+pBVe0A/gjc4Kfc94EfAW3uQvbGqv3VjM9KYnxWstehGGPMsBgwQYhI\nFFCkqsuAHGCaqi5Q1SMuzj0GONZnu8TZ1/f8FwLjVPVvg8Rxu4hsFJGNVVVVLi59brV1dvPWgWou\nt+olY0wEGTBBOKOm73NeN6tq47m6sJN8fg7822BlVfVBVS1S1aKcnOH/kn7nYA1tnT0smZY77Nc2\nxhivuKliek1Evioi40Qks/fh4rhSfOtZ9xrr7OuVCswC3hCRw8A84IVgbKheuaeShNgo5k/M8joU\nY4wZNm7mVPq483xHn30KTBzkuA1AoYhMwJcYbgE+efIEqvXAyRZfZ0rxr6rqRhcxDRtV5R97K7l0\nUratHmeMiShuRlJPOJMTq2qXiNwJvAxEA4+o6k4RuR/YqKovnMl5h9uBqmaOnWjl9kWTvA7FGGOG\nVUBnZVXVFcCKfvu+c5qyiwMZy5lauacSgCVTrYHaGBNZbEjwIFburWRKXgpjRyR5HYoxxgwrSxAD\naGzrZP2hE9Z7yRgTkVxVMYnIGGB83/KquipQQQWLNfur6epRlk61BGGMiTxu1qT+Eb6eTLuAbme3\nAmGfIFburSQ1IYYLx4/wOhRjjBl2bu4gbgSmqmp7oIMJJj09ysq9VSyakmOztxpjIpKrJUeB2EAH\nEmx2ljVQ1dhu1UvGmIjl5g6iBdgqIq8DJ+8iVPXugEUVBFburUQELrfurcaYCOUmQbzgPCLK63sq\nOW9sBtkp8V6HYowxnnAzkvpxEYkDpji79qpqZ2DD8lZFQxvvHqvjq1dNGbywMcaEKTe9mBYDjwOH\nAQHGichnwrmb66u7KgC4auZIjyMxxhjvuKli+hlwlaruBRCRKcBTwEWBDMxLr+6qoCAricLcFK9D\nMcYYz7jpxRTbmxwAVHUfYdyrqbGtk7cOVHPljDxExOtwjDHGM27uIDaKyMPAk872p4CgmpL7XHpz\nXxWd3WrVS8aYiOcmQfwrvrUgeru1rgZ+E7CIPPbKzgqykuO4MN9GTxtjIpubXkzt+JYG/Xngw/FW\nR1cPK/dUcu3skURHWfWSMSZXRli/AAAO9ElEQVSy2RwSfaw7VENjexdXzbDqJWOMsQTRxys7K0iM\njWZBYfbghY0xJswNmCBEJFpEfjpcwXhJVXl1VwWLptja08YYA4MkCFXtBhYMUyye2l5az/GGNqte\nMsYYh5teTFtE5AXgT0Bz705V/UvAovLASzuOEx0lLLXV44wxBnCXIBKAGmBpn30KhE2CUFVWbC/n\n0klZjEiO8zocY4wJCm66ud42HIF4aVd5A4drWvjC5ZO8DsUYY4LGoL2YRGSsiDwnIpXO488iMnY4\nghsuK7aXEx0lXG2jp40x5iQ33VwfxbcexGjnsdzZFxZ81UvHmT8xi0yrXjLGmJPcJIgcVX1UVbuc\nx2NA2Cyztru8kUPVzXxg9iivQzHGmKDiJkHUiMitzpiIaBG5FV+jdVh4r3opz+tQjDEmqLhJEP8C\n3AwcB8qBm4CwaLju7b00b2ImWba0qDHGnGLAXkwiEg18RFWvH6Z4htWe440crG7mswsneB2KMcYE\nHTcjqT8xTLEMuxXby4kSrPeSMcb44Wag3FoR+RXwNKeOpN4csKiGgaqy/N0y5k3MItuql4wx5n3c\nJIg5zvP9ffYpp46sDjnvltRzuKaFLy2e7HUoxhgTlAZrg4gCHlDVZ4YpnmHz1y2lxMVEcc1sq14y\nxhh/BmuD6AHuO9OTi8g1IrJXRIpF5Bt+3r9XRHaJyDYReV1Exp/ptYaiq7uHF7eVsWx6LmkJscNx\nSWOMCTluurm+JiJfFZFxIpLZ+xjsIKcH1K+Ba4EZwCdEZEa/YluAIlU9D3gW+PEQ4z8ja4qrqW7q\n4IY5Y4bjcsYYE5LctEF83Hm+o88+BSYOctwlQLGqHgQQkT8CNwC7Tp5EdWWf8u8At7qI56w9v7WM\ntIQYFk8NmwHhxhhzzrmZzfVMBwmMAY712S4B5g5Q/rPA3/29ISK3A7cD5Ofnn2E4Pi0dXby88zg3\nzBlNfIytHGeMMadz2iomEbmvz+uP9Xvvf53LIJzpO4qAn/h7X1UfVNUiVS3KyTm7v/pf3VVBS0e3\nVS8ZY8wgBmqDuKXP62/2e+8aF+cuBcb12R7r7DuFiCwD/h24XlXbXZz3rPx1Symj0xO4pGDQZhRj\njIloAyUIOc1rf9v+bAAKRWSCiMThSzgvnHISkQuA3+JLDpUuznlWqhrbWbW/muvnjCEqys2PYIwx\nkWugBKGnee1v+/0Hq3YBdwIvA7uBZ1R1p4jcLyK9czv9BEgB/iQiW521rwPmr1tK6e5RbrrIqpeM\nMWYwAzVSny8iDfjuFhKd1zjbCW5OrqorgBX99n2nz+tlQwv3zKkqT288xoX5GUzOTR2uyxpjTMg6\nbYJQ1bDq4rPlWB3FlU388COzvQ7FGGNCgpuBcmHhTxuPkRgbzXXn2cpxxhjjRkQkiJaOLpa/W851\n540i1abWMMYYVyIiQazYfpym9i5uLho3eGFjjDFAhCSIZzYeoyAriYsLRngdijHGhIywTxCHq5tZ\nf+gEHysah4iNfTDGGLfCPkE8teEo0VHCTReN9ToUY4wJKWGdINo6u3lmwzGunJ5HXpqroRvGGGMc\nYZ0gVmwvp7alk0/PH5Z1iIwxJqyEdYJ44p0jTMxJZv6kLK9DMcaYkBO2CWJHaT1bjtbxT/PGW+O0\nMcacgbBNEE++c4TE2Gg+cqE1ThtjzJkIywRR39rJX7eWcuMFo0lPtJHTxhhzJsIyQfx5UwltnT3c\nOs8ap40x5kyFXYLo7lEef/swF+ZnMHN0utfhGGNMyAq7BPHqrgqO1LTwuYUTvQ7FGGNCWtgliIdX\nH2RcZiJXzxzpdSjGGBPSwipBbDlay8Yjtdx26QSibc1pY4w5K2GVIB5ec4jUhBhuvtim9TbGmLMV\nNgni2IkW/r69nE/OzSclfqClto0xxrgRNgnisbcOEyXCP19a4HUoxhgTFsIiQdS3dPL0hmNcd94o\nRqUneh2OMcaEhbBIEI+9dZim9i6+ePkkr0MxxpiwEfIJoqm9i0fWHuLKGXlMH5XmdTjGGBM2Qj5B\nPPH2EepbO7lzyWSvQzHGmLAS0gmitaObh1cfZNGUHM4fl+F1OMYYE1ZCOkE8tf4oNc0d3LXU7h6M\nMeZcC9kE0d7VzW9XHWDexEwuLsj0OhxjjAk7IZsgnt9aRkVDO3ctLfQ6FGOMCUshO+T4xjljSEuI\n4VJbb9oYYwIiZBNEXEwU18wa5XUYxhgTtkK2iskYY0xgBTRBiMg1IrJXRIpF5Bt+3o8Xkaed99eJ\nSEEg4zHGGONewBKEiEQDvwauBWYAnxCRGf2KfRaoVdXJwP8BfhSoeIwxxgxNIO8gLgGKVfWgqnYA\nfwRu6FfmBuBx5/WzwBUiYiv9GGNMEAhkghgDHOuzXeLs81tGVbuAeuB93ZJE5HYR2SgiG6uqqgIU\nrjHGmL5CopFaVR9U1SJVLcrJyfE6HGOMiQiBTBClQN+1P8c6+/yWEZEYIB2oCWBMxhhjXArkOIgN\nQKGITMCXCG4BPtmvzAvAZ4C3gZuAf6iqDnTSTZs2NYnI3gDEG4qygWqvgwgS9lm8xz6L99hn8Z6p\nQz0gYAlCVbtE5E7gZSAaeERVd4rI/cBGVX0B+B3whIgUAyfwJZHB7FXVokDFHUpEZKN9Fj72WbzH\nPov32GfxHhHZONRjAjqSWlVXACv67ftOn9dtwMcCGYMxxpgzExKN1MYYY4ZfKCaIB70OIIjYZ/Ee\n+yzeY5/Fe+yzeM+QPwsZpE3YGGNMhArFOwhjjDHDwBKEMcYYv0IqQQw2O2ykEJFxIrJSRHaJyE4R\n+bLXMXlJRKJFZIuIvOh1LF4TkQwReVZE9ojIbhGZ73VMXhCRrzi/GztE5CkRSfA6puEkIo+ISKWI\n7OizL1NEXhWR/c7ziMHOEzIJwuXssJGiC/g3VZ0BzAPuiODPAuDLwG6vgwgS/xd4SVWnAecTgZ+L\niIwB7gaKVHUWvnFYbsZYhZPHgGv67fsG8LqqFgKvO9sDCpkEgbvZYSOCqpar6mbndSO+L4H+EyFG\nBBEZC1wHPOx1LF4TkXRgEb4BqKhqh6rWeRuVZ2KARGcKnySgzON4hpWqrsI3+LivvrNnPw7cONh5\nQilBuJkdNuI4iyxdAKzzNhLP/AK4D+jxOpAgMAGoAh51qtweFpFkr4MabqpaCvwUOAqUA/Wq+oq3\nUQWFPFUtd14fB/IGOyCUEoTpR0RSgD8D96hqg9fxDDcR+SBQqaqbvI4lSMQAFwIPqOoFQDMuqhHC\njVO3fgO+hDkaSBaRW72NKrg4c94NOsYhlBKEm9lhI4aIxOJLDn9Q1b94HY9HLgOuF5HD+Kocl4rI\nk96G5KkSoERVe+8mn8WXMCLNMuCQqlapaifwF+BSj2MKBhUiMgrAea4c7IBQShAnZ4cVkTh8jU4v\neByTJ5xV934H7FbVn3sdj1dU9ZuqOlZVC/D9f/iHqkbsX4qqehw4JiK9s3ZeAezyMCSvHAXmiUiS\n87tyBRHYWO9H7+zZOM/PD3ZAQCfrO5dONzusx2F55TLgn4DtIrLV2fctZ3JEE9nuAv7g/BF1ELjN\n43iGnaquE5Fngc34evxtIcKm3BCRp4DFQLaIlADfBX4IPCMinwWOADcPeh6basMYY4w/oVTFZIwx\nZhhZgjDGGOOXJQhjjDF+WYIwxhjjlyUIY4wxflmCMMYY45clCGOMMX5ZgjCuiUi3iGzt84i4eX7c\nEpHvichXz7acs77Dl85tdKe91luDvD9ssZjgEDIjqU1QaFXVOad705nWQFS1x9+22+PMKTKALwG/\nCfSFVHWw+YqGLRYTHOwOwpwVESlwVvn7PbADWNhve5yI3Ous7LVDRO45zXHj+p33fcf0OW63iDzk\nrBj2iogk9js2WUT+JiLvOsd/3Nl/q4isd+5+fussQoWIfFpEtjnlnxgohoGuLyL/LiL7RGQNMJXT\nOF05EfmriGxyznu7s/uHwCQn5p8MUM7fv8seEfmDE++zIpI0yGfbNMjP+L5Y+l3zfBFZJb6VDntE\nREXk/tN9DiYEqKo97OHqAXQDW/s8Pg4U4FuLYZ5Tpv/2RcB2IBlIAXbiW7/ilHL9ruP3mD7n7wLm\nONvPALf2O/6jwEN9ttOB6cByINbZ9xvg08BMYB+Q7ezPdBH3+67fp3wSkAYUA18d4Gd7X7k+107E\nlzSznOvt6HeO95Xzc50CfNM5X+ZsPwJ8dZDPtmmgz9hfLH2ulwDsAS5xtr8P/ARnOh97hObD7iDM\nULSq6pw+j6ed/UdU9Z0+5fpuLwCeU9VmVW3CN/XywtMch4tjwDeVc+8khZvwfXH1tR24UkR+JCIL\nVbUe34yeFwEbnAkOrwAmAkuBP6lqNYCq9q7CNVAM/q6/0Cnfor61OU430/BA5e4WkXeBd/DdURWe\n5hxuyx1T1bXO6yedn2mwz7bXYJ9xf8uAzaq63tnehi+R2WRvIcwShDkXmgfZdnucW+19XnfTry1N\nVffhWwdhO/BfIvIdQIDH+yS3qar6vUBc/0yIyGJ8X7LzVfV8fDOQJpxpOUf/L+ehfFkP9Weche/z\n7nUhvtlUEZHbRORa8Xmkf5WgCV6WIEygrQZuFN/c/MnAh5195/qYk0RkNNCiqk/iq+a4EN8i7TeJ\nSK5TJlNExgP/AD4mIlm9+88whlVO+UQRSQU+NMRy6UCtqraIyDRgnrO/EUjtc/zpyvmTLyLzndef\nBNacwc/VV/9Y+qoBzgMQkSnAR/At4gS+n3kB8FngaVVtdXk94zHrxWSGIlHeW38C4CXgvwc6QFU3\ni8hjQG/Vw8OqukV8a2kP6ZghxDkb+ImI9ACdwL+q6i4R+TbwiohEOfvvUNV3ROQHwJsi0o3vL/J/\nHmrcTvmngXfxrdS1YYjlXgK+KCK7gb34qo9Q1RoRWSsiO4C/A9/2V+409gJ3iMgj+BYOesBJLO/7\nuQY4R9/YT4lFVb/W5+2n8K3utwOoBj6hqjXOcQdE5EIgXVUfdnMtExxsPQhjwpCTyF5U1VkehwKA\niDwPfF5VB13m0gQPq2IyxgSMiKSLyC/xtf9YcggxdgdhjDHGL7uDMMYY45clCGOMMX5ZgjDGGOOX\nJQhjjDF+WYIwxhjjlyUIY4wxflmCMMYY45clCGOMMX79f+3lkPVvHPq+AAAAAElFTkSuQmCC\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x10c2106a0>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plot(Sigma_y, results)\n",
"xlabel(r\"Error on second data point $\\sigma_y$\")\n",
"ylabel(r\"Error on recovered gradient\")\n",
"xlim(0,10)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.4"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment