{"metadata":{"language_info":{"name":"python","version":"3.7.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kernelspec":{"name":"python3","display_name":"Python 3 (ipykernel)","language":"python"}},"nbformat_minor":5,"nbformat":4,"cells":[{"cell_type":"markdown","source":"**How to Use**\n\nFirst, you need to load the pyemcee library, and define your function. For example: ","metadata":{},"id":"f1be8d57-d006-4024-877e-3eabc37ffae7"},{"cell_type":"code","source":"import pyemcee\nimport numpy as np\n\ndef myfunc21(input1):\n   result1 = np.sum(input1)\n   result2 = input1[1] ** input1[0]\n   return [result1, result2]","metadata":{"trusted":true},"execution_count":1,"outputs":[],"id":"5f544ead-f366-49f7-8ef8-572f88e59ba7"},{"cell_type":"markdown","source":"Then, specify the upper and lower uncertainties of the prior parameters:","metadata":{},"id":"8d134940-9e02-4b62-84de-62e798148221"},{"cell_type":"code","source":"input1 = np.array([1., 2.])\ninput1_err = np.array([0.2, 0.5])\ninput1_err_p = input1_err\ninput1_err_m = -input1_err\noutput1 = myfunc21(input1)\noutput1_num = len(output1)","metadata":{"trusted":true},"execution_count":2,"outputs":[],"id":"98891d4b-d6fc-489a-8dc2-053b9506f448"},{"cell_type":"markdown","source":"Choose the appropriate uncertainty distribution. For example, for a uniform distribution, use_gaussian=0, and for a Gaussian distribution, use_gaussian=1. Then, specify the number of walkers and the number of iterations, e.g. walk_num=30 and iteration_num=100. ou can then create the MCMC sample and propagate the uncertainties of the input parameters into your defined functions as follows::","metadata":{},"id":"d5ffb8e3-f2d0-4104-a333-f46879360d83"},{"cell_type":"code","source":"use_gaussian=0 # uniform distribution from min value to max value\nwalk_num=30 # number of walkers\niteration_num=100 # number of samplers\nmcmc_sim = pyemcee.hammer(myfunc21, input1, input1_err_m, input1_err_p, \n                          output1, walk_num, iteration_num, \n                          use_gaussian, print_progress=1)","metadata":{"trusted":true},"execution_count":3,"outputs":[{"name":"stdout","text":"Progress: 100% \n\n","output_type":"stream"}],"id":"bfe18fce-502c-482b-8955-e113570af668"},{"cell_type":"markdown","source":"To determine the upper and lower errors of the function outputs, you need to run with the chosen appropriate confidence level. For example, a 1.645-sigma standard deviation can be specified with clevel=0.90. For a 1-sigma standard deviation, we have clevel=0.682:","metadata":{},"id":"972852c4-9dbf-40a8-851c-66ebce26d5a7"},{"cell_type":"code","source":"clevel=0.68268949 # 1-sigma\noutput1_error = pyemcee.find_errors(output1, mcmc_sim, clevel, do_plot=1)","metadata":{"trusted":true},"execution_count":4,"outputs":[{"output_type":"display_data","data":{"text/plain":"<Figure size 600x600 with 1 Axes>","image/png":"iVBORw0KGgoAAAANSUhEUgAAAgQAAAH5CAYAAAD+5ibMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAiX0lEQVR4nO3dfZBW1WE/8O8CsqujrFJkeXETbBLfangJFGa1GSGuUnRI/aMTqxlhSLTVwgy60zbQKtTaiiaK2IrBlyBNWwvRiaYtVkMQZFRa5WWnxFETX2GMizI2LGKzJLvP7w9+rC7uwj7LvsDu5zNzZ3junvvc8xwO+3w559x7SwqFQiEAQJ/Wr6crAAD0PIEAABAIAACBAACIQAAARCAAACIQAABJBvR0Bdqjqakpv/jFL3LSSSelpKSkp6sDAMeMQqGQPXv2ZMSIEenXr+1xgGMiEPziF79IZWVlT1cDAI5ZO3bsyGmnndbmz4+JQHDSSScl2f9hBg0a1MO1AYBjR319fSorK5u/S9tyTASCA9MEgwYNEggAoAMON+VuUSEAIBAAAAIBABCBAACIQAAARCAAACIQAAARCACACAQAQAQCACACAQAQgQAAiEAAAEQgAAAiEAAA6UAg2LBhQ6ZPn54RI0akpKQkjz/+eLuPfe655zJgwICMHTu22NMCAF2o6ECwd+/ejBkzJkuXLi3quF/+8peZMWNGLrzwwmJPCQB0sQHFHjBt2rRMmzat6BNde+21ufLKK9O/f//Djio0NDSkoaGh+XV9fX3R5wMA2q9b1hA89NBDeeONN7Jw4cJ2lV+0aFHKy8ubt8rKyi6uIQD0bV0eCH7+859n3rx5+ed//ucMGNC+AYn58+dn9+7dzduOHTu6uJbQu02Zsn8DaEvRUwbFaGxszJVXXpmbb745Z5xxRruPKy0tTWlpaRfWDAD4pC4NBHv27MmmTZuydevWzJkzJ0nS1NSUQqGQAQMG5Mc//nG+8pWvdGUVAIB26NJAMGjQoGzbtq3FvnvvvTdPP/10Hn300Zx++uldeXoAoJ2KDgQffvhhXnvttebXb775ZmprazN48OB85jOfyfz58/POO+/k+9//fvr165dzzz23xfFDhw5NWVnZp/YDAD2n6ECwadOmTPnE6qSampokycyZM7NixYq8++672b59e+fVEADociWFQqHQ05U4nPr6+pSXl2f37t0ZNGhQT1cHjjkHMvy6dT1bD6D7tfc71LMMgCPikkboHQQCAEAgAAAEAgAgAgEAEIEAAIhAAABEIAAAIhAAABEIAIAIBNAnubsgcDCBAAAQCAAAgQAAiEAAAEQgAACSDOjpCgBdx5UEQHsZIQAABAIAQCAAACIQAAARCACACAQAQAQCACACAQAQgQAAiEAAAEQgAAAiEAAAEQigT5syxQOQgP0EAgBAIAAABAIAIAIBABCBAACIQAAARCAAACIQAAARCACACAQAQAQCACACAQAQgQAAiEAAAEQgAAAiEAAASQb0dAWAo8eUKR//ed261n928H6gdzBCAAAIBACAQAAARCAAACIQAADpQCDYsGFDpk+fnhEjRqSkpCSPP/74Icv/8Ic/zEUXXZRTTz01gwYNSlVVVZ566qmO1hcA6AJFB4K9e/dmzJgxWbp0abvKb9iwIRdddFGeeOKJbN68OVOmTMn06dOzdevWoisLAHSNou9DMG3atEybNq3d5ZcsWdLi9a233pof/ehH+fd///eMGzeu1WMaGhrS0NDQ/Lq+vr7YagIARej2NQRNTU3Zs2dPBg8e3GaZRYsWpby8vHmrrKzsxhoCQN/T7YHgjjvuyIcffpivfe1rbZaZP39+du/e3bzt2LGjG2sIAH1Pt966+OGHH87NN9+cH/3oRxk6dGib5UpLS1NaWtqNNQOAvq3bAsHKlStz9dVX55FHHkl1dXV3nRYAaIduCQT/+q//mm984xtZuXJlLr300u44JVCETz7UCOibig4EH374YV577bXm12+++WZqa2szePDgfOYzn8n8+fPzzjvv5Pvf/36S/dMEM2fOzN13351Jkyalrq4uSXL88cenvLy8kz4GAHAkil5UuGnTpowbN675ksGampqMGzcuCxYsSJK8++672b59e3P5+++/P7/5zW8ye/bsDB8+vHmbO3duJ30EAOBIFT1CMHny5BQKhTZ/vmLFihav169fX+wpgGPYJ6cf1q3ruXoAxfEsAwBAIAAABAIAIN18YyKg93LpIhzbjBAAAAIBAGDKAPoUw/pAW4wQAAACAQAgEABFmjLF1AP0RgIBACAQAAACAQAQgQAAiEAAAEQgAADiToXQKx1tlwUeqM+6dT1bD6BtRggAAIEAABAIAIAIBABABAIAIK4yADroaLuSATgyRggAAIEAABAIAIAIBEAbpkyxTgD6EoEAABAIAACXHQKHYdoA+gYjBACAQAAACAQAQAQCACACAQAQgQAAiEAAAEQgAAAiEAAAEQiAXsLDmODICAQAgEAAAAgEAEAEAgAgAgEAEIEAOEa4igC6lkAAAAgEAIBAAABEIAAAIhAAABEIAIB0IBBs2LAh06dPz4gRI1JSUpLHH3/8sMesX78+X/rSl1JaWprPf/7zWbFiRQeqCgB0laIDwd69ezNmzJgsXbq0XeXffPPNXHrppZkyZUpqa2tz/fXX5+qrr85TTz1VdGUBgK4xoNgDpk2blmnTprW7/LJly3L66afnzjvvTJKcffbZefbZZ3PXXXdl6tSprR7T0NCQhoaG5tf19fXFVhMAKEKXryHYuHFjqqurW+ybOnVqNm7c2OYxixYtSnl5efNWWVnZ1dUEjiIH7kp4qDsTunMhdK4uDwR1dXWpqKhosa+ioiL19fX5v//7v1aPmT9/fnbv3t287dixo6urCQB9WtFTBt2htLQ0paWlPV0NAOgzujwQDBs2LDt37myxb+fOnRk0aFCOP/74rj490MuYJoCu0eVTBlVVVVm7dm2LfWvWrElVVVVXnxoAaKeiA8GHH36Y2tra1NbWJtl/WWFtbW22b9+eZP/8/4wZM5rLX3vttXnjjTfyF3/xF3nllVdy77335gc/+EFuuOGGzvkEAMARKzoQbNq0KePGjcu4ceOSJDU1NRk3blwWLFiQJHn33Xebw0GSnH766Vm9enXWrFmTMWPG5M4778yDDz7Y5iWHAED3K3oNweTJk1MoFNr8eWt3IZw8eXK2bt1a7KkAgG7iWQYAgEAAAAgEAEAEAgAgAgEAEIEAOAocbQ8qOtrqA91BIAAABAIAQCAAACIQAAARCACACAQAQDrwcCOA7uTyP+geRggAAIEAABAIAIAIBABABAIAIAIBABCBAACIQAAARCAAACIQAN1oyhR3HoSjlUAAAAgEAICHGwG91CenJtat67l6wLHCCAEAIBAAAAIBABCBAACIQAAAxFUGQBc6Wm5CdKAerjaAthkhAAAEAgBAIAAAIhAAPcBDjuDoIxAAAAIBACAQAAARCACACAQAQNypEOhBnXGlgasVoHMYIQAABAIAQCAAAGINAXAUsR4Aeo4RAgBAIAAABAIAIAIBABCBAACIQAAApIOBYOnSpRk1alTKysoyadKkvPDCC4csv2TJkpx55pk5/vjjU1lZmRtuuCG/+tWvOlRhAKDzFR0IVq1alZqamixcuDBbtmzJmDFjMnXq1Lz33nutln/44Yczb968LFy4MC+//HK+973vZdWqVfnLv/zLI648ANA5ig4EixcvzjXXXJNZs2blnHPOybJly3LCCSdk+fLlrZZ//vnnc/755+fKK6/MqFGjcvHFF+eKK6445KhCQ0ND6uvrW2wAQNcpKhDs27cvmzdvTnV19cdv0K9fqqurs3HjxlaPOe+887J58+bmAPDGG2/kiSeeyCWXXNLmeRYtWpTy8vLmrbKysphqAgBFKurWxbt27UpjY2MqKipa7K+oqMgrr7zS6jFXXnlldu3ald/7vd9LoVDIb37zm1x77bWHnDKYP39+ampqml/X19cLBQDQhbr8KoP169fn1ltvzb333pstW7bkhz/8YVavXp1bbrmlzWNKS0szaNCgFhsA0HWKGiEYMmRI+vfvn507d7bYv3PnzgwbNqzVY2666aZcddVVufrqq5MkX/ziF7N379788R//cf7qr/4q/fq58hHoHgcenrRuXeuvoS8r6tt44MCBGT9+fNauXdu8r6mpKWvXrk1VVVWrx3z00Uef+tLv379/kqRQKBRbXwCgCxT9+OOamprMnDkzEyZMyMSJE7NkyZLs3bs3s2bNSpLMmDEjI0eOzKJFi5Ik06dPz+LFizNu3LhMmjQpr732Wm666aZMnz69ORgAAD2r6EBw+eWX5/3338+CBQtSV1eXsWPH5sknn2xeaLh9+/YWIwI33nhjSkpKcuONN+add97JqaeemunTp+fv/u7vOu9TQB9m2LvzHGhL6IuKDgRJMmfOnMyZM6fVn61fv77lCQYMyMKFC7Nw4cKOnAoA6AZW9AEAAgEAIBAAABEIAIAIBABABAIAIAIBABCBAACIQAAARCAAACIQAADp4LMMgJ7nQTxAZzJCAAAIBACAKQPoNUwhAEfCCAEAIBAAAAIBABBrCADa9Ml1GevW9Vw9oDsYIQAABAIAQCAAACIQAAARCACAuMoA6GXcsRE6xggBACAQAAACAQAQgQAAiEAAAEQgAAAiEAAAEQgAgAgEAEDcqRDogw6+m2Fn3N3wwHusW3fk7wU9wQgBACAQAACmDACOiIcp0VsYIQAABAIAQCAAACIQAAARCACACAQAQAQCACACAQAQgQAAiEAAAEQgAAAiEAAAEQgAgAgEAO0yZYonG9K7dSgQLF26NKNGjUpZWVkmTZqUF1544ZDlf/nLX2b27NkZPnx4SktLc8YZZ+SJJ57oUIUBgM43oNgDVq1alZqamixbtiyTJk3KkiVLMnXq1Lz66qsZOnTop8rv27cvF110UYYOHZpHH300I0eOzNtvv52TTz65M+oPAHSCogPB4sWLc80112TWrFlJkmXLlmX16tVZvnx55s2b96nyy5cvzwcffJDnn38+xx13XJJk1KhRhzxHQ0NDGhoaml/X19cXW02ALmX6gN6mqCmDffv2ZfPmzamurv74Dfr1S3V1dTZu3NjqMf/2b/+WqqqqzJ49OxUVFTn33HNz6623prGxsc3zLFq0KOXl5c1bZWVlMdUEAIpUVCDYtWtXGhsbU1FR0WJ/RUVF6urqWj3mjTfeyKOPPprGxsY88cQTuemmm3LnnXfmb//2b9s8z/z587N79+7mbceOHcVUEwAoUtFTBsVqamrK0KFDc//996d///4ZP3583nnnnXznO9/JwoULWz2mtLQ0paWlXV01AOD/KyoQDBkyJP3798/OnTtb7N+5c2eGDRvW6jHDhw/Pcccdl/79+zfvO/vss1NXV5d9+/Zl4MCBHag2ANCZipoyGDhwYMaPH5+1a9c272tqasratWtTVVXV6jHnn39+XnvttTQ1NTXv+9nPfpbhw4cLAwBwlCj6PgQ1NTV54IEH8o//+I95+eWXc91112Xv3r3NVx3MmDEj8+fPby5/3XXX5YMPPsjcuXPzs5/9LKtXr86tt96a2bNnd96nAOgmblBEb1X0GoLLL78877//fhYsWJC6urqMHTs2Tz75ZPNCw+3bt6dfv49zRmVlZZ566qnccMMNGT16dEaOHJm5c+fmW9/6Vud9CgDgiJQUCoVCT1ficOrr61NeXp7du3dn0KBBPV0dOCr4X+rRad26nq4BtNTe71DPMgAABAIAQCAAACIQAAARCOCY4XI3oCsJBACAQAAAdMPDjYDOZdoA6ApGCAAAgQAAEAgAgFhDAEc9awaA7mCEAAAQCAAAgQAAiEAAAEQgAAAiEECP8KCi3u9Qf8f+/jkaCQQAgEAAAAgEAEAEAgAgAgEAEIEAAIhAAABEIAAAIhAAABEIAIAIBABABAIAIMmAnq4A0DoPvzk2+XvjWGWEAAAQCAAAgQAAiEAAAEQgAAAiEAAAEQgAgAgEAEAEAgAgAgEAEIEAAIhAAABEIAAAIhAAABEIAIAIBABABAKAo8aUKfs36AkCAQAgEAAAAgEAEIEAAEgHA8HSpUszatSolJWVZdKkSXnhhRfaddzKlStTUlKSyy67rCOnBQC6SNGBYNWqVampqcnChQuzZcuWjBkzJlOnTs177713yOPeeuut/Nmf/Vm+/OUvd7iyAEDXKDoQLF68ONdcc01mzZqVc845J8uWLcsJJ5yQ5cuXt3lMY2Njvv71r+fmm2/Ob//2bx/2HA0NDamvr2+xAQBdp6hAsG/fvmzevDnV1dUfv0G/fqmurs7GjRvbPO5v/uZvMnTo0Hzzm99s13kWLVqU8vLy5q2ysrKYagIARSoqEOzatSuNjY2pqKhosb+ioiJ1dXWtHvPss8/me9/7Xh544IF2n2f+/PnZvXt387Zjx45iqgkAFGlAV775nj17ctVVV+WBBx7IkCFD2n1caWlpSktLu7BmcPRypzqgJxQVCIYMGZL+/ftn586dLfbv3Lkzw4YN+1T5119/PW+99VamT5/evK+pqWn/iQcMyKuvvprPfe5zHak3ANCJipoyGDhwYMaPH5+1a9c272tqasratWtTVVX1qfJnnXVWtm3bltra2ubtq1/9aqZMmZLa2lprAwDgKFH0lEFNTU1mzpyZCRMmZOLEiVmyZEn27t2bWbNmJUlmzJiRkSNHZtGiRSkrK8u5557b4viTTz45ST61H6C3Mx3E0azoQHD55Zfn/fffz4IFC1JXV5exY8fmySefbF5ouH379vTr5waIAHAs6dCiwjlz5mTOnDmt/mz9+vWHPHbFihUdOSUA0IX8Vx4AEAgAAIEAAEgX35gI2O/A6vJ163q2HnQ/VxZwrDBCAAAIBACAQAAARCAAOGpNmWINAt1HIAAABAIAQCAAACIQAAARCACAuFMhHBWsJKej3AWTzmKEAAAQCAAAUwbQo0wVAEcLIwQAgEAAAAgEAECsIYAuZY0A7aGfcDQwQgAACAQAgEAAAEQgAAAiEAAAEQgAgAgEAEAEAgAgAgEAEHcqBOgxXXmHwgPvvW5d152D3sUIAQAgEAAApgygS3hYDXCsMUIAAAgEAIBAAHBMmTLFlBRdQyAAAAQCAMBVBtCtDPXSHu3pJ22V0cfoKCMEAIBAAAAIBABABAL6KJdugX8HtCQQAAACAQAgEAAAEQgAgAgEAEDcqRCO2CdXaa9b13P1oG9ztQBHyggBACAQAAAdDARLly7NqFGjUlZWlkmTJuWFF15os+wDDzyQL3/5yznllFNyyimnpLq6+pDlAYDuV/QaglWrVqWmpibLli3LpEmTsmTJkkydOjWvvvpqhg4d+qny69evzxVXXJHzzjsvZWVluf3223PxxRfnpZdeysiRIzvlQwD0Zp25PsBaA9pS9AjB4sWLc80112TWrFk555xzsmzZspxwwglZvnx5q+X/5V/+JX/6p3+asWPH5qyzzsqDDz6YpqamrF27ts1zNDQ0pL6+vsUGAHSdogLBvn37snnz5lRXV3/8Bv36pbq6Ohs3bmzXe3z00Uf59a9/ncGDB7dZZtGiRSkvL2/eKisri6kmAFCkogLBrl270tjYmIqKihb7KyoqUldX1673+Na3vpURI0a0CBUHmz9/fnbv3t287dixo5hqQrPufniLh8UAx6puvQ/BbbfdlpUrV2b9+vUpKytrs1xpaWlKS0u7sWYA0LcVFQiGDBmS/v37Z+fOnS3279y5M8OGDTvksXfccUduu+22/OQnP8no0aOLrykA0GWKmjIYOHBgxo8f32JB4IEFglVVVW0e9+1vfzu33HJLnnzyyUyYMKHjtYVuYNgf6IuKnjKoqanJzJkzM2HChEycODFLlizJ3r17M2vWrCTJjBkzMnLkyCxatChJcvvtt2fBggV5+OGHM2rUqOa1BieeeGJOPPHETvwoAEBHFR0ILr/88rz//vtZsGBB6urqMnbs2Dz55JPNCw23b9+efv0+Hnj47ne/m3379uUP//APW7zPwoUL89d//ddHVnsAoFN0aFHhnDlzMmfOnFZ/tn79+hav33rrrY6cAgDoRp5lAAAIBACAQAAARCCgj+jMSwldlsixpJj+qm/3bQIBACAQAAACAQAQgQAAiEAAAKSbH38Mx5IDq63XrWv5Gvqag/8t0DsZIQAABAIAQCAAAGINAX1MZ64DsKaA3upI1gxYb3DsMkIAAAgEAIBAAABEIAAAIhAAAHGVATRz1QC9WUf698HHuIKgdzNCAAAIBACAKQP6uPYMo5pKgCN3qIeFmYI4OhghAAAEAgBAIAAAIhAAABEIAIAIBABAXHbIMeRQd0lzaSB0H5cM9k5GCAAAgQAAMGXAUcyDVODYYdru2GeEAAAQCAAAgQAAiDUEHAWsFYCe1ZXz/9YWHDuMEAAAAgEAYMqAY5ihSOh5/h32HkYIAACBAAAwZUA3cSUB9C1tTSW0tr+t3w+H+71RzEOW2vM76OAyfe33lhECAEAgAABMGdAO3T1sdvCQ4uFeA31DR34X+X3RfkYIAACBAAAQCACAWEPQ6x1qzq0z1gYc7jKdtub/D1UfoG9q76WKh7p0saPn6CzF/l49mi5tNEIAAAgEAEAHpwyWLl2a73znO6mrq8uYMWPyD//wD5k4cWKb5R955JHcdNNNeeutt/KFL3wht99+ey655JIOV7qzdPTuWN2lI3fWOnh/R85zuGG79rTL4c5vegA4WrX3UufO/I4o5q6LXaXoEYJVq1alpqYmCxcuzJYtWzJmzJhMnTo17733Xqvln3/++VxxxRX55je/ma1bt+ayyy7LZZddlp/+9KdHXHkAoHMUPUKwePHiXHPNNZk1a1aSZNmyZVm9enWWL1+eefPmfar83Xffnd///d/Pn//5nydJbrnllqxZsyb33HNPli1b1uo5Ghoa0tDQ0Px69+7dSZL6+vpiq3tIv/lN/v/7tm9/d2tPPQ73GQ5o7T3aKnPw/oN98r0OPv/hjgXoLQ71e7W93x+t/c7s7O+eA9+dhULh0AULRWhoaCj079+/8Nhjj7XYP2PGjMJXv/rVVo+prKws3HXXXS32LViwoDB69Og2z7Nw4cJCEpvNZrPZbJ207dix45Df8UWNEOzatSuNjY2pqKhosb+ioiKvvPJKq8fU1dW1Wr6urq7N88yfPz81NTXNr5uamvLBBx/kt37rt1JSUlJMlXtEfX19Kisrs2PHjgwaNKinq9OjtMV+2uFj2mI/7bCfdvhYV7VFoVDInj17MmLEiEOWOyrvQ1BaWprS0tIW+04++eSeqcwRGDRoUJ/v4Adoi/20w8e0xX7aYT/t8LGuaIvy8vLDlilqUeGQIUPSv3//7Ny5s8X+nTt3ZtiwYa0eM2zYsKLKAwDdr6hAMHDgwIwfPz5r165t3tfU1JS1a9emqqqq1WOqqqpalE+SNWvWtFkeAOh+RU8Z1NTUZObMmZkwYUImTpyYJUuWZO/evc1XHcyYMSMjR47MokWLkiRz587NBRdckDvvvDOXXnppVq5cmU2bNuX+++/v3E9yFCktLc3ChQs/Ne3RF2mL/bTDx7TFftphP+3wsZ5ui5JC4XDXIXzaPffc03xjorFjx+bv//7vM2nSpCTJ5MmTM2rUqKxYsaK5/COPPJIbb7yx+cZE3/72t4+KGxMBAPt1KBAAAL2LZxkAAAIBACAQAAARCACACAQdsmHDhkyfPj0jRoxISUlJHn/88UOWX79+fUpKSj61Her2zceCRYsW5Xd/93dz0kknZejQobnsssvy6quvHva4Rx55JGeddVbKysryxS9+MU888UQ31LbrdKQdVqxY8an+UFZW1k017jrf/e53M3r06OY7rVVVVeU///M/D3lMb+sPSfHt0Fv7w8Fuu+22lJSU5Prrrz9kud7YJw7Wnrbo7n4hEHTA3r17M2bMmCxdurSo41599dW8++67zdvQoUO7qIbd45lnnsns2bPzX//1X1mzZk1+/etf5+KLL87evXvbPKY3Pg67I+2Q7L896Sf7w9tvv91NNe46p512Wm677bZs3rw5mzZtyle+8pX8wR/8QV566aVWy/fG/pAU3w5J7+wPn/Tiiy/mvvvuy+jRow9Zrrf2iU9qb1sk3dwvDveEQw4tyaee/niwdevWFZIU/vd//7db6tRT3nvvvUKSwjPPPNNmma997WuFSy+9tMW+SZMmFf7kT/6kq6vXbdrTDg899FChvLy8+yrVg0455ZTCgw8+2OrP+kJ/OOBQ7dDb+8OePXsKX/jCFwpr1qwpXHDBBYW5c+e2Wba394li2qK7+4URgm40duzYDB8+PBdddFGee+65nq5Op9u9e3eSZPDgwW2W2bhxY6qrq1vsmzp1ajZu3NildetO7WmHJPnwww/z2c9+NpWVlYf93+OxqLGxMStXrszevXvbvFV5X+gP7WmHpHf3h9mzZ+fSSy/91N91a3p7nyimLZLu7RdH5dMOe5vhw4dn2bJlmTBhQhoaGvLggw9m8uTJ+e///u986Utf6unqdYqmpqZcf/31Of/883Puuee2Wa4jj8M+lrS3Hc4888wsX748o0ePzu7du3PHHXfkvPPOy0svvZTTTjutG2vc+bZt25aqqqr86le/yoknnpjHHnss55xzTqtle3N/KKYdenN/WLlyZbZs2ZIXX3yxXeV7c58oti26u18IBN3gzDPPzJlnntn8+rzzzsvrr7+eu+66K//0T//UgzXrPLNnz85Pf/rTPPvssz1dlR7V3naoqqpq8b/F8847L2effXbuu+++3HLLLV1dzS515plnpra2Nrt3786jjz6amTNn5plnnmnzy7C3KqYdemt/2LFjR+bOnZs1a9b0ykWSxehIW3R3vxAIesjEiRN7zZfnnDlz8h//8R/ZsGHDYVNrb34cdjHtcLDjjjsu48aNy2uvvdZFtes+AwcOzOc///kkyfjx4/Piiy/m7rvvzn333fepsr25PxTTDgfrLf1h8+bNee+991qMhDY2NmbDhg2555570tDQkP79+7c4prf2iY60xcG6ul9YQ9BDamtrM3z48J6uxhEpFAqZM2dOHnvssTz99NM5/fTTD3tMb3wcdkfa4WCNjY3Ztm3bMd8nWtPU1JSGhoZWf9Yb+0NbDtUOB+st/eHCCy/Mtm3bUltb27xNmDAhX//611NbW9vqF2Bv7RMdaYuDdXm/6Lbli73Inj17Clu3bi1s3bq1kKSwePHiwtatWwtvv/12oVAoFObNm1e46qqrmsvfddddhccff7zw85//vLBt27bC3LlzC/369Sv85Cc/6amP0Cmuu+66Qnl5eWH9+vWFd999t3n76KOPmstcddVVhXnz5jW/fu655woDBgwo3HHHHYWXX365sHDhwsJxxx1X2LZtW098hE7RkXa4+eabC0899VTh9ddfL2zevLnwR3/0R4WysrLCSy+91BMfodPMmzev8MwzzxTefPPNwv/8z/8U5s2bVygpKSn8+Mc/LhQKfaM/FArFt0Nv7Q+tOXhlfV/pE605XFt0d78QCDrgwGWEB28zZ84sFAqFwsyZMwsXXHBBc/nbb7+98LnPfa5QVlZWGDx4cGHy5MmFp59+umcq34laa4MkhYceeqi5zAUXXNDcLgf84Ac/KJxxxhmFgQMHFn7nd36nsHr16u6teCfrSDtcf/31hc985jOFgQMHFioqKgqXXHJJYcuWLd1f+U72jW98o/DZz362MHDgwMKpp55auPDCC5u/BAuFvtEfCoXi26G39ofWHPwl2Ff6RGsO1xbd3S88/hgAsIYAABAIAIAIBABABAIAIAIBABCBAACIQAAARCAAACIQAAARCACACAQAQJL/B4HCRI0jf/ThAAAAAElFTkSuQmCC\n"},"metadata":{}},{"output_type":"display_data","data":{"text/plain":"<Figure size 600x600 with 1 Axes>","image/png":"iVBORw0KGgoAAAANSUhEUgAAAgQAAAH5CAYAAAD+5ibMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAgC0lEQVR4nO3df3DX9WHH8VcIJdhJohQJgrHYdWqt5Uex5jLr1rRRRl06/9jNU08YV93Zoz01102zKpTZGtdViztRqqtlvR2F6lXdpsOxDMo58RQod3ZXba06mJIo85pAuoY2yf5wpkYC5AtJvvnxeNx9rvf98Pl8v+/0a5Jn3t/Pj5Kenp6eAADj2oRiDwAAKD5BAAAIAgBAEAAAEQQAQAQBABBBAAAkmVjsAQxEd3d3XnvttUyZMiUlJSXFHg4AjBo9PT3Zv39/Zs6cmQkTDj8PMCqC4LXXXktVVVWxhwEAo9aePXty2mmnHfbfR0UQTJkyJclbX0x5eXmRRwMAo0d7e3uqqqp6f5cezqgIgrc/JigvLxcEAHAMjvaRu4MKAQBBAAAIAgAgggAAiCAAACIIAIAIAgAgggAAiCAAACIIAIAIAgAgggAAiCAAACIIAIAIAgAgggAAiCAAACIIAIAIAgAggoARoLb2rQWA4hEEAIAgAAAEAQAQQQAARBAAABEEAEAEAQAQQQAARBAAABEEAEAEAQAQQQAARBAAABEEjBHumAhwfAQBACAIAABBAABEEAAAEQQAQAQBABBBAABEEAAAEQQAQAQBAJBjCIKtW7emvr4+M2fOTElJSR555JEB7/sf//EfmThxYubNm1foywIAQ6jgIOjo6MjcuXOzevXqgvb7+c9/nsWLF+dTn/pUoS8JAAyxiYXusGjRoixatKjgF7r22mtzxRVXpLS0tKBZBQBg6A3LMQTf/va389JLL2XFihUD2r6zszPt7e19FgBg6Ax5EPz0pz/NTTfdlH/4h3/IxIkDm5BoampKRUVF71JVVTXEowSA8W1Ig6CrqytXXHFFVq5cmTPPPHPA+zU2Nqatra132bNnzxCOEgAo+BiCQuzfvz/bt2/PD3/4w3z+859PknR3d6enpycTJ07Mv/7rv+aTn/zkIfuVlZWlrKxsKIcGALzDkAZBeXl5nnvuuT7r7rnnnvz7v/97HnrooZxxxhlD+fIAwAAVHAQHDhzIiy++2Pv45Zdfzq5duzJ16tScfvrpaWxszKuvvprvfOc7mTBhQs4999w++0+fPj2TJ08+ZD0AUDwFB8H27dtTW1vb+7ihoSFJsmTJkqxduzZ79+7N7t27B2+EAMCQK+np6ekp9iCOpr29PRUVFWlra0t5eXmxh8Mge7svN28u7nMAjEUD/R3qXgYAgCAAAAQBABBBAABEEAAAEQSMQLW1vzlrAIDhIQgAAEHA2GJ2AeDYCAIAQBAAAIIAAIggAAAiCACACAIAIIIAAIggAAAiCACACAIAIIIAAIggAAAiCBjB3KgIYPgIAgBAEAAAgoBRwEcHAENPEAAAggAAEAQAQAQBABBBAABEEAAAEQQAQAQBABBBAABEEAAAEQQAQAQBABBBAABEEAAAEQQAQAQBABBBAABEEAAAEQQAQAQBABBBAABEEAAAEQQAQAQBAJBjCIKtW7emvr4+M2fOTElJSR555JEjbv/9738/F110UU455ZSUl5enpqYmTzzxxLGOFwAYAgUHQUdHR+bOnZvVq1cPaPutW7fmoosuyuOPP54dO3aktrY29fX1+eEPf1jwYAGAoTGx0B0WLVqURYsWDXj7VatW9Xl822235dFHH80//dM/Zf78+f3u09nZmc7Ozt7H7e3thQ6TcaC2ttgjABg7hv0Ygu7u7uzfvz9Tp0497DZNTU2pqKjoXaqqqoZxhAAw/gx7EHz961/PgQMH8id/8ieH3aaxsTFtbW29y549e4ZxhAAw/hT8kcHxWLduXVauXJlHH30006dPP+x2ZWVlKSsrG8aRAcD4NmxBsH79+lx99dV58MEHU1dXN1wvCwAMwLB8ZPDd7343S5cuzXe/+91ccsklw/GSAEABCp4hOHDgQF588cXexy+//HJ27dqVqVOn5vTTT09jY2NeffXVfOc730ny1scES5YsyV133ZXq6uq0tLQkSU444YRUVFQM0pcBAByPgmcItm/fnvnz5/eeMtjQ0JD58+dn+fLlSZK9e/dm9+7dvdvfd999+fWvf51ly5bl1FNP7V2uu+66QfoSAIDjVfAMwSc+8Yn09PQc9t/Xrl3b5/GWLVsKfQkAYJgN61kGcCSDeaGht59r8+bBe06AsczNjQAAQQAACAIAII4hYBRyUyOAwWeGAAAwQ8DoYWYAYOiYIQAABAEAIAgAgAgCACCCgDGuttbBiAADIQgAAEEAAAgCACCCAACIIAAAIggAgAgCACCCAACIIAAAIggAgAgCACCCAACIIAAAIggAgAgCACCCAACIIAAAIggAgAgCACCCAACIIAAAIggAgAgCACCCAACIIAAAIggAgAgCxona2rcWAPonCAAAQQAACAIAIIIAAIggAAAiCACACAIAIIIAAMgxBMHWrVtTX1+fmTNnpqSkJI888shR99myZUs++tGPpqysLB/84Aezdu3aYxgqADBUCg6Cjo6OzJ07N6tXrx7Q9i+//HIuueSS1NbWZteuXbn++utz9dVX54knnih4sADA0JhY6A6LFi3KokWLBrz9mjVrcsYZZ+SOO+5IknzoQx/Kk08+mW984xtZuHBhoS8PAAyBIT+GYNu2bamrq+uzbuHChdm2bdth9+ns7Ex7e3ufBQAYOgXPEBSqpaUllZWVfdZVVlamvb09//u//5sTTjjhkH2ampqycuXKoR4aReRGQwAjy4g8y6CxsTFtbW29y549e4o9JAAY04Z8hmDGjBlpbW3ts661tTXl5eX9zg4kSVlZWcrKyoZ6aADA/xvyGYKampo0Nzf3Wbdp06bU1NQM9UsDAANUcBAcOHAgu3btyq5du5K8dVrhrl27snv37iRvTfcvXry4d/trr702L730Uv7iL/4izz//fO65555873vfyw033DA4XwEAcNwKDoLt27dn/vz5mT9/fpKkoaEh8+fPz/Lly5Mke/fu7Y2DJDnjjDPy2GOPZdOmTZk7d27uuOOO/N3f/Z1TDgFgBCnp6enpKfYgjqa9vT0VFRVpa2tLeXl5sYfDICjWWQabNxfndQGKZaC/Q0fkWQYAwPASBACAIIBC1da6sBIw9ggCAEAQAACCAACIIAAAIggAgAgCACCCAACIIAAAIggAgAgCACCCAI7KpYqB8UAQAACCAAAQBABABAEAEEEAAEQQAABJJhZ7AIwvTt8DGJnMEAAAgoDxxUWGAPonCAAAQQAACAIAIIIAAIggAAAiCACACAIAIIIAAIgggEMM9OJFLnIEjCWCAAAQBACAIAAAIggAgAgCACCCAABIMrHYA4CRwimEwHhmhgAAMEMAA2UGARjLzBAAAIIAABAEAEAEAQAQQQAARBAAABEEAECOMQhWr16d2bNnZ/Lkyamurs4zzzxzxO1XrVqVs846KyeccEKqqqpyww035Je//OUxDRgAGHwFB8GGDRvS0NCQFStWZOfOnZk7d24WLlyY119/vd/t161bl5tuuikrVqzIj3/843zrW9/Khg0b8pd/+ZfHPXgAYHAUHAR33nlnrrnmmixdujTnnHNO1qxZk/e+97154IEH+t3+qaeeygUXXJArrrgis2fPzsUXX5zLL7/8qLMKAMDwKSgIDh48mB07dqSuru43TzBhQurq6rJt27Z+9/nd3/3d7NixozcAXnrppTz++OP59Kc/fdjX6ezsTHt7e58FABg6Bd3LYN++fenq6kplZWWf9ZWVlXn++ef73eeKK67Ivn378vGPfzw9PT359a9/nWuvvfaIHxk0NTVl5cqVhQwNADgOQ36WwZYtW3Lbbbflnnvuyc6dO/P9738/jz32WG699dbD7tPY2Ji2trbeZc+ePUM9TAAY1wqaIZg2bVpKS0vT2traZ31ra2tmzJjR7z633HJLrrrqqlx99dVJko985CPp6OjIn/3Zn+VLX/pSJkw4tEnKyspSVlZWyNAAgONQ0AzBpEmTsmDBgjQ3N/eu6+7uTnNzc2pqavrd5xe/+MUhv/RLS0uTJD09PYWOFwAYAgXNECRJQ0NDlixZkvPOOy/nn39+Vq1alY6OjixdujRJsnjx4syaNStNTU1Jkvr6+tx5552ZP39+qqur8+KLL+aWW25JfX19bxgwdtXWvvW/mzcXdxwAHFnBQXDZZZfljTfeyPLly9PS0pJ58+Zl48aNvQca7t69u8+MwM0335ySkpLcfPPNefXVV3PKKaekvr4+X/3qVwfvqwAAjktJzyiYt29vb09FRUXa2tpSXl5e7OFQgHfPELz9uNj6m7F499gGOmazH8BINtDfoe5lAAAIAgBAEAAAEQQAQI7hLAMYS0bKQY4AxWaGAAAQBACAIAAAIggAgAgCACCCAACI0w4Zp5xuCNCXGQIAwAwBHI5ZBGA8MUMAAAgCAEAQAAARBABABAEAEEEAAEQQAAARBABABAEAEEEAAEQQAAARBABABAEAEEEAAEQQAAARBABABAEAEEEAAEQQAAARBABABAEAEEEAAEQQAAARBABABAEAEEEAAEQQAAARBABABAEMidratxaA0UIQAACCAAAQBABAkonFHgDjg8/TAUY2MwQAwLEFwerVqzN79uxMnjw51dXVeeaZZ464/c9//vMsW7Ysp556asrKynLmmWfm8ccfP6YBAwCDr+CPDDZs2JCGhoasWbMm1dXVWbVqVRYuXJgXXngh06dPP2T7gwcP5qKLLsr06dPz0EMPZdasWfmv//qvnHTSSYMxfhjV3v4oZfPm4o4DoOAguPPOO3PNNddk6dKlSZI1a9bkscceywMPPJCbbrrpkO0feOCBvPnmm3nqqafynve8J0kye/bs4xs1ADCoCvrI4ODBg9mxY0fq6up+8wQTJqSuri7btm3rd59//Md/TE1NTZYtW5bKysqce+65ue2229LV1XXY1+ns7Ex7e3ufBUaqty9CdKQDJ12oCBjpCgqCffv2paurK5WVlX3WV1ZWpqWlpd99XnrppTz00EPp6urK448/nltuuSV33HFHvvKVrxz2dZqamlJRUdG7VFVVFTJMAKBAQ36WQXd3d6ZPn5777rsvCxYsyGWXXZYvfelLWbNmzWH3aWxsTFtbW++yZ8+eoR4mDAszBcBIVdAxBNOmTUtpaWlaW1v7rG9tbc2MGTP63efUU0/Ne97znpSWlvau+9CHPpSWlpYcPHgwkyZNOmSfsrKylJWVFTI0AOA4FDRDMGnSpCxYsCDNzc2967q7u9Pc3Jyampp+97ngggvy4osvpru7u3fdT37yk5x66qn9xgAAMPwK/sigoaEh999/f/7+7/8+P/7xj/O5z30uHR0dvWcdLF68OI2Njb3bf+5zn8ubb76Z6667Lj/5yU/y2GOP5bbbbsuyZcsG76sAAI5LwacdXnbZZXnjjTeyfPnytLS0ZN68edm4cWPvgYa7d+/OhAm/6Yyqqqo88cQTueGGGzJnzpzMmjUr1113XW688cbB+yoAgONS0tPT01PsQRxNe3t7Kioq0tbWlvLy8mIPhwKMtwPo3r7A0LsvOHS0/x9cmAgYKgP9HepeBgCAIAAABAEAEEEAAEQQAAA5htMOoT+FHlUPwMhihgAAEAQAgCAAACIIAIA4qBAGlYMpgdHKDAEAIAgAAEEAAEQQAAARBABABAEAEEEAAEQQAAARBABABAEAEEEAAEQQAABxcyMGmZv7AIxOZggAAEEAAAgCACCCAACIgwphSDnIEhgtzBAAAIIAABAEAEAEAQAQQQAARBAAAHHaIYwI7zw9cfPm4o0DGL/MEAAAggAAEAQAQBxDwHFyaV6AscEMAQAgCAAAQQAARBAAABEEAEAEAQAQQQAARBAAADnGIFi9enVmz56dyZMnp7q6Os8888yA9lu/fn1KSkpy6aWXHsvLAgBDpOAg2LBhQxoaGrJixYrs3Lkzc+fOzcKFC/P6668fcb9XXnklX/ziF3PhhRce82ABgKFRcBDceeedueaaa7J06dKcc845WbNmTd773vfmgQceOOw+XV1dufLKK7Ny5cp84AMfOK4BAwCDr6AgOHjwYHbs2JG6urrfPMGECamrq8u2bdsOu99f/dVfZfr06fnsZz87oNfp7OxMe3t7nwUAGDoFBcG+ffvS1dWVysrKPusrKyvT0tLS7z5PPvlkvvWtb+X+++8f8Os0NTWloqKid6mqqipkmABAgYb0bof79+/PVVddlfvvvz/Tpk0b8H6NjY1paGjofdze3i4KRhB3OAQYewoKgmnTpqW0tDStra191re2tmbGjBmHbP+zn/0sr7zySurr63vXdXd3v/XCEyfmhRdeyG//9m8fsl9ZWVnKysoKGRoAcBwK+shg0qRJWbBgQZqbm3vXdXd3p7m5OTU1NYdsf/bZZ+e5557Lrl27epfPfOYzqa2tza5du/zVDwAjRMEfGTQ0NGTJkiU577zzcv7552fVqlXp6OjI0qVLkySLFy/OrFmz0tTUlMmTJ+fcc8/ts/9JJ52UJIesBwCKp+AguOyyy/LGG29k+fLlaWlpybx587Jx48beAw13796dCRNcABGG2tvHcmzeXNxxAGNDSU9PT0+xB3E07e3tqaioSFtbW8rLy4s9nHHPQYVDa6C/4AUBMBAD/R3qT3kAQBAAAIIAAMgQX5gIKNy7j9FwjAAwHMwQAACCAAAQBABABAEAEEEAAMRZBjDquFIkMBTMEAAAggAA8JEBBTBVDTB2mSEAAAQBACAIAIAIAgAgggAAiCCAEa+21hkewNATBACAIAAABAEAEEEAAEQQAAARBABABAEAEEEAAMTtjxkAF8UZGQ73Pry9fvPm4RsLMPaYIQAABAEAIAgAgAgCACAOKoQx450HHTrAECiUGQIAQBAAAIIAAIgggHGlttaFpoD+CQIAwFkG49nRLnnrL0mA8cMMAQAgCAAAQQBjmoMIgYESBACAgwo5lL8oRz/vIVAoMwQAgCDA58wACAIAIIIAAMgxBsHq1asze/bsTJ48OdXV1XnmmWcOu+3999+fCy+8MCeffHJOPvnk1NXVHXF7AGD4FRwEGzZsSENDQ1asWJGdO3dm7ty5WbhwYV5//fV+t9+yZUsuv/zybN68Odu2bUtVVVUuvvjivPrqq8c9eABgcJT09PT0FLJDdXV1Pvaxj+Xuu+9OknR3d6eqqipf+MIXctNNNx11/66urpx88sm5++67s3jx4n636ezsTGdnZ+/j9vb2VFVVpa2tLeXl5YUMl3dx8OD49Pb9Ko52/wpg7Glvb09FRcVRf4cWNENw8ODB7NixI3V1db95ggkTUldXl23btg3oOX7xi1/kV7/6VaZOnXrYbZqamlJRUdG7VFVVFTJMAKBABQXBvn370tXVlcrKyj7rKysr09LSMqDnuPHGGzNz5sw+UfFujY2NaWtr61327NlTyDABgAIN65UKb7/99qxfvz5btmzJ5MmTD7tdWVlZysrKhnFkADC+FRQE06ZNS2lpaVpbW/usb21tzYwZM46479e//vXcfvvt+bd/+7fMmTOn8JECg+bdx5I4pgAo6CODSZMmZcGCBWlubu5d193dnebm5tTU1Bx2v6997Wu59dZbs3Hjxpx33nnHPloAYEgU/JFBQ0NDlixZkvPOOy/nn39+Vq1alY6OjixdujRJsnjx4syaNStNTU1Jkr/+67/O8uXLs27dusyePbv3WIMTTzwxJ5544iB+KcCxcvYBUHAQXHbZZXnjjTeyfPnytLS0ZN68edm4cWPvgYa7d+/OhAm/mXi49957c/DgwfzxH/9xn+dZsWJFvvzlLx/f6AGAQVHwdQiKYaDnUHJ0rkMwPr37OgRH2w4YO4bkOgSMHu5gCEAhBAEAIAiAwph9grFJEAAAw3ulQmD0MisAY5sZAgBAEAAAggAAiCAAAOKgQuAd3nngoKsWwvhihgAAMEMA48GxnDJ4tDsgukMijC1mCAAAQQAACAIAIIIAAIggAAAiCACACIJR72j3pnfvegAGQhAAAIJgrDATwFjgv2MoHkEAALh08VhzuL+u/NXFUBnqSxi7RDIMDzMEAIAgAAAEAQAQQQAAxEGFI5qDqRhNjue/Vwe9QvGZIQAABAEw/FyACEYeQQAAOIYAOLLj+UvecTAwepghAADMEIwG/spiNClkRsGltmHkMEMAAAgCAMBHBiPCsXwkYEqV0WQw/nt99/fJ0R4X8lyDMR4Y7cwQAABmCEYSf3HA4HnnrMRAv6fePZPxzv18fzLWmSEAAMwQjCaOG4Cx9X1g1oGRxAwBAGCGYKi4FSyMHIf7fjza9+lwfy8e7vVGwwyC2Y7RzwwBACAIAAAfGQy6d0/59TcFaEoNBs9AvueO9m9DeUfHYzn9cTTxEefYYYYAABjfMwQDPQimv8J30CCMTwOdZejvZ8Ox/twYyCzD8VzKudBxDNfloRlexzRDsHr16syePTuTJ09OdXV1nnnmmSNu/+CDD+bss8/O5MmT85GPfCSPP/74MQ0WABgaBc8QbNiwIQ0NDVmzZk2qq6uzatWqLFy4MC+88EKmT59+yPZPPfVULr/88jQ1NeUP//APs27dulx66aXZuXNnzj333EH5IgbbsdzPXQUDxVDIMRQD+ffjNdDnH+7jq47nJnIjbVxDpeAguPPOO3PNNddk6dKlSZI1a9bkscceywMPPJCbbrrpkO3vuuuu/MEf/EH+/M//PEly6623ZtOmTbn77ruzZs2afl+js7MznZ2dvY/b2tqSJO3t7YUO94h+/ev8//P2v74/b2/77m0Otx4Yn470M+Hd/1bIz49Cfl4dy/YD8c7nLPT1j7TfIP+I7/f1CnmNY9mnUMPxGm//7uzp6Tnyhj0F6Ozs7CktLe15+OGH+6xfvHhxz2c+85l+96mqqur5xje+0Wfd8uXLe+bMmXPY11mxYkVPEovFYrFYLIO07Nmz54i/4wuaIdi3b1+6urpSWVnZZ31lZWWef/75fvdpaWnpd/uWlpbDvk5jY2MaGhp6H3d3d+fNN9/M+973vpSUlBQy5BGlvb09VVVV2bNnT8rLy4s9nHHL+1B83oPi8x4U33C9Bz09Pdm/f39mzpx5xO1G5FkGZWVlKSsr67PupJNOKs5ghkB5eblvwBHA+1B83oPi8x4U33C8BxUVFUfdpqCzDKZNm5bS0tK0trb2Wd/a2poZM2b0u8+MGTMK2h4AGH4FBcGkSZOyYMGCNDc3967r7u5Oc3Nzampq+t2npqamz/ZJsmnTpsNuDwAMv4I/MmhoaMiSJUty3nnn5fzzz8+qVavS0dHRe9bB4sWLM2vWrDQ1NSVJrrvuuvz+7/9+7rjjjlxyySVZv359tm/fnvvuu29wv5JRoKysLCtWrDjk4xCGl/eh+LwHxec9KL6R9h6U9PQc7TyEQ9199935m7/5m7S0tGTevHn527/921RXVydJPvGJT2T27NlZu3Zt7/YPPvhgbr755rzyyiv5nd/5nXzta1/Lpz/96UH7IgCA43NMQQAAjC1ubgQACAIAQBAAABEEAEAEwbDYunVr6uvrM3PmzJSUlOSRRx4p9pDGnaampnzsYx/LlClTMn369Fx66aV54YUXij2scefee+/NnDlzeq/MVlNTk3/5l38p9rDGrdtvvz0lJSW5/vrriz2UceXLX/5ySkpK+ixnn312sYclCIZDR0dH5s6dm9WrVxd7KOPWD37wgyxbtixPP/10Nm3alF/96le5+OKL09HRUeyhjSunnXZabr/99uzYsSPbt2/PJz/5yfzRH/1R/vM//7PYQxt3nn322Xzzm9/MnDlzij2UcenDH/5w9u7d27s8+eSTxR7SyLyXwVizaNGiLFq0qNjDGNc2btzY5/HatWszffr07NixI7/3e79XpFGNP/X19X0ef/WrX829996bp59+Oh/+8IeLNKrx58CBA7nyyitz//335ytf+UqxhzMuTZw4ccRdwt8MAeNSW1tbkmTq1KlFHsn41dXVlfXr16ejo8OlzIfZsmXLcskll6Surq7YQxm3fvrTn2bmzJn5wAc+kCuvvDK7d+8u9pDMEDD+dHd35/rrr88FF1yQc889t9jDGXeee+651NTU5Je//GVOPPHEPPzwwznnnHOKPaxxY/369dm5c2eeffbZYg9l3Kqurs7atWtz1llnZe/evVm5cmUuvPDC/OhHP8qUKVOKNi5BwLizbNmy/OhHPxoRn9mNR2eddVZ27dqVtra2PPTQQ1myZEl+8IMfiIJhsGfPnlx33XXZtGlTJk+eXOzhjFvv/Ah5zpw5qa6uzvvf//5873vfy2c/+9mijUsQMK58/vOfzz//8z9n69atOe2004o9nHFp0qRJ+eAHP5gkWbBgQZ599tncdddd+eY3v1nkkY19O3bsyOuvv56PfvSjveu6urqydevW3H333ens7ExpaWkRRzg+nXTSSTnzzDPz4osvFnUcgoBxoaenJ1/4whfy8MMPZ8uWLTnjjDOKPST+X3d3dzo7O4s9jHHhU5/6VJ577rk+65YuXZqzzz47N954oxgokgMHDuRnP/tZrrrqqqKOQxAMgwMHDvQpv5dffjm7du3K1KlTc/rppxdxZOPHsmXLsm7dujz66KOZMmVKWlpakiQVFRU54YQTijy68aOxsTGLFi3K6aefnv3792fdunXZsmVLnnjiiWIPbVyYMmXKIcfN/NZv/Vbe9773OZ5mGH3xi19MfX193v/+9+e1117LihUrUlpamssvv7yo4xIEw2D79u2pra3tfdzQ0JAkWbJkSZ/bRDN07r333iRv3Z77nb797W/nT//0T4d/QOPU66+/nsWLF2fv3r2pqKjInDlz8sQTT+Siiy4q9tBg2Pz3f/93Lr/88vzP//xPTjnllHz84x/P008/nVNOOaWo43L7YwDAdQgAAEEAAEQQAAARBABABAEAEEEAAEQQAAARBABABAEAEEEAAEQQAABJ/g+ojLxqko2algAAAABJRU5ErkJggg==\n"},"metadata":{}}],"id":"996f516e-cfcd-45a0-8089-c415d11b1ba1"},{"cell_type":"markdown","source":"To prevent plotting, you should set do_plot=None. To print the results:","metadata":{},"id":"1aac1992-fec9-446f-9f97-b6c2f02b2b95"},{"cell_type":"code","source":"for i in range(0, output1_num):\n   print(output1[i], output1_error[i,:])","metadata":{"trusted":true},"execution_count":5,"outputs":[{"name":"stdout","text":"3.0 [-0.31710491  0.37039269]\n2.0 [-0.34198337  0.37857234]\n","output_type":"stream"}],"id":"e70b6e38-9959-4f0f-9ced-3c1d43ae1c5b"},{"cell_type":"markdown","source":"Now defined a function having two arguments:","metadata":{},"id":"d4c8febc-cce5-490a-971e-421e03a18de1"},{"cell_type":"code","source":"def myfunc22(input1, functargs=None):\n   result1 = functargs['scale1']*np.sum(input1)\n   result2 = functargs['scale2']*input1[1] ** input1[0]\n   return [result1, result2]\n\nfcnargs = {'scale1':0.0, 'scale2':0.0}\ninput1 = np.array([1., 2.])\ninput1_err = np.array([0.2, 0.5])\ninput1_err_p = input1_err\ninput1_err_m = -input1_err\nscale1=2.\nscale2=3.\n#fcnargs = fcnStruct(scale1, scale2)\nfcnargs['scale1']=scale1\nfcnargs['scale2']=scale2\noutput1 = myfunc22(input1, functargs=fcnargs)\noutput1_num = len(output1)","metadata":{"trusted":true},"execution_count":6,"outputs":[],"id":"be38644d-e08c-4f32-9f98-7dcacdf01d80"},{"cell_type":"code","source":"use_gaussian=0 # uniform distribution from min value to max value\nwalk_num=30 # number of walkers\niteration_num=100 # number of samplers\nmcmc_sim = pyemcee.hammer(myfunc22, input1, input1_err_m, input1_err_p, \n                          output1, walk_num, iteration_num, \n                          use_gaussian, print_progress=1, functargs=fcnargs)","metadata":{"trusted":true},"execution_count":7,"outputs":[{"name":"stdout","text":"Progress: 100% \n\n","output_type":"stream"}],"id":"16ce1c2c-5b5d-4754-b83a-95d7d1d700a9"},{"cell_type":"code","source":"output1_error = pyemcee.find_errors(output1, mcmc_sim, clevel)","metadata":{"trusted":true},"execution_count":8,"outputs":[],"id":"31e5d85b-4aa0-4176-87a4-8de8b83519d2"},{"cell_type":"code","source":"for i in range(0, output1_num):\n   print(output1[i], output1_error[i,:])","metadata":{"trusted":true},"execution_count":9,"outputs":[{"name":"stdout","text":"6.0 [-0.87131245  0.8385104 ]\n6.0 [-1.36297413  1.47110334]\n","output_type":"stream"}],"id":"5b1d3518-24cf-4069-a8b1-ae6283d87b56"},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[],"id":"65bc80d5-c091-4e64-aa07-78134b5d4bb5"}]}