pyerrors/examples/01_basic_example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic pyerrors example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import numpy, matplotlib and pyerrors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import pyerrors as pe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.style.use('./base_style.mplstyle')\n",
    "import shutil\n",
    "usetex = shutil.which('latex') != ''\n",
    "plt.rc('text', usetex=usetex)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use numpy to generate some fake data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_sample1 = np.random.normal(2.0, 0.5, 1000)\n",
    "test_sample2 = np.random.normal(1.0, 0.1, 1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From this we can construct `Obs`, which are the basic object of `pyerrors`. For each sample we give to the obs, we also have to specify an ensemble/replica name. In this example we assume that both datasets originate from the same gauge field ensemble labeled 'ensemble1'. **For correct error propagation it is crucial that observations from the same Monte Carlo chain are initialized with the same name.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "obs1 = pe.Obs([test_sample1], ['ensemble1'])\n",
    "obs2 = pe.Obs([test_sample2], ['ensemble1'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now combine these two observables into a third one:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "obs3 = np.log(obs1 ** 2 / obs2 ** 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`pyerrors` overloads all basic math operations, the user can work with these `Obs` as if they were real numbers. The proper error propagation is performed in the background via automatic differentiation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we are now interested in the error of `obs3`, we can use the `gamma_method` to compute it and then print the object to the notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.367(20)\n"
     ]
    }
   ],
   "source": [
    "obs3.gamma_method()\n",
    "print(obs3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the method `details` we can take a look at the integrated autocorrelation time estimated by the automatic windowing procedure as well as the detailed content of the `Obs` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Result\t 1.36706932e+00 +/- 2.04253682e-02 +/- 1.02126841e-03 (1.494%)\n",
      " t_int\t 5.01998002e-01 +/- 4.47213595e-02 S = 2.00\n",
      "1000 samples in 1 ensemble:\n",
      "  · Ensemble 'ensemble1' : 1000 configurations (from 1 to 1000)\n"
     ]
    }
   ],
   "source": [
    "obs3.details()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected the random data from numpy exhibits no autocorrelation ($\\tau_\\text{int}\\,\\approx0.5$). It can still be interesting to have a look at the window size dependence of the integrated autocorrelation time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlkAAAGJCAYAAAC5Lib1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAss0lEQVR4nO3dbWxc133n8d/hsx7IGVG2bMWyaw2TIHHVbDKW4gZJsUg83A2QpovsSlYXzr4oUJObAt0mQEEu/cbImxDUiz5ssciSKbAvGqOQxK2xTQt4V7SLtClSxzabTZWkSKxR4iiRS1vUzFAPHJLDsy94Lz0kZ6ghOXfOffh+gAF5n4Z/jcdzf3POuecaa60AAADQXG2uCwAAAIgjQhYAAEAACFkAAAABIGQBAAAEgJAFAAAQAEIWAABAAAhZAAAAAehwXYArxhgj6T2SFlzXAgAAIqVX0i/sPSYbTWzI0lrAuua6CAAAEEnHJP18ux2SHLIWJOlnP/uZ+vr6XNeCJnu7tKgLr13TUyeP6f6+HtflAABiolQq6eGHH5Ya6AlLcsiSJPX19RGyYmhRXeo5cFC9fX3qI2QBABwIzcB3Y0zaGDNkjLm0y+N3dRziqbuzXSceSqm7s911KQCAhApFS5YxJivppKS0pP5dHH9aUq7JZSHCUvs6NfjYA67LAAAkWChasqy1s9baKUn5nR5rjElrF8EM8bZcWdU7t8parqy6LgUAkFChCFl79JSkC/fayRjTbYzp8x9au/wSMXXz9pL+7Ns/1c3bS65LAQAkVKRDljEmJ2mmwd3HJBWrHkzfAAAAAhPpkCUpba1ttItxXFKq6nEssKoAAEDihWLg+24YY4a8cVwNsdaWJZWrjg+kLgAAACmiLVne1Yivua4D4dbeRpAGALgT1ZasfklZb0yWJA1IkjFmRFLeWjvtrDKEwpG+Hv2XJ9/nugwAQIKFLWTVnIrBGJORdNpae06SrLUzqhrw7rVsDfnbAQDRMlda1NxCecv6I73dOsJdGxBRoQhZfoiSdFZrLVQTkl6tapHKSRqWtCVEeRORnvV+n5B0yQthSLAbt8p68ftv6dO//KAOH+x2XQ6Ae3j+lTf1xy/9eMv633vyffrS4PsdVATsnbHWuq7BCW+urGKxWOTehTE0V1rU86+8qaefeIRvwUAE+C1Zb8zd0hfPf1d/dPbDeu+Rg7RkIXRKpZJSqZQkpay1pe32DUVLFgAg2Y709WwIU+89clAnHko5rCjcZmZmNDo6qv7+fl26VPvWvbOzsxodHVU+n9eVK1daXKF07ty7nU83btzQ8PCwpqenNTIy0vJaXCFkAQBCobJq9b1rBUnS964V9MGjfVwlXEcul9PY2JjGx8fr7pPNZjU6Oqrh4eEdPXehUNCFCxd08eLFugHuXoaHhzU8PKxsNru+7syZM7t6rs388OYHx8nJyYaOSafTktb+fa0KepGcwgEAEC8vXr6uT0y8rGdfuCxJevaFy/rExMt68fJ1x5WFlx8attPfv7Nb+87OzurChQsqFAqan5/fZWXShQsXNgQsSfra17626+fzjY6OamRkRCMjI+vhanBwcNtj/FA2NDSkoaEhZbPZHQfP3SJkIZb69nXqMx86qr59na5LAXAPL16+ri98fVbXi4sb1r9VXNQXvj5L0GqhbDaroaEhZTKZPT1PoVBQPr/xhizpdFqnTp3a03POzs6qUCisrxseHtbMzMyWv1VtfHxcQ0ND68u5XE5TUw3PZb4ndBcilno62/X+B7gHOBB2lVWrL3/jB6p1CZaVZCR9+Rs/0OBjD7as6/DcuXPKZDLK5/PKZDI6ffr0+hgoaa1FJp/PK5/P68aNG5qYmFg/dmpqSplMZj1kpNPp9RN8o887Pz+v119/XZOTk5qamlJ/f7/Onz+vsbGxLa1DkjQ9vXYh/vz8fMNdYbVqaaZsNqvBwUFNTk4ql8utr99rN91rr72mfD6//jr4YbA6eFXL5/MqFAo1W/1mZmY21BYEQhZi6XZ5Rf/81oI+8GCvDnTzNgfC6jtX57e0YFWzkq4XF/Wdq/P62MDhwOs5c+aMzp49ux46BgcHlclklMvlNDExoeHhYc3Pz69vHxgY0NmzZ5XNZtfDjn/izufzmpmZ2fXzjo6ObghwzzzzjF5//fUN9c7OziqXy62HiKmpKQ0PD287TqleLbUC3G5dvHhRg4OD6115uVxOo6Ojewo16XRaN2/e3LDOf33rtbzVa+FKp9N1g1kzcfZBLN0ur+hvf/S2Hj60j5AFhNjcQv2AtZv99iKfz2t6eloXL15cX3fmzBlNTk5qcnJS/f39yufzG4KC3xrkB5SLFy/qqaeeUjqdViaT0cmTJ3f9vNWy2WzNwJDNZje00gwNDckYo9HR0ZrB4161NEsmk9GVK1c0MzOjS5cuaWZmRoODg7p48eJ6uGt0XNTjjz++obuv2vj4uCYnJxsan1atv79/T2POGsXZBwDgzJHexubAanS/vZiZmVE6nV5vHZHWrmCrDjebg0s6nV4/WZ8+fVqTk5M6dOiQstmszp49q5GREU1NTe3qeQcGBnb178hkMpqdna0Zshr5NzZTLpdbD4+jo6N65pln1kPWXkPd6Oiozp49WzeAbacVAUsiZAEAHPro8X4dTfXoreJizXFZRtKDqR599PjOrpLbjUKhsN6F59tp99alS5c0OzurmZmZDSFir8+7E9sFiGb8G++lUChoZmZmyziviYkJnTt3ru4YqZ2Ynp7WwMDAPQNWvW5E/3UIGiELAOBMe5vRc599TF/4+qyMtCFo+cPcn/vsYy0Z9J7NZmvOO9VoKJiamlqfIsC/Su/JJ5/UxMTEnp53pwqFQt3xVXv9Nzbq1VdfrTmYPpPJrP+d3XYX+q1w/jp/uolaocn/e/4A/2pBD3qXmMIBMdXV0abM/QfU1cFbHAi7T584qq9+PqsHUxu7BB9M9eirn8/q0yeOtqSOXC6nkydPrg9g9124cKHuMdWDpwuFwpapAfxWo70873b8q+d8ftCr10rTaC31WsPy+XxD0x9MTU1t6JKUtKV1yx8Hdq9HdcCanZ3V7Ozs+hg1vx5/PrB8Pr9hpnlJGhsb21DL9PT0rroYd4N7F3LvQgAIhcqq1flX39SzL1zWVz53QmdPPeJkxvfR0VENDAysn7hPnz6t2dlZjY+Pa3p6WhMTExoZGdG5c+c0Pj6uTCajsbGx9WBSfcIfGhpab7nZ6fOePHly/dY5/j4jIyMaGxtTOp3W7Ozs+rQN0tYpHKqfe2RkZMOVirVq8Wuenp7W+fPnNTs7q5GREZ06dWp9+9TUlEZHR3X16tW6LV/+jPH+2LBqe5nCoVAo6Pjx4zUDqJ9lpqenNTo6uuU2Qv6UFdJaK1v1a7FTO7l3ISGLkBVLlVWr8kpF3R3t3JYDiJDLPy/q1//kW/qr3/0E9y4MKX+QfCvGNIXRTkIWfSmIpRu3ypr8Zl43bpVdlwIAsVLvykVsxcB3AIBzc6VFzS2U9cbcLUla/3mkt1tH+oKfvgGNa9X0B3FAyAIAOPf8K2/qj1/68fryF89/V5L0e0++T18afL+jqrBZPp/XyZMnXZcRGYQsAIBzTz/xiAYfe2DL+iO93Q6qQT10E+4MIQsA4NyRvh66BRE7hCzE0n0Hu/U7nxxQZxvXdgAA3CBkIZba2oy629pdlwEASDC+5iOWbt5e0l/MXtPN20uuSwEAJBQhC7G0XFnVT2/c0XJl1XUpAICEImQBAAAEgJAFAAAQAEIWAABAAAhZiKWDPR365AeO6GAPF9ACANzgDIRY2t/VoQ8/nHZdBgAgwWjJQiwtLlf0w+slLS5XXJcCAEgoQhZiqXR3WS9efkulu8uuSwEAJBQhCwAAIACELAAAgAAQsgAAAAJAyEIsdbS36WiqRx3tvMUBAG6EZgoHY0xa0lOSzlhrBxs8ZsT7dUCSrLXDwVSHqOk/0KXf/OgjrssAACRYKEKWMSYr6aSktKT+Bo+ZsNaOVi1PGmMuNRrQAAAAghSKvhRr7ay1dkpSvpH9vVavrPfTNykpZ4zJNL9CRM1caVF/eOl
      "text/plain": [
       "<Figure size 640x395.55 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "obs3.plot_tauint()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This figure shows the windowsize dependence of the integrated autocorrelation time. The vertical line signalizes the window chosen by the automatic windowing procedure with $S=2.0$.\n",
    "Choosing a larger windowsize would not significantly alter $\\tau_\\text{int}$, so everything seems to be correct here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Correlated data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now generate fake data with a given covariance matrix and integrated autocorrelation times:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "cov = np.array([[0.5, -0.2], [-0.2, 0.3]])  # Covariance matrix\n",
    "tau = [4, 8]  # Autocorrelation times\n",
    "c_obs1, c_obs2 = pe.misc.gen_correlated_data([2.8, 2.1], cov, 'ens1', tau)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and once again combine the two `Obs` to a new one with arbitrary math operations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Result\t 3.27194697e-01 +/- 1.53249111e+00 +/- 2.49471479e-01 (468.373%)\n",
      " t_int\t 4.75187177e+00 +/- 1.33949719e+00 S = 2.00\n",
      "1000 samples in 1 ensemble:\n",
      "  · Ensemble 'ens1' : 1000 configurations (from 1 to 1000)\n"
     ]
    }
   ],
   "source": [
    "c_obs3 = np.sin(c_obs1 / c_obs2 - 1)\n",
    "c_obs3.gm()  # gm is a short alias for gamma_method\n",
    "c_obs3.details()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This time we see a significant autocorrelation so it is worth to have a closer look at the normalized autocorrelation function (rho) and the integrated autocorrelation time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmYAAAGfCAYAAAD1WR7GAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAA9hAAAPYQGoP6dpAABCXklEQVR4nO3dfXAc52Hn+d8DvgAg8TIEKZAwKUoERVmEKVkCrbeYsXwW4Cgbs3K+o8ykpFzqyisy8l3O0lWWtLR1x2ivSgqYeNc++1Yhld3UlqVyKDHHxPLGzImS7axtWZIJMRJMxSYJyhQpiLAEDgC+gC/Ac390NzQYzMvTwDSme/D9VE21Zrp75sE0NfOb59VYawUAAIDyqyp3AQAAAOAhmAEAAMQEwQwAACAmCGYAAAAxQTADAACICYIZAABATBDMAAAAYoJgBgAAEBMEMwAAgJiYW+4CAEgeY0xKUpekDkmtkrol/SzjkCZJe6y1e7POa/fP+4Sk7dba3TNS4DIwxnRI2mqtvbcEz9UqabukdMbDu+S9/westb3TfQ3/ddolbZaUkn9drbXbHc9Nybu2B/3zF7ueC+BDhiWZAEyVMWaTpOckdVprD2Tte0FSk7V2fY7zzqgCg1lGOJG80NRrre2c5nO2ynuP77bWpjMe3yXpC5LWlyKYBaEsM0wVuoY5zj/mlyXt39/kP9+0gykwm9CUCWA60gX23Sup3Q8Q2QaiKU55WWvT1tqt1tqt8moRS6FL0hOZocx/ra0lev7A1hw1XFvlXcMthU40xmyTV7uWzijfXkkdfuAD4IhgBiASGV/SHeUsRwXoUP4AXMoaxy3ZITqjJq5Yrd9mSa/leLxXXH8gFIIZgEj4zXpS6WqOZqteef3Lctmj0tU+HpB0bIrntit3eOxV8VAHIAOd/wFEpUveF/MD+Q7w+yE1yessfqukB7Kb7Pw+Vlv1YWhYHaZTeUaH9mOSVks6FvRtyxqMEJQzb3n8Dv2t8sJQcJystTtdyzMFT0h6zu/DtV1eZ/+0/7olC725+sL5772UuzYs+5hcAXFAXmgD4IhgBqAU7s34gl4tL7y8VqQfVKe8AQC9kmSMeU5eSBo/xw9OT2V2PjfGtGZ3NM/HD1K7rLWrMx47aIyRtXa3H2w6/cEInZK68pUnCIiZndn954+0Rshau9cYs90vy3P+6/b6f1eUgVDy/vb0NF+nqVSFAWYDghmAUnguGJXpB5hd8mqcCsoaTfiaMkJZ8Lz6cJTj+DnGmG5lhbg8dmWfL68GqksT+2cNyBt9WKg8rf4tsywHMgJpZKy1O40xu+WNwuyU12+ryxizOSu07souo8Nz5wyWfijeJqnoiMwCCGVASAQzACXlh5tOY8wZY8yuArVm2c1j6cw7fjBoldf3KdsLKhLMMs7/WdaubuUOLwXLE4Qwv3btWUkvWGv3ztSUH37t4G7/FoyE7DLGbAnKUOKRmk/JmwalWHNp0ISZL4RV5AhcICp0/gcQlQPyanjySRc5v1DfJUlKZQwwKHR+hzFmS3CTV9uUK8AUK48krZIXjDrk9fs6E/V0EH4/vEn85sUDiqAp1a95eyJ7bro85Uj7/5nKsTslr58hAEfUmAGIUsoYkyrWFyyP4Au9SZNDU0qaEAoKnX+gFJ3kjTGtfm3gdknbMyaTfUrTa+4rplPS3jz7XlBGk/EUmjLT2RPA+jVxL2Su2mCMaS/yHnZLWpzj8SZ5I0cBOCKYAYhKUNMVLNmkMCHNWtttjEnLq53Kbi68VfnDSvb5n1DWlB3GmA6X2qAs7X5A2es/f1rSVr9pM0pfMMZsz/O+rZYXzuSXaVpNmX7tXDp7KS1516BQMNuj3DV37SowKhfAZDRlApiO9qxtpiAwdEgTpsYIpLKOz74veasHbM1ssvSbDl2/8O/Vh7Vbwfkp5Q4RLuV5JPOO/1z5Al5KefpdGWNS/ujQgjPqZ3gqu8nUv/+JUvVx858vGIG6JePWlXVcUPbMJtbdklozB0L4f9vuUtRWArMJa2UCCC2jGS+d8XBK3vQX6YzjtsgLR8/5D/1M3pf/FnlNjXuttdv9L/9N8mrX9srr3xTUsgWBIZjHbLFyLFFUoKyZ56clKWseM6fy+P8dzF0WvHZKXvjI/Ju7/MeD0LXbP368zP77d1xeM2vBtSSDART+exlM+5GSV7NVskXC/Zq/VJ7d9wa1aBlln7DWacZC6yxiDkwDwQwAyiRzRCUASDRlAkA5pcpdAADxQjADgDLwm1HpfwVgAoIZAJTHVEaGAqhw9DEDAACICWrMAAAAYoIJZkMwxhhJH5E0XO6yAACARKmX9K4t0lRJMAvnI5JOlrsQAAAgkVZIOlXoAIJZOMOS9M4776ihoaHcZQEQd8PD0qFD0s03S/X15S4NgDIZGhrS1VdfLTm0uBHMpqChoYFgBqA4Y6SFC6WGBoIZACd0/geAqMybJ7W0eFsAcECNGQBEpaZG+uhHy10KAAlCjRkARGVsTDp3ztsCgAOCGQBE5dw56bXXvC0AOCCYAQAAxATBDAAAICYIZgAAADHBqMxpGB2zevX4gPqHR9RcX6PbVjVpTpUpd7EAxIUxUlWVtwUAB4kOZsaYlKQvSLrXWtvpeM42SWn/bspau3Mqr72/p0+PPX9YfYMj44+1NNZox8Y23bOuZSpPCaDS1NVJn/pUuUsBIEES25RpjGmXF8pSkpocz9kmSdba3dba3ZK6jTG7wr72C4ff04NPd08IZZL03uCIHny6W/t7+sI+JQAAQHKDmbW22w9XvSFOe0TS7oznOCBpS9jX/rPv/YtyLQ0fPPbY84c1OlZw8XgAs8H589LPfuZtAcBBYoNZWMaYVnlNl+kc+zrynFNtjGkIbpLqJen00MW8r2Ml9Q2O6NXjAyUpN4AEGx2Vzp71tgDgYNYEM0mteR5Py2sOzeURSYMZt5OuL9Y/PFL8IAAAgAyzKZjlM6D8fdSekNSYcVvh+qTN9TXTLxkAAJhVEj0qs0TyDhyw1l6UNN5uafwh70sbqvX+ReXsZ2YkLWv0ps4AAAAIYzbVmOUbJJAqsC+nr/z2DZK8EJYpuL9jYxvzmQGQamulj33M2wKAg1kTzKy1vZLS/iCA7H0HwjxXZ9syPXl/u5Y1TmyuXNZYoyfvb2ceMwCeuXOlq67ytgDgoBI+LXK2GfoBbFPWBLJPSOqQP2WGMWaTMqbPCOOedS3qbFvGzP8A8rt0STp9Wlq6VJo/v9ylAZAAiQ1mQfCStFlSuzGmS9Jr1tq9/iEdkrZKGg9m1tqdxphtfiCTpFuttVunWoY5VUZ3rl481dMBVLqLF6Vjx6RUimAGwImxlolQXflzmQ0ODg6qoaGh3MUBEHfDw9LBg9L69VJ9fblLA6BMhoaG1NjYKEmN1tqhQsfOmj5mAAAAcUcwAwAAiAmCGQBEZe5cafFiRmUCcManBQBEpbZWuvHGcpcCQIJQYwYAUbFWunzZ2wKAA4IZAETl7Fnpxz/2tgDggGAGAAAQEwQzAACAmCCYAQAAxATBDAAAICaYLgMAolJXJ23YIM2ZU+6SAEgIghkARMUYJpcFEApNmQAQlQsXpDfe8LYA4IBgBgBRuXJFGhjwtgDggGAGAAAQEwQzAACAmCCYAQAAxATBDACiUlMjrVnjbQHAAeO4ASAq8+ZJy5eXuxQAEoQaMwCIyuXL0unT3hYAHBDMACAqIyPSW295WwBwQDADAACICYIZAABATBDMAAAAYoJgBgBRmTNHamjwtgDggOkyACAqCxZI7e3lLgWABKHGDAAAICYIZgAQleFh6Qc/8LYA4IBgBgAAEBMEMwAAgJggmAEAAMQEwQwAACAmmC4DAKKycKF0++1SdXW5SwIgIQhmABCVqiqptrbcpQCQIDRlAkBURkakt97ytgDggBqziPUPjah/+OKkx5vrq9XcUFOGEgGYMZcvS6dPSytWSDX8/w6guMQHM2PMNklp/27KWrvT4ZwtklL+easlPWGtTRc4ZcqeeeWEvv7ikUmPf/nuNXq48/ooXhIAACRUooOZH8pkrd3t3+8wxuyy1m4tcs7uIIgZY1KSnpJ
      "text/plain": [
       "<Figure size 640x395.55 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlQAAAGJCAYAAABM0K1FAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAyeElEQVR4nO3dXXBb533n8d/Dd0kkAVE2LcWyGlG2m6hK69CyXU+zs5uY7LrTajvpUFJn7N2ZThuq6U6aZKYVK99ofBOW2pkmmU7HpdKrrT07krj17jqdais63bTZdS3baJoySppYkK3IpkNbFADqhW/gsxc4oAHwkDggDnBwgO9nBgPx4CH5Fw8J/PC8HWOtFQAAADavKegCAAAAwo5ABQAAUCYCFQAAQJkIVAAAAGUiUAEAAJSJQAUAAFAmAhUAAECZWoIuoFqMMUbSRyTNBV0LAAAIlS5J79oNNu9smEClTJi6FnQRAAAglHZLeme9BxspUM1J0k9+8hN1d3cHXQuAGvd+al5nX7+mIwd36+7ujqDLARCQVCql++67TyoywtVIgUqS1N3dTaACUNS82tSxrVNd3d3qJlABKIJJ6QDgor21WQfujai9tTnoUgCEQMP1UAGAF5EtrRrcf0/QZQAICXqoAMDFUnpFH9xc0FJ6JehSAIQAgQoAXNy4tai/fOVt3bi1GHQpAEKgZgKVMSZqjBk2xlzw0LZoGwAAgGqpiTlUxph+SQclRSX1FGk7JGmgCmUBAAB4UhOBylobkxRzwtK6jDFRFQlcAAAA1VYzQ34eHZF0NugiADSG5iYTdAkAQqImeqi8MMYMSJosoX27pPacQ12+FwWgbvV2d+j3n3gg6DIAhERoApWkqLU27gz7eXFC0skK1gMAqGMzqXnNzC2sOd7b1a5eds9HgVAEKmPMsLX2dImfNirpT3I+7hIXRwbg0fWbCzr//ff05M/t1I7O9uKfgLrzwqtX9fWXf7zm+BefeEBfHnwwgIpQy2o+UDkrAF8v9fOstQuSVt9aGMNcCADepVesZlILSq/YoEtBQJ56bI8G99+jN2du6ktnvquvHX1I9/d2qreLgI21aj5QKbOqr9+ZQyVJ+yTJGHNcUtxaOxFYZQCAutXb3ZE3tHd/b6cO3BsJsKLaderUqdV/X79+XceOHdPExISOHz8eYFXVVWuBas2WCNbaSeVMRnd6rIattacK2wIA4Kf0itX3riUkSd+7ltDHd3XX1erPRCKhs2fP6ty5c7pwYXN7Zh87dkzHjh1Tf3//6rHDhw/7Ul82qF2+fFmSND4+7ulzotGopMz/r1qhria2TTDG9Dk9TseU6Y0ac9uTyjl2wvn3WE6vFQAAvjo/Na1PjX1Lz7w4JUl65sUpfWrsWzo/NR1wZf6IxWI6e/asEomEZmdnN/11zp49mxemJOkb3/hGueVpZGREx48f1/Hjx1eD1ODg4Iafkw1gw8PDGh4eVn9/v44dO1Z2LV4YaxtjfoAxpltSMplMqru7O+hyANS4+aW0rs7e1p6erepobQ66HFTZ+alpff75mApfIbN9U8893a8nD+yqdlkVMTExodHRUb3xxhub+nxjjC5fvqy+vr6846dOndp071AikdDhw4d17ty51d6mWCymhx9+2PV7ZW3fvl1XrlxZ/ZxsfeVknVQqpUgkIkkRa21qvXa1NuQHADWho7VZD97D9nWNKL1i9exLl9aEKUmyyoSqZ1+6pMH9O6s2/Hfq1Cn19fUpHo+rr69PQ0NDmpyc1MjIiKRMj1A8Hlc8Htf169c1Nja2+rmnT59WX1+fEomE4vG4otGohoeHfautv79fg4ODGh8f18DAhwNH5Q61vf7664rH46u9X9kQlUgkXNvH43ElEom8MJU1OTmZV1sl1MSQHwDUmlsLy3rj7Ru6tbAcdCmosotXZjWdnF/3cStpOjmvi1c2P0xWisOHD6+GqOzwVywW08DAgMbGxlaH7LKPT0xMKBaLScr0PknSwMCAhoaGNDS04RXeNuXcuXOSMsNxxhgNDg5qctLzPtyuotGobty4kTeUmP2a6/VOxePxdb/WeiHMT/RQAYCLWwvL+vsfva/7tm/RtnaeKhvJzNz6YWoz7coRj8c1MTGxGlqkTMAaHx/X+Pi4enp6FI/H83pfsj1Z2TBy7tw5HTlyRNFoVH19fTp48KCvNfb19eny5cuanJzUhQsXNDk5qcHBQZ07d241wHmdx/Twww+v23s2Ojqq8fFx1x6ojfT09JQ1R8wrniUAAMjR2+VtF3Sv7coxOTmpaDSa1+Nz+fLlvN6Ywh6baDS6GiCGhoY0Pj6u7du3q7+/X0ePHq3YqreBgYHVYDcyMqLPfe5zq4HKy+q8jYyMjOjo0aObGqqsRpiSCFQAAOR5dG+PdkU69F5y3nUelZG0M9KhR/eu2enHd4lEQn19fXk9UKXOBbpw4YJisZgmJydXg41foSqRSGhycnLNUOLY2JhOnTq17pymUkxMTGjfvn1Fw9R6Q4HZn2GlEagAAMjR3GR08tB+ff75mIyUF6qyU9BPHtpflQnp/f39Gh0dXXPca1A5ffr06vYB/f39Gh4e1hNPPOFrL9Vrr73mOjerr69vtcbNDvlle+ayx7LzxdwCUvb7ZSfu56r0hHSJQAUArtpamtR39za1tbB2pxE9eWCXnnu6X8++dClvgvrOSIdOHtpftS0TBgYGdPDgQU1MTOSFlrNnz67bY5M7ATuRSKyGqqzCsLHekFg8Htfk5GTRnqHTp09rcHAwL7QU9lptZsgvFospFotpaGhodYhzYmJitZ7s/LLccHjixIm8mnPbVxr7UAEAsI70itWZ167qmRen9JXPHtDRR/YEslP6yMiI9u3bp56ezDDj0NCQYrGYRkdHNTExobGxMR0/flynTp3S6Oio+vr6dOLEidWwlP28eDyu4eHh1Z6ciYkJnTlzRrFYTMePH9cjjzyyGoROnz6tkZGRNfs65crutN7X17e6sjCrnF6wRCKhvXv3uq7Oy+aWiYkJjYyMrO6inpXdYkLK9J7lbiGxGV73oSJQAYCL9IrVwnJa7S3NdXWpEZRu6p2kfu1Pv6NvfuFTDXctv2zPUDXmINUqr4GKvmwAcHH95oLGvx3X9ZsLQZcCBCYWizV0mCoFc6gAAHAxk5rXzNyC3py5KUmr971d7ertrvyWCbWgWlsO1AMCFQAALl549aq+/vKPVz/+0pnvSpK++MQD+vLggwFVVT3xeNz3TUDrGYEKAAAXTz22R4P771lzvLerPYBqqo+hvtIQqAAAcNHb3dEwQ3soH4EKAFzc1dmu3/v0PrU2sXYHQHEEKgBw0dRk1N7UHHQZAEKCt14A4OLGrUX9VeyabtxaDLoUACFAoAIAF0vpFb19/baW0itBlwIgBAhUAAAAZSJQAQAAlIlABQAAUCYCFQC46Oxo0ac/1qvODhZDAyiOZwoAcLG1rUUP3RcNugwAIUEPFQC4mF9K6wfTKc0vpYMuBUAIEKgAwEXqzpLOT72n1J2loEsBEAIEKgAAgDIRqAAAAMpEoAIAACgTgQoAXLQ0N2lXpEMtzTxNAiiObRMAwEXPtjb95qN7gi4DQEjw1gsAAKBMBCoAcDGTmtdXL/xIM6n5oEsBEAIEKgAAgDLVzBwqY0xU0hFJh621gy6PH3f+uU+SrLXHqlcdAADA+moiUBlj+iUdlBSV1OPy+Ji1diTn43FjzAW34AUAAFBtNTHkZ62NWWtPS4oXPub0XPU791njkgaMMX3VqRAAAGB9NRGoPDgoKTc8ZYNXtPqlAGgEPdva9Fu/9FH1bGsLuhQAIVATQ34bsdYmJG0vODzg3K/p0coyxrRLas851OVvZQDqWUtzk6JbCVMAvAlLD1WhE5KOOWFrozbJnNu1KtQFoE4kby/p/NS0kreXgi4FQAiELlAZY8YknXHmXG1kVFIk57a70rUBqB8Ly2n9YHpOC8vpoEsBEAI1P+SXyxgzJOmyhzAla+2CpIWcz61kaQAAoIGFpofKGDMgSdkwZYyJssoPAADUgloLVGv2oJJW96nqlxQzxvQ5QWpY0mw1iwMAAHB
      "text/plain": [
       "<Figure size 640x395.55 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "c_obs3.plot_rho()\n",
    "c_obs3.plot_tauint()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now redo the error analysis and alter the value of S or attach a tail to the autocorrelation function via the parameter `tau_exp` to take into account long range autocorrelations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Result\t 3.27194697e-01 +/- 1.67779862e+00 +/- 2.08884244e-01 (512.783%)\n",
      " t_int\t 5.69571763e+00 +/- 2.09295390e+00 tau_exp = 20.00,  N_sigma = 1\n",
      "1000 samples in 1 ensemble:\n",
      "  · Ensemble 'ens1' : 1000 configurations (from 1 to 1000)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlQAAAGJCAYAAABM0K1FAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAA9hAAAPYQGoP6dpAABEZUlEQVR4nO3deXxcV33///eZ0b7NWLZleY0tOZswWWQnISwtSaw2FBwKeGm/CfRBCzYppUDbr01o+zUpFD/klq2Upnb4lv5o0jaxgW8IoS5WCFsIsR0lpIqz2XIWO3JkWxqt1jZzfn/MkpE00lx5rmbRvJ6Phx6y7r2aOb6S5r7nnM8511hrBQAAgAvnyXQDAAAAch2BCgAAIEUEKgAAgBQRqAAAAFJEoAIAAEgRgQoAACBFBCoAAIAUFWS6AelijDGSlkjqy3RbAABATqmU9JqdZvHOvAlUCoepk5luBAAAyEnLJJ2aamc+Bao+SXr11VdVVVWV6bYgh5zpHdL9R05q87plWlhVkunmAADSqLe3V8uXL5eSjHDlU6CSJFVVVRGoMCNDKlJJeYUqq6pURaACACRAUTqQRHGhV2uW+lRc6M10UwAAWSrveqiAmfKVFqqpYVGmmwEAyGL0UAFJjAZDOts/rNFgKNNNAQBkKQIVkET3wIj+7bGX1T0wkummAACyVNYEKmOM3xiz1Rhz0MGxSY8BAABIl6yooTLGNEpaJ8kvqTrJsRslrU9DswAAABzJikBlrW2V1BoJS1MyxviVJHABAACkW9YM+Tm0WdL9mW4E8o/XYzLdBABAFsuKHionjDHrJbXM4PhiScVxmypdbxTyQk1Vif70posz3QwAQBbLmUAlyW+tbY8M+zlxh6Sds9geAECO6OwdUmff8KTtNZXFquEOCHBBTgQqY8xWa+3eGX7bLklfjvu6UtwcGRfgXP+wDjxzWje/qVbzK4qTfwOArHPv46/oaw+/OGn7J2+6WJ9uuiQDLcJck/WBKjID8MhMv89aOywp9nbEGGpgcGGCIavO3mEFQzbTTQFwgW69boWaGhbpWGe/PnXfU/rqlqu0uqZCNZW8SYI7sj5QKTyrrzFSQyVJ9ZJkjNkuqd1auz9jLQMA5ISaqpJxQ3urayq0Zqkvgy1yXyAQ0P333699+/bp4MG5uVzj7t27JUnHjx+XJO3Zs2fSfr/fLyl8PrZv3562tmVboJq0JIK1tkVxxeiRHqut1trd6WwYACC3BUNWT58MSJKePhnQ5Yur5swM3tbWVh05ckSBQEBdXV2Zbs6s2LFjh5qbm2Nfb9u2TU1NTbHwGA1bW7dulSS1tLRo27Ztk0LXbDHWZn4YwxhTJ2mjpC2SGiXtlnR4Yu9TZJ2qLZFjd0s6GAlcTp6jSlJPT0+Pqqqq3Gw+5rjO3iHd+/gruvW6FRSvAjnqQFuH7nzwqDp6hmLbFvtKtHNDg25esziDLXPX/v37tWvXLj3xxBOZboqrAoGANm3apH379sV6oFpbW7V27VodP35cdXV1mjdvnk6cOBHbL4XLfVLNOb29vfL5fJLks9b2TnVcVvRQWWvbFQ5I0/Y6RQIWQ3xIq6rSQr37isWqKi3MdFMAXIADbR26/Z5WTbysnu4Z0u33tOqu2xrTGqp2796turo6tbe3q66uThs3blRLS4t27NghSbr77rvV3t6u9vZ2nTt3blyvzN69e1VXV6dAIKD29nb5/f5Yj4wbor06XV1dCgQC48LJxLCSbkeOHFF7e7saGxslSXV1dZIUOxcT2xvV0tKi9etn/wYrWRGogGxWUujVJYtYxgzIRcGQ1Z0PHp0UpiTJSjKS7nzwqJoaatMy/Ldp0yZt2bJFGzeGbwzS1NSkuro6rV+/Xs3NzbEwE91fX1+vLVu2qLGxUfv3h/sTouGgvb1dLS2Ol2dMKhAI6ODBg7H6pN27d6e1Bmk6fr9f3d3d47ZF/+91dXU6ciTx3DW/369AIDDbzZNEoAKSGhge03On+3RZbaXKi/mTAXLJoRNd44b5JrKSOnqGdOhEl66vnz+rbWlvb9f+/fu1b9++2LZNmzZpz5492rNnj6qrq9Xe3j6uNyXakxXtldm3b582b94sv9+vuro6rVu3ztX2RXvDZjuEbNu2zdFxa9eunbIHbteuXdqzZ8+0vWbV1dVpqynj6gAkMTA8pp+9cEbL55USqIAc09k3dZi6kONS0dLSIr/fP65X6fjx42pvb499HR3GivL7/bFAsHHjRu3Zs0fz5s1TY2OjtmzZ4moPUjS0Rdsa/7XbUi0U37Fjh7Zs2ZJ0uDOdBfpcHQAAc1ZNpbOJJE6PS0UgEIgN70XNtLbn4MGDam1tVUtLSyyUzMaw3H333ae7777b9cd1w/79+1VfXz8uTE0MolHRc54OBCoAwJx17apqLfaV6HTPUMI6KiOp1leia1dNWrXHdY2Njdq1a9ek7VMVU0+0d+9ebd26VY2NjWpsbNTWrVt10003zUqgam1tndUC9Asd8ov27kW3RZeJqKurk9/vjxX6x0tHQbpEoAIAzGFej9HODQ26/Z5WGWlcqIqWoO/c0JCWgvT169dr3bp12r9/f6zoXJLuv//+KYeu4muZAoFALFRFTQwPUw1xRQvYncwIjM4enGjHjh265pprYnVeXV1d2rZtm5qbm7V+/XqtXbtWzc3N8vv92rRpk5qbm1VdXa2DBw9qx44d49p6IUN+ra2tam1t1caNG2PDpPv374/9n+64445x/8f4femQFetQpQPrUOFCBQZH9NMXzug3L1kof1lRppsD4AJk0zpUO3bsUH19vaqrw71iGzduVGtrq3bt2qX9+/erublZ27dv1+7du7Vr1y7V1dXpjjvuiIWl6Pe1t7dr69atsZ6Z/fv367777lNra6u2b9+ua665Jhbc9u7dqx07djha+mD//v1qb28f1/O1d+/ecSuPRxfUbG9v17Zt27Rt2zatX78+9thr167Vww8/HGtbU1NTbPbghQgEAlq1alXCYvn4HBNdkkKSDh8+PG7JiQvldB0qAhUAIC8EQ1b3HX5Fn/1em774vjXacs2KObNSuhPRXp0LqSnatm2b/H6/mpqaJIVDXbRofffu3Tp8+PC42Ytr164dt7hookU3c4XTQOVJX5OA3BQMWQ2OjHFzZCDHeT1GVyzzS5KuWObPqzAlhYfMLrRAe+3atZLCw5bxNUnRmqW6urrYOllR0d6kQCCg6urqnAxTM0ENFZDEuf5hbj0DIOelsoTA1q1btXv3bu3du1fV1dWqq6vT3r171dzcrCeeeEKNjY1au3aturq6YnVL999/v6qrq3X48OE5e7PmeAQqAMCc19k7pM6+YR3r7Jek2OeayuK8eKPU3t6e8iKgE2cTRmcaSolXMo/uiy/An8sIVACAOe/ex1/R1x5+Mfb1p+57SpL0yZsu1qebLslQq9InXWsxxXO6HMRcQaACAMx5t163Qk0NiyZtr6kszkBr5rbW1la1t7dr7969WXMvwHQgUAEA5ryaqpK8GNrLBo2NjZOG//IBgQpIYkFFsf74hnoVepgUCwBIjEAFJOHxGBV7vJluBgAgi/GWG0iie2BE3209qe6BkUw3BQCQpQhUQBKjwZBePjeo0WAo000BAGQpAhUAAECKCFQAAAApIlABAACkiEAFJFFRUqAbLqtRRQmTYgEAiXGFAJIoKyrQVcv9mW4GACCL0UMFJDE0GtSzHb0aGg1muikAgCxFoAKS6D0/qgNtp9V7fjTTTQEAZCkCFQAAQIoIVAAAACkiUAEAAKSIQAUkUeD1aLGvRAVe/lwAAImxbAKQRHV5kX7v2hWZbgYAIIvxlhsAACBFBCogic7eIX3l4Avq7B3KdFMAAFmKQAUAAJCirKmhMsb4JW2WtMla25Rg//bIP+slyVq7LX2tAwAAmFpWBCpjTKOkdZL8kqoT7G+21u6I+3qPMeZgouAFAACQblkx5GetbbXW7pXUPnFfpOeqMfI5ao+k9caYuvS0EAAAYGpZEagcWCcpPjxFg5c//U1BvqkuL9KH37ZS1eVFmW4KACBLZcWQ33SstQFJ8yZsXh/5PKl
      "text/plain": [
       "<Figure size 640x395.55 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "c_obs3.gamma_method(tau_exp=20)\n",
    "c_obs3.details()\n",
    "c_obs3.plot_tauint()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}