random_multivariate.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "d6b63136",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sandy\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "from scipy.stats import norm\n",
    "from scipy.stats import lognorm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 451,
   "id": "e342afe8",
   "metadata": {},
   "outputs": [],
   "source": [
    "dummy_cov = np.asarray([[5/100,-4/70],[-4/70,3/49]])\n",
    "dummy_mean = np.array([10,7])\n",
    "dummy_center = np.ones(2)\n",
    "dummy_std = np.array([np.sqrt(5/100), np.sqrt(3/49)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 452,
   "id": "5676dc82",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.00180555,  0.11303004])"
      ]
     },
     "execution_count": 452,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.eig(dummy_cov)[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7bf2eadc",
   "metadata": {},
   "source": [
    "## Normal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 453,
   "id": "1605f7d6",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\fgrimald\\Anaconda3\\envs\\sandy-devel\\lib\\site-packages\\ipykernel_launcher.py:1: RuntimeWarning: covariance is not positive-semidefinite.\n",
      "  \"\"\"Entry point for launching an IPython kernel.\n"
     ]
    }
   ],
   "source": [
    "df = pd.DataFrame(np.random.multivariate_normal(dummy_center, dummy_cov, size=100000))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 454,
   "id": "53a16cb6",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\fgrimald\\Anaconda3\\envs\\sandy-devel\\lib\\site-packages\\seaborn\\distributions.py:2619: FutureWarning: `distplot` is a deprecated function and will be removed in a future version. Please adapt your code to use either `displot` (a figure-level function with similar flexibility) or `histplot` (an axes-level function for histograms).\n",
      "  warnings.warn(msg, FutureWarning)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "covariance:\n",
      "             0            1\n",
      "0  5.20672e-02 -5.53618e-02\n",
      "1 -5.53618e-02  6.28059e-02\n",
      "\n",
      "\n",
      "mean: [0.9998048670863129, 1.0001423669985525]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAApXklEQVR4nO3de1iU95338feXAeTkGVRUEBVQERVhRCzGQw5qbGNMn7aJadOm7dakbdrmuXpId/ts0t1un6fd7LPdZptuatI0aY26aZuDSTxgEhMPUWRQ8EQwqIigBlREAaPA/PYPBksIyCAD98w939d1zeXMfZj5MNf48fbHPb9bjDEopZSyrxCrAyillOpbWvRKKWVzWvRKKWVzWvRKKWVzWvRKKWVzoVYH6ExsbKxJSkqyOoZSSgWMwsLCs8aYuM7W+WXRJyUl4XK5rI6hlFIBQ0ROdLVOh26UUsrmtOiVUsrmtOiVUsrmtOiVUsrmtOiVUsrmuj3rRkSeBT4DVBtj0jtZ/0Pgi+2ebwoQZ4w5LyLlwCWgBWg2xjh9FVwppZR3vDmifw5Y0tVKY8zjxpgMY0wG8PfAu8aY8+02WehZryWvlFIW6LbojTHbgPPdbeexAljbq0RK2VRLSwtut9vqGCoI+ewLUyISReuR/0PtFhsgT0QM8DtjzKrr7L8SWAmQmJjoq1hKWer48eM8//zzbNq0ibq6OoYMGcL8+fN54IEHGD9+vNXxVJDw5Tdj7wB2dhi2yTXGnBKREcAWEXnf8z+ET/D8I7AKwOl06tVQVMB75ZVX+Kd/+ieMMSxatIjExEQqKip488032bRpEz/4wQ/40pe+ZHVMFQR8WfT30GHYxhhzyvNntYi8DGQDnRa9Unby3HPP8ctf/pLs7Gwef/xxRowYcW1ddXU1jz32GD//+c9paGjggQcesDCpCgY+KXoRGQzMB77Ublk0EGKMueS5vwj4Z1+8nlL+Zk1+xbX7B3dv5eWnfskU503c9sDf8+bxj+B4xce2z/3yjzh7JYT/+I//4IMLhswFS7l3tg5Zqr7hzemVa4EFQKyIVAKPAWEAxpinPJvdBeQZYxra7ToSeFlE2l5njTFmk++iK+V/zp2p4rVnf0ViajrLV/6IEIej0+1CQhzc+Y0f0HjxAptf+C/GTJwCWvSqj4g/Xhzc6XQanb1SBZI1+RW43S08/3+/z9lTJ3ng579j0NDYbverr6vl6Ue/RUR0DG9tfI3w8PB+SKvsSEQKuzqNXb8Zq5SPuN56jcqyEpbc922vSh4gZvBQPvPVhzl7qoLVq1f3cUIVrLTolfKByw2X2PbKasZPnUl6zsIe7ZuSMZuUjByefPJJqqur+yihCmZa9Er5wI7X1nK5sZ5b7/4Gnt9L9cjiex+kqamJ3/72t32QTgU7LXqleqmmpoaCN9czI/dWRiVOvKHnGDoins997nO89NJLVFVV+TihCnZa9Er10vPPP4+7pYW5d6zo1fN84xvfAODpp5/2RSylrvHLa8YqFSjq6upYt24dabPnMWzkmF4919aKJqbftJg//+WvjJ6znJghwwD0/HrVa3pEr1QvvPjiizQ0NJD76bt98nyzF30Wt7uFwq1v+OT5lAIteqVuWEtLC+vWrSMnJ4eRCRN88pzDR40hZfosCre+TnPTVZ88p1Ja9ErdoG3btnHq1ClWrOjd2HxH2bctp+HiBQ7tedenz6uClxa9UjdozZo1jBgxgptvvtmnzzt+aibDRo1l37sbffq8Knhp0St1A06dOsXOnTv5/Oc/T2iob89pEBEyblrEySOHOHem0qfPrYKTFr1SN2D9+vUYY7jrrrv65Pln5N6GhIRQtG1znzy/Ci5a9Er1kDGGV155hezsbMaM6d0plV2JGTKMlBnZ7N/5Js3NzX3yGip4aNEr1UNFRUWcOHGC5cuX9+nrzJi7iPq68+Tn5/fp6yj706JXqodef/11IiIiWLRoUZ++TvK0WQyIjOKNN/ScetU7WvRK9YDb7SYvL4958+YRHR3dp68VGh7O5KxctmzZwpUrV/r0tZS9adEr1QP79u3j7NmzfX403yY9ZyH19fW8+66eU69unM51o1QP/Ocf/4ojNIzq6Ikfu05sX0maksGwYcPYtGlTv/3jouxHj+iV8pIxhhLXDiZOy2JAZFS/vGaIw8HChQvZtm0bV6/qlAjqxmjRK+WlAwcOcPF8DVOcN/Xr695yyy00NDSwZ8+efn1dZR9a9Ep5afPmzYQ4QknNyOnX150zZw6RkZG89dZb/fq6yj66LXoReVZEqkXkYBfrF4hInYgUeW6Ptlu3RERKRaRMRH7sy+BK9SdjDHl5eUyYOpOI6Jh+fe2IiAjmzp3L22+/jdvt7tfXVvbgzRH9c8CSbrbZbozJ8Nz+GUBEHMCTwO1AGrBCRNJ6E1Ypq5SUlFBZWcnkrLmWvP4tt9xCdXU1hw4dsuT1VWDrtuiNMduA8zfw3NlAmTHmmDHmKrAOuPMGnkcpy+Xl5eFwOJiUOceS158/fz4Oh0OHb9QN8dUY/RwRKRaRjSIy1bNsDHCy3TaVnmWdEpGVIuISEVdNTY2PYinlG9u2bSMjI4OogYP7/bXX5FewofQiCanpvPT6pn45rVPZiy+Kfi8wzhgzA/hP4BXPculkW9PVkxhjVhljnMYYZ1xcnA9iKeUb1dXVlJSUMG/ePEtzTMr8FDVVJzj/YZWlOVTg6XXRG2MuGmPqPfc3AGEiEkvrEXxCu03HAqd6+3pK9bcdO3YAWF70qTNbh41K975naQ4VeHpd9CIySkTEcz/b85zngAIgRUTGi0g4cA+wvrevp1R/2759OyNGjGDSpEmW5hgSO5IRCeMp219gaQ4VeLqdAkFE1gILgFgRqQQeA8IAjDFPAZ8DvikizcBl4B5jjAGaReQhYDPgAJ41xugpAyqgNDc3895773HrrbfiOZ6x1MR0J/l5L9PQ0NDnk6op++i26I0x173ysTHmN8Bvuli3AdhwY9GUsl5xcTEXL160fNimzcRpWeza+GcKCgpYsGCB1XFUgNBJzZTqRNuZLW//5XUkJISq8ES/ONslIXkqoeED2Llzpxa98ppOgaDUdZTtLyAhZSoRUf4xTBIaHk7S5OnXfkGslDe06JXqwqUL5/iw4ijJ02ZZHeVjJqRnUV5eTlWVnmapvKNFr1QXjh5wATBxutPiJB83MT0LQI/qlde06JXqwtEDhcQMGcbIhAlWR/mY4fEJxMfHs3PnTqujqAChRa9UJ4zbTXlJEePTZvrFaZXtiQi5ubns3r2b5uZmq+OoAKBFr1QnqqvKabxUx/i0mVZH6dTcuXO5dOkS+/fvtzqKCgBa9Ep14vjhIgCSpsywNkgXcnJyCAkJ0XF65RUteqU6UX54H8NGjmHw8BFWR+nU4MGDSU9PZ9euXVZHUQFAi16pDpqamjhReoDxaRlWR7munJwcDh48SENDg9VRlJ/Toleqg4MHD3L1o8sk+en4fJvZs2fT3NzM3r17rY6i/JwWvVId7Nq1C0T8dny+zcyZMwkLC2P37t1WR1F+Tue6UaqD3bt3MypxIlExg6yO0qW2eXfix09i49YdjJl3NwD3zk60MpbyU3pEr1Q7ly9fpqioiPFTMqyO4pWkKRmcKS/jo4Z6q6MoP6ZFr1Q7e/fupampiSQ//0Vsm6QpMzDGzYkjB6yOovyYFr1S7ezatYvQ0FASJ02zOopXxkycTGhYOCdKiq2OovyYFr1S7RQUFDB9+nTCB0RYHcUroWHhJKSkcbykyOooyo9p0Svl0djYyOHDh3E6/Wu2yu4kTcmg+uRxGi/VWR1F+SkteqU8iouLaW5uJisry+ooPdJ2GuiJ93XeG9U5LXqlPFwuFyEhIcyc6d9flOooPimV8IhIynWcXnVBi14pj8LCQiZNmsTAgQOtjtIjjtBQElKmcqJUj+hV57otehF5VkSqReRgF+u/KCL7Pbf3RGRGu3XlInJARIpExOXL4Er5ypr8Cv608yiFe/cxcOwkv7gIeE8lTppGTdUJamtrrY6i/JA3R/TPAUuus/44MN8YMx34GbCqw/qFxpgMY0xg/YZLBZXT5R/Q3HSVxNR0q6PckHGe00ELCwstTqL8UbdFb4zZBpy/zvr3jDFthxG7gbE+yqZUvzlR2vqFo0A5f76j+KQUQsPCtehVp3w9Rv91YGO7xwbIE5FCEVl5vR1FZKWIuETEVVNT4+NYSl3fySMHGT5qLNGDhlgd5YaEhoUzZuJkXC4dIVWf5LOiF5GFtBb9I+0W5xpjMoHbgW+LyLyu9jfGrDLGOI0xzri4OF/FUqpbbncLFUcOBezRfJvE1HQOHz6s89OrT/BJ0YvIdOAZ4E5jzLm25caYU54/q4GXgWxfvJ5SvlRdWc6Vyw0BOz7fJnHSNNxuN/v27bM6ivIzvS56EUkEXgLuM8Ycabc8WkQGtt0HFgGdnrmjlJUqSls/loF+RD82OY3Q0FAdvlGf0O189CKyFlgAxIpIJfAYEAZgjHkKeBQYDvxWRACaPWfYjARe9iwLBdYYYzb1wc+gVK9UHDnAoGFxDIkdaXWUXgkfEEFaWpoWvfqEboveGLOim/V/B/xdJ8uPAf59iR4V9IwxVJQeZPzUwPo2bFecTid/+tOfuHLlCgMGDLA6jvIT+s1YFdROnDhBw8XagB+fb+N0OmlqamL/fv2WrPobLXoV1NqGOQJ9fL5NZmYmIqLDN+pjtOhVUCssLCRq4GBi4xOsjuITgwcPJjU1VYtefYwWvQpqLpeLhNR0PCcN2ILT6aSoqIjm5maroyg/oUWvgtaZM2eorKy0zfh8m6ysLBobGykpKbE6ivITWvQqaLXNCzPOJuPzbdounKLDN6qNFr0KWi6Xi6ioKEYmTLA6ik+NGDGCcePGadGra7ToVdByuVxkZmYS4nBYHcXnnE4nhYWFuN1uq6MoP6BFr4JSbW0tZWVlAXd9WG85nU7q6uooKyuzOoryA91+M1YpO2qb+MvpdHKkxeIwPtR2daza0NEAPPWXLThvieDe2YlWxlIW0yN6FZRcLhdhYWFMm2avX8S2GRI7koFDY6k4ovMIKi16FaQKCwuZPn26beeDERESU9OpKD2AMcbqOMpiWvQq6DQ0NHDo0CGcTntfxjgxNZ1LF85xoeaM1VGUxbToVdApLi6mpaXF9kWf4PkimA7fKC16FXRcLhchISFkZGRYHaVPjRgzjojoGC16pUWvgk9hYSFTpkwhJibG6ih9SkJCSEhJ16JXWvQquFy9epXi4mLbnj/fUWJqOufPVHL27FmroygL6Xn0Kmisya/g5AeHuHLlCpcHJ10759zO2iZsKywsZPHixRanUVbRI3oVVNouBJ6QMtXiJP0jPimZ0PAB1yZwU8FJi14FlYojB4kdnUj0oCFWR+kXjtAwxk6crBOcBblui15EnhWRahHp9Dc60uoJESkTkf0iktlu3RIRKfWs+7EvgyvVU253Cyc/OGi7+ee7k5g6jdLSUurr662OoizizRH9c8CS66y/HUjx3FYC/wUgIg7gSc/6NGCFiKT1JqxSvVF9spwrlxtJTLXntAddSZyUjtvtvja/jwo+3Ra9MWYbcP46m9wJ/NG02g0MEZF4IBsoM8YcM8ZcBdZ5tlXKEhVHDgCtxRdMxkycQmhoqI7TBzFfjNGPAU62e1zpWdbV8k6JyEoRcYmIq6amxgexlPq4itKDDB4+ksHDR1gdpV+FD4ggLS1Nx+mDmC+KvrOrKpvrLO+UMWaVMcZpjHHGxcX5IJZSf2OMoeLIgaAbn2+TlZXF/v37uXLlitVRlAV8UfSVQEK7x2OBU9dZrlS/Ky8vp+HihaAbtmmTlZVFU1MTBw4csDqKsoAvin498GXP2Tc5QJ0x5jRQAKSIyHgRCQfu8WyrVL9rG7YI1iP6zMzWk+F0nD44dfvNWBFZCywAYkWkEngMCAMwxjwFbACWAmVAI/BVz7pmEXkI2Aw4gGeNMYf64GdQqluFhYVEDRzM8PiE7je2oaFDh5KcnKxFH6S6LXpjzIpu1hvg212s20DrPwRKWcrlcpGYmo5IZ786Cg5ZWVm8/vrrtLS04LDhBdFV1/Sbscr2zpw5Q1VVFYmTguv8+Y6cTicNDQ28//77VkdR/UyLXtlesI/Pt2mbsVOHb4KPFr2yPZfLRUxMDCMTJ1gdxVLx8fGMHj1aiz4IadEr2ysoKCAzM5OQEB2XzsrKorCwUC8YHmS06JWtnT17lmPHjjFr1iyro/gFp9PJuXPnKC8vtzqK6kda9MrW2sbng73o1+RXsCa/gjNhowH47Z+3BMWFV1QrLXplawUFBURFRZGWphOnAgyPTyBq4OBrE7yp4KBFr2zN5XKRkZFBWFiY1VH8goiQkKoXDA82WvTKtmprazly5EjQD9t0lJiazoWaM1ys1QuGBwstemVbbacRatF/XNv3Cdqun6vsT4te2VZBQQEDBgxg2rTg/kZsR6MSJxIeEclJHb4JGlr0yrbaxufDw8OtjuJXQhwOxian6Th9ENGiV7azJr+CZ7eWUFLyPgNGpVw7tVD9TWJqOtWVx6mtrbU6iuoHWvTKlk4eOYgxbsZN1mGbziRNmQG0Dm8p+9OiV7Z0ovQAjtAwxkycYnUUvzR6fCph4QPIz8+3OorqB1r0ypZOlO5n9ITWMlOf5AgNIyE1XYs+SGjRK9u5+tFlTpd/wLhJ062O4teSpszg6NGjnD2r59PbnRa9sp2TZYcxbnfQX2ikO23j9Hv27LE4ieprWvTKdk68vx8JCSEhWee3uZ74cSlER0dr0QcBLXplOxWlBxidlEp4RKTVUfxaiMOB0+nUcfogoEWvbKWxsZGqY6Uk6mmVXpk9ezbl5eV8+OGHVkdRfciroheRJSJSKiJlIvLjTtb/UESKPLeDItIiIsM868pF5IBnncvXP4BS7e3btw93SzPjp2RYHSUgZGdnA+hRvc11W/Qi4gCeBG4H0oAVIvKxwU9jzOPGmAxjTAbw98C7xpjz7TZZ6Fnv9F10pT5p165dhDhCSQjyC4F7a/LkyQwaNEiL3ua8OaLPBsqMMceMMVeBdcCd19l+BbDWF+GU6qn8/HzGTpxM+IAIq6MEBIfDwaxZs/QXsjbnTdGPAU62e1zpWfYJIhIFLAH+2m6xAfJEpFBEVnb1IiKyUkRcIuKqqanxIpZSH1dXV8ehQ4dISptpdZSAkp2dTWVlJVVVVVZHUX3Em6KXTpZ1dQn5O4CdHYZtco0xmbQO/XxbROZ1tqMxZpUxxmmMccbFxXkRS6mPKygowBjDeM/54co7s2fPBvR8ejvzpugrgYR2j8cCp7rY9h46DNsYY055/qwGXqZ1KEgpn9u9ezeRkZGMmTjZ6igBJSUlhaFDh+o4vY15U/QFQIqIjBeRcFrLfH3HjURkMDAfeLXdsmgRGdh2H1gE6CTYqk/s3r2brKwsHKF6fVhvrcmvYF1BJaMmprN1+3u8sPuE1ZFUH+i26I0xzcBDwGagBHjRGHNIRB4UkQfbbXoXkGeMaWi3bCSwQ0SKgT3AG8aYTb6Lr1Sr6upqjh49em0YQvVM0pQZXDxfQ23NaaujqD4Q6s1GxpgNwIYOy57q8Pg54LkOy44BOmCq+lzb+PKcOXMorrc4TAAa5/m9RvnhIrgjx9owyuf0m7HKFnbt2sWgQYOYPFnH529EbHwCMUOGcfzwPqujqD6gRa8CnjGG/Px8srOzcTgcVscJSCLChKmZHD+0j5aWFqvjKB/TolcB78SJE1RVVZGTo0MOvTFhahaXGy5RUlJidRTlY16N0Svlj9ou+F3w5msA1A6coBcB74UJ6ZkA7Nixg/R0nULCTvSIXgW8owcKGRoXz7CRnX5hW3kpetAQRo1LZufOnVZHUT6mRa8CWktzE+XvFzNxms6X5wsTpmZSVFREQ0ND9xurgKFFrwLayQ8O0XTlIyZMy7I6ii1MSM+iublZp0OwGS16FdCOHigkxOEgabJ+XcMXElLSiIyMZMeOHVZHUT6kRa8C2tGDhSQkT2VAZJTVUWwhNCycWbNm6Ti9zWjRq4BVf+E8H1Yc1fF5H8vNzeXEiRNUVlZaHUX5iBa9ClhHDxUC6Pi8j+Xm5gLoUb2NaNGrgHXsQGHrKYEJE6yOYisTJkwgPj6e7du3Wx1F+YgWvQpIbrebY4f2MmFqFhKiH2NfEhHmzZvHrl27uHLlitVxlA/o3xAVkA4ePEjjpTodn+8jCxYsoLGxkYKCAqujKB/QolcB6e2330ZCQkierkXfF3JycoiIiOCdd96xOoryAZ3rRgWkd955h8SUdCJjBlkdxXba5gtKmJzBG5vfZPyt9yEi3Ds70eJk6kbpEb0KOFVVVZSWlpKSoVeT6kupGbO5cPZDak7p5QUDnRa9CjhtwwmpM3Va4r6UPCMbgA+K9KLhgU6LXgWcd955h6SkJIaPGmt1FFsbNDSWUeOS+aBYiz7QadGrgNLQ0EB+fj4LFy60OkpQSJmRTeUHJTTWX7Q6iuoFr4peRJaISKmIlInIjztZv0BE6kSkyHN71Nt9leqJnTt30tTUxIIFC6yOEhRSMmZjjJuj+11WR1G90G3Ri4gDeBK4HUgDVohIWiebbjfGZHhu/9zDfZXyytatWxk8eDCZmZlWRwkKo5NSiR40lCNFu6yOonrBmyP6bKDMGHPMGHMVWAfc6eXz92ZfpT6mpaWFd999l3nz5hEaqmcG9wcJCWFS5hzK9hfw0UcfWR1H3SBvin4McLLd40rPso7miEixiGwUkak93Fep61qTX8G/vrCJ2tpawhKm6bVh+9EU51yufnRZJzkLYN4UvXSyzHR4vBcYZ4yZAfwn8EoP9m3dUGSliLhExFVTU+NFLBVsSvZsIzQsnJTp2VZHCSrjJs8gMnogeXl5VkdRN8iboq8EEto9Hgucar+BMeaiMabec38DECYisd7s2+45VhljnMYYZ1xcXA9+BBUM3O4WSgp3kDx9FuERkVbHCSqO0FBSM+ewdetWrl69anUcdQO8KfoCIEVExotIOHAPsL79BiIySkTEcz/b87znvNlXKW+c/OAw9RfOM2XWPKujBKUpzpu4dOkSu3bpL2UDUbdFb4xpBh4CNgMlwIvGmEMi8qCIPOjZ7HPAQREpBp4A7jGtOt23L34QZW9twzapOu2BJcanZTBwoA7fBCoxptMhc0s5nU7jcul5u6qV2+1mdu5NjJk4hS9859Hud1B9ouivv2Hbtm1s376dsLAwq+OoDkSk0BjT6XSu+s1Y5ff27t1L/YXzpM26yeooQW3x4sXU1dWxZ88eq6OoHtKiV35v8+bNrWfbzNBhGyvl5uYSFRXF5s2brY6iekiLXvk1t9tNXl4eE6fPYkBklNVxglpERAQ333wzeXl5evZNgNGiV37N5XJRXV2twzZ+YE1+BQNTc6irq+Nfnvkra/Ir9ItrAUKLXvm1V199lejoaFJnzrE6igImTM0ketBQ9r/3ltVRVA9o0Su/1djYyKZNm1i8eDHhAyKsjqOAEIeD9DkL+aAon8s6dXHA0KJXfuvNN9+ksbGR5cuXWx1FtTPtU7fgbmnm0J53rY6ivKRFr/zWyy+/zNixY8nKyrI6impnVOJE4sYmcUCHbwKGFr3yS6dPnyY/P5/ly5cTEqIfU38iIkz/1C1UlpVw/sMqq+MoL+jfIOWX1q9fjzGGO+/Uyxf4o/SchSDCgffetjqK8oIWvfI7xhheeeUVnE4nY8fqBcD90aBhcYxPy6B4xxZaWlqsjqO6oUWv/ErbBUbKy8uJnz5Pz9X2Y5kLPk3duQ/Ztm2b1VFUN7Told9xvfUaAyKjSMvWKYn92aSZc4gZMox169ZZHUV1Q4te+ZX6C+c5XLCdGXMX6QVG/JwjNJSZ829n+/btnDx5svsdlGW06JVf2fvuRtwtzThvucPqKMoLmfNvJyQkhP/+7/+2Ooq6Di165TeamprYu/UNJqZnMXyU/hI2EAwaFsfNN9/MSy+9xJUrV6yOo7qgRa/8xttvv82lC+dw3rrM6iiqB+655x5qa2t1+mI/pkWv/MYLL7zAkNiRJE+fZXUU1QM5OTkkJSXxwgsv4I9XrFNa9MpPlJaWUlBQQNbNdxAS4rA6juqBdQWVTJ77afbv388v/viGng7rh7TolV94+umniYqKYua8xVZHUTcgY95iogcNYefreqqlP9KiV5Y7fvw4GzduZMWKFUTGDLI6jroBYeEDyF50F0cPFnK6/AOr46gOvCp6EVkiIqUiUiYiP+5k/RdFZL/n9p6IzGi3rlxEDohIkYi4fBle2cMzzzxDWFgY999/v9VRVC84b76DAZFR7HxDT7X0N90WvYg4gCeB24E0YIWIpHXY7Dgw3xgzHfgZsKrD+oXGmAxjjNMHmZWNVFVVsX79er7whS8QGxtrdRzVCxFR0WTdfAclrh0cP37c6jiqHW+O6LOBMmPMMWPMVWAd8LEpBY0x7xljaj0PdwN6ErTyytNPP42I8LWvfc3qKMoHZi+6i9DQMH7/+99bHUW1403RjwHaf7+50rOsK18HNrZ7bIA8ESkUkZVd7SQiK0XEJSKumpoaL2KpQLYmv4KnNu3lL399iWm5t/H2iat6toYNxAweysx5S3j11Vf1qN6PeFP00smyTk+WFZGFtBb9I+0W5xpjMmkd+vm2iHQ6U5UxZpUxxmmMccbFxXkRSwW6bS//CWMMn/r03VZHUT40d9m9hIeH86tf/crqKMrDm6KvBBLaPR4LnOq4kYhMB54B7jTGnGtbbow55fmzGniZ1qEgFeTOVBxl3/bNzLp1GUPjRlkdR/lQzOChfP3rX2fLli3s3bvX6jgK74q+AEgRkfEiEg7cA6xvv4GIJAIvAfcZY460Wx4tIgPb7gOLgIO+Cq8CkzGGLWtXERkVw03L7rU6juoD999/P3FxcTz++OP6bVk/0G3RG2OagYeAzUAJ8KIx5pCIPCgiD3o2exQYDvy2w2mUI4EdIlIM7AHeMMZs8vlPoQLKO++8Q3lJEfOW30dk9ECr46g+EBUVxXe/+12KiorIy8uzOk7QE3/819bpdBqXS0+5t6OmpiaWLVvGpSstPPCzp3CEhlodSfURt7uFVY9+i+arV3jgX35HWPgA7p2daHUs2xKRwq5OYddvxqp+tXr1asrLy7nt7m9oydtcSIiDxfd+k9rq02x79QWr4wQ1LXrVb06ePMkTTzzBggULSJ6hv5MPBuPTMsiYt4RdG/+sUyNYSIte9QtjDD/96U8JCQnh0UcfRaSzs3aVHd129zeIHjSE1579FU1NTVbHCUpa9KpfrF27lvfee4/vf//7xMfHWx1H9aOI6Bhuv+87fFhxlGeffdbqOEFJi171uWPHjvH4448zd+5cVqxYYXUcZYHJWZ9iyqybePLJJzl4UM+w7m9a9KpPPb/tCF994CEkNBzn//oWa/ec1KkOgtTSr3yX2NhYvve971FbW9v9DspntOhVnzHGsHH1k1RXlbN85Y8YOHS41ZGUhaJiBvHrX/+ampoaHnnkEVpaWqyOFDS06FWfWb16NcXb85j7mRV6HVgFwLRp0/jJT37C9u3beeqpp6yOEzT0RGbVJ959911+8YtfMCnzUyy46z6r4yg/sSa/AjMuh+m5t/Kb3/yGyqvRTPvUzfpFqj6mRa98bt++fTz88MNMmjSJZSt/hITofxzV34gIn/7K97h4voZXn3mcAZFRoEXfp/RvoPKp4uJiHnzwQUaOHMnTTz9NeESk1ZGUHwoND+cL3/0p8eOS+cuT/0J+fr7VkWxNi175zM//8Cr3feV+HBHRLPvOz9hc1mB1JOXHBkRGseL7P2fYyNF861vfYufOnVZHsi0teuUTW7ZsYe2//yND4+L5yj/8f4bEjrQ6kgoAUTGD+OIP/x8JCQk8+OCDvPTSS1ZHsiUtetUrbreb3/3udzz88MOMGpfMl3/8rwwcoqdRKu8NHDKc1atXk52dzU9+8hOeeOIJ3G631bFsRacpVjestraWRx55hO3bt7N06VKmLVupY/LqhrU0N7Ph+Sco2r6ZCVMzWfZ3P+CBJTOtjhUwdJpi5XNbt27ls5/9LLt37+axxx7j3/7t37TkVa84QkP5zNf+N0u//B0qPjjE7/7PA3rREh/RI3rVI5WVlXzrh//IB0W7iR2dyPKVPyI+KcXqWMpmzp4+ySur/pXTx48wf/58fvjDHzJx4kSrY/m16x3Ra9Err5w+fZo//OEPvPjii7gR5i//Etm33aUXD1F9pqW5mfy8l9nx2hquXvmIzAVLyf303XxzaaddFvS06NUNMcZw+PBh1q5dy/r16zHGsGzZMhLmfpZBw+KsjqeCRMPFC2x7dTWFW98A4PYlS/jKV77C9OnTLU7mX7ToVY9UVlayZcsWnlvzZ6orjxMaFk7GvMXMuf3zetqkssyFmjMUvPkqB3fmUV9fT3JyMkuXLuX2228nKSnJ6niW06JX13X27FmKi4spKChg+/btHDt2DIAxE6cwI/dW0mbPJzJ6oMUplWq1bOow1q9fz4YNGygsLARg/PjxzJkzh5ycHGbOnElsbKzFKftfr4teRJYAvwYcwDPGmF90WC+e9UuBRuB+Y8xeb/btjBa97xljuHDhAuXl5Rw/fvzarbS0lMrKSgAcoWEkTZ5O8oxskqfPYtjIMRanVur66s5VU+LawfFDezlVdojGxkYARo0axdSpU0lOTiYxMfHaLS4uzraXsexV0YuIAzgC3AZUAgXACmPM4XbbLAW+Q2vRzwZ+bYyZ7c2+nfHnom97v7z983r7ebNtc3MzTU1NNDU1cfXq1U7vNzY2Ul9fT319PZcuXaKhoYFLly5RW1tLTU3NtVv763U6QsMYOiKeuNHjGDNxMmOTpzBqXDJh4QN892Yp1Y9ampuoOlbKqeNHOH38CKdPlFFbfQp3u3nvQ8PCiR40lOjBQ5g0bjTDhw9n8ODBREdHEx0dTUxMzLX7ERERhIWFERYWRmho6LX77R87HA5CQkIQkU/c2i/vD9crem9OmcgGyowxxzxPtg64E2hf1ncCfzStDbVbRIaISDyQ5MW+PpObm8vly5eB3hXw9bYPBKHhA4iIjCYiOoaBQ4YxPGkK4zJyGTh4GENHjiY2PoEhsaMIcTisjqqUzzhCw0hMTScxNf3aMndLC3XnqjlffYrzZ6q4cPYMDXUXaLhYy+GyEzTsLeZywyVamvv2ouWd/UPQ8QYwfPhwtmzZ4vPX96boxwAn2z2upPWovbttxni5LwAishJY6XlYLyKlXmQLFLHAWatDWEzfg1b6Puh7ANd5D3rxP4BxXa3wpug7e9WOh7ldbePNvq0LjVkFrPIiT8AREVdX/6UKFvoetNL3Qd8D6P/3wJuirwQS2j0eC5zycptwL/ZVSinVh7yZ66YASBGR8SISDtwDrO+wzXrgy9IqB6gzxpz2cl+llFJ9qNsjemNMs4g8BGym9RTJZ40xh0TkQc/6p4ANtJ5xU0br6ZVfvd6+ffKT+DdbDkn1kL4HrfR90PcA+vk98MsvTCmllPIdnaZYKaVsToteKaVsToveh0RkiYiUikiZiPy4k/UiIk941u8XkUwrcvYlL96DBSJSJyJFntujVuTsSyLyrIhUi8jBLtYHw+egu/fA1p8DEUkQka0iUiIih0Tke51s03+fg7av4uutdzdaf9l8FJhA62mlxUBah22WAhtp/X5BDpBvdW4L3oMFwOtWZ+3j92EekAkc7GK9rT8HXr4Htv4cAPFApuf+QFqngrGsD/SI3neuTRVhjLkKtE330N61qSKMMbuBtqki7MKb98D2jDHbgPPX2cTunwNv3gNbM8acNp6JHY0xl4ASWmcKaK/fPgda9L7T1TQQPd0mkHn7880RkWIR2SgiU/snml+x++fAW0HxORCRJGAmkN9hVb99DvQ6cL7Tm6ki7MKbn28vMM4YU++Z9fQVINguOmv3z4E3guJzICIxwF+Bh40xFzuu7mSXPvkc6BG97/Rmqgi76PbnM8ZcNMbUe+5vAMJEJNiuEmH3z0G3guFzICJhtJb8C8aYlzrZpN8+B1r0vtObqSLsotv3QERGiWd6PhHJpvUzeK7fk1rL7p+Dbtn9c+D52X4PlBhj/r2Lzfrtc6BDNz5iejFVhF14+R58DvimiDQDl4F7jOcUBLsQkbW0nlUSKyKVwGNAGATH5wC8eg/s/jnIBe4DDohIkWfZPwCJ0P+fA50CQSmlbE6HbpRSyua06JVSyua06JVSyua06JVSyua06JVSyua06JVSyua06JVSyub+BwEXbQV2ttLlAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.distplot(x=df.iloc[:,0], fit=norm, kde=False)\n",
    "print(f\"\"\"covariance:\n",
    "{df.cov()}\n",
    "\n",
    "\n",
    "mean: {[df.iloc[:,0].mean(), df.iloc[:,1].mean()]}\"\"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 455,
   "id": "c4f63b08",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\fgrimald\\Anaconda3\\envs\\sandy-devel\\lib\\site-packages\\seaborn\\distributions.py:2619: FutureWarning: `distplot` is a deprecated function and will be removed in a future version. Please adapt your code to use either `displot` (a figure-level function with similar flexibility) or `histplot` (an axes-level function for histograms).\n",
      "  warnings.warn(msg, FutureWarning)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<AxesSubplot:>"
      ]
     },
     "execution_count": 455,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD4CAYAAADlwTGnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAAsk0lEQVR4nO3deXiU9d3v8fc3CWEJayAECAmBEEJYIwkBwQooKFBbtJtgH0tbe6yKrdbHnsce2/O0x56rXq0+dvPo40LVVsStClUsWgRZJCEJS8IWCIGEkJBE1hAC2b7nj0x0jIHMkEnuzMz3dV25ZuZePzPXZL4zv/t3/25RVYwxxgSfEKcDGGOMcYYVAGOMCVJWAIwxJkhZATDGmCBlBcAYY4JUmNMBvDFo0CCNj493OoYxxviVnJycT1Q1quV0vyoA8fHxZGdnOx3DGGP8iogUtTbdmoCMMSZIWQEwxpggZQXAGGOClBUAY4wJUlYAjDEmSHlUAERkvojki0iBiDzUyvyxIrJVRC6KyINu05NEZKfb31kRud8175cicsxt3kKfPStjjDFtarMbqIiEAk8C84ASIEtEVqvqXrfFTgI/Bm52X1dV84EUt+0cA95yW+QJVX2sHfmNMcZcIU9+AaQDBapaqKq1wEpgkfsCqlqhqllA3WW2cz1wSFVb7Y9qTLBSVerr652OYYKQJwUgBjjq9rjENc1bi4FXWky7V0RyRWS5iAxobSURuVNEskUku7Ky8gp2a0zXc+HCBV566SVuuukmJk6cSEpKCkuXLuXdd9/FrtFhOosnBUBamebVO1REwoGvAq+7TX4KSKCpiagMeLy1dVX1GVVNU9W0qKgvnMlsjN8pKipi8eLF/OY3v6FPnz7ccccdLF26lE8++YQHH3yQu+++m7Nnzzod0wQBT4aCKAFi3R4PB0q93M8CYLuqljdPcL8vIs8C73i5TWP8zuHDh7nttttQVZ5++mlmzZr16bwHHniAv/3tbzz++OPccccdPP/88/Tt29fBtCbQeVIAsoBEERlJ00HcxcBtXu5nCS2af0RkqKqWuR7eAuz2cpvGdHkrMos/vX/u9En+8uv7qW1QvvvwExzrEfO5+QDdxs7ha/f05PU/P8I3b7+D2//jt/zbjJGdHdsEiTabgFS1HrgXWAvsA15T1T0icpeI3AUgIkNEpAR4APi5iJSISF/XvF409SD6e4tN/1ZE8kQkF5gD/MRnz8qYLkZVWfXcY5w7e5olP3mEgUMufRhtzFXT+codD1B8YDcfvf3XTkxpgo1Ho4Gq6hpgTYtpT7vdP05T01Br654HBrYy/Xavkhrjx3LWv0Ph7hwWfOdeho1KanP5STOup2jfLja/s5KcW79MampqJ6Q0wcbOBDamg1WdOsG/Xn2OURNSSZ1zk8fr3fhv99AvMopHHnnEuomaDmEFwJgOtuHvL9LY0MDC7/wIkdY61bUuvHsP5i35Ifn5+axcubIDE5pgZQXAmA5UfrSQnZvfZ+rcrzJg8FCv1x+bOpOrr76aJ598kurq6g5IaIKZFQBjOtBHb/2VHr0iuOYrS65ofRHh/vvv5/Tp06xYscLH6UywswJgTAc5ePAg+ds/Jn3uzfSM6HPF25k0aRLXXnsty5cvt18BxqesABjTQZ577jm6de/B1HmL2l74MlZkFjN69tc5ffo0P3v8WVZkFn/h/AFjroQVAGM6wPHjx3n33XdJnb2QXr3bfzZvzKixxI4ZT9a/VtHY2OCDhMZYATCmQ7z66qs0NjYydW77vv27S593C6crj3NwZ6bPtmmCmxUAY3ystraWN954g1mzZtE/aojPtjt2ygz6DRzMtg/e9tk2TXCzAmCMj33wwQd88skn3Habt0NmXV5IaChTZi/kyL5dnCw/5tNtm+BkBcAYH3vjjTcYPnw4M2fO9Pm2J10zD5EQdm563+fbNsHHCoAxPlRaWkpmZiY333wzISG+//fqO2AQoyelkbv5AxsewrSbFQBjfGj16tWoKosW+e7gb0sp186n6vQJtmzZ0mH7MMHBCoAxPqKqvP3220ydOpXhw1sdHNcnEidPo2fvvrzzjl1DybSPFQBjfGTPnj0UFRV16Ld/gNCwMJLTrmHdunWcP3++Q/dlApsVAGN8ZO3atYSFhTF37twO39eE6XOoqalh/fr1Hb4vE7isABjjA6rK+++/z/Tp0+nXr1+H7y9uzASio6N59913O3xfJnB5dEUwY8zlPfHGBoqLi5l43S2dMk6PhISwYMECXn75ZaqqqujT58oHmzPBy34BGOMD+7I2IxJC0pQZnbbPuXPnUldXx6ZNmzptnyawWAEwxgf252xmxNhJRPTt32n7TElJITIyknXr1nXaPk1gsQJgTDsVFBTwSWkxyWnXdOp+Q0NDmT17Nhs3bqS2trZT920Cg0cFQETmi0i+iBSIyEOtzB8rIltF5KKIPNhi3hERyRORnSKS7TY9UkQ+EJGDrtsB7X86xnS+tWvXgghJqZ3X/NNs7ty5nDt3jqysrE7ft/F/bRYAEQkFngQWAOOAJSIyrsViJ4EfA49dYjNzVDVFVdPcpj0ErFPVRGCd67ExfueDDz4gNnEcffoP7PR9X3311fTs2dOagcwV8eQXQDpQoKqFqloLrAQ+d6aLqlaoahZQ58W+FwEvuu6/CNzsxbrGdAlFRUXk5+eTnNq5zT/NevTowcyZM/nwww9RVUcyGP/lSQGIAY66PS5xTfOUAu+LSI6I3Ok2PVpVywBct4NbW1lE7hSRbBHJrqys9GK3xnS8jRs3AjDmqqs7fd/Nl4bsNTKF8vJyHn91nV0q0njFkwIgrUzz5qvGTFWdQlMT0jIRudaLdVHVZ1Q1TVXToqKivFnVmA63adMm4uPjGTB4qGMZEidPQ0JCyN/+sWMZjH/ypACUALFuj4cDpZ7uQFVLXbcVwFs0NSkBlIvIUADXbYWn2zSmK6ipqSEzM5Nrr/XqO43P9erdl7gxE6wAGK95UgCygEQRGSki4cBiYLUnGxeRCBHp03wfuAHY7Zq9Gljqur8UWOVNcGOctm3bNmprax0vANDUBFV5rIgzJ+x7lPFcmwVAVeuBe4G1wD7gNVXdIyJ3ichdACIyRERKgAeAn4tIiYj0BaKBzSKyC9gGvKuq/3Rt+lFgnogcBOa5HhvjNzZt2kTPnj1JS0tre+EOljAhFYBDedltLGnMZzwaC0hV1wBrWkx72u3+cZqahlo6C0y+xDZPANd7nNSYLmbTpk2kp6fTvXt3p6MwaFgcfSMHUbh7u9NRjB+xweCM8dKKzGJOHD9GcXEx42d9tUv0vBERRo1PZX/OFurr6wkLs39t0zYbCsKYK1CQuw2AhElTHU7ymVETU7lw/hx5eXlORzF+wgqAMVfgUG4WA4cMZ0DUEKejfGrUuKsQCbFrBRuPWQEwxkt1Fy9wZH8uo7vQt3+Anr37MmzkGDZv3ux0FOMnrAAY46Wi/Dwa6utImOh875+WRk2YQl5eHmfOnHE6ivEDVgCM8dKRvTsJDetG3JgJTkf5goSJaTQ2NpKRkeF0FOMHrAAY46XD+3YyfHQy3br3cDrKF8SMGkufPn2sGch4xAqAMV44deoUx4sPEZ+c4nSUVoWEhjJ9+nS2bNlio4OaNlkBMMYL27ZtA1VGjktxOsolzZgxg7KyMo4cOeJ0FNPFWQEwxgtbt24lvEdPho1McjrKJU2bNg1wFStjLsMKgDFeyMjIYETSREK78Jm28fHxREdHk5mZ6XQU08VZATDGQ2VlZRQVFRE/7iqno1yWiJCens62bdvsOIC5LCsAxniouWtlV27/h6axiohK5MSJE/zhrU1dYqwi0zVZATDGQ1u3biUyMpLBMfFOR2lTfHLTILxF+3Y5nMR0ZVYAjPGAqpKRkcG0aU2XX+zq+kcNof+gaA5bATCX0fXfycZ0AYWFhVRWVjJ9+nSno3gsPjmFov270MZGp6OYLsoKgDEeaO5S2dzF0h+MSJ7MhepzHD9a6HQU00VZATDGAzk5OURFRREXF+d0FI/ZcQDTFisAxrRBVcnOziY1NRURcTqOx/oOGETkkOEcsQJgLsEKgDFtOHbsGOXl5V3i4u/eGpk8maL8POrr652OYrogKwDGtCE7OxvALwvAiLGTqL1wnv379zsdxXRBHhUAEZkvIvkiUiAiD7Uyf6yIbBWRiyLyoNv0WBFZLyL7RGSPiNznNu+XInJMRHa6/hb65ikZ4xsrMotZkVnMq+99RI9evck+2d3vTqqKS5oIfFbEjHHXZgEQkVDgSWABMA5YIiLjWix2Evgx8FiL6fXAv6tqMjAdWNZi3SdUNcX1t+ZKn4QxHan4wG5iE8f7Rf//lvr0H0hk9DArAKZVnryj04ECVS1U1VpgJbDIfQFVrVDVLKCuxfQyVd3uul8F7ANifJLcmE5w7vRJTh4v+fSbtD+KGzORnJwcGu18ANOCJwUgBjjq9riEK/gQF5F44CrAfYjCe0UkV0SWi8iAS6x3p4hki0h2ZWWlt7s1pl2KD+4B6JKXf/RUXNJETp8+zaFDh5yOYroYTwpAa/3evBpiUER6A28C96vqWdfkp4AEIAUoAx5vbV1VfUZV01Q1LSoqypvdGtNuxfl5dAvvztD40U5HuWLNxcuagUxLnhSAEiDW7fFwoNTTHYhIN5o+/F9W1b83T1fVclVtUNVG4FmampqM6VKKD+wmJiGZ0LBuTke5Yv2jhhAdHW0FwHyBJwUgC0gUkZEiEg4sBlZ7snFpOmvmeWCfqv5Xi3lD3R7eAuz2LLIxneNC9TnKjxb6dfs/NF0fIC0tjezsbLs+gPmcNguAqtYD9wJraTqI+5qq7hGRu0TkLgARGSIiJcADwM9FpERE+gIzgduB61rp7vlbEckTkVxgDvAT3z89Y67c0YK9oOrX7f/N0tLSqKiooKSkxOkopgvx6Lp2ri6aa1pMe9rt/nGamoZa2kzrxxBQ1ds9j2lM5yvOzyMkNJThCWOdjtJuzSexZWdnExsb28bSJlj4X8dmYzpJ8YHdDI0fQ7fuPZyO0m4JCQn079/fjgOYz7ECYEwrLly4QOnhA8Ql+X/zDzQdB0hNTbUCYD7HCoAxrcjNzaWxoT4g2v+bpaWlUVxcTEVFhdNRTBdhBcCYVmRlZYEIsYnjnY7iM1OnTgXsfADzGSsAxrQiJyeHwcPj6RnRx+koPpOUlESvXr2sAJhPWQEwpoW6ujp27tzJiDH+3f+/pbCwMKZMmWIFwHzKCoAxLezbt4+amhpiA+QAsLu0tDQOHjzIqVOnnI5iugArAMa00PwNOZAOADdrPh9g+/btDicxXYEVAGNayMnJIS4ujj79BzodxecmTpxIeHi4NQMZwMMzgY0JFo2NjeTk5HD99dc7HcWn3K9kNmRkEu9/9DGxs4u5bVqcg6mM0+wXgDFuCgoKOHPmjF9e/9dTcWMmUFZUQO2FGqejGIdZATDGTU5ODuCfF4D3VNyYCWhjIyUFe52OYhxmBcAYN9nZ2QwePJjhw1sb2zAwDB+djEgIxQdsBPZgZwXAGBdVJScnh7S0NJouZRGYuveMIDpulBUAYwXAmGYlJSWUl5eTmprqdJQONyJpIscO7ae2ttbpKMZBVgCMcQmG9v9msWMmUF9Xy549e5yOYhxkBcAYl6ysLPr168fo0f57AXhPNZ/k1lz0THCy8wBM0GvuI79hSybRo5JZmRX4l02M6NufgUOGk52dzQ9+8AOn4xiH2C8AY4Bzp09ysvxYQA7/cClxYyawY8cOGhsbnY5iHGIFwBj4tEdMXFJgjQB6ObFJEzh79iwHDx50OopxiEcFQETmi0i+iBSIyEOtzB8rIltF5KKIPOjJuiISKSIfiMhB1+2A9j8dY65McX4e3cK7MyQu8Nv/mzUPd23jAgWvNguAiIQCTwILgHHAEhEZ12Kxk8CPgce8WPchYJ2qJgLrXI+NcUTxgd0MTxxHaFjwHBbrNyiaIUOG2IHgIObJL4B0oEBVC1W1FlgJLHJfQFUrVDULqPNi3UXAi677LwI3X9lTMKZ9LlSfo7zkcFC1/8NnF4rPyclBVZ2OYxzgSQGIAY66PS5xTfPE5daNVtUyANft4NY2ICJ3iki2iGRXVlZ6uFtjPHf04B5QJS7ArgDmidTUVCoqKjh69GjbC5uA40kBaO2ceE+/LrRn3aaFVZ9R1TRVTYuKivJmVWM8UnQgj5DQMGISxjodpdM1n/VszUDByZMCUALEuj0eDpR6uP3LrVsuIkMBXLcVHm7TGJ86emA3w0aOoVt4d6ejdLrRo0fTr18/OxAcpDwpAFlAooiMFJFwYDGw2sPtX27d1cBS1/2lwCrPYxvjGzU1NZQePhB07f/NQkJCmDJliv0CCFJtFgBVrQfuBdYC+4DXVHWPiNwlIncBiMgQESkBHgB+LiIlItL3Uuu6Nv0oME9EDgLzXI+N6VS5ubk0NjQEVf//ltLS0igqKsKOsQUfj/q8qeoaYE2LaU+73T9OU/OOR+u6pp8AAuu6e8bvZGdngwixo1v2bA4ezccBtm/fzo033uhwGtOZ7ExgE9Sys7OJjh1Fj4jeTkdxzLhx4+jZs6cdBwhCVgBM0Kqrq2PXrl1B2/7frFu3bkyePNmOAwQhKwAmaO3Zs4eamhpGBHH7f7PU1FT2799PVVWV01FMJ7ICYIJWVlYWEFwDwF1KamoqqsqOHTucjmI6kRUAE7SysrJISEggom9/p6M4bvLkyYSFhVkzUJAJnpGvjHFTX1/P9u3buemmm5yO4qjmi+EARMeN5r0NW4ie8XVumxbnYCrTWewXgAlK+/fvp7q6mqlTpzodpcuIGzOB0sID1NuF4oOGFQATlJq7PFoB+Exs0gQa6us4djjf6Simk1gBMEEpKyuLuLg4Bg9udRDaoBSbOB5oujiOCQ5WAEzQaWxsJCcnx779t9Crd1+iYkZw1HV5TBP4rACYoHPgwAHOnDljBaAVcWMmcLRgHw0NDU5HMZ3ACoAJOtb+f2lxSROpvXCe/fv3Ox3FdALrBmqCRnOXxzfXfkS/gdFsOFoPR4vbWCu4NA+LkZWVxfjx4x1OYzqa/QIwQUVVKc7PY8RYO/u3NX0jo4iMHkZmZqbTUUwnsAJggsonpcWcrzrDiKRJTkfpskaMnUx2djb19fVORzEdzAqACSpF+bkAxNkvgEuKT57MuXPn2Ldvn9NRTAezAmCCSnH+bvoMGMSAqKFOR+my4sdOBrBmoCBgBcAEDVWlKD+XuKQJiIjTcbqs3v0jSUhIYNu2bU5HMR3MCoAJGifLSzl3+qS1/3sgPT2dnJwc6urqnI5iOpAVABM0il3t/3YBmLZNmzaN8+fPs3u3nRUcyKwAmKBxZH8uEX37M3BorNNRurz09HQAawYKcB4VABGZLyL5IlIgIg+1Ml9E5I+u+bkiMsU1PUlEdrr9nRWR+13zfikix9zmLfTpMzPGjapyZN9O4pNTrP3fAwMGDGDMmDFkZGQ4HcV0oDYLgIiEAk8CC4BxwBIRGddisQVAouvvTuApAFXNV9UUVU0BUoHzwFtu6z3RPF9V17T3yRhzKYWFhZw7fZKR41KcjuI3pk2bxo4dO6i16wMELE9+AaQDBapaqKq1wEpgUYtlFgEvaZMMoL+ItOxndz1wSFWL2p3aGC81f5ONT05xNogfSU9P5+LFi+zatcvpKKaDeFIAYoCjbo9LXNO8XWYx8EqLafe6moyWi8iA1nYuIneKSLaIZFdWVnoQ15gvysjIoP+gaAYMtv7/npo6dSoiYscBApgnBaC1BlP1ZhkRCQe+CrzuNv8pIAFIAcqAx1vbuao+o6ppqpoWFRXlQVxjPq+hoYGsrCz79u+lfv36MXbsWDshLIB5UgBKAPduE8OBUi+XWQBsV9Xy5gmqWq6qDaraCDxLU1OTMT63f/9+zpw5Q7y1/3tt+vTp7Ny5kwsXLjgdxXQATwpAFpAoIiNd3+QXA6tbLLMa+I6rN9B04IyqlrnNX0KL5p8WxwhuAazDsekQn7X/T3Y4if9YkVnMisxiavqNpK6ujsdfef/T4bRN4GjzegCqWi8i9wJrgVBguaruEZG7XPOfBtYAC4ECmnr6fK95fRHpBcwDfthi078VkRSamoqOtDLfGJ/IyMggISGBPv0HOh3F78QlTUBCQjiyb5f1oApAHl0QxtVFc02LaU+73Vdg2SXWPQ984T9PVW/3KqkxV6C2tpacnBy+9rWvOR3FL3XvGcHQ+ESO7NsJLHU6jvExOxPYBLTc3FxqamqYNm2a01H8VnxyCqWH87lYU+10FONjVgBMQMvMzEREPh3awHgvYcIUGhsaOLI/1+koxsesAJiAtmXLFsaPH0+/fv2cjuK3ho8eR7fw7hTuznE6ivExuyi8CUgrMou5UH2Onbt2MfPLt1oPlnYI6xbOiLGTrQAEIPsFYALW4b070MZGEiakOh3F7yVMTOVkeSlHjx5te2HjN6wAmIB1aHcO3Xv2IiYh2ekofm/U+CkAfPzxxw4nMb5kBcAEJFWlcHcO8eOuIjTMWjrba+DQWPpGRrFlyxanoxgfsgJgAtKJsqOcOVFhzT8+IiKMmpBKRkYG9fX1TscxPmIFwASkQ64DllYAfCdhQipVVVXk5eU5HcX4iBUAE5AKd+cwcMhw+kcNcTpKwBg5rulqatYMFDisAJiAc/HiRY7szyVhYprTUQJKz959mThxIps3b3Y6ivERKwAm4GRnZ1Nfe5FR1vzjczNnziQvL49Tp045HcX4gBUAE3C2bNlCaFg3Royd5HSUgDNr1iwaGxvtV0CAsAJgAs7GjRuJS5pAePceTkcJOBMnTmTgwIFs2LDB6SjGB6wAmIBSVFTEoUOHGJMy3ekoASkkJIRZs2axefNm6urqnI5j2skKgAko69evByDRCkCHWJFZTMiw8Zw9e5bHVqy1MZb8nBUAE1A2bNhAYmIiA6z7Z4cZNX4KIaFhHNxpF4v3d1YATMA4e/YsOTk5zJkzx+koAa17z17Ej53EgV1WAPydFQATMDZt2kR9fT2zZ892OkrAS0yZxomyo5wsP+Z0FNMOVgBMwFi/fj2RkZFMmmTdPzva6MlNV1g7uGubw0lMe3hUAERkvojki0iBiDzUynwRkT+65ueKyBS3eUdEJE9EdopIttv0SBH5QEQOum4H+OYpmWBUV1fHpk2bmDVrFqGhoU7HCXiRg4cxaFgcB60ZyK+1WQBEJBR4ElgAjAOWiMi4FostABJdf3cCT7WYP0dVU1TV/dz8h4B1qpoIrHM9NuaK7Nixg7Nnz1r7fydKTJlG0f48qqqqnI5irpAnvwDSgQJVLVTVWmAlsKjFMouAl7RJBtBfRIa2sd1FwIuu+y8CN3se25jPW79+Pd26dWPGjBlORwkaY6fMoLGh3k4K82OeFIAYwP06cCWuaZ4uo8D7IpIjIne6LROtqmUArtvB3gQ3Bpr6pb+cUcSqd9cyIjmFVbtPWN/0ThIzaix9Bgxi7dq1TkcxV8iTAiCtTFMvlpmpqlNoaiZaJiLXepEPEblTRLJFJLuystKbVU2QOF5UwKnKMpLTvuR0lKAiISEkp13D5s2bqa6udjqOuQKeFIASINbt8XCg1NNlVLX5tgJ4i6YmJYDy5mYi121FaztX1WdUNU1V06KiojyIa4LN3m0bCQkNJWnK1U5HCTpj067h4sWLbNy40eko5gp4UgCygEQRGSki4cBiYHWLZVYD33H1BpoOnFHVMhGJEJE+ACISAdwA7HZbZ6nr/lJgVTufiwlCqsrerI2MTE6hZ+++TscJOrGJ4xg0aBDvv/++01HMFWizAKhqPXAvsBbYB7ymqntE5C4Rucu12BqgECgAngXucU2PBjaLyC5gG/Cuqv7TNe9RYJ6IHATmuR4b45XjRQWcrjxO8lSvWhaNj4SEhDJ37lw++ugjampqnI5jvBTmyUKquoamD3n3aU+73VdgWSvrFQKTL7HNE8D13oQ1pqW92zYiISEkpVrvH6fccMMNrFy5ks2bNzNv3jyn4xgv2JnAxm+pKvuyNjFy3FX0suYfx0ydOpX+/ftbM5AfsgJg/NbevXs5VVnGOGv+cVRYWBhz585l/fr1XLx40ek4xgtWAIzfWrt2rTX/dBHz58+nurraTgrzMx4dAzCmq2lsbOS9996z5p8uYEVmMY0hw+jdP5L/9+JKTvRPBuC2aXEOJzNtsV8Axi/l5ORQUlLCxBnWj6ArCAkJZeL06yjIzeJ81Rmn4xgPWQEwfuntt9+mV69ejE2d6XQU4zJx5vU0NjSwJ/Mjp6MYD1kBMH7n/PnzrF27lvnz5xPevYfTcYxLdOwoBseOJPfjfzkdxXjICoDxO//617+orq5m0aKWg9Iap02aMZfSwnw+KTva9sLGcVYAjN95++23GT58OGlpaW0vbDrVhKvnIBJC3scfOh3FeMAKgPErZWVlZGRksGjRIkJC7O3b1fTpP5CR469i99Z1NDY2Oh3HtMH+g4xf+cc//oGqWvNPFzZpxlxOf1JOZqZdLrKrswJg/MKKzGJe3nqEF195nbgxE9hUqnbhly4qOe0aekb0YeXKlU5HMW2wAmD8xuG9Ozh5vISrZi90Ooq5jLDwcFKunc+6deuoqGj1Mh+mi7ACYPxG1rrV9OrTj3FT7cpfXd2UOQtpaGjg9ddfdzqKuQwrAMYvnK48zoGdmUyZvZCwbuFOxzFtiBw8jC996Uu89tpr1NXVOR3HXIIVAOMXste/g4gwxZp//MaSJUuoqKhg/fr1Tkcxl2AFwHR5Fy5cYOfGtSRNmUG/gYOdjmM8dO211zJ06FBeeeUVp6OYS7ACYLq89957j5pzZ5l6/VecjmK8EBoayq233kpGRgaHDh1yOo5phRUA06WpKn/7298YNCyOEWNbvbqo6aJWZBYTMupqwsK78/Cjf2BFZrF13e1irACYLu3jjz9m7969TLvhFkTE6TjGSxF9+3PVtfPJ2/ohZ05Yl9CuxgqA6dKeeuophgwZwqSZc52OYq7Q1Qu+AUDGP990OIlpyaMCICLzRSRfRApE5KFW5ouI/NE1P1dEprimx4rIehHZJyJ7ROQ+t3V+KSLHRGSn68+6d5jPycrKIicnhzvuuMO6fvqxfgMHM2H6dWz/6D2qz552Oo5x02YBEJFQ4ElgATAOWCIi41ostgBIdP3dCTzlml4P/LuqJgPTgWUt1n1CVVNcf2va91RMoHnqqacYOHAg3/jGN5yOYtppxpe/RX1dLds+WOV0FOPGk18A6UCBqhaqai2wEmg5Etci4CVtkgH0F5GhqlqmqtsBVLUK2AfE+DC/CVA7d+5k69atfO9736NHD7voi7+LGhZH0pQZZK9bzblz55yOY1w8KQAxgPvVHUr44od4m8uISDxwFeA+ROC9riaj5SIyoLWdi8idIpItItmVlZUexDX+bkVmMf/70d/TM6IPoaOvsZ4jAeKamxZz4fw5XnjhBaejGBdPCkBrXS/Um2VEpDfwJnC/qp51TX4KSABSgDLg8dZ2rqrPqGqaqqZFRUV5ENf4u+IDuzm4K5NpN36N8B49nY5jfGTYyDEkT/0Sy5cvt0HiughPCkAJEOv2eDhQ6ukyItKNpg//l1X1780LqGq5qjaoaiPwLE1NTSbINTY28v4r/02fAYOYduPXnI5jfOy6b36f+vp6/vznPzsdxeBZAcgCEkVkpIiEA4uB1S2WWQ18x9UbaDpwRlXLpKnj9vPAPlX9L/cVRGSo28NbgN1X/CxMwHjnnXcoO3yAOd/4rl3wPQBFDh7G4sWLefPNNzl48KDTcYJemwVAVeuBe4G1NB3EfU1V94jIXSJyl2uxNUAhUEDTt/l7XNNnArcD17XS3fO3IpInIrnAHOAnPntWxi/V1NTwxBNPMHTkGCZdfb3TcUwHufvuu4mIiODxx1tt9TWdKMyThVxdNNe0mPa0230FlrWy3mZaPz6Aqt7uVVIT8P7yl79w/Phxlv7sQcSu9xuw3jtQRfqCb7Hutef59fNvMWpCKgC3TYtzOFnwsf8y0yUcO3aM5557jhtvvJG4pIlOxzEdLH3uzURGx7DmxT9Re/GC03GClhUA4zhV5Ve/+hUiwk9/+lOn45hOEBYezk3f/wmnKsv46O8vOh0naFkBMI5755132LRpE/fddx8xMXaeYLAYkTSR1OtuIvP9tzl2aL/TcYKSFQDjqOPHj/PrX/+ayZMn8+1vf9vpOKaTXf/N79NnQCSrn3+c2tpap+MEHSsAxjGNjY387Gc/o66ujkcffZTQ0FCnI5lO1r1nBAuX3scnpcU88cQTTscJOh71AjKmIyz7xW/JyMjgy9+9j4/LQ/i43IZ8CEaJk9NJu/4rvPDCC6SkpHDjjTc6HSlo2C8A44jMzEw+euuvjJ8+m6tmLXA6jnHYDUt+yKRJk3j44Yc5fPiw03GChhUA0+mKi4u5//77iRwSw03fvd+u9GUIDevG73//e8LDw/nRj35EdXW105GCghUA06nOnDnDPffcg6qy+P5f2WBv5lPri+tY+IP/SWHhYRZ//x7+uqXQ6UgBzwqA6TQ1NTXcfffdFBcX88c//pHIaOvyaT5v1PgpLFz6Iwpyt7Hq2d/R0NDgdKSAZgXAdIrz58+zbNkydu3axWOPPUZ6ug3+alo3ZfZCrv/WHezJ3MAjjzxC00gzpiNYLyDT4Z7/cA8rn/gFxw7l85U7HuCTfmPtIi/msmYs/BY11ed49dVXCQkJ4eGHH7Zuwh3ACoDpUJWVlbz06E85UVbC15c9THLaNU5HMn7ium98j7HRESxfvpzy8nJ+97vf0atXL6djBRRrAjIdZteuXdx6662cqihj8U/+j334G680jw31i1/8gg0bNrB06VLssrC+ZQXA+Jyq8tJLL3H77bcTEhLCdx76HaPGT3E6lvFDKzKLIeEavvmj/yT/wEEW3PRVfvnfrzkdK2BYATA+VVpayrJly/jNb37Dl770Jd58802GjRzjdCzj58ZcNZ3v/+8/EtFvAK/+/j/51a9+RU1NjdOx/J740xH2tLQ0zc7OdjqGaUVtbS0/eeQPbFr9MgrM+fp3mXbDLXaSl/Gp+rpa1r/5Ahn/fJOhQ4fywAMPsHDhQkLsAkKXJSI5qpr2helWAEx71NbWsnr1ap599lmKi4tJSp3JjbfdRb+Bg52OZgJYcX4ea1c8zfGiAoaNHMN13/w+8ckpfHv6CKejdUlWAIxPnTp1ilWrVvHiiy9y/Phxxo8fz6QblzB60lSno5kgoY2N5G5dx/o3XqDq1CcMjh3J/Xf/D7785S8THh7udLwuxQqAabfq6mq2bNnCP/7xDz5cv4HGhnrixkzgmq8sYdSEVGvuMY6or60lb+uHZL7/dyqPFdG3b1/mzp3LwoULmTZtGmFh1tvdCoDx2sWLF9m7dy87duxg8+bNZGdnU1dXx6BBgxidNotJM+cRHTvS6ZjGAE29z0Y2HmPVqlWsW7eO6upq+vbtS3p6OldffTVTp05l1KhRQXlCWbsKgIjMB/4AhALPqeqjLeaLa/5C4DzwXVXdfrl1RSQSeBWIB44A31LVU5fLYQWgY1y4cIGSkhIKCws5cuQIhYWFHDp0iPz8fOrq6gCIihnB6EnpjJ6URtyYiYQE4T+R8R91tRcpyM2iIDeLw3u2c+ZEBQC9evUiOTmZpKQkRowYQVxcHCNGjCAmJiagm42uuACISChwAJgHlABZwBJV3eu2zELgRzQVgGnAH1R12uXWFZHfAidV9VEReQgYoKr/cbksXbkANL+Ont56u17L5RsaGqirq6Ouro7a2tpW71+4cIHq6mqqqqo4d+4c586do6qqiqqqKiorKz/9q6qq+tw++kYOYuCQ4QyNH0NMQjLDE8bSu3+kD14lYzqfqnKqopSjBXs5fqSA0sMHqCg5Qu2F858uIxJCrz79iOjXn9GxQxk4cCCRkZFERETQu3dvIiIiiIiIoFevXkRERNCtW7fP/YWFhX3utlu3bogIISEhiMinzaPuj92nd7RLFQBPGsfSgQJVLXRtaCWwCNjrtswi4CVt+nTKEJH+IjKUpm/3l1p3ETDbtf6LwAbgsgXgSj366KO8/vrr7fpwbmtZfyAhIfToGUGPXr2J6DeA3pHDGDtyAr0HDKRf5GAGDYslMjqG7j3tdHsTOESEyOgYIqNjmDxzHtD0/3u+6gynKko5WV7KyYpSzp06QXXVaYqOn2TvgUNUV52h7uKFTsnXWmFoWSD+9Kc/MXPmTJ/u25MCEAMcdXtcQtO3/LaWiWlj3WhVLQNQ1TIRabXfoIjcCdzpenhORPI9yOxPBgGfOB3CYfYa2GsA9hrAZV6Da65p11AqrfaP9aQAtPYbpeXX30st48m6l6WqzwDPeLOOPxGR7NZ+mgUTew3sNQB7DaDzXwNPTp8rAWLdHg8HSj1c5nLrlruaiXDdVnge2xhjTHt5UgCygEQRGSki4cBiYHWLZVYD35Em04Ezruady627Gljqur8UWNXO52KMMcYLbTYBqWq9iNwLrKWpK+dyVd0jIne55j8NrKGpB1ABTd1Av3e5dV2bfhR4TUTuAIqBb/r0mfmPgG3e8oK9BvYagL0G0MmvgV+dCGaMMcZ3bAg9Y4wJUlYAjDEmSFkBcIiIzBeRfBEpcJ0JHXRE5IiI5InIThHpmqd4dwARWS4iFSKy221apIh8ICIHXbcDnMzY0S7xGvxSRI653g87XSMMBCQRiRWR9SKyT0T2iMh9rumd+j6wAuAA1xAZTwILgHHAEhEZ52wqx8xR1ZQg6//9AjC/xbSHgHWqmgiscz0OZC/wxdcA4AnX+yFFVdd0cqbOVA/8u6omA9OBZa7PgE59H1gBcManw2uoai3QPESGCQKquhE42WLyIpqGRMF1e3NnZupsl3gNgoaqljUPmKmqVcA+mkZO6NT3gRUAZ1xq6Ixgo8D7IpLjGvIjmH1uaBQgWC+pdq+I5LqaiAK6GayZiMQDVwGZdPL7wAqAM9o9REaAmKmqU2hqClsmItc6Hcg46ikgAUgByoDHHU3TCUSkN/AmcL+qnu3s/VsBcIYnw2sEPFUtdd1WAG/R1DQWrIJ+aBRVLVfVBlVtBJ4lwN8PItKNpg//l1X1767Jnfo+sALgDE+G1whoIhIhIn2a7wM3ALsvv1ZAC/qhUZo/+FxuIYDfD66LaD0P7FPV/3Kb1anvAzsT2CGuLm6/57MhMv6vs4k6l4iMoulbPzQNSbIiWF4DEXmFpmthDALKgf8E3gZeA+JwDY2iqgF7kPQSr8Fsmpp/lKarBP6wuT080IjINcAmIA9odE3+XzQdB+i094EVAGOMCVLWBGSMMUHKCoAxxgQpKwDGGBOkrAAYY0yQsgJgjDFBygqAMcYEKSsAxhgTpP4/O/h7A3GlfQoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "samples = df.iloc[:,0]*dummy_mean[0]\n",
    "sns.distplot(x=samples, fit=norm, kde=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63fcd173",
   "metadata": {},
   "source": [
    "## Lognormal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 456,
   "id": "a2eb140e",
   "metadata": {},
   "outputs": [],
   "source": [
    "normal_equivalent_cov = (np.exp(dummy_cov)-1)*np.prod(dummy_center) # vale per una matrice 2x2 va bene la divisione per 1 ??\n",
    "normal_equivalent_mean = np.log(dummy_center)-dummy_std**2/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 457,
   "id": "4a3f2f2e",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.DataFrame(np.random.multivariate_normal(normal_equivalent_mean, normal_equivalent_cov, size=100000)).apply(np.exp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 458,
   "id": "e5e83021",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\fgrimald\\Anaconda3\\envs\\sandy-devel\\lib\\site-packages\\seaborn\\distributions.py:2619: FutureWarning: `distplot` is a deprecated function and will be removed in a future version. Please adapt your code to use either `displot` (a figure-level function with similar flexibility) or `histplot` (an axes-level function for histograms).\n",
      "  warnings.warn(msg, FutureWarning)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "covariance:\n",
      "             0            1\n",
      "0  5.27697e-02 -5.41199e-02\n",
      "1 -5.41199e-02  6.52275e-02\n",
      "\n",
      "\n",
      "mean: [1.0011624892185653, 1.0002198407308642]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAApI0lEQVR4nO3deXiU533u8e9PowUQAoEkFiFAQoh9ESB2GePgUHAc48SmgXhJUyccfOz2pPVJ4iZxcly3OTlX0jSOl7q0iXfwsWOHuDGLgzHGGBAIELIEzAi0oAWQxC6xaHv6h0bOWEhoQDN6Zvl9rmsuRu8yc2s8vnl5l+cVYwxKKaVCV4TtAEoppfxLi14ppUKcFr1SSoU4LXqllApxWvRKKRXiIm0H6EhiYqJJTU21HUMppYLGvn37ao0xSR3NC8iiT01NJTc313YMpZQKGiJS1tk83XWjlFIhToteKaVCnBa9UkqFOC16pZQKcVr0SikV4rTolVIqxGnRK6VUiNOiV0qpEKdFr5RSIS4gr4xVgWFtznFqKssoc+YTHdOb5FFjSRw6nK/PHmE7mlLqBmjRqw5VVVXx8k8f47ir4M8TRRg3fR63p/2UQYMG2QunlLohWvTqGgUFBXz729/m0pUGFq/8H4yZPo/mpkYKdm1l96a3WbFiBS+88AJjxoyxHVUp5QUJxHvGZmVlGR3UzI5Tp06xfPlyoqOjWfY3/0jCkJTPzT9RWsQfnnuSlpYW3nzzTZKTky0lVUp5EpF9xpisjubpwVj1mebmZv72b/+W+vp6nn/++WtKHmBoagYvvvgiV65c4dFHH+Xy5csWkiqlboQWvfrMG2+8QX5+Pk8++eR1d8ukp6fzi1/8gsOHD/P000/3YEKl1M3QolcA1NbW8vTTTzN37ly+9KUvdbn8woUL+drXvsarr75Kfn5+DyRUSt0sLXoFwJo1a7h8+TJPPPEEIuLVOo899hhJSUn85Cc/oaWlxc8JlVI3S4tecfr0ad566y3uuusu0tLSulx+bc5x1uYc578OnWXeV77JkSNH+MHTL/dAUqXUzejy9EoR+S1wJ1BtjJnUwfzvAvd5vN54IMkYc0ZESoGLQDPQ1NkRYWXXq6++ypWrVxk88w7W5hy/oXUnzrqVXRt/x0fvvEzDI/cRHR3tp5RKqZvlzRb9S8CSzmYaY35ujMk0xmQC/wB8ZIw547HIbe75WvIB6NKlS6xdu5bxM7JJHDr8hteXiAi+cO83OVd7irffftsPCZVS3dVl0RtjtgNnulrObSWwrluJVI/avHkzFy9eZNYXl930a4yaNINho8bx4osv0tTU5MN0Silf8Nk+ehHpQ+uWv+dmnQHeF5F9IrLKV++lfOett94iLS2N4WOu2SvnNRFh7h3LKS8v509/+pMP0ymlfMGXB2O/DHzSbrfNfGPMdGAp8IiILOhsZRFZJSK5IpJbU1Pjw1iqM0VFRRw4cIB7773X6zNtOjN2+lxSU1P5zW9+QyBeba1UOPNl0a+g3W4bY0yV+89q4PfArM5WNsasMcZkGWOykpKSfBhLdWb9+vVERkZy9913d/u1IiIcPPjggxQWFup59UoFGJ8UvYj0B24F/uAxLVZE4tqeA4uBgo5fQfU0YwwbN25k3rx5DBw40Cevedddd9GnTx/WrdPDNEoFki6LXkTWAbuAsSJSISIPichqEVntsdhXgPeNMfUe0wYDO0TkILAHeM8Ys8mX4dWNazsH/hfr/sSJEyeIHzPrhk+p7ExsbCzLli1j48aNnD171ievqZTqvi7PozfGrPRimZdoPQ3Tc1oxMPVmgyn/OrTnIxyRUYydNtdnr7k25zj9Ji6kYd06fvz0i8xdeq/epESpAKBXxoYhYwyH9n7MqEkz6BXb16evPSgllWHp4zm44309KKtUgNCiD0Mny45y4UwN42fM98vrT83+IjWVZZwoLfLL6yulbowWfRgqOpgDIoye2ulJUN0ycfZCIqOiOfjx+355faXUjdGiD0NFB/cwbNQ4YvvF++X1e/WJZdyM+RTkfEhDQ4Nf3kMp5T0t+jBTd/4sVcVOMvy0Nd9m0pzbuFJfx86dO/36PkqprmnRh5mj+XsAyJg626/vM2rSdHrF9mXjxo1+fR+lVNe06MPMsU/30Td+IINHjPLr+zgioxg3fT4ffPABV69e9et7KaWuT4s+jLS0tFB6OI+0CdO6PbaNNybMvpX6+np27Njh9/dSSnVOiz6MOJ1OLl08T9rEaT3yfqnjphIfH6+7b5SyTIs+jOzevRuAtAk9U/SOyEgWL17Mhx9+yOXLl3vkPZVS19KiDyO7du0iYehw+g1I7LH3XLp0KZcuXWL79u099p5Kqc/Tog8TDQ0N5ObmMqqHtubbzJw5k8TERN19o5RFWvRhorCwkMuXLzNyfM+OM+dwOFi8eDEfffQRly5d6tH3Vkq10qIPE7m5uQCM6MYtA2/W7bffzpUrV9i1a1ePv7dSSos+bOTm5jJq1Ci/DXtwPVlZWcTFxbF169Yef2+llBZ9WGhubmb//v1kZWVZef+oqCgWLFjAtm3baG5utpJBqXCmRR8GXC4XdXV1Voq+7Y5WMSOmcObMGX6x7n2f3dFKKeUdLfow0LZ/fsaMGdYypE+eSYTDgfOA7qdXqqdp0YeB/fv3M3ToUJKTk61l6NUnlpHjpuDar0WvVE/z5ubgvxWRahEp6GT+QhE5LyJ57sePPeYtERGniBwVkcd9GVx5Ly8vj2nTevb8+Y6MnTaX0ycrqD1RbjuKUmHFmy36l4AlXSzzsTEm0/34RwARcQDPAUuBCcBKEZnQnbDqxp08eZKTJ0+SmZlpOwoZmXMAcB3YbTmJUuGly6I3xmwHztzEa88Cjhpjio0xDcAbwLKbeB3VDQcPHgRg6tSevVCqI/GJgxk8Ih1Xnu6+Uaon+Wof/VwROSgiG0VkonvaMMDz3+gV7mkdEpFVIpIrIrk1NTU+iqXy8vKIiYlh3LhxtqMAMGbaHMqLDnH27FnbUZQKG74o+v3ASGPMVOAZYL17ekcDnpvOXsQYs8YYk2WMyUpKSvJBLAWtRT9x4kSio6NtRwFovYWhMTpGvVI9qNtFb4y5YIypcz/fAESJSCKtW/DDPRZNAaq6+37Kew0NDRQWFgbE/vk2yalj6BPXn48//th2FKXCRmR3X0BEhgCnjDFGRGbR+pfHaeAckCEiaUAlsAL4enffT3lnbc5xqoqdNDY2cqF3csBcpCQREYyaNIMdO3bQ0tJCRISe4auUv3lzeuU6YBcwVkQqROQhEVktIqvdi9wLFIjIQeDXwArTqgl4FNgMHAbeNMYU+ufXUB2pKnUBMDRtjOUknzd68kzOnj1LQUGHZ+wqpXysyy16Y8zKLuY/CzzbybwNwIabi6a660TpUfrE9ad/wiDbUT4nffIMRISPP/6YKVOm2I6jVMjTfzeHsBMlLoamZvTIjcBvRJ+4/kyePFn30yvVQ7ToQ1Rjw1WqK0sZmpphO0qHbrnlFvLz8/U0S6V6gBZ9iDp1vBjT0sLQ1MDaP98mOzsbYww7d+60HUWpkKdFH6JOlBYBkJwWmFv0kydPJj4+Xm8arlQP0KIPUSdKXcT2G0DcgETbUTrkcDiYP3/+Z6dZKqX8R4s+RFWVFgXkgVhPt9xyC2fOnOHQoUO2oygV0rToQ9ClS5eorTzO0ADdbQOtF3RV90kD4Nm1//XZnaiUUr6nRR+CnE4nxrSQHKAHYtvE9otnaNoYjn2613YUpUKaFn0IarviNFBPrfQ0enIWlcecXK67YDuKUiFLiz4EFRQU0Dd+IHEDEmxH6dLoKTMxpoVjBfttR1EqZGnRh6DCwsKAPX++veRRY+kdG6e7b5TyIy36EFNfX09xcXHAnj/fXkSEg7SJ0yku2I8xnd6uQCnVDVr0Iab1QKxhyMjRtqN4bfTkLOrOn+FUebHtKEqFJC36EON0OgEYPHyU5STeGzVpBgDFBfssJ1EqNGnRhxin00m/fv3oNzB4bscYNyCBQSlpHPs013YUpUKSFn2IKSoqIiMjsK+I7Uj6pBkcdxVSX19vO4pSIUeLPoQYY3A6nYwdO9Z2lBs2avIMWpqb2LNnj+0oSoUcLfoQUlVVRX19fVAW/YiMSURGx/DJJ5/YjqJUyNGiDyFtB2IzMoLj1EpPkdHRpI6bwo4dO2xHUSrkeHNz8N+KSLWIdHgnZxG5T0Ty3Y+dIjLVY16piHwqInkiokfa/Mzlar0Z+JgxwXGxVHvpk7MoKyujoqLCdhSlQoo3W/QvAUuuM78EuNUYMwV4CljTbv5txphMY0zWzUVU3nI6nQwfPpzY2FjbUW5K22mWulWvlG91WfTGmO3AmevM32mMabvx524gxUfZ1A1yuVxBuzUPkDAkheTkZC16pXzM1/voHwI2evxsgPdFZJ+IrLreiiKySkRyRSS3pqbGx7FC39WrVyktLQ3qohcRsrOz2b17N42NjbbjKBUyfFb0InIbrUX/fY/J840x04GlwCMisqCz9Y0xa4wxWcaYrKSk4LnYx7a2G3Y8u771lnzVEQlBfQOP+fPnU19fT15enu0oSoUMnxS9iEwB/hNYZow53TbdGFPl/rMa+D0wyxfvp65VXV4CwKCUVLtBumnu3Lk4HA49zVIpH+p20YvICOAd4AFjjMtjeqyIxLU9BxYDHZ65o7qvuqKEyKhoBg5Oth2lW+Li4pg6darup1fKh7w5vXIdsAsYKyIVIvKQiKwWkdXuRX4MJADPtzuNcjCwQ0QOAnuA94wxm/zwOyjgVEUJSSmpREQ4bEfptuzsbA4dOsSZM52eA6CUugGRXS1gjFnZxfxvAd/qYHoxMPXaNZQ/VJeXkDF1tu0YPpGdnc2vf/1rdu7cyZ133mk7jlJBT6+MDQF1589Sf+Fc0O+fbzNhwgTi4+N1941SPtLlFr0KfNUV7gOxQTQGfWfazhgaNjaTLR9u5/VdpUhEBF+fPcJyMqWCl27Rh4BQOePGU/qkGdRfOMsp9++mlLp5WvQhoLqihL79BxLbL952FJ9pGw7hWIEOkaRUd2nRh4BT5SUMGp5mO4ZPtd11Sm8vqFT3adEHuZbmZmoqyxicElpFD62jWR53FdJw5bLtKEoFNS36IHf6VCXNTY0kDU+1HcXn0t13nSo9ctB2FKWCmhZ9kGs742ZwSvCfcdPe8IyJREXHUPyp7r5Rqju06INcdXkJEhFBYvJw21F8LjIqmpHjpnJM99Mr1S1a9EHuVHkJCUNSiIyKth3FL9Inz+DMqUrKy8ttR1EqaGnRB7maihIGh8CFUp3Ru04p1X1a9EHs4sWLnKs9FXKnVnpKGJJC/4TBOmyxUt2gRR/EioqKgNC6IrY9ESF98gy965RS3aBFH8Rcrtbh/0N51w20nk+vd51S6uZp0Qcxp9NJTO9Y+g0M7Vsvpo7P1LtOKdUNWvRBzOVyMWh4KiJiO4pf9eoTS2Zmph6QVeomadEHKWMMLpcrJC+U6sj8+fMpLCzk9OnTXS+slPocLfogVVVVRV1dXUifceMpOzsbgJ07d1pOolTw0aIPUn8+EBseRa93nVLq5nlzc/Dfiki1iBR0Ml9E5NciclRE8kVkuse8JSLidM973JfBw53T6QQgaViq3SA9xOFwMH/+fHbu3ElLS4vtOEoFFW+26F8Cllxn/lIgw/1YBfwbgIg4gOfc8ycAK0VkQnfCqj9zuVykpKQQ07uP7Sg9Jjs7m9ra2s/+klNKeafLojfGbAfOXGeRZcArptVuIF5EhgKzgKPGmGJjTAPwhntZ5QMul4uxY8fajtFj1uYc51TvkQA88/q7n91bVinVNV/sox8GeI44VeGe1tn0DonIKhHJFZHcmpoaH8QKXVevXqW0tJSMjAzbUXpUXHwCQ9PGUHRwj+0oSgUVXxR9Rydxm+tM75AxZo0xJssYk5WUFNoXAHXXsWPHaG5uDqst+jYZU2dTcewwly6etx1FqaDhi6KvADwHQ08Bqq4zXXVT2xk34Vj0YzJngzEczd9rO4pSQcMXRf8u8KD77Js5wHljzAlgL5AhImkiEg2scC+rusnpdBITE8OIESNsR+lxQ0ak0zd+IEUHc2xHUSpoRHa1gIisAxYCiSJSAfwEiAIwxrwAbADuAI4Cl4Bvuuc1icijwGbAAfzWGFPoh98h7LhcLkaPHo3D4bAdpcdJRAQZU2ZxaO92GhsbiYqKsh1JqYDXZdEbY1Z2Md8Aj3QybwOtfxEoH3K5XCxYsMB2DGsyMmdzYPsm9u3bx5w5c2zHUSrg6ZWxQeb06dPU1taG3Rk3ntImTMMRGcW2bdtsR1EqKGjRB5lwPhDbJrpXb1LHT9WiV8pLWvRBZG3Ocdb9qfUgZEFdbFhfNJQxdTZlZWWUlJTYjqJUwNOiDzLV5cXE9htAbL9421GsysicDcBHH31kOYlSgU+LPshUV5SGzYiV1xOfOJgxY8bo7hulvKBFH0RampupqSwLmzHou7Jw4UL27dvHhQsXbEdRKqBp0QeRM6cqaWpsYFCKFj20Fn1TU5PeS1apLmjRB5FTFa0HHnWLvtWUKVOIj4/X3TdKdUGLPohUl5ciEREkDQ2/oQ864nA4WLBgAdu3b6epqcl2HKUClhZ9EKmuKCZhSAqR0dG2owSMRYsWce7cOfbv3287ilIBS4s+iFSXl+j++Xays7OJiYlhy5YttqMoFbC06INEXV0d52pP6amV7fTp04f58+ezZcsWWoddUkq1p0UfJIqKigA9EOtpbc5x1uYcp2/aNE6cOMG//P8PwvpqYaU6o0UfJNpuiK27bq6VMW0OEhHBkX16mqVSHdGiDxIul4uY3n3onzDIdpSA06dvP0aOnYJz/07bUZQKSFr0QcLlcjFoeBoiHd2KV42bMY/aquPUVumuG6Xa06IPAsYYnE6n7ra5jrHT5wFwRLfqlbqGFn0QqKyspK6ujsHDR9mOErD6DUwiOW0MTt1Pr9Q1vCp6EVkiIk4ROSoij3cw/7sikud+FIhIs4gMdM8rFZFP3fNyff0LhIMjR44ArTfGVp0bNyObqhIXlZWVtqMoFVC6LHoRcQDPAUuBCcBKEZnguYwx5ufGmExjTCbwD8BHxpgzHovc5p6f5bvo4cPpdCIiJKWk2o4S0MbPugWAzZs3W06iVGDxZot+FnDUGFNsjGkA3gCWXWf5lcA6X4RTrY4cOcLIkSOJjullO0pAGzgomaGpGWzcuNF2FKUCijdFPwwo9/i5wj3tGiLSB1gCvO0x2QDvi8g+EVl1s0HDmdPpZNy4cbZjBIUJs26loKCA8vLyrhdWKkx4U/Qdnc/X2bXmXwY+abfbZr4xZjqtu34eEZEFHb6JyCoRyRWR3JqaGi9ihYe6ujrKy8vD+mbgN2LCrNav16ZNmywnUSpweFP0FcBwj59TgKpOll1Bu902xpgq95/VwO9p3RV0DWPMGmNMljEmKykpyYtY4aHtiljdovdOfOJgpkyZortvlPLgTdHvBTJEJE1Eomkt83fbLyQi/YFbgT94TIsVkbi258BioMAXwcNF2xk3WvTeW7p0KYcPH6a0tNR2FKUCQpdFb4xpAh4FNgOHgTeNMYUislpEVnss+hXgfWNMvce0wcAOETkI7AHeM8bov6lvgNPppH///gwePNh2lKCxZMkSAN2qV8ot0puFjDEbgA3tpr3Q7ueXgJfaTSsGpnYrYZhrOxCrQx94b8iQIUyfPp1Nmzbx8MMP246jlHV6ZWwAa25uxuVy6YHYm3DHHXfgcrlwuVy2oyhlnRZ9ACsrK+PKlSuMHz/edpSgs2TJEhwOB+++e83hJKXCjle7bpQdbQdidYv+xrTdfGTUpBm8+c4fGDL/HiIiHHx9tt5UXYUn3aIPYE6nk8jISNLTdYybmzF53iIunq2l7Ei+7ShKWaVFH8COHDnCqFGjiI6Oth0lKI2ZNpeY3n3I/+QD21GUskqLPoAdOXJEz5/vhqjoGMbPXMCRfTtovHrFdhylrNGiD1C1tbVUV1frgdhumjJvEQ1XLusNSVRY06IPUIWFhQBMnDjRcpLgNmLMJPonDOLTnbr7RoUvPesmAK3NOc72zTtBhIL6OIpy9D6oN0siIpg09wvsfO9NqqurGTRIb66uwo9u0QeoE6VFJAweRkzvPrajBL2p2YsxpoX169fbjqKUFVr0AepkWRFDUzNsxwgJCUOGMXLsFH73u9/R0tJiO45SPU6LPgDVnT/LhTO1DE0bYztKyJh261LKy8vJycmxHUWpHqdFH4BOlBYB6Ba9D43PyqZ///689dZbtqMo1eO06APQidIiEGHICL0i1lcio6O566672LJlC2fPnrUdR6kepUUfgPRArH8sX76cxsZG/vCHP3S9sFIhRIs+AOmBWP/IyMggMzOTN998E2M6u+2xUqFHiz7A1NbW6oFYP1q+fDklJSXs3bvXdhSleowWfYBpuyJWt+j9Y+nSpfTv35/XXnvNdhSleoxeGRtgCgsL9UCsn7SNUz8xewlbNrzF83/cQ3zSEB2nXoU8r7boRWSJiDhF5KiIPN7B/IUicl5E8tyPH3u7rvq8Q4cO6YFYP8v6wp0gkLv1j7ajKNUjuix6EXEAzwFLgQnAShGZ0MGiHxtjMt2Pf7zBdZVbYWGh7rbxs/4Jgxg3Yz4HPtpIgw5frMKAN1v0s4CjxphiY0wD8AawzMvX7866Yae6upqTJ0+SrAdi/W7WF+/myqU6HdVShQVvin4YUO7xc4V7WntzReSgiGwUkbaxdb1dFxFZJSK5IpJbU1PjRazQk5eXB0BKhv6jx9+GZ0xkyMjR7PnTej3VUoU8b4peOpjW/v+M/cBIY8xU4Blg/Q2s2zrRmDXGmCxjTFZSUpIXsULPgQMHiImJ0QOxPUBEmPXFu6mtOs6OHTtsx1HKr7wp+gpguMfPKUCV5wLGmAvGmDr38w1AlIgkerOu+rO8vDwmTZqEIzLKdpSwMGnOQuIGJPIf//EftqMo5VfeFP1eIENE0kQkGlgBvOu5gIgMERFxP5/lft3T3qyrWl29epXCwkIyMzNtRwkbjsgo5iy5h71793LgwAHbcZTymy6L3hjTBDwKbAYOA28aYwpFZLWIrHYvdi9QICIHgV8DK0yrDtf1xy8S7A4dOkRjYyPTpk2zHSWsTL91KfHx8axZs8Z2FKX8xqsLpty7Yza0m/aCx/NngWe9XVddq22LMjMzk81H6y2nCR/RvXpz//338+yzz+J0Ohk7dqztSEr5nA6BECDy8vIYMWIECQkJtqOEnfvuu48+ffrovnoVsrToA4AxhgMHDuj+eUvi4+NZsWIFGzdupKyszHYcpXxOiz4AVFZWUltbq0Vvydqc4/SbcjsRjki+++TPPxsTR6lQoUVv0dqc46zNOc5zb20B4ETkEC0ZS+LiE5h5+zI+3b2V6opS23GU8ikt+gBQcewQ0b16Mygl1XaUsDbvS39JTK/ebHvnZdtRlPIpLfoAUHH0MMPSxxER4bAdJaz16duPOUvuxbl/JwcPHrQdRymf0aK37OrlS5w6XkxK+njbURQwe/FX6BPXn1/96le2oyjlM1r0lpUXFWJMCyPGTrYdRQExvfuQfecKdu/ezSeffGI7jlI+oUVvWenhg0Q4Ihk+WkesDBQzbruTlJQUfvazn9HU1GQ7jlLdpkVvWemRg6SkjyMqppftKMotMjqa73//+xw9epQ33njDdhyluk2L3qIrl+o5WXqUkeOn2o6i2lm0aBHz5s3jmWee4cyZM7bjKNUtWvQWHXd+ijEtpI7Tog806/aUM+XOb1JXX8/f/PCfP7vmQalgpEVvUemRgzgio/SMmwCVlDyCmYvuYv9HGzlRWmQ7jlI3TYveouKCfYwYM5HI6GjbUVQnFiy7n9i4eP744q9o1gOzKkhp0Vty8uRJairLGDUpy3YUdR29Yvuy9MFHOVl2lF2bfmc7jlI3RYvekrZztNMnz7CcRHVlfFY247NuYfv61zh27JjtOErdMC16S3bs2EHf+IEMSkmzHUV5YckDjxAd04sf/ehHNDc3246j1A3RoregubmZXbt2MWriDNy32lUBrm//ASy+72Hy8vJ45ZVXbMdR6oZo0VuQn5/P+fPnSZ+s++eDyeS5X+D222/nX//1XykoKLAdRymveVX0IrJERJwiclREHu9g/n0iku9+7BSRqR7zSkXkUxHJE5FcX4YPVh988AGRkZGM1qIPKiLCU089RWJiIn//939PXV2d7UhKeaXLohcRB/AcsBSYAKwUkfYDs5QAtxpjpgBPAWvazb/NGJNpjNFmA7Zu3crMmTPpFdvXdhR1gzY4L/AXf/1dKior+atHv8vru8v0QioV8LzZop8FHDXGFBtjGoA3gGWeCxhjdhpjzrp/3A2k+DZm6CgpKaGkpIRFixbZjqJu0vCMiSz8yoMU5mzjwEcbbcdRqkveFP0woNzj5wr3tM48BHh++w3wvojsE5FVna0kIqtEJFdEcmtqaryIFZy2bt0KwG233WY5ieqOeV/6S9InzWDjq89x3Pmp7ThKXZc3Rd/RaSGmwwVFbqO16L/vMXm+MWY6rbt+HhGRBR2ta4xZY4zJMsZkJSUleRErOL3//vtMmDCB5ORk21FUN0REOPjqwz8gPmkIbz37FJWVlbYjKdUpb4q+Ahju8XMKUNV+IRGZAvwnsMwYc7ptujGmyv1nNfB7WncFhaWKigry8/NZunSp7SjKB3rF9mXFd56kpbmZhx9+mPr6etuRlOqQN0W/F8gQkTQRiQZWAO96LiAiI4B3gAeMMS6P6bEiEtf2HFgMhO15aRs2bADQog8hCUNSuOeRH1JcXMx3vvMdGhoabEdS6hpdFr0xpgl4FNgMHAbeNMYUishqEVntXuzHQALwfLvTKAcDO0TkILAHeM8Ys8nnv0WAaxvi9vXfrSdl9AQ+qmjWMzVCyKiJ03nyySfZsWMH3/ve9/SuVCrgRHqzkDFmA7Ch3bQXPJ5/C/hWB+sVAzrYOlBdWUp1eQl/cd//tB1F+cE999xDXV0dP/vZz4iNjeWpp54iIkKvR1SBwauiV92Xt30zEY5IJs6+1XYU5Qdrc44TNe42Fiyr4J13XqPsXAOvPPP/tOxVQNBvYQ9oamwg/5MtjJ0+l9h+8bbjKD9acPf9zL1jOfu2/pHHH3+cxsZG25GU0i36nuA8sIvLdReYdqsehA11IsKi5Q/Rq3cs//X2S9TX1/PLX/6SmJgY29FUGNMt+h6w/8P36J8wiFETptmOonqAiJD95ZU88cQTbN26lYceekhvMK6s0qL3s8OHD1N6+CBZi76M6P7a8JKezVcf/gcO5n/KHXd9lV++9aGebaWs0ObxsxdffJHoXr2ZfusdtqMoCybOXsg3fvAvtLQ089I//R2FOdtsR1JhSIvej06cOMHGjRvJvOUvdKTKMJacNoZv/eQZBo9I551/+7/88Ic/1KtoVY/Sovejf//3f0dEmL34K7ajKMv6xg/kwcd/TvaXV/L73/+er371q+Tn59uOpcKEFr2flJWV8fbbb7N8+XLik4bYjqMCgCMyktvu+StefvllGhsbWblyJT/96U/1BibK77To/eSZZ54hMjKS1atXd72wCitFLYO5/0fPMn3hl3j1tdf4wheX8KNnXsOYDgeFVarbtOh9bG3OcX760ru89957zFz8Vf5UfNl2JBWAesX2ZemDj/LXP/oVfeL68/bz/8yKFSvYs2eP7WgqBGnR+1hTYwMbXn6GAUlDmX/nCttxVIAblj6Ob/2fZ7nzm3/HqVOn+MY3vsGqVavIzc3VLXzlM1r0Pvbh2y9x+mQFSx54hKhovRpSdS3C4WDarUvYtGkTjz32GAUFBTzwwAN8/etfZ8uWLTQ3N9uOqIKcBOJWQ1ZWlsnNze16wQCzbds2Hn74YbK+8GWWPvio7TgqSDVevULex5vZveltztWeYsiQIdx7773cc889DBmiB/ZVx0RknzEmq8N5WvS+UVRUxP3330/v/on89RNPExkdbTuSCnItzc248naz/8P3OFa4H4ARYyYxYeYCfrDqa4TyLTfVjdOi97OKigruv/9+WlpaWPH9fyE+cbDtSCrEnK05Sf4nWzi8dzs1lWWICDNmzGDRokXMmzePjIwMRDq6vbMKF1r0fnTkyBG+/e1v09DQwCuvvMK+c71tR1IhrqayjMgT+WzatImjR48CkJCQwJw5c5gzZw7Tpk0jLS1Nx8IPM1r0fvD67jLytm9m8+vP07tvHCsf+2cGDUu1HUuFmfOnqyk9lEfJoTxKDh2g7nzrKJmxsbFMmjSJyZMnM2HCBEaPHs3IkSOJ1l2KIUuL3scKCwv530/8E6WH80gdP5Vlq75HvwGJtmOpMGeM4fSJcipLnFQVtz5OHi+mpbn1HraRkZGMHDmS9PR0hg8fzrBhwz57JCcn07u3/ms0mHW76EVkCfA04AD+0xjzs3bzxT3/DuAS8FfGmP3erNuRQCz6c+fOsW3bNtavX09OTg69+/bj1rvvZ8YX7iQiwmE7nlIdampsoLaqnJqqMmoqy6itOk5N1XHO156iuenzd7+K7RdP3/gExoxMJiEhgYSEBBITE0lISGDAgAHExcV99ujXr5/eTCXAXK/ou7zDlIg4gOeALwIVwF4RedcYc8hjsaVAhvsxG/g3YLaX6/YoY8xnj5aWls/+vHTpEhcvXqSuro76+npqamooKyujrKwMp9OJy+XCGMOwYcN47LHHiBg1T0ekVAEvMiqaISPTGTIy/XPTTUsLdefPcq72JOdqT3G+9hTnak5Sd/4MRytOcbDwCPUXzl3zl4Gn6Oho+vXrR9++fenduze9e/emV69eHT7a5sXExBAZGUlkZCRRUVFERkbicDiumeb5s8PhwOFwEBERQUREBCLyueftf76ReZ4HsNuedzQt2HlzK8FZwFFjTDGAiLwBLAM8y3oZ8Ipp/efBbhGJF5GhQKoX6/rM/PnzuXz58udKvH2p36h+AxNJTB7JLcvuJ2PKTIamjQmZ//gqfElEBHEDEogbkMDwjIkdLmOM4eqleurOn+Fy/UWuXKrnyqU6rl6u//Nz95+NV69y5vJVmi6cpfHqVRobrtLUeNXjeUMP/4b+09lfBN7+RXG99RMSEtiyZYtP84J3RT8MKPf4uYLWrfaulhnm5boAiMgqYJX7xzoRcXqRrYfs8GahRKDWz0GCgX4O+hm00c+h1Q19Dt3YkBzZ2Qxvir6jd22/adzZMt6s2zrRmDXAGi/yBCQRye1s/1g40c9BP4M2+jm0CoTPwZuirwCGe/ycAlR5uUy0F+sqpZTyI2+uqNgLZIhImohEAyuAd9st8y7woLSaA5w3xpzwcl2llFJ+1OUWvTGmSUQeBTbTeorkb40xhSKy2j3/BWADradWHqX19MpvXm9dv/wm9gXtbicf089BP4M2+jm0sv45BOQFU0oppXxHB8NQSqkQp0WvlFIhTov+BojIEhFxishREXm8g/kLReS8iOS5Hz+2kdPfROS3IlItIgWdzBcR+bX7c8oXkek9ndHfvPgMwuW7MFxEPhSRwyJSKCL/q4NlwuH74M3nYO874Xn1qD46f9B6MPkYMIrW00YPAhPaLbMQ+KPtrD3wWSwApgMFncy/A9hI63UUc4Ac25ktfAbh8l0YCkx3P48DXB38fxEO3wdvPgdr3wndovfeZ0NBGGMagLbhHMKOMWY7cOY6i3w2JIYxZjfQNiRGyPDiMwgLxpgTxj2AoTHmInCY1iviPYXD98Gbz8EaLXrvdTbMQ3tzReSgiGwUkY4HEQl93n5WoS6svgsikgpMA3LazQqr78N1Pgew9J3w5spY1cqb4Rz2AyONMXUicgewntYRPcON10NfhLCw+i6ISF/gbeA7xpgL7Wd3sEpIfh+6+BysfSd0i957XQ4FYYy5YIypcz/fAESJSDjekcSbYTNCWjh9F0QkitZye90Y804Hi4TF96Grz8Hmd0KL3ntdDucgIkPEPfSciMyi9fM93eNJ7etsSIywES7fBffv+BvgsDHml50sFvLfB28+B5vfCd114yXj3VAQ9wIPi0gTcBlYYdyH20OJiKyj9QyCRBGpAH4CRMH1h8QIJV58BmHxXQDmAw8An4pInnvaD4ARED7fB7z7HKx9J3QIBKWUCnG660YppUKcFr1SSoU4LXqllApxWvRKKRXitOiVUirEadErpVSI06JXSqkQ999W4GElq3xiJQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.distplot(x=df.iloc[:,0], fit=lognorm, kde=False)\n",
    "print(f\"\"\"covariance:\n",
    "{df.cov()}\n",
    "\n",
    "\n",
    "mean: {[df.iloc[:,0].mean(), df.iloc[:,1].mean()]}\"\"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 459,
   "id": "b42e30ad",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\fgrimald\\Anaconda3\\envs\\sandy-devel\\lib\\site-packages\\seaborn\\distributions.py:2619: FutureWarning: `distplot` is a deprecated function and will be removed in a future version. Please adapt your code to use either `displot` (a figure-level function with similar flexibility) or `histplot` (an axes-level function for histograms).\n",
      "  warnings.warn(msg, FutureWarning)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<AxesSubplot:>"
      ]
     },
     "execution_count": 459,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD4CAYAAADlwTGnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAAsKElEQVR4nO3de3RU9bn/8feTkDuESwgkQLiFcA8ECaCAigpWsKfoOdUD3qjL/ii19GpbPaftObar9eexVntc9if1QrWnRqutHrGg1OKNmyHBQCCQGwGSQEiiUEIIkNvz+2MmdIwJmVz3XJ7XWrNmZl9mPzOZNZ/s7/7u/RVVxRhjTPAJcboAY4wxzrAAMMaYIGUBYIwxQcoCwBhjgpQFgDHGBKl+ThfQGUOHDtWxY8c6XYYxxviV3bt3f6Kq8a2n+1UAjB07luzsbKfLMMYYvyIiR9uabk1AxhgTpCwAjDEmSFkAGGNMkLIAMMaYIGUBYIwxQcoCwBhjgpQFgDHGBCkLAGOMCVIWAMYYE6T86kxg4zvWv3uAor27aGy4wMC4YYydmsYdV4xzuixjTCdYAJhOaWpq4sknn+SZZ5+jqbHh4vQhCaNIevDHXH311Q5WZ4zpDAsA47X6+nruvfdetm/fzvQrrmXu4uX0HzSE8kMH2frGi3z961/n/vvvZ9WqVU6XaozxggWA8Yqq8h//8R9s376dBx98kKaxV1ycNzBuGBNnzmP3n57k4YcfJjIykn/91391sFpjjDfsILDxyksvvcQbb7zB2rVr2/xxD4uI5PHHH+fKK6/k5z//OVlZWQ5UaYzpDAsA06Gqqioee+wxFixYwNe//vV2lwsNDeXRRx9l1KhR3H///Zw9e7YPqzTGdJYFgOnQI488QkNDAz/5yU8ICbn0VyY2NpaHHnqIEydO8Pjjj/dRhcaYrrAAMJeUl5fHxo0bueeeexgzZoxX68yaNYs77riDjIwM9u/f38sVGmO6ygLAXNIzzzzDgAEDuPvuuztcNiOz9OItcf7NRPWP5Yf/+Ys+qNIY0xVeBYCI3CAiBSJSLCIPtDF/sojsFJELIvJ9j+mTRGSPx61GRL7jnvegiBzzmLesx96V6RElJSVs/utfmbnoi7x54NRnfuA7EhEVw8J/WsnhA3vYvn17H1RrjOmsDgNAREKB3wBLganAShGZ2mqxk8C3gEc9J6pqgaqmqWoaMBuoA173WOTxlvmquqnrb8P0hueff55+YeHMXXJTl9affc2NDBo6nCeeeAJV7dnijDHd5s0ewFygWFVLVLUeeBlY7rmAqlapahbQ0NYLuF0HHFLVNgcnNr7l7NmzbNy4kWnzriYmdlCXXqNfWDjzl91Kbm6udQs1xgd5EwAjgTKP5+XuaZ21Anip1bS1IpIrIutFZHAXXtP0ko0bN1JXV8dlVy/t1uvMWLiEuLg4nn322R6qzBjTU7wJAGljWqf250UkHPgS8KrH5KeAZCANqAB+1c66q0UkW0Syq6urO7NZ0w1/+tOfmDBhAiOTp3TrdcLCI7jjjjvYunUrhYWFPVSdMaYneBMA5UCSx/NRwPFObmcp8LGqVrZMUNVKVW1S1WbgGVxNTZ+jqk+rarqqpsfHx3dys6YriouL2bdvH1/+8pcRaSv/O2fFihVERESQkZHRA9UZY3qKNwGQBaSIyDj3f/IrgA2d3M5KWjX/iEiix9ObAesw7iPefvttRIRly3qmY9agQYO48cYbefPNN6mtre2R1zTGdF+HF4NT1UYRWQtsBkKB9aqaJyJr3PPXiUgCkA3EAs3urp5TVbVGRKKBJcDXWr30IyKShqs56Ugb800faunaqaq8/NoGRk9K5Z2Scz32+itXruS1117jjTfe4Pbbb++x1zXGdJ1XVwN1d9Hc1GraOo/HJ3A1DbW1bh0Q18b0OztVqekTVeVH+LSijLlLlne8sJdc4RJL4riJ/Pb5P0DyQkSE2+aN7rFtGGM6z84ENp9xYNeHiIQwJf3KHn/tWVfdQHX5ESqOFPX4axtjOs8CwHxGYc5ORk+c3uW+/5cybd4i+oWFs3frX3v8tY0xnWcBYC46/WkVVeWHSUmb1yuvHxkdw+TZC9if+R6N9fW9sg1jjPcsAMxFRXszAUiZ2WaP3B4xc+ESzp+tpdC9LWOMcywAzEVFe3YxeFgicYlJHS/cRWOnphETO5gDmR/02jaMMd6xADAANFw4z5GDe0iZOa9HTv5qT0hIKFPmLKQod5eNGGaMwywADABHC/bR2FDPhBlzen1b0+ZeTWP9Bd5///1e35Yxpn0WAAaAwwdyCO0XxuiJ03t9W0kp0xgweChvvfVWr2/LGNM+CwADwOG8HJJSphIWEdnr25KQEKbOuZIPP/zQLg1hjIMsAAwnT56ksqyEsVPS+mybU+deRUNDA1u2bOmzbRpjPssCwPDRRx8BMH7aZX22zZHJU0hMTLRmIGMcZAFg2LlzJxFR0SSOTemzbYoIS5cuZceOHZw+fbrPtmuM+QcLAENWVhajJ6YSEhrap9tdunSpNQMZ4yALgCBXXV3N0aNHGTM5tc+3PW3aNBITE3n33Xf7fNvGGAuAoLd7924ARk/s+wAQEa655hp27NjB+fPn+3z7xgQ7C4Agl52dTVRUFAljJjiy/WuvvZZz586xc+dOR7ZvTDDzakAYE7h2795NWloaof36/quQkVlKEwmER0bzzMsbqIhOBrCBYozpI7YHEMRqamooKChg9uzZjtUQ2i+MCTPSKdyTiTY3O1aHMcHIAiCI5eTkoKqOBgDAxFlXcLbmFMdKChytw5hg41UAiMgNIlIgIsUi8kAb8yeLyE4RuSAi328174iI7BORPSKS7TF9iIi8IyJF7vvB3X87pjP27t1LSEgIqal9fwDY04QZc5CQEApy7DiAMX2pwwAQkVDgN8BSYCqwUkSmtlrsJPAt4NF2XuYaVU1T1XSPaQ8AW1Q1Bdjifm760J49e5g0aRIxMTGO1hEVM4Axk1IptAAwpk95swcwFyhW1RJVrQdeBpZ7LqCqVaqaBTR0YtvLgRfcj18AburEuqabmpqayM3NZebMmU6XAriagT45XsrJymNOl2JM0PAmAEYCZR7Py93TvKXAX0Vkt4is9pg+XFUrANz3w9paWURWi0i2iGRXV1d3YrPmUg4dOsTZs2dJS0tzuhTAFQAAhTkfOVyJMcHDmwBoa3go7cQ2FqjqZbiakL4hIld1Yl1U9WlVTVfV9Pj4+M6sai5hz549AD4TAIPjExg6YjTF+7KcLsWYoOFNAJQDnoPEjgKOe7sBVT3uvq8CXsfVpARQKSKJAO77Km9f03Tfnj17GDx4MKNH+06f+wmp6ZQW7LehIo3pI94EQBaQIiLjRCQcWAFs8ObFRSRGRAa0PAauB/a7Z28AVrkfrwLe6EzhpmsyMkvJyCxla+Zuho5O4aVdZWRkljpdFgDJM+bQ1NjArl27nC7FmKDQYQCoaiOwFtgMHAReUdU8EVkjImsARCRBRMqB7wE/FpFyEYkFhgPbRGQvsAvYqKpvu1/6YWCJiBQBS9zPTR+oP3+OTyrKSBw3yelSPmP0xOmERUTy4YcfOl2KMUHBq/P/VXUTsKnVtHUej0/gahpqrQZos5uJqn4KXOd1pabHnCg9BKqM6MPr/3ujX1g4Y6eksXXrVlQVkbYOPxljeoqdCRyEKg4XAvTpADDemjAjnWPHjnH48GGnSzEm4FkABKHjR4oYMHgo/QcNcbqUz5mQOgfAmoGM6QMWAEGo4nChT/73DzAoPoHx48ezdetWp0sxJuBZAASZC+fO8mnlMUaMm+h0Ke268sorycrKoq6uzulSjAloFgBB5sRR1wFgX90DAFcANDQ0kJmZ6XQpxgQ0C4Agc9yHDwC3mDNnDlFRUWzbts3pUowJaBYAQabiSBED44YREzvI6VLa9aecE4yaOIO33nmPFz866jMnqhkTaCwAgkzFkSKf/u+/xYQZczhVXWFXBzWmF1kABJGamhpOVh7ziwBInuHqDlqcaxeHM6a3WAAEkQMHDgCQ6MM9gFoMjk8gLmEUhywAjOk1FgBBZP9+13X4/GEPAFzNQEfyc2m4cN7pUowJSBYAQSQvL49BQ4cT3T/W6VK8kpyaTlNjA0cL9jldijEByQIgiBw4cIAEP/nvH2D0pFT6hYVzaF+206UYE5AsAILE2bNnKSsrY3jSeKdL8VpYeASjJ6VyaP9up0sxJiBZAASJ4uJiVJVho8Y5XUqnJE+fzacVZRw/7vUgdMYYL1kABInCQtcZwMOT/CsAxqfOBrCzgo3pBRYAQaKwsJDo6GgGDR3udCmdEj9iDAMGD2X79u1Ol2JMwLEACBIFBQWkpKQgIf71JxcRklPT2blzJ42NjU6XY0xA8a9fA9MlqkphYSGTJvnWGMDeSp4+mzNnzrBvn3UHNaYneRUAInKDiBSISLGIPNDG/MkislNELojI9z2mJ4nIeyJyUETyROTbHvMeFJFjIrLHfVvWM2/JtFZVVcXp06eZONH3zwBuy7hpswgJCbHjAMb0sA4DQERCgd8AS4GpwEoRmdpqsZPAt4BHW01vBO5T1SnA5cA3Wq37uKqmuW+bML2ioKAAwG/3AKJiBpCammoBYEwP82YPYC5QrKolqloPvAws91xAVatUNQtoaDW9QlU/dj8+AxwERvZI5cZrLT2AUlL85ySw1hYuXMi+ffs4deqU06UYEzC8CYCRQJnH83K68CMuImOBWYDnME9rRSRXRNaLyOB21lstItkikl1dXd3ZzRpcAZCYmMjAgQOdLqXLFixYgKqyc+dOp0sxJmB4EwDSxjTtzEZEpD/wZ+A7qlrjnvwUkAykARXAr9paV1WfVtV0VU2Pj4/vzGaDXkZmKRmZpWTm7Kf/sCS/HlglNTWV2NhY6w5qTA/yJgDKgSSP56MAr0/LFJEwXD/+L6rqay3TVbVSVZtUtRl4BldTk+lhTY0NfFJRSryfnQHcWr9+/Zg/fz7btm1DtVP/fxhj2uFNAGQBKSIyTkTCgRXABm9eXEQEeA44qKqPtZqX6PH0ZmC/dyWbzvikopzmpia/OwO4LQsWLKCqqoqioiKnSzEmIPTraAFVbRSRtcBmIBRYr6p5IrLGPX+diCQA2UAs0Cwi38HVY2gGcCewT0T2uF/y3909fh4RkTRczUlHgK/14PsyblVlhwEYFgABsHDhQgC2b9/ut11ajfElHQYAgPsHe1Oraes8Hp/A1TTU2jbaPoaAqt7pfZmmqyrLSwgJ7Ufc8Lb+PP4lISGB5ORktm7dyt133+10Ocb4PTsTOMBVlx8hfsRoQvt5lfU+b+HChezevZtz5845XYoxfi8wfhVMuyrLDjN2ykyny+i2lh5MDUMnUl9fzy//sIkJM+Zw27zRDldmjP+yPYAAdq62hjOnPvG7MQAuxUYJM6bnWAAEsKryI4D/jQFwKWHhEYyxUcKM6REWAAGsMoB6AHkan5rOpxVlnP60yulSjPFrFgABrKqshKj+sfQfOMTpUnpU8nTXKGHWDGRM91gABLCq8iMMTxqH63y8wDF0xGhihwy1ZiBjuskCIEA1NzdTVX44oA4AtxARxk9P53Bejo0SZkw3WAAEqLKyMhrqLwRkAICrGejCubPk5uY6XYoxfssCIEC1DAITSD2API2bNgsRGyXMmO6wAAhQRUVFIEL8yDFOl9IromIGMDJ5kl0e2phusAAIUAUFBQwZPpKwiEinS+k146en2yhhxnSDBUCAKigoYNiosU6X0auSp8+2UcKM6QYLgABUV1dHWVlZwLb/txgxfiIDBw60ZiBjusgCIAAVFxejqgwbNd7pUnpVSEgoV1xxhY0SZkwXWQAEoJYeQMOSxjpbSB+wUcKM6ToLgABUVFREdHQ0g4cmOF1Kr2sZJcy6gxrTeRYAAaigoICUlBQkJPD/vAkJCUyYMMECwJguCPxfiCCjqhQUFATVmLk2SpgxXeNVAIjIDSJSICLFIvJAG/Mni8hOEbkgIt/3Zl0RGSIi74hIkft+cPffjqmurub06dNMmjTJ6VL6zMKFC6mvrycrK8vpUozxKx0GgIiEAr8BlgJTgZUiMrXVYieBbwGPdmLdB4AtqpoCbHE/N93UcgA4mPYAZs+eTUREhDUDGdNJ3uwBzAWKVbVEVeuBl4HlnguoapWqZgENnVh3OfCC+/ELwE1dewvGU7AFQEZmKa/trWLUxFQ2/vVdMjJLL44fbIy5NG8CYCRQ5vG83D3NG5dad7iqVgC474e19QIislpEskUku7q62svNBq/CwkISEhIYOHCg06X0qZSZ8zhZeYxPT5Q7XYoxfsObAGhrNBFvz7rpzrquhVWfVtV0VU2Pj4/vzKpBqbCwMGj++/eUMnMuAEV7dzlciTH+w5sAKAeSPJ6PAo57+fqXWrdSRBIB3Pc2wGs31dfXU1JSEpQBMCg+gfhRYyna85HTpRjjN7wJgCwgRUTGiUg4sALY4OXrX2rdDcAq9+NVwBvel23acuTIERoaGoKqB5CnlJnzKC3cz/mztU6XYoxf6DAAVLURWAtsBg4Cr6hqnoisEZE1ACKSICLlwPeAH4tIuYjEtreu+6UfBpaISBGwxP3cdENhYSEQPAeAW5uYNo/mpiYO5dlYwcZ4o583C6nqJmBTq2nrPB6fwNW849W67umfAtd1plhzaQUFBYSFhTFuXGBfBbQ9I5MnExUzgKI9mcCdTpdjjM+zM4EDSEFBAePHjycsLMzpUhwREhLKhBlzKc7NoqmpyelyjPF5FgABpKioKGibf1qkpM3lXG2NDRZvjBe8agIyvi0js5RztTWcOHGCuqjhQX0iVPL0dEJCQ3n//feZNWuW0+UY49NsDyBAVJUfAQj4YSA7EhnTn9ETp/P+++87XYoxPs8CIEBUlR8GYHhSYI8C5o2UmfMoLCzk2LFjTpdijE+zAAgQlWWHieofS/9BQ5wuxXEpafMA+OCDDxyuxBjfZgEQIKrKDjNs1FhE2rr6RnCJSxjFmDFjrBnImA5YAAQAbW6m6tgRa/7xsGjRIjIzM6mrq3O6FGN8lgVAADj1yQkaLpxn2KjgPAGsLYsWLaK+vp4dO3Y4XYoxPssCIABUlbkOAA9LGutsIT5k9uzZxMbGsmXLFqdLMcZnWQAEgMqywyBC/MixTpfiM8LCwli0aBHvvfceDQ2txykyxoAFQECoKj/MkGEjCI+IdLoUn7J48WJOnz5Ndna206UY45MsAAJASw8g8w8ZmaUcjxhDv/AI/t+Lrwf12dHGtMcCwM/V1dVxsuo4w6wH0OeERUSSPH02BR/vRJubnS7HGJ9jAeDnDh06BKoMT7IeQG2ZPHsBZ059wvEjhU6XYozPsQDwcwUFBYBdA6g9KWnzCAkNJT97u9OlGONzLAD8XEFBAWERkQyOT3S6FJ8UFTOAMZNnkP/xDlTV6XKM8SkWAH4uPz+fYaPGISH2p2zP5MsWcPJEuau5zBhzkVe/GiJyg4gUiEixiDzQxnwRkSfc83NF5DL39EkissfjViMi33HPe1BEjnnMW9aj7ywIqCoFBQXW/t+BSZfNBxE2b97sdCnG+JQOA0BEQoHfAEuBqcBKEZnaarGlQIr7thp4CkBVC1Q1TVXTgNlAHfC6x3qPt8x3jx1sOuH48eOcOXOGhNHJTpfi0wYMjmN0yjTefvttp0sxxqd4swcwFyhW1RJVrQdeBpa3WmY58Ht1+QgYJCKtG6WvAw6p6tFuV22AfxwAHj7auoB2ZOrcqyguLqaoqMjpUozxGd4EwEigzON5uXtaZ5dZAbzUatpad5PRehEZ7EUtxkN+fj4iYheB88KU9CsJCQmxvQBjPHgTAG1dYL51d4pLLiMi4cCXgFc95j8FJANpQAXwqzY3LrJaRLJFJLu6utqLcoNHfn4+SUlJhEdGOV2Kz+s/aAhz5szhrbfest5Axrh5EwDlQJLH81HA8U4usxT4WFUrWyaoaqWqNqlqM/AMrqamz1HVp1U1XVXT4+PjvSg3eBQUFDB58mSny/AbS5cu5fDhwxQW2klhxoB3AZAFpIjIOPd/8iuADa2W2QDc5e4NdDlwWlUrPOavpFXzT6tjBDcD+ztdfRA7e/YspaWlFgCdsGTJEkJDQ3nrrbecLsUYn9BhAKhqI7AW2AwcBF5R1TwRWSMia9yLbQJKgGJc/83f27K+iEQDS4DXWr30IyKyT0RygWuA73b3zQSTlv9iLQC8N2TIEObNm2fNQMa49fNmIXcXzU2tpq3zeKzAN9pZtw6Ia2P6nZ2q1HxGfn4+AJMmTaKirNHhavzH0qVL+clPfsKBAweYNm2a0+UY4yg7fdRP5efnExsbS2KiXQKiMxYvXkxYWBgbN250uhRjHGcB4Kfy8/OZPHkyIm11wDLtGTRoEFdddRV/+ctfaGpqcrocYxxlAeCHmpqaKCoqYtKkSU6X4lcyMkvJyCxl8OQrqK6u5qHnN9hAMSaoWQD4odLSUs6dO2cHgLsoJW0ekdH92bfdBow3wc0CwA+1XALC9gC6pl9YOFPnXkX+7m3Unz/ndDnGOMYCwA/l5+cTGhrKhAkTnC7Fb6XOv46G+gvk77aBYkzwsgDwQwcOHCA5OZmIiAinS/FbSSnTGBSfwL4d1gxkgpcFgJ9RVfLy8pg+fbrTpfg1ESH1ims5fGAPVVVVTpdjjCMsAPxIRmYp697azcmTJ6mLTrzYq8V0Ter8xag28+abbzpdijGOsADwMxVHigFIHJvicCX+Ly5hJEkp0/jzn/9sl4YwQckCwM9UHClCQkJsEJgeMutq1xVCs7OznS7FmD5nAeBnKo4UEj9iDGHhdgC4J0ydcyUDBgzg1Vdf7XhhYwKMBYAfUVUqjhSROM6af3pKWEQkX/ziF9m8eTN///vfnS7HmD5lAeBHak5WU3fmNIljLAB60i233EJ9fb0dDDZBxwLAj9gB4N4xZcoUpk+fzquvvmoHg01QsQDwI3YAuPfccsstFBUVsXfvXqdLMabPWAD4kRNHi4gfaQeAe8ONN95IdHQ0f/zjH50uxZg+YwHgJ1SV44cLrfmnl8TExLB8+XI2btzIp59+6nQ5xvQJCwA/ceLECTsA3EtazqgelHotDQ0N/PixZ+wMaxMUvAoAEblBRApEpFhEHmhjvojIE+75uSJymce8I+7B3/eISLbH9CEi8o6IFLnvB/fMWwpMeXl5gB0A7k1DR4wmefpsst99k6ZGG2fZBL4OA0BEQoHfAEuBqcBKEZnaarGlQIr7thp4qtX8a1Q1TVXTPaY9AGxR1RRgi/u5aUdeXp4dAO4Dc5bcRO3fT3Iwe6vTpRjT67zZA5gLFKtqiarWAy8Dy1stsxz4vbp8BAwSkY5GK18OvOB+/AJwk/dlB599+/YxbORYOwDcyyakpjNk+Eh2vfOG06UY0+u8CYCRQJnH83L3NG+XUeCvIrJbRFZ7LDNcVSsA3PfD2tq4iKwWkWwRya6urvai3MDT3NzM3r17GZXSesfL9DQJCWHO4uUcO3SQ3Nxcp8sxpld5EwDSxrTWZ8tcapkFqnoZrmaib4jIVZ2oD1V9WlXTVTU9Pj6+M6sGjOLiYmpraxmVPMXpUoLCzIWLCY+M5vnnn3e6FGN6lTcBUA4keTwfBRz3dhlVbbmvAl7H1aQEUNnSTOS+t1E52rFnzx4A2wPoIxFRMcy+5kY2b97M0aNHnS7HmF7jTQBkASkiMk5EwoEVwIZWy2wA7nL3BrocOK2qFSISIyIDAEQkBrge2O+xzir341WANbq2Iycnh7i4OAbHd3RYxfSUeV+4mX79+rF+/XqnSzGm13QYAKraCKwFNgMHgVdUNU9E1ojIGvdim4ASoBh4BrjXPX04sE1E9gK7gI2q+rZ73sPAEhEpApa4n5s25OTkkJaWhkhbLW2mNwwYFMfNN9/M66+/bkNGmoAl/nTxq/T0dA22gTtOnTrF/Pnzue+++4hOvd7pcoLKlSOEG264gVWrVvHDH/7Q6XKM6TIR2d2qGz5gZwL7vJb2/1mzZjlbSBBKSkrixhtv5I9//COnTp1yuhxjepwFgI/LyckhLCyMadOmOV1KUPrqV79KXV0dL774otOlGNPj+jldgLm0nJwcpkyZQmRkpNOlBB3X9YAimXTZfJ5d/zsipywiun8st80b7XRpxvQI2wPwURmZpfzP9hL25O4jKiHZLk7moEX/vIoL58+xY9MrTpdiTI+yAPBhlWUlNNZfYNQE6//vpGGjxpJ6+bVk/W0DZ07ZpaJN4LAA8GHlxQcAGDXBzgB22lU330FzUyNb38xwuhRjeowFgA8rKzpA7JB4YocE5yUwfMmQYSOYdfVScj54i7Kyso5XMMYPWAD4KFXlaH4uoyelOl2Kcbvyn24jJLQfTzzxhNOlGNMjLAB81CcVZZytOcXYyTOdLsW4DRgcx7zrb+Yvf/kLOTk5TpdjTLdZAPioowf3AjB2igWAL1n4xRUMGzaMX/ziFzQ3NztdjjHdYgHgo44c3MvAuGEMik9wuhTjITwyih/84Afk5eXx2muvOV2OMd1iJ4L5oObmZo4W5DJhxly7AJwPOj00laSUafzfRx6lZshUImP628lhxi/ZHoAPKioqou7MaWv+8VEiwhfuuJe62ho++N//cbocY7rMAsAHbd++HYBxU9OcLcS0K3HMBGYvWkbW3zZwvKTA6XKM6RILAB+0bds2ho4Ybf3/fdy1t9xD/0FD2PDcr6ivr3e6HGM6zQLAx5w7d47du3eTnPq5S3cbHxMZHcONq75F9bGj/Pa3v3W6HGM6zQLAx2RlZVFfX0/y9NlOl2K8kJI2j9QrruPpp5+moMCagox/sQDwMdu2bSMiIsLOAPYj19++hoEDB/Jv//Zv1hRk/IoFgA9RVT788EPmzJlDWHiE0+UYL0X3j+WnP/0pBw8e5LHHHnO6HGO85lUAiMgNIlIgIsUi8kAb80VEnnDPzxWRy9zTk0TkPRE5KCJ5IvJtj3UeFJFjIrLHfVvWc2/LP5WUlHD06FGuueYap0sxnXTddddx++2388ILL/Dee+85XY4xXukwAEQkFPgNsBSYCqwUkdYXqF8KpLhvq4Gn3NMbgftUdQpwOfCNVus+rqpp7tum7r0V//fuu+8CcO211zpciemsjMxSRl+9guGjk7nvhw+w7u2PbRAf4/O82QOYCxSraomq1gMvA8tbLbMc+L26fAQMEpFEVa1Q1Y8BVPUMcBAY2YP1B5QtW7Ywbdo0EhLs8g/+qF94OP/89X+jsaGe19c9TFNjo9MlGXNJ3gTASMDzAujlfP5HvMNlRGQsMAvI9Ji81t1ktF5EBre1cRFZLSLZIpJdXV3tRbn+qbq6mtzcXPvv388NTUzixlXfprRgH5sznup4BWMc5E0AtHUxGu3MMiLSH/gz8B1VrXFPfgpIBtKACuBXbW1cVZ9W1XRVTY+PD9wTo7Zs2YKqsnjxYqdLMd2UOv9a5i+7ld3v/oWMDBtBzPgubwKgHEjyeD4KOO7tMiIShuvH/0VVvXj5RFWtVNUmVW0GnsHV1BS0Nm3aRHJyMikpKU6XYnrANV/+Cilpl/PQQw+xY8cOp8sxpk3eBEAWkCIi40QkHFgBbGi1zAbgLndvoMuB06paIa5LWT4HHFTVz/SPE5FEj6c3A/u7/C78XGVlJdnZ2Sxbtsyu/hkgQkJCuflr9zN+/Hi++93vkp+f73RJxnxOhwGgqo3AWmAzroO4r6hqnoisEZE17sU2ASVAMa7/5u91T18A3Alc20Z3z0dEZJ+I5ALXAN/tsXflRzIyS/nFb19CVWkamUZGZqn1HgkQEVHRrFu3jpiYGL761a9y+PBhp0sy5jO8Gg/A3UVzU6tp6zweK/CNNtbbRtvHB1DVOztVaQDb/9F7JIyZQFzCKKdLMT1sxIgRPPfcc9xxxx3cc889/OEPf2DEiBFOl2UMYGcCO66yrISKw4XMXLjE6VJML8jILGVnVSj/8u2fc/LvNdxy2108tSnb6bKMASwAHJfzwduE9gsj9Qrr/hnIEsYkc9t9P+dszSleeOg+aw4yPsECwEHnz59n344tTE5fQFT/WKfLMb1s1ISp3PnAL2mov8Cdd97JwYMHnS7JBDkLAAe9/fbbnK+rZdZVNzhdiukjiWMmsOrff0V4eDh33XUX27Ztc7okE8QsAByiqrzwwgsMHTGasVPSnC7H9KGhiUnc8v1HiB4cz+rVX+NbP32cFz86ar2/TJ+zAHDIRx99RH5+Ppd/4V+s738QGhg3jLt//GsmzZ7POy8/zYZnH6XhwnmnyzJBxgLAIb/73e+Ii4uzg79BLDwiki/f+yOuvvlOcnds4Zn//AYHDhxwuiwTRCwAHJCTk8PWrVu566676Bce7nQ5xkESEsJVy+/gjh88TP35c6xYsYL169fT1NTkdGkmCFgA9DFV5de//jVxcXHcfvvtTpdjfMS4qWms/vk6rr76an75y19y22232d6A6XUWAH1sx44d7Nq1izVr1hATE+N0OcaHRPeP5fLbv89NX7ufQ0dK+fKXb+Er3/kRNTU1Ha9sTBeI6yoO/iE9PV2zs/3zLMqMzFIa6+v57U/WoNrMml88Tb8wa/4xbTt39gzv/el37H5/EwNjY1mzZg0rV64kIsLGijadJyK7VTW99XTbA+hDO956hZOVx1h61zftx99cUlTMAJat+hZfffBJpk+fzn/913+xdOlSXnrpJc6ft95CpmdYAPSRiiNFbN3wEtPmLSJ5+mynyzF+InHMBJ599lnWr19PfHw8P/vZz1i8eDHr1q3j1KlTTpdn/Jw1AfWB2tparr/xSzQ1NvB/fvYU0XbZB9MFqkppwT52bHqF4twswsLCuP7667n11luZM2eOnU9i2tVeE5BXl4M2XdfQ0MD3vvc9/v5JJXc98Ev78TddJiKMmTyDMZNnUHXsCB+/t4m/vfs3Nm7cyMC44UydcyVT5l7J91cssTAwXrE9gF7U1NTEj370I9544w1u/Mq3uWzRso5XMqYTGi6c52D2Ng7s+oBD+z+muamRESNGsGTJEhYuXMjs2bOJiopyukzjsPb2ACwAesmFCxf4wQ9+wDvvvMM3v/lNBs3+otMlmQB3/mwtBTk7OV2UxY4dO2hoaCAsLIxZs2ZxxRVXMHv2bKZOnWrdj4OQBUAfycgspfp4Ka+ve5jK0kNcv/JrzPvCPztdlgkyDRfOU1q4n8N5OZQcyKGy9BAAISEhTJgwgdTUVFJTU0lJSSE5OZmBAwc6XLHpTRYAfaCmpobv/uwxMje/RnhkFF+65z4mzrrc6bKMoe7MaY4dyufY4QKOlxRw/HAh52r/cYJZfHw8ycnJjBs3jpEjR37mNnjwYDum4Oe6dRBYRG4A/hsIBZ5V1YdbzRf3/GVAHfAVVf34UuuKyBDgj8BY4Ahwq6r6Xb+2xsZGPv74YzZt2sSbb75JXV0dqfOv47pb72HAoDinyzMGgOgBA0lJm0dK2jzA1aPo9CeVVB8vpfrYUT45fpTSY0fZs28/58/WfmbdsPAIYuOG0T92MFPGjWDo0KHExcVdvMXGxhIbG0v//v2JjY0lJibGAsNPdLgHICKhQCGwBCgHsoCVqnrAY5llwDdxBcA84L9Vdd6l1hWRR4CTqvqwiDwADFbV+y9VS2/vAagqzc3NF+/BNWpXbW3txdvp06cpLS2ltLSUkpIScnNzqaurIyIigmXLljE0bQkJY5J7rUZjetv5urOc/rSSv39SyelPXPc1n1ZztuYUtadPcbbmFBfO1bW7fkhICAMGDGDAgAFER0cTGRlJVFQUkZGRn7lFRUURERFx8T4sLIywsDD69et3yZvnMqGhoYSEhBASEoKIXHzs+VxEPjOv9eO25nkGWFuP/S3gurMHMBcoVtUS9wu9DCwHPK9UtRz4vbrS5CMRGSQiibj+u29v3eXAIvf6LwDvA5cMgK56+OGHefXVVz/3497yvOXWGRFR0QxJGMXUy69j7JSZJM+YQ3hEZG+Ub0yfioyOITJ6PMOTxre7TGN9PbU1p6g7c5rzdbWcrzvLhbqzHo9rOVdXS8OF89TWX+DUqVoa60/SUH+BhgsXaGy4cPGxanMfvrve1VZAdCZALrX+k08+yYIFC3q0Xm8CYCRQ5vG8HNd/+R0tM7KDdYeragWAqlaIyLC2Ni4iq4HV7qe1IlLgRc195GNvFhoKfNLLhfgD+xxc7HOwz6BFpz6HhQsXdmdbY9qa6E0AtBVVrf9dbm8Zb9a9JFV9Gni6M+v4EhHJbmvXK9jY5+Bin4N9Bi184XPw5lpA5UCSx/NRwHEvl7nUupXuZiLc91Xel22MMaa7vAmALCBFRMaJSDiwAtjQapkNwF3icjlw2t28c6l1NwCr3I9XAW90870YY4zphA6bgFS1UUTWAptxdeVcr6p5IrLGPX8dsAlXD6BiXN1A777Uuu6Xfhh4RUTuAUqBW3r0nfkOv22+6mH2ObjY52CfQQvHPwe/OhHMGGNMz7HxAIwxJkhZABhjTJCyAOhFInJERPaJyB4R8d2LGPUwEVkvIlUist9j2hAReUdEitz3g52ssbe18xk8KCLH3N+HPe4z6AOaiCSJyHsiclBE8kTk2+7pwfZ9aO9zcPQ7YccAepGIHAHSVTWoTnoRkauAWlxnh093T+v0pT/8WTufwYNArao+6mRtfcndxTtRVT8WkQHAbuAm4CsE1/ehvc/hVhz8TtgegOlxqvohcLLV5OW4LvmB+/6mvqypr7XzGQQdVa1ouTCkqp4BDuK6QkCwfR/a+xwcZQHQuxT4q4jsdl/SIph95tIfQJuX/ggCa0Uk191EFNDNHq2JyFhgFpBJEH8fWn0O4OB3wgKgdy1Q1cuApcA33M0CJng9BSQDaUAF8CtHq+lDItIf+DPwHVWt6Wj5QNXG5+Dod8ICoBep6nH3fRXwOq4rqwaroL/0h6pWqmqTui5/+QxB8n0QkTBcP3ovqupr7slB931o63Nw+jthAdBLRCTGfbAHEYkBrgf2X3qtgBb0l/5o+cFzu5kg+D6I61rGzwEHVfUxj1lB9X1o73Nw+jthvYB6iYiMx/VfP7guuZGhqr9wsKQ+IyIv4RrrYShQCfwn8L/AK8Bo3Jf+UNWAPUjazmewCNeuvuIaBe9rLe3ggUpEFgJbgX1Ay4X//x1X+3cwfR/a+xxW4uB3wgLAGGOClDUBGWNMkLIAMMaYIGUBYIwxQcoCwBhjgpQFgDHGBCkLAGOMCVIWAMYYE6T+P6GHbycUgwiqAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "samples = df.iloc[:,0]*dummy_mean[0]\n",
    "sns.distplot(x=samples, fit=lognorm, kde=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 460,
   "id": "fe416d00",
   "metadata": {},
   "outputs": [],
   "source": [
    "dd = pd.concat([df.iloc[:,0]*dummy_mean[0], df.iloc[:,1]*dummy_mean[1]], ignore_index=True, axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 461,
   "id": "e35db894",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>5.27697e+00</td>\n",
       "      <td>-3.78839e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>-3.78839e+00</td>\n",
       "      <td>3.19615e+00</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             0            1\n",
       "0  5.27697e+00 -3.78839e+00\n",
       "1 -3.78839e+00  3.19615e+00"
      ]
     },
     "execution_count": 461,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dd.cov()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c039c323",
   "metadata": {},
   "source": [
    "# Test case: Eu155"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d43ab36b",
   "metadata": {},
   "source": [
    "$$ cov(x_i, x_j) = (e^{cov(ln(x_i),ln(x_j))}-1)\\cdot<x_i><x_j> $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 462,
   "id": "5595bbc3",
   "metadata": {},
   "outputs": [],
   "source": [
    "endf6 = sandy.get_endf6_file('jeff_33','xs', 631550)\n",
    "# endf6 = sandy.get_endf6_file('jeff_33','xs', 922350)  # This was tested. It gives no error, but its smaller std gives less evident results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 463,
   "id": "cdf509cf",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:  found duplicate energy for MAT6337/MF3/MT1 at 1.00000e+05 MeV, keep only first value\n",
      "WARNING:  found duplicate energy for MAT6337/MF3/MT2 at 1.00000e+05 MeV, keep only first value\n",
      "WARNING:  found duplicate energy for MAT6337/MF3/MT102 at 1.00000e+05 MeV, keep only first value\n"
     ]
    }
   ],
   "source": [
    "Cov = sandy.XsCov.from_endf6(endf6)\n",
    "BE = sandy.Xs.from_endf6(endf6)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9bc8011",
   "metadata": {},
   "source": [
    "## Normal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 464,
   "id": "3ab4d12e",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\fgrimald\\Anaconda3\\envs\\sandy-devel\\lib\\site-packages\\ipykernel_launcher.py:1: RuntimeWarning: covariance is not positive-semidefinite.\n",
      "  \"\"\"Entry point for launching an IPython kernel.\n"
     ]
    }
   ],
   "source": [
    "df = pd.DataFrame(np.random.multivariate_normal(np.ones(len(Cov)), Cov, size=1000000))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 465,
   "id": "8f2828d4",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\fgrimald\\Anaconda3\\envs\\sandy-devel\\lib\\site-packages\\seaborn\\distributions.py:2619: FutureWarning: `distplot` is a deprecated function and will be removed in a future version. Please adapt your code to use either `displot` (a figure-level function with similar flexibility) or `histplot` (an axes-level function for histograms).\n",
      "  warnings.warn(msg, FutureWarning)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "mean: [0.9636815054506004, 0.9830082955835906]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAD8CAYAAABekO4JAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAAZ+0lEQVR4nO3deXQVZb7u8e+PkAAJCUFGCYOkERDFMU6INjdXo5fjQSaRSVTEXHAe0NZuGzwQ9doKDQ4rMmoTAg6A2EbwoLQubVtDB0Vl8NgqKBCRQCIJApKQ9/6RSKMGspHs1Jvs57PWXtlJFXs/u9bLsypV9VbMOYeIiPirQdABRETkyFTUIiKeU1GLiHhORS0i4jkVtYiI51TUIiKeC6mozSzRzBaZ2admtsHMzg93MBERqdAwxPWmA6855wabWQwQG8ZMIiJyCKtuwouZJQAfAclOs2NERGpdKHvUyUAB8IyZnQasBm5zzn1/6Epmlg6kA8TFxZ3VvXv3ms4qIlJvrV69eodzrlVVy0LZo04B3gcucM7lmtl0oNg598fD/ZuUlBSXl5d3LJlFRCKKma12zqVUtSyUk4lbgC3OudzK7xcBZ9ZUOBERObJqi9o5tw3YbGbdKn/0v4H1YU0lIiIHhXrVxy1AduUVH18C14UvkoiIHCqkonbOrQGqPHYiIiLhpZmJIiKeU1GLiHhORS0i4jkVtYiI51TUIiKeC/XyPJGIsyD365DXHX5uxzAmkUinPWoREc+pqEVEPKeiFhHxnIpaRMRzKmoREc+pqEVEPKeiFhHxnIpaRMRzKmoREc+pqEVEPKeiFhHxnIpaRMRzKmoREc+pqEVEPKeiFhHxnIpaRMRzKmoREc+pqEVEPKeiFhHxnIpaRMRzKmoREc+F9FfIzWwTUAIcAMqccynhDCUiIv8WUlFX+l/OuR1hSyIiIlXSoQ8REc+FWtQOWGFmq80sPZyBRETkp0I99HGBcy7fzFoDr5vZp865tw9dobLA0wE6duxYwzFFRCJXSHvUzrn8yq/bgZeAc6pYZ6ZzLsU5l9KqVauaTSkiEsGqLWozizOz+B+fA2nA2nAHExGRCqEc+mgDvGRmP66/wDn3WlhTiYjIQdUWtXPuS+C0WsgiIiJV0OV5IiKeU1GLiHhORS0i4jkVtYiI51TUIiKeU1GLiHhORS0i4jkVtYiI51TUIiKeU1GLiHhORS0i4jkVtYiI51TUIiKeU1GLiHhORS0i4jkVtYiI51TUIiKeU1GLiHhORS0i4jkVtYiI51TUIiKeU1GLiHhORS0i4jkVtYiI51TUIiKeU1GLiHhORS0i4jkVtYiI50IuajOLMrMPzSwnnIFEROSnjmaP+jZgQ7iCiIhI1UIqajNrD/wHMDu8cURE5OdC3aOeBtwDlB9uBTNLN7M8M8srKCioiWwiIkIIRW1mlwPbnXOrj7Sec26mcy7FOZfSqlWrGgsoIhLpQtmjvgDoZ2abgOeAVDObH9ZUIiJyULVF7Zy7zznX3jl3AjAU+JtzbmTYk4mICKDrqEVEvNfwaFZ2zr0FvBWWJCIiUiXtUYuIeE5FLSLiORW1iIjnVNQiIp5TUYuIeE5FLSLiORW1iIjnVNQiIp5TUYuIeE5FLSLiORW1iIjnVNQiIp5TUYuIeE5FLSLiORW1iIjnVNQiIp5TUYuIeE5FLSLiORW1iIjnVNQiIp5TUYuIeE5FLSLiORW1iIjnVNQiIp5TUYuIeE5FLSLiORW1iIjnqi1qM2tsZqvM7CMzW2dm/1UbwUREpELDENb5AUh1zu02s2jg72a23Dn3fpiziYgIIRS1c84Buyu/ja58uHCGEhGRfwvpGLWZRZnZGmA78LpzLreKddLNLM/M8goKCmo4pohI5AqpqJ1zB5xzpwPtgXPM7JQq1pnpnEtxzqW0atWqhmOKiESuo7rqwzn3HfAWcFk4woiIyC+FctVHKzNLrHzeBLgY+DTMuUREpFIoV30cD/zFzKKoKPYXnHM54Y0lIiI/CuWqj4+BM2ohi4iIVEEzE0VEPKeiFhHxnIpaRMRzKmoREc+pqEVEPKeiFhHxnIpaRMRzKmoREc+pqEVEPKeiFhHxnIpaRMRzKmoREc+pqEVEPKeiFhHxnIpaRMRzKmoREc+pqEVEPKeiFhHxnIpaRMRzKmoREc+pqEVEPKeiFhHxnIpaRMRzKmoREc+pqEVEPKeiFhHxnIpaRMRz1Ra1mXUwszfNbIOZrTOz22ojmIiIVGgYwjplwF3OuQ/MLB5YbWavO+fWhzmbiIgQwh61c+4b59wHlc9LgA1AUriDiYhIhaM6Rm1mJwBnALlVLEs3szwzyysoKKiheCIiEnJRm1lTYDFwu3Ou+OfLnXMznXMpzrmUVq1a1WRGEZGIFlJRm1k0FSWd7ZxbEt5IIiJyqFCu+jBgDrDBOTc1/JFERORQoexRXwBcDaSa2ZrKR98w5xIRkUrVXp7nnPs7YLWQRUREqqCZiSIinlNRi4h4TkUtIuI5FbWIiOdU1CIinlNRi4h4TkUtIuI5FbWIiOdCuR+1SL1TXl5OQUEB+/btO/izZs2a0axZMyrumiDiDxW11Hvl5eV8+umnrFq1ilWrVrFx40a2bt1KaWnpL9Zt2rQp7du3p1u3briWXeh00mkktmwTQGqRf1NRS721YcMGXnrpJV599VUKCwsBOOGEE+jWrRupqam0b9+euLg4AJxzFBUVsWXLFjZv3szbb79NUdHLALRqfwKnXXAJPc9PpWnicYF9HolcKmqpV8rLy1mxYgWzZs1i/fr1REdHk5qaSmpqKueeey5t2oS2d1xeXs70Je+waf2HrFv1Nm88P4uVL86h2xnn07vfcI7v1CXMn0Tk31TUUi+Ul5ezfPlyMjMz+eKLL0hOTub++++nb9++NG/e/Khfr0GDBrTp0Jk2HTpz7qUD2ZH/NR/9/XU+eGsZn65+lxNPO4eL+l9Nu85dw/BpRH7KnHM1/qIpKSkuLy+vxl9XpCpr165l8uTJfPzxx5x44omMHTuWSy+9lKioqF+suyD362N6r317vuefb7xM7n8vYe/3JZx+4aWkXjmaGy459ZheV8TMVjvnUqpapj1qqbOKi4uZMmUKL774Ii1atODhhx+mX79+NGgQvqtOG8fGcWG/4ZxzSX/e+esCclcsYUPe34ndfgfDhg0L63tL5NKokjrp3Xff5YorrmDx4sWMGjWKZcuW0b9//1orykZNYrn4qjH838lP0y65KxkZGYwePZr8/PxaeX+JLCpqqVP27dtHRkYGY8aMITY2loULF3LvvfcSHx8fSJ6W7ToyYvzDTJo0iU8++YQrrriCpUuXEo5DihK5VNRSZ2zatImhQ4eSnZ3N1VdfzeLFi+nZs2fQsTAzrrzySl5++WW6d+/Offfdxx/+8Af27t0bdDSpJ1TUUiesWLGCK6+8km3btjFjxgx+//vf07hx46Bj/UT79u159tlnufHGG1m6dClDhw5l48aNQceSekBFLV4rLy9n2rRp3HbbbSQnJ7N48WIuuuiioGMdVlRUFLfccgszZ85k+/btDBkyhHfeeSfoWFLHqajFW3v37uXOO+9kxowZDB48mKysLJKSkoKOFZLevXuzePFiOnTowNixY8nOzg46ktRhKmrx0vbt2xk1ahQrVqzgnnvuYdKkScTExAQd66i0a9eOrKws+vTpQ0ZGBpMnT6asrCzoWFIHqajFOxs2bOCqq67iyy+/5Mknn+S6666rs3e0i4uL4/HHH2f06NEsWLCAsWPHUlxcHHQsqWNU1OKV3NxcRo4cCcD8+fNJTU0NONGxi4qK4u6772by5Mnk5uZy9dVXs3379qBjSR2iohZvvPHGG6Snp9OuXTuef/55TjrppKAj1ajBgwczY8YMtmzZwvDhw9m0aVPQkaSOUFGLF5YuXcptt91G9+7dycrKonXr1kFHCotevXrx7LPPsmfPHkaMGMG6deuCjiR1gIpaAjdv3jzuu+8+zj33XObOnUtiYmLQkcKqZ8+eZGdn06hRI6655hpyc3ODjiSeq/bueWY2F7gc2O6cOyWUF9Xd8yQUzjmeeOIJMjMzSUtL49FHH/1VV3Yc6x3xatPwczsefP7tt98yZswYvvrqKx577DHS0tICTCZBO9Ld80LZo34WuKxGE0nEKy8v58EHHyQzM5OBAwcyZcqUOnf53bFq06YNWVlZnHzyydxxxx288MILQUcST1Vb1M65t4HCWsgiEaK0tJR7772X7OxsRo8eTUZGBg0bRuYddxMTE5k7dy4XXHABEydOZNasWUFHEg/pGLXUqn379nHrrbfyyiuvcMcddzB+/Pg6e410TWnSpAlPPvkkffv2ZerUqTz22GO6+578RI3txphZOpAO0LFjx2rWlkhUUlLCjTfeyOrVq5k4cSJDhw4NOpI3YmJi+NOf/kR8fDxz5syhuLiYiRMnVvlXaiTy1FhRO+dmAjOh4mRiTb2u1A+FhYXccMMNfPbZZzz22GP07ds36EjeiYqKYuLEiSQmJjJjxgxKSkp45JFHIu7YvfxSZB4YlFqVn5/PmDFj+Oabb3jqqae8vvtd0MyM22+/nYSEBB599FF2797N9OnTiY2NDTqaBKjaY9RmthB4D+hmZlvM7Prwx5L6YuPGjYwcOZIdO3Ywe/ZslXSIRo8ezeTJk/nHP/7BmDFj2LVrV9CRJEChXPUxzDl3vHMu2jnX3jk3pzaCSd23bt06RowYwf79+5k3bx5nnXVW0JHqlMGDB/PnP/+ZtWvXcs0111BQUBB0JAmIrvqQsHjvvfcYNWoUTZo0Yf78+XTv3j3oSHVSWloamZmZbN68mZEjR7J169agI0kAqp2Z+GtoZmJkW758Ob/73e844YQTmDVrFm3atAnr+9WlmYlH49BZjGvWrGHs2LE0atSIOXPm0KVLlwCTSTgc68xEkZDNnz+fu+66i1NPPZX58+eHvaQjxemnn868efMoLy9n5MiRfPLJJ0FHklqkopYa4Zxj2rRpPPjgg6SmpjJ79mwSEhKCjlWvdO3alQULFhAfH8+1117L+++/H3QkqSUqajlmZWVl3H///cyYMYMhQ4Ywbdo07/5CeH3RoUMHsrOzSUpKIj09nZycnKAjSS1QUcsx+f7777nllltYsmQJ48aN44EHHojY+3bUltatWzNv3jxOP/107r77bjIzMzXlvJ7T/yj51fLz87npppv417/+xYQJExg2bFjQkeqV6k6SXpw+gX0Np/H444+zZcsWHnjgAaKjo2spndQmFbX8Kh9//DE33XQT+/bt4+mnn6Z3795BR4o4DaNjuCL9bhJbt2XJkmzy8/OZPn26zg3UQzr0IUdt+fLljBo1isaNG7Nw4UKVdIDMjD4DRvHwww+Tl5fHVVddxRdffBF0LKlhKmoJ2YEDB5g2bRp33nknJ598Ms8//7yu5/VE//79mTt3LiUlJQwZMoQ33ngj6EhSg1TUEpLCwkLGjBnDjBkzGDhwIM888wzHHXdc0LHkEGeffTaLFi3iN7/5DbfccgvTp0+nvLw86FhSA1TUUq0PP/yQQYMG8eGHH5KRkcGDDz6oW296qm3btmRlZTFw4ECefvppxo0bR1FRUdCx5BipqOWwnHNkZWUxatQooqOjWbBgAYMGDQo6llSjUaNGZGRkMGHCBN577z369++vyTF1nIpaqrR9+3bS09N56KGHuPDCC1m0aBE9evQIOpaEyMwYNmwYzz33HHFxcYwePZqpU6dSWloadDT5FVTU8gvLli2jX79+5OXl8cc//pGnnnpKl3zVUT169GDRokUMHjyYWbNmMWLECDZt2hR0LDlKKmo5aNeuXYwfP5677rqLTp06sWTJEoYPHx7xf3y2rouNjWXSpElMnz6dr7/+mgEDBjB79mztXdchmvByjI7mFpuH3raypl73aF7zcJxz5OTk8Mgjj7Br1y5uvfVWbrjhhkCngtfXW5cGKS0tjdNOO42MjAymTJnCq6++yqRJk+jZs2fQ0aQa2qOOcBs3buT666/nnnvuISkpiRdeeIFx48bpfh31VJs2bXjiiSd4/PHHKSwsZOjQoTz88MPs3r076GhyBCrqCFVUVMRDDz1Ev379WLt2LRMmTGDBggWcdNJJQUeTWnDJJZeQk5PDkCFDmDdvHpdeeinPP/88ZWVlQUeTKqioI8yePXuYM2cOl112GdnZ2QwaNIhly5YxbNgwoqKigo4ntSg+Pp6JEyfy4osv0rlzZx544AEGDBjAypUrdTc+z+j32wixd+9ennvuOebMmcPOnTu58MILGT9+PF27dg06mgTslFNOISsri9dff50pU6Zw880306NHD26++Wb69Omjk8keUFHXczt37mTBggUsXLiQoqIiOp98Bv859g90OPFk8oog7zAn7WriJGVVdJLQT2ZGWloaqampvPLKK2RmZnLjjTfSrVs3rr32Wvr27avZqAFSUddDzjnWrFnDokWLyMnJYf/+/fTp04dO519Ox246wy+H17BhQwYMGMDll19OTk4OzzzzDPfddx9Tp07lyiuvZODAgSQlJQUdM+KoqOuRgoICXn75ZZYsWcLGjRuJjY2lf//+XHPNNSQnJ2tvVkIWHR3NgAED6N+/P++++y5/+ctfyMzMJDMzk/POO49BgwZx8cUX06hRo6CjRgQVdR23bds2Vq5cycqVK1m1ahUHDhzgzDPP5Prrr+eyyy4jLi4u6IhSh5kZvXv3pnfv3mzdupWlS5eyZMkSxo8fT3x8PL/97W+5+OKL6d27t8ZaGKmo65iy0v3kb/yMTRs+4rM17zN542cAJCcnc/3119O/f386d+4ccEqpj5KSkrjpppsYN24cubm5vPLKK7z55pvk5OQQExNDr169uOiiizjnnHNITk7WScgaZOG4DCclJcXl5eXV+Ov6KNyHE/buLmbb11+w5fMNfPXpR2z+fANl+38AICm5O93O7EW3M8+nZbvwnPyTuiGok7/lBw6w+V/rIP8TVq5cydatWwFo0aIFZ599NmeffTY9e/aka9euBw+ThGs2b11nZqudcylVLdMetSdK9/9A4bf57Ny2hZ3fbGbbV5+z7avP+W7HtxUrmNGmfWfO/O3/oVP30+jY7RRim+pGSRKsBlFRdOp+KsOvuZx7772XzZs3s2rVqoOP1157Dag4SdmlSxd69OjB7satadG2PS3atqdZy9Y0aKDr96sTUlGb2WXAdCAKmO2c+39hTVXPlJXuZ+/uEnbvKqS4aAclhTsoLiqgpGgnxYU7KNr+DbsKt8Mhv90c16Yd7ZK7cVbq5bTt1IXjO3WhiYpZPGZmdOzYkY4dOzJ48GCcc+Tn57Nu3TrWr1/P+vXreeuttygsLDz4b6IaRtO89fEktmxLfPMWJBzXkoTmrYg/riUJzVtSVBRPQkJCxE/GqraozSwKeAq4BNgC/NPM/uqcWx/ucEfinPvJ48c/OfTzn/94aKeqn//4OHDgAKWlpZSVlVFWVnbweWlpaZU//+GHH9i7dy979uzh/c++ofSHfez/YW/F13372LenhD27S9j3fQl7dhdT+sO+X+RvEBVFfGIL4pu3pP2JJ3Fa27TKvYwkWrRtT0zjJrW6PUVqmpmRlJREUlISaWlpQMX/w9lvfFL5m+OWiq/btlBcWMA3mz7j++LvfvIaT1d+TUhIIDExkcTERJo1a0bTpk2JjY09+GjSpAlNmjQ5+DwmJoaGDRsSExNDdHQ0DRs2JDo6usrnDRo0oEGDBpgZZvaT76v6+uPz2hTKHvU5wOfOuS8BzOw54Aqgxou6V69e7Nmzp9qC9U1Uw2hiGjUmpnETohs1pnFsUxKat6RNh2SaNI2nSVw8TZrG07RZCxKatyD+uJbEJSTqVz6JOGZGXEIicQmJdOx6yi+Wl5Xup+S7QkoKCygu2sFJxzXgu+++Y9euXQe/FhUVsXnz5oM7S3v27OHAgQOBfJZDCx6gZcuW/O1vf6vx9wqlqJOAzYd8vwU49+crmVk6kF757W4z+59jj/ertQR2BPj+dYG2UfXqzDYaEczbHtw+Ab2/l362t300Y6jT4RaEUtRV7eP/YrfWOTcTmBlioLAys7zDnT2VCtpG1dM2OjJtn+rV1DYK5e55W4AOh3zfHsg/1jcWEZHQhFLU/wRONLPOZhYDDAX+Gt5YIiLyo2oPfTjnyszsZuC/qbg8b65zbl3Ykx0bLw7BeE7bqHraRkem7VO9GtlGYZmZKCIiNUd/4UVExHMqahERz9Xpojazy8zsf8zsczO7t4rlZmaPVy7/2MzODCJnkELYRn3MbJeZral8TAgiZ1DMbK6ZbTeztYdZHtFjKITtE9HjB8DMOpjZm2a2wczWmdltVaxzbOPoSFOrfX5QcWLzCyAZiAE+Anr8bJ2+wHIqrgU/D8gNOreH26gPkBN01gC30UXAmcDawyyP9DFU3faJ6PFTuQ2OB86sfB4PfFbTXVSX96gPTm13zu0HfpzafqgrgHmuwvtAopkdX9tBAxTKNopozrm3gcIjrBLRYyiE7RPxnHPfOOc+qHxeAmygYkb3oY5pHNXloq5qavvPN04o69RnoX7+883sIzNbbmYn1060OiPSx1AoNH4qmdkJwBlA7s8WHdM4qsv3ow5lantI09/rsVA+/wdAJ+fcbjPrCywFTgx3sDok0sdQdTR+KplZU2AxcLtzrvjni6v4JyGPo7q8Rx3K1PZIn/5e7ed3zhU753ZXPl8GRJtZy9qL6L1IH0NHpPFTwcyiqSjpbOfckipWOaZxVJeLOpSp7X8FRlWecT0P2OWc+6a2gwao2m1kZm2t8nZfZnYOFWNiZ60n9Vekj6Ej0vipuKIDmANscM5NPcxqxzSO6uyhD3eYqe1mNrZy+dPAMirOtn4O7AGuCypvEELcRoOBcWZWBuwFhrrK09SRwMwWUnHlQksz2wJMBKJBYwhC2j4RPX4qXQBcDXxiZmsqf/Z7oCPUzDjSFHIREc/V5UMfIiIRQUUtIuI5FbWIiOdU1CIinlNRi4h4TkUtIuI5FbWIiOf+P2feKIuL0gkqAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.distplot(x=df.iloc[0,:], bins= 20, fit=norm, kde=False)\n",
    "print(f\"\"\"\n",
    "\n",
    "mean: {[df.iloc[0,:].mean(), df.iloc[-1,:].mean()]}\"\"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 466,
   "id": "49905b10",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-3.888276182763171"
      ]
     },
     "execution_count": 466,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.values.min()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a6b3728b",
   "metadata": {},
   "source": [
    "## Lognormal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 467,
   "id": "2de5cde3",
   "metadata": {},
   "outputs": [],
   "source": [
    "COV = Cov.apply(np.exp) - 1 * 1 * 1  # Target mean to 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 468,
   "id": "0b791b4e",
   "metadata": {},
   "outputs": [],
   "source": [
    "std = Cov.std(axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 469,
   "id": "db2ef620",
   "metadata": {},
   "outputs": [],
   "source": [
    "normal_equivalent_mean = np.log(np.ones(len(Cov))) - std**2/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 470,
   "id": "dac8601e",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\fgrimald\\Anaconda3\\envs\\sandy-devel\\lib\\site-packages\\ipykernel_launcher.py:1: RuntimeWarning: covariance is not positive-semidefinite.\n",
      "  \"\"\"Entry point for launching an IPython kernel.\n"
     ]
    }
   ],
   "source": [
    "df = pd.DataFrame(np.random.multivariate_normal(normal_equivalent_mean, COV, size=1000000)).apply(np.exp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 471,
   "id": "a86b3051",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\fgrimald\\Anaconda3\\envs\\sandy-devel\\lib\\site-packages\\seaborn\\distributions.py:2619: FutureWarning: `distplot` is a deprecated function and will be removed in a future version. Please adapt your code to use either `displot` (a figure-level function with similar flexibility) or `histplot` (an axes-level function for histograms).\n",
      "  warnings.warn(msg, FutureWarning)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "mean: [1.088601700903344, 1.039112103246743]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAAjVklEQVR4nO3deXQUZd728e8vIexhkSSsQhAQEERkx5VFEBiU0VFHHUBUiCLwKI+OII6OcHAU4QVBdCCyCC64jKAygg4+igIOyK6AihFZIkgCSCAQ1tzvH2mZTEhIh3RSSeX6nNMn3V13V19piouiuhZzziEiIsVfmNcBREQkNFToIiI+oUIXEfEJFbqIiE+o0EVEfKKUV28cFRXlYmNjvXp7EZFiae3atfucc9HZTfOs0GNjY1mzZo1Xby8iUiyZ2Y6cpmmTi4iIT6jQRUR8QoUuIuITKnQREZ9QoYuI+IQKXUTEJ1ToIiI+oUIXEfEJFbqIiE94dqSoSDDeWLUzqHF3tq9bwElEij6toYuI+ESuhW5mZc3sKzPbaGabzWx0NmM6mVmKmW0I3J4smLgiIpKTYDa5HAe6OOdSzSwCWG5mi51zK7OMW+ac6x36iCIiEoxcC91lXEU6NfAwInDTlaVFRIqYoLahm1m4mW0AkoAlzrlV2QzrGNgss9jMmuUwnzgzW2Nma5KTk88/tYiInCWoQnfOnXbOtQTqAO3MrHmWIeuAes65y4AXgPdymE+8c66Nc65NdHS252cXEZHzlKe9XJxzB4GlQI8szx9yzqUG7i8CIswsKkQZRUQkCMHs5RJtZlUC98sB1wHfZRlTw8wscL9dYL77Q55WRERyFMxeLjWBOWYWTkZRv+2c+6eZ3Q/gnJsG3AIMNrNTQBpwe+DLVBERKSTB7OXyNXB5Ns9Py3R/KjA1tNFERCQvdKSoiIhPqNBFRHxChS4i4hMqdBERn1Chi4j4hApdRMQnVOgiIj6hQhcR8QkVuoiIT6jQRUR8QoUuIuITKnQREZ9QoYuI+IQKXUTEJ1ToIiI+oUIXEfEJFbqIiE+o0EVEfEKFLiLiEyp0ERGfyLXQzaysmX1lZhvNbLOZjc5mjJnZFDNLMLOvzaxVwcQVEZGclApizHGgi3Mu1cwigOVmttg5tzLTmJ5Ao8CtPfD3wE8RESkkua6huwypgYcRgZvLMqwPMDcwdiVQxcxqhjaqiIicS1Db0M0s3Mw2AEnAEufcqixDagO7Mj1ODDyXdT5xZrbGzNYkJyefZ2QREclOUIXunDvtnGsJ1AHamVnzLEMsu5dlM59451wb51yb6OjoPIcVEZGc5WkvF+fcQWAp0CPLpETgwkyP6wC78xNMRETyJpi9XKLNrErgfjngOuC7LMM+APoH9nbpAKQ45/aEOqyIiOQsmL1cagJzzCycjH8A3nbO/dPM7gdwzk0DFgG9gATgKHB3AeUVEZEc5FrozrmvgcuzeX5apvsOGBLaaCIikhc6UlRExCdU6CIiPqFCFxHxCRW6iIhPqNBFRHxChS4i4hMqdBERn1Chi4j4hApdRMQnVOgiIj6hQhcR8QkVuoiIT6jQRUR8QoUuIuITKnQREZ9QoYuI+EQwVywSKfLeWLUzqHF3tq9bwElEvKM1dBERn1Chi4j4hApdRMQnci10M7vQzD4zs2/NbLOZPZjNmE5mlmJmGwK3JwsmroiI5CSYL0VPAQ8759aZWSSw1syWOOe2ZBm3zDnXO/QRRUQkGLmuoTvn9jjn1gXuHwa+BWoXdDAREcmbPG1DN7NY4HJgVTaTO5rZRjNbbGbNcnh9nJmtMbM1ycnJeU8rIiI5CrrQzawi8C7wkHPuUJbJ64B6zrnLgBeA97Kbh3Mu3jnXxjnXJjo6+jwji4hIdoIqdDOLIKPMX3fOzc863Tl3yDmXGri/CIgws6iQJhURkXMKZi8XA2YC3zrnJuYwpkZgHGbWLjDf/aEMKiIi5xbMXi5XAv2Ab8xsQ+C5UUBdAOfcNOAWYLCZnQLSgNudcy70cUVEJCe5FrpzbjlguYyZCkwNVSgREck7HSkqIuITKnQREZ9QoYuI+IQKXUTEJ1ToIiI+oUIXEfEJFbqIiE+o0EVEfEKFLiLiEyp0ERGfUKGLiPiECl1ExCdU6CIiPqFCFxHxCRW6iIhPqNBFRHxChS4i4hMqdBERn1Chi4j4RK6FbmYXmtlnZvatmW02swezGWNmNsXMEszsazNrVTBxRUQkJ7leJBo4BTzsnFtnZpHAWjNb4pzbkmlMT6BR4NYe+Hvgp4iIFJJc19Cdc3ucc+sC9w8D3wK1swzrA8x1GVYCVcysZsjTiohIjvK0Dd3MYoHLgVVZJtUGdmV6nMjZpS8iIgUo6EI3s4rAu8BDzrlDWSdn8xKXzTzizGyNma1JTk7OW1IRETmnoArdzCLIKPPXnXPzsxmSCFyY6XEdYHfWQc65eOdcG+dcm+jo6PPJKyIiOQhmLxcDZgLfOucm5jDsA6B/YG+XDkCKc25PCHOKiEgugtnL5UqgH/CNmW0IPDcKqAvgnJsGLAJ6AQnAUeDukCcVEZFzyrXQnXPLyX4beeYxDhgSqlAiIpJ3OlJURMQnVOgiIj6hQhcR8QkVuoiIT6jQRUR8QoUuIuITweyHLhK0N1btDGrcne3rFnASkZJHa+giIj6hQhcR8QkVuoiIT6jQRUR8QoUuIuITKnQREZ9QoYuI+IQKXUTEJ1ToIiI+oSNFpcg5ePAg69evJyEhgeXfbOP4saOEh5eiXMVKVI2uQY16Dal+4UWEhYd7HVWkSFGhS5Fw5MgRPvzwQ95//302bNhAeno6AGXKladMufKcPnWao6kpuDPPV6DBpW249IouNLy0rcpdBBW6eCwtLY25c+cye/ZsUlJSaNiwIYMHD6Z9+/Y0adKEhVt+PTM2Pf00B/ftZfe2rfy0eR1bN6xky1efU7lada664XYuu6ob4aUiPPxtRLylQhfPfPHFF4wePZrdu3fTqVMn4uLiaNmyJWaZL2H7n0IPCwvngphaXBBTi+YdOnH61El+2PgVKz58iw9fmcyyD97gmj59aXl1dyxMXw9JyZNroZvZLKA3kOSca57N9E7A+8BPgafmO+fGhDCj+MypEycYM2YM8+bNo0GDBsyZM4d27drleT7hpSJo0vpKGre6gm2b1/H5glf55+xJfL1iCb0GPEh0LZ3RUUqWYNbQXwGmAnPPMWaZc653SBKJrx0+uJ93pozh523fcddddzF8+HDKlCmTr3maGQ2at+aiZq3YuOxfLHkrnvgnBnNNnz9xZe8/Eham7etSMuRa6M65L8wsthCyiM/t/yWR18c/xtHUQ0yePJnu3buHdP5mRstrrqdRy/Z8/PrfWTp/Dju+28jv7xtJxcpVQ/peIkVRqDY0djSzjWa22Mya5TTIzOLMbI2ZrUlOTg7RW0txsHfXT7zy9P9y8sRx7npsQsjLPLMKlapw8+DHuOHe/2VXwre8/OQDbP92Y4G9n0hREYpCXwfUc85dBrwAvJfTQOdcvHOujXOuTXR0dAjeWoqD5J938NpzIwgvFcGAxydSM7ZRobxvy6uv594np1CmfAVeGz+SdUsXFcr7ingl34XunDvknEsN3F8ERJhZVL6TiS8cPrifN/7fKMzC6DdiHNVq1CnU94+pE8u9T06hQfPWfPjKZMaPH39mH3cRv8l3oZtZDQvsZ2Zm7QLz3J/f+Urxd/L4Md6a/BRpqYe54+GxhV7mvylTrjx/fHA0rbv0ZtasWQwfPpwTJ054kkWkIAWz2+I8oBMQZWaJwF+BCADn3DTgFmCwmZ0C0oDbnXOuwBJLseDS03kv/jn2bP+B2/7nr9Ss19DTPGHh4fTsN5Tr2jZj3LhxpKam8sILL1C+fHlPc4mEUjB7udyRy/SpZOzWKHLG0gVz+W7tCq774yAaX97R6zhAxl4wAwYMoHLlyvzlL39h4MCBTJs2jUqVKnkdTSQkdDidhNyPm9ayfOE8Wl59PR16/MHrOGe56aabmDRpEps2bWLAgAEcPHjQ60giIaFCl5BKTfmV9+PHE127Hj36DclyGH/R0b17d1588UUSEhIYOHAghw4d8jqSSL6p0CVk0tPT+eDl8RxPO8LNg0cRUTp/R4AWtKuvvpopU6awdetWBg0aRGpqqteRRPJFhS4hM2fOHH7ctJZud8QRUyfW6zhB6dSpExMnTmTLli3cf//9HDlyxOtIIudNhS4hsW3bNiZNmkTjVlfQunPxOq3Pddddx/jx41m/fj0PPPAAaWlpXkcSOS8qdMm39PR0nnjiCcqXL0+vu/6nyG43P5cePXowbtw4Vq9ezYMPPsjJkye9jiSSZyp0ybd58+axbt06Ro4s3ifB6t27N0899RTLli3j8ccf1xGlUuyo0CVfdu/ezcSJE7nyyivp06eP13Hy7bbbbuPBBx9k4cKFPPfcc+gYOSlOdMUiOW/OOcaMybiWyejRo4vlppbs3Hfffezfv585c+ZQrVo1Bg0a5HUkkaCo0OW8ffbZZ3z++eeMGDGC2rVrex0nZMyMxx57jAMHDjBx4kSqVavGzTff7HUskVyp0OW8HD9+nGeffZYGDRrwpz/9yes4IRcWFsYzzzzDwYMHefLJJ6latSqdO3f2OpbIOWkbupyX2bNns2vXLh5//HEiIiK8jlMgSpcuzZQpU2jatCnDhw9n7dq1XkcSOScVuuTZnj17iI+Pp3v37nTsWDROvFVQKlSowPTp06lZsyZDhgwhISHB60giOVKhS54999xzpKen8+ijj3odpVBccMEFvPzyy0RERBAXF8fevXu9jiSSLRW65Mn69ev56KOPGDhwoK++CM1NnTp1mD59OikpKcTFxXH48GGvI4mcRV+KlnBvrNoZ1Lg729fFOcf48eOJiorinnvuKeBkBSMvv29Wl1xyCVOmTOH+++9n2LBhxMfHU7p06VBHFDlvWkOXoC1ZsoT169czbNiwEnulnyuvvJKxY8eyatUqHnvsMR1NKkWK1tAlKCdPnmTixIk0aNCgxO+T3adPH5KSkpg4cSLVq1cvMd8lSNGnQpegvP322+zYsYOXXnqJUqW02AwcOJC9e/cye/ZsYmJiGDBggNeRRFTokrvjaUeZ8dJLtG3blk6dOnkdp0j47WjSpKQkxo0bR0xMDL169fI6lpRwuW5DN7NZZpZkZptymG5mNsXMEszsazNrFfqY4qWvlrzHgQMHePjhh31zvpZQCA8P57nnnqN169aMHDmSVatWeR1JSrhgvhR9Behxjuk9gUaBWxzw9/zHkqLi2JFUVn70Lp07d+ayyy7zOk6RU7ZsWaZOnUrdunUZOnQoW7du9TqSlGC5Frpz7gvgwDmG9AHmugwrgSpmVjNUAcVbKz9+l2NHUxk2bJjXUYqsKlWqEB8fT/ny5Rk0aBC7d+/2OpKUUKHYbbE2sCvT48TAc2cxszgzW2Nma5KTk0Pw1lKQjh5OYdW/FtC07dU0bdrU6zhFWq1atYiPj+fo0aPExcWRkpLidSQpgUJR6NltVM32qgDOuXjnXBvnXJvo6OgQvLUUpC8XvcPJ48e59qb+XkcpFho3bszUqVPZuXMnQ4YM4fjx415HkhImFIWeCFyY6XEdQP/nLOZSDx5g9f99QPOOnYmudfZRk5K99u3b88wzz7B27VpGjBjB6dOnvY4kJUgoCv0DoH9gb5cOQIpzbk8I5iseWvHhW5w+dZJr+vjvXOcF7Xe/+x0jRozg448/5tlnn9Vl7KTQ5LofupnNAzoBUWaWCPwViABwzk0DFgG9gATgKHB3QYWVwpGyP4m1n31Iy6u6c0H1knMCrlAaMGAAv/zyC3PmzKFGjRrce++9XkeSEiDXQnfO3ZHLdAcMCVki8dzyhfMAuFpr5/ny6KOPsnfvXiZMmEB0dDQ33nij15HE53SkqPyXA0m72bDsY1p3+h2Vq8V4HadYCwsLY9y4cRw4cIDHH3+cyMhIXcZOCpTOtij/Zdn7rxMWFs6VN9zudRRfKF26NFOnTqVx48Y89NBDrFy50utI4mMqdDlj3+6dfPPlp7TteiORVap5Hcc3IiMjefnll6lXrx5Dhgxh/fr1XkcSn1Khyxmfv/caEWXK0LHXrV5H8Z2qVasyc+ZMoqKiuO+++9i8ebPXkcSHVOgCwN5d29jy1ee063YTFSpV8TqOL0VHRzN79mwqVqzIoEGDdMFpCTkVugCwdP5cypSrQIceJfviFQWtVq1azJo1i/DwcO655x527gzukngiwVChC7u3fc/W9f+mY89bKFch0us4vhcbG8vMmTM5efIkd999Nz///LPXkcQnVOjC0vlzKFexEu26/d7rKCXGxRdfzIwZM0hNTaVfv34kJiZ6HUl8QIVewu3cuokfN63lil63UaZcybzws1eaNWvGrFmzOHr0KP3799fmF8k3FXoJ5pzj03/MomKVC2jb9Qav45RIzZo1Y/bs2aSlpdG/f3+2b9/udSQpxlToJdjy5cvZtXUzV994JxFlynodp8Rq2rQpc+bM4eTJk/Tv35+ffvrJ60hSTKnQSyjnHM8//zxVomtw+TXnusKgFIaLL76YV155hfT0dPr27cuWLVu8jiTFkAq9hFqyZAlbtmzh2t/3I7xUhNdxBGjUqBGvvvoqZcqUoX///nz11VdeR5JiRoVeAp0+fZrJkyfToEEDmnfUyaKKkvr16/PGG29Qo0YNBg0axKeffup1JClGdLbFEmjhwoVs27aNyZMnsy8sPKjXvLEqtHtghHp+fpD5M/n98GeZN/EvDB02jBvu+V8uu6rbmWl3tg/uClLBfsZ+mZ9oDb3EOXHiBC+88ALNmjWjW7duub9APFG+YiX6PTqO2KaX8cGMCSxfOE9XPpJcqdBLmHfeeYfdu3czfPhwzLK7vrcUFaXLluP2h8bQvENnPnv3FRbOmsjpUye9jiVFmDa5lCBHjhxh2rRptG3bliuuuMLrOBKEUhGl+f19I7igem2+eP81Uvbt5cbm06lcubLX0aQI0hp6CTJz5kz27dvHI488orXzYsTMuPamfvSJe5RdP2zhjjvu0FGlki0Vegmxd+9eZs+eTa9evWjRooXXceQ8tLiiK33//AwHDhzgtttuY8WKFV5HkiImqEI3sx5m9r2ZJZjZyGymdzKzFDPbELg9Gfqokh9Tpkzh9OnTDB8+3Osokg91G1/KW2+9RUxMDHFxccTHx+vLUjkj10I3s3DgRaAncAlwh5ldks3QZc65loHbmBDnlHz47rvvWLBgAX379qVOnTpex5F8qlevHm+++SY9e/Zk0qRJDBs2jNTUVK9jSREQzBp6OyDBObfNOXcCeBPoU7CxJFScc4wfP55KlSpx3333eR1HQqR8+fKMHz+ekSNHsnTpUm699Va2bt3qdSzxWDCFXhvYlelxYuC5rDqa2UYzW2xmzbKbkZnFmdkaM1uTnJx8HnElrz755BO+/PJLHnjgAe0Z4TNmxl133cXs2bNJTU3l1ltv5fXXX9cmmBIsmN0Ws9sdIusSsw6o55xLNbNewHtAo7Ne5Fw8EA/Qpk0bLXUFLC0tjWeeeYZGjRpx5513eh2nWAnlkawFfaRj27ZtWbBgAaNGjWLs2LEsX76cp59+ukDfU4qmYNbQE4ELMz2uA+zOPMA5d8g5lxq4vwiIMLOokKWU8/Lyyy+zZ88ennjiCUqV0iEHfhYVFcX06dMZNWoUK1asoE+fPmzbtNbrWFLIgin01UAjM6tvZqWB24EPMg8wsxoW2LHZzNoF5rs/1GEleDt37mTmzJn07t2btm3beh1HCoGZ0a9fP95++20qV67M6xNG8eErkzmedsTraFJIci1059wpYCjwMfAt8LZzbrOZ3W9m9weG3QJsMrONwBTgdqcNeZ5xzvG3v/2NUqVK8ec//9nrOFLImjRpwjvvvEOHHn9g/ecfMe3xOH7YqFPxlgRB/T88sBllUZbnpmW6PxWYGtpocr4WLVrE559/zogRI4iJifE6jnigXLlydLs9jkvaXcPCmRN5c9ITXNqxK93vvI/ykfpy3K90pKjPHDhwgKeffpoWLVrQr18/r+OIx2pf1IRBo1/kmj592fzVUl4aeS9rPl1Ievppr6NJAVCh+8zYsWNJTU1l7NixhIcHd65z8bfwUhFce1M/Bo15iep1L2Lx3KnM+OtQdnz/jdfRJMRU6D6yZMkSFi9ezODBg2nU6Ky9RqWEi6kdS99Hx/GHIX/h2NFU5j7zCPP//gyJiYleR5MQUaH7xP79+xkzZgxNmzZl4MCBXseRIsrMuKTt1Qz+28tc06cv36/7kl69ejF27Fj27dvndTzJJ+2cXETk53Jc6enpjBgxgsOHDzNjxgwiInTRZzm3iDJlufamflzeqSe/rPyAN998k/nz59O/f3/uueceKlWqdN7z1uUFvaM1dB+YOXMmK1asYNSoUTRu3NjrOFKMVKoaxVNPPcWHH35Ily5dmD59Ol27dmXixIlaYy+GVOjF3Lp165g8eTI9evTg1ltv9TqOFFP16tVjwoQJLFiwgKuuuooZM2bQtWtXRo8erW3sxYg2uRRjycnJPPLII9SqVYsxY8boKkSSb02aNGHSpEls376dWbNm8e677/LOO+/QpUsXarS6jnpNLtNyVoSp0IuptLQ0hgwZwsGDB3n11VeJjIz0OpL4SGxsLGPGjGHo0KHMnTuXf/zjH6QsWUJUrbq06XIDLa7sSplyFbyOKVlok0sxdPr0aR599FE2bdrEhAkTaNYs27MVi+RbTEwMjzzyCEuXLuXGex8mokxZPnrtRSY9dCfvvzyen7ZswKWnex1TArSGXgxNmDCBTz75hFGjRtGlSxev40gJULZsWS67ujuXXd2dn7d9x/qli9my+gu+XvEJlavFcOkVXWneoTPRtet5HbVEU6EXI845nn/+eV555RX69u2rQ/vFE7UvakLti5pwfd8H+H7dl3y9fAkr/vkWyxfOI6pWXZq0voqmba+i+oUXaXt7IVOhFxPOOf7vrRn8+6N/cMsttzBy5FnX6hYpVBGly9C8Q2ead+jM4YP7+W7NCr5ds4wV/3yT5QvfoGpMTRq2aEeDS9sQ26QFEWXKeh3Z91ToxYBLT+fjN6ax+pP3adP1Bpr0upe31vx8ztcU9FVyJHgl4UCbyCrVaHvdjbS97kaOHDrI9+u+5Pt1X7L+88Ws/uR9wktFUK9JCxpc2oYGzVsTVauu1t4LgAq9iEtLPcR78eNJ+PorOvT4A9f9cZD+IkiRVqFSFVp16kWrTr04eeI4O7du4sdv1vDjN2tYMm86S4DykZW58OLm1G3UjLoXN6d63QaE66pa+aZPsAj7ZcePvDN1DIcO7KNn/6G07txbZS7FSkTpMjRo3poGzVvDHfdxcN9eftq8jp1bN7Fr62a+X7vizLjaDZpS66LG1KzXkBqxDakaXVPLex6p0Iug06dOsvLj+Xzx3muUq1iJux6bQJ2GTb2OJZJvVaKqc/m1Pbn82p4AHP51P7t+2MyuHzazc+smVn70LumnTwFQplwFasY25OeOrWnUqBENGzakfv36VKig/d9zokIvYnZ8/w2L5kxh3+6dNG51Bb3u+h8qVq7qdSyRAhFZtRqXtLuGS9pdA8CpkydI/nkHe7b/wC87Etiz/Qdee+01Tpw4ceY1NWvWpGHDhlx00UXExsZy4YUXUqdOHWrWrEnp0qW9+lWKBBV6EeCc46ctG/j3orf5cdNaqkRV548Pjebilh28jiZSqEpFlKZmbCNqxv7nfP63tqpJYmIiCQkJ/Pjjj2zbto2EhARWr17NsWPHzowLCwujevXqZwq+Tp061KhRg5iYGGJiYqhevTqRkZG+3oyjQvdQUlISS5YsYf78+WzZsoUKlarS+Za7ad/t99rFSyQgIiKC+vXrU79+fbp163bm+fT0dJKSkkhMTCQxMZFdu3axa9cuEhMTWbZsGcnJyWfNq2zZsmcKPjo6mqpVq3LBBRdQpUqVM/erVq1K1apVqVKlSrFb4w+q0M2sBzAZCAdmOOeezTLdAtN7AUeBAc65dSHOWuylpqayceNGNmzYwL///W/WrVuHc45GjRrR++6HuLRjV0oVswVIxCthYWHUqFGDGjVq0KZNm7OmHzt2jKSkJJKSkti7dy/Jycns3buXvXv3kpSUxObNm0lJSSElJSXH96hQoQKRkZFUrFiRihUrEhkZme1zFStWpEKFCpQrV46yZcv+18/f7pctW5awsII920quhW5m4cCLQDcgEVhtZh8457ZkGtYTaBS4tQf+HvhZZDnnzvz87Zb5ceZxOY2FjLWEY8eOcfToUdLS0jh69CiHDx9m//797N+/n19++YUdO3awfft2du/ejXOOsLAwGjduzNChQ7n++utp0KBBidhXWaQwlS1blrp161K37rmPyTh58iQpKSn8+uuvZ90OHjxIamoqhw8fJjU1lYMHD5KYmMjhw4c5cuQIaWlpecpUpkwZypUrR//+/Rk8eHB+fr1sBbOG3g5IcM5tAzCzN4E+QOZC7wPMdRlNt9LMqphZTefcnlAHXrJkCSNHjsx3EReWyMhI6tWrx+WXX87NN99My5YtadGiBRUrViz0LCJytoiICKKiooiKisrza0+ePMmRI0dITU09U/DHjh0jLS3tv+5nfa6grvkbTKHXBnZlepzI2Wvf2Y2pDfxXoZtZHBAXeJhqZt/nKW1wooAidamV1atXZ35YKPn+dP4vLXKfXxYlPl8+/mwhm3z5nF+oFfV84P0ymOMZ0IIp9Oy+Es66qhvMGJxz8UB8EO953sxsjXPu7A1qRYTy5Y/y5Y/y5V9RzhjMFvpE4MJMj+sAu89jjIiIFKBgCn010MjM6ptZaeB24IMsYz4A+luGDkBKQWw/FxGRnOW6ycU5d8rMhgIfk7Hb4izn3GYzuz8wfRqwiIxdFhPI2G3x7oKLnKsC3aQTAsqXP8qXP8qXf0U2o3mx54eIiISerikqIuITKnQREZ8oNoVuZj3M7HszSzCzs66/ZmZ/MrOvA7cvzeyyTNO2m9k3ZrbBzNZ4lK+TmaUEMmwwsyeDfW0hZvxzpnybzOy0mV0QmFagn6GZzTKzJDPblMN0M7Mpgexfm1mrYH+vQsrn9fKXWz5Pl78g8nm27AXe40Iz+8zMvjWzzWb2YDZjPF0Gg5L5KMqieiPjy9gfgYuA0sBG4JIsY64Aqgbu9wRWZZq2HYjyOF8n4J/n89rCyphl/A3Ap4X4GV4DtAI25TC9F7CYjGMeOvz251uIn19u+Txb/oLM5/Xyd858Xi57gfeoCbQK3I8Etmbzd9jTZTCYW3FZQz9z+gHn3Angt9MPnOGc+9I592vg4Uoy9oUvMvkK6LUFmfEOYF4B5MiWc+4L4MA5hpw5vYRzbiVQxcxqUkifX275PF7+gvn8clIkPr8sCnXZA3DO7XGBEwo65w4D35JxtHtmni6DwSguhZ7TqQVyci8Z/5L+xgH/MrO1lnH6Aa/ydTSzjWa22Mya5fG1hZURMysP9ADezfR0QX+Guckpf2F9fnlR2MtfsLxc/oJSFJY9M4sFLgdWZZlU5JfB4nI+9KBOLQBgZp3J+At1Vaanr3TO7TazGGCJmX0XWGMozHzrgHrOuVQz6wW8R8bZKYP+3fIpL+9zA7DCOZd5jaqgP8Pc5JS/sD6/oHi0/AXD6+UvWJ4ue2ZWkYx/TB5yzh3KOjmblxSpZbC4rKEHdWoBM2sBzAD6OOf2//a8c2534GcSsICM/yIVaj7n3CHnXGrg/iIgwsyignltYWXM5Hay/Je3ED7D3OSUv8icdsLD5S9XRWD5C5Zny56ZRZBR5q875+ZnM6TIL4OFvtH+fG5k/E9iG1Cf/3zp0CzLmLpkHKl6RZbnKwCRme5/CfTwIF8N/nMgVztgJxn/suf62sLKGBhXmYxtnRUK8zMMzDuWnL/U+x3//YXUV3n5vQohn2fLX5D5PF3+cstXBJY9A+YCz59jjOfLYG63YrHJxQV3+oEngWrAS5ZxzcBTLuOMaNWBBYHnSgFvOOc+8iDfLcBgMzsFpAG3u4ylIdvXhjJfHjIC3AT8yzl3JNPLC/wzNLN5ZOyJEWVmicBfgYhM2bI9vUROv1coswWZz7PlL8h8ni5/QeQDj5a9gCuBfsA3ZrYh8NwoMv6hLhLLYDB06L+IiE8Ul23oIiKSCxW6iIhPqNBFRHxChS4i4hMqdBERn1Chi4j4hApdRMQn/j860B1WKCufAgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.distplot(x=df.iloc[0,:], bins= 20, fit=lognorm, kde=False)\n",
    "print(f\"\"\"\n",
    "\n",
    "mean: {[df.iloc[0,:].mean(), df.iloc[-1,:].mean()]}\"\"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 472,
   "id": "aa2cc96e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.002691316318290994"
      ]
     },
     "execution_count": 472,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.values.min()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2e01a91",
   "metadata": {},
   "source": [
    "## Comparison with what is implemented in CategoryCov"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40330c0d",
   "metadata": {},
   "source": [
    "Consider the negative eigenvalues are set to 0 in CategoryCov sampling method, but not in our proposed procedure.\n",
    "Eigenvalues could be set to zero in the latter as well, if one wants."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 473,
   "id": "b4e410fe",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:  skip MAT125/MF3/MT1 because interpolation schme is not lin-lin\n",
      "WARNING:  skip MAT125/MF3/MT102 because interpolation schme is not lin-lin\n"
     ]
    }
   ],
   "source": [
    "zam = 10010\n",
    "endf6 = sandy.get_endf6_file('jeff_33','xs', zam)\n",
    "Cov = sandy.XsCov.from_endf6(endf6)\n",
    "BE = sandy.Xs.from_endf6(endf6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 474,
   "id": "55884b69",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Cov.reset_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 475,
   "id": "41e369e7",
   "metadata": {},
   "outputs": [],
   "source": [
    "CatCov = sandy.CategoryCov(Cov)\n",
    "CatCov_semp = pd.DataFrame(np.cov(CatCov.sampling(nsmp=10000, tolerance=0).data.T))\n",
    "# CatCov_semp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 476,
   "id": "08e4134d",
   "metadata": {},
   "outputs": [],
   "source": [
    "Cov_ = sandy.CategoryCov(Cov)\n",
    "E, V = Cov_.eig(tolerance = 0)\n",
    "Cov_ = V.dot(np.diag(E)).dot(np.linalg.inv(V))\n",
    "df = pd.DataFrame(np.random.multivariate_normal(np.ones(len(Cov)), Cov_, size=10000))\n",
    "# df.cov()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 477,
   "id": "bc70114e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Maximum absolute difference in the matrix: 0.0018864399152059563\n",
      "Maximum absolute difference in the eigenvalues: 1.0681724236465893\n"
     ]
    }
   ],
   "source": [
    "Delta_CategoryCov = abs(CatCov_semp.values-Cov.values)\n",
    "Delta_E_CategoryCov = abs(np.linalg.eig(CatCov_semp)[0]-np.linalg.eig(Cov)[0])\n",
    "print(f\"\"\"Maximum absolute difference in the matrix: {Delta_CategoryCov.max()}\n",
    "Maximum absolute difference in the eigenvalues: {Delta_E_CategoryCov.max()}\"\"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 478,
   "id": "f6473742",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Maximum absolute difference in the matrix: 0.0010094720053583917\n",
      "Maximum absolute difference in the eigenvalues: 1.058357502846415\n"
     ]
    }
   ],
   "source": [
    "Delta_df = abs(df.cov().values-Cov.values)\n",
    "Delta_E_df = abs(np.linalg.eig(Delta_df)[0].max()-np.linalg.eig(Cov)[0].max())\n",
    "print(f\"\"\"Maximum absolute difference in the matrix: {Delta_df.max()}\n",
    "Maximum absolute difference in the eigenvalues: {Delta_E_df}\"\"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 479,
   "id": "901f784f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# BEST = {}\n",
    "# BEST_E = {}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 480,
   "id": "1d44ec52",
   "metadata": {},
   "outputs": [],
   "source": [
    "if Delta_df.max() < Delta_CategoryCov.max():\n",
    "    BEST[zam] = 1 # if the newly peroposed method gives better covariance representation accoring to this measure\n",
    "else:\n",
    "    BEST[zam] = 0\n",
    "    \n",
    "if Delta_E_df.max() < Delta_E_CategoryCov.max():\n",
    "    BEST_E[zam] = 1 # if the newly peroposed method gives better covariance representation accoring to this measure\n",
    "else:\n",
    "    BEST_E[zam] = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 481,
   "id": "01713e6d",
   "metadata": {},
   "outputs": [],
   "source": [
    "b = pd.Series(BEST)\n",
    "e = pd.Series(BEST_E)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 482,
   "id": "66d61fe6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.4166666666666667"
      ]
     },
     "execution_count": 482,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(b[b==1])/len(BEST)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 483,
   "id": "93e9bbea",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5833333333333334"
      ]
     },
     "execution_count": 483,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(e[e==1])/len(BEST_E)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 484,
   "id": "8d96d3a8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{10010: 1,\n",
       " 922350: 0,\n",
       " 942390: 1,\n",
       " 942410: 0,\n",
       " 922330: 0,\n",
       " 10020: 1,\n",
       " 30060: 1,\n",
       " 40090: 1,\n",
       " 100200: 0,\n",
       " 441020: 0,\n",
       " 631550: 0,\n",
       " 812050: 0}"
      ]
     },
     "execution_count": 484,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "BEST"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 485,
   "id": "47d9cae6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{10010: 1,\n",
       " 922350: 0,\n",
       " 942390: 0,\n",
       " 942410: 1,\n",
       " 922330: 1,\n",
       " 10020: 1,\n",
       " 30060: 1,\n",
       " 40090: 0,\n",
       " 100200: 1,\n",
       " 441020: 0,\n",
       " 631550: 0,\n",
       " 812050: 1}"
      ]
     },
     "execution_count": 485,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "BEST_E"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 486,
   "id": "4b021f1c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 486,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import os\n",
    "os.system('jupyter nbconvert --to html random_multivariate.ipynb')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python3 (sandy-devel)",
   "language": "python",
   "name": "sandy-devel"
  },
  "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.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}