{ "cells": [ { "cell_type": "code", "execution_count": 2, "id": "8d31132f", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import pandas as pd\n", "import seaborn as sns\n", "from fitter import Fitter, get_common_distributions, get_distributions\n", "from scipy import io" ] }, { "cell_type": "code", "execution_count": 3, "id": "945db1eb", "metadata": {}, "outputs": [], "source": [ "suv_filename = \"/home/marija/maki/suv_percentilesSLOthenUWM.mat\"\n", "flags_filename=\"/home/marija/maki/flags_combined.mat\"\n", "normal_range_filename =\"/home/marija/maki/normal_range.mat\"" ] }, { "cell_type": "code", "execution_count": 4, "id": "e3ac38f3", "metadata": {}, "outputs": [], "source": [ "\"\"\"\n", "Loading all data\n", "\"\"\"\n", "suv_dict = io.loadmat(suv_filename)\n", "flags_dict = io.loadmat(flags_filename)\n", "normal_range_dict = io.loadmat(normal_range_filename)" ] }, { "cell_type": "code", "execution_count": 5, "id": "b8dadd66", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(58, 17, 100)\n", "(58,)\n" ] } ], "source": [ "\"\"\"\n", " Extracting only data associated to lungs\n", "\"\"\"\n", "lung_suv = suv_dict['lung_SUVperc_COMBINED'][0:58,:]\n", "flags_lung = flags_dict['flags'][0:58,3] # 0 == NC, 1 == AE\n", "\n", "print(lung_suv.shape)\n", "print(flags_lung.shape)\n" ] }, { "cell_type": "code", "execution_count": 6, "id": "f92743a3", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0\n", " 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", "NC_lung: 53\n", "AE_lung: 5\n" ] } ], "source": [ "print(flags_lung)\n", "print(\"NC_lung:\", np.count_nonzero(flags_lung==0))\n", "print(\"AE_lung:\", np.count_nonzero(flags_lung==1))" ] }, { "cell_type": "code", "execution_count": 7, "id": "ee2874bd", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(58, 17, 100)\n", "(58,)\n" ] } ], "source": [ "\"\"\"\n", " Extracting only data associated to bowel\n", "\"\"\"\n", "bowel_suv = suv_dict['bowel_SUVperc_COMBINED'][0:58,:]\n", "flags_bowel = flags_dict['flags'][0:58,1] # 0 == NC, 1 == AE\n", "\n", "print(bowel_suv.shape)\n", "print(flags_bowel.shape)" ] }, { "cell_type": "code", "execution_count": 8, "id": "7f427125", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0\n", " 0 0 1 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0]\n", "NC_bowel: 52\n", "AE_bowel: 6\n" ] } ], "source": [ "print(flags_bowel)\n", "print(\"NC_bowel:\", np.count_nonzero(flags_bowel==0))\n", "print(\"AE_bowel:\", np.count_nonzero(flags_bowel==1))" ] }, { "cell_type": "code", "execution_count": 9, "id": "07fe6dcd", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(58, 17, 100)\n", "(58,)\n" ] } ], "source": [ "\"\"\"\n", " Extracting only data associated to thyroid\n", "\"\"\"\n", "thyroid_suv = suv_dict['thyroid_SUVperc_COMBINED'][0:58,:]\n", "flags_thyroid = flags_dict['flags'][0:58,5] # 0 == NC, 1 == AE\n", "\n", "print(thyroid_suv.shape)\n", "print(flags_thyroid.shape)" ] }, { "cell_type": "code", "execution_count": 10, "id": "9bbea04e", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 1 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0\n", " 0 1 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0]\n", "NC_thyroid: 49\n", "AE_thyroid: 9\n" ] } ], "source": [ "print(flags_thyroid)\n", "print(\"NC_thyroid:\", np.count_nonzero(flags_thyroid==0))\n", "print(\"AE_thyroid:\", np.count_nonzero(flags_thyroid==1))" ] }, { "cell_type": "code", "execution_count": 11, "id": "8e0492f1", "metadata": {}, "outputs": [], "source": [ "p = 95 # percentile of interest\n", "X_lung = np.nanmax(lung_suv[:,:,p], axis=1) # taking SUV max and ignoring NaNs\n", "\n", "X_bowel = np.nanmax(bowel_suv[:,:,p], axis=1)\n", "\n", "X_thyroid = np.nanmax(thyroid_suv[:,:,75], axis=1)" ] }, { "cell_type": "code", "execution_count": 12, "id": "55666bc8", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(53,)" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "X_lung_NC = X_lung[flags_lung == 0]\n", "X_lung_AE = X_lung[flags_lung == 1]\n", "\n", "X_lung_NC.shape" ] }, { "cell_type": "code", "execution_count": 13, "id": "bfb51a60", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(52,)" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "X_bowel_NC = X_bowel[flags_bowel == 0]\n", "X_bowel_AE = X_bowel[flags_bowel == 1]\n", "\n", "X_bowel_NC.shape" ] }, { "cell_type": "code", "execution_count": 14, "id": "32430668", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(49,)" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "X_thyroid_NC = X_thyroid[flags_thyroid == 0]\n", "X_thyroid_AE = X_thyroid[flags_thyroid == 1]\n", "\n", "X_thyroid_NC.shape" ] }, { "cell_type": "code", "execution_count": 15, "id": "dc4a8459", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "numpy.ndarray" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "type(X_bowel_NC)" ] }, { "cell_type": "code", "execution_count": 16, "id": "16fba333", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "numpy.ndarray" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "type(X_lung_NC)" ] }, { "cell_type": "code", "execution_count": 17, "id": "a63a15cd", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "numpy.ndarray" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "type(X_thyroid_NC)" ] }, { "cell_type": "code", "execution_count": 18, "id": "36ca3761", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 2.37582922 8.32889462 2.49695492 2.62231517 3.1779449 2.48935366\n", " 2.64083982 3.21895695 2.65162754 2.81365871 13.13622665 4.00094509\n", " 6.43130398 7.14293861 2.63764477 4.91679096 3.0133276 4.00735998\n", " 2.82554626 2.3756032 3.05244803 5.04009581 2.75375438 2.96612239\n", " 2.84024048 3.02956247 3.56154656 2.63551474 3.54600501 3.10588789\n", " 5.65029192 3.35723162 2.60515237 2.84188676 2.41808271 2.51086831\n", " 2.64900827 2.2847507 3.04761338 2.94424272 2.26703072 2.83462191\n", " 3.44818497 2.7778244 5.96339226 3.04745531 2.66029692 5.41591644\n", " 2.07756186 2.52008891 1.9738239 7.01769114]\n" ] } ], "source": [ "print(X_bowel_NC)" ] }, { "cell_type": "code", "execution_count": 19, "id": "412eeb4e", "metadata": {}, "outputs": [], "source": [ "X_NC = np.concatenate((X_lung_NC,X_bowel_NC,X_thyroid_NC))" ] }, { "cell_type": "code", "execution_count": 20, "id": "527906e5", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(154,)" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "X_NC.shape" ] }, { "cell_type": "code", "execution_count": 21, "id": "97acc5a0", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(array([0.43247677, 0.29007588, 0.03164464, 0.02109643, 0.01054821,\n", " 0.01582232, 0.00527411, 0. , 0. , 0.00527411]),\n", " array([ 0.82417732, 2.05538226, 3.28658719, 4.51779212, 5.74899706,\n", " 6.98020199, 8.21140692, 9.44261186, 10.67381679, 11.90502172,\n", " 13.13622665]),\n", " )" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAMvklEQVR4nO3df4zfB13H8efL1kWZEtSdv9pqqzbMhgwh55xiNIpLuo3QEU3sRBgCaZY4fhiNKzHxHxIzojGQOGmaUYdxoSFjaMOKY5km/AEzvSGZbGPQjEmPDXeAgoixK7z9474zx+3a7+d+tJ/e2+cjufQ+P/r9vj/b9ZnPfb7fz12qCknS5vcdYw8gSdoYBl2SmjDoktSEQZekJgy6JDWxdawnvuyyy2rnzp1jPb0kbUoPPvjgl6pqZqVtowV9586dzM3NjfX0krQpJfnXs23zkoskNWHQJakJgy5JTRh0SWrCoEtSEwZdkpow6JLUhEGXpCYMuiQ1Mdqdouux8+A9oz33E7deN9pzS9K5eIYuSU0YdElqwqBLUhMGXZKaMOiS1IRBl6QmDLokNWHQJakJgy5JTRh0SWrCoEtSE4OCnmRvkseSnExy8Bz7/WySbyb5jY0bUZI0xNSgJ9kC3AZcA+wBbkiy5yz7vQO4d6OHlCRNN+QM/UrgZFU9XlWngaPAvhX2exPwAeDpDZxPkjTQkKBvA04tWZ6frPs/SbYBrwIOneuBkhxIMpdkbmFhYbWzSpLOYUjQs8K6Wrb8TuCWqvrmuR6oqg5X1WxVzc7MzAwcUZI0xJBfcDEP7FiyvB14ctk+s8DRJACXAdcmOVNVf7sRQ0qSphsS9BPA7iS7gC8A+4HfWrpDVe169vMkdwAfMuaSdGFNDXpVnUlyM4vvXtkCHKmqh5PcNNl+zuvmkqQLY9DvFK2q48DxZetWDHlVvW79Y0mSVss7RSWpCYMuSU0YdElqwqBLUhMGXZKaMOiS1IRBl6QmDLokNWHQJakJgy5JTRh0SWrCoEtSEwZdkpow6JLUhEGXpCYMuiQ1YdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJgy6JDVh0CWpCYMuSU0YdElqwqBLUhMGXZKaMOiS1IRBl6QmDLokNWHQJakJgy5JTRh0SWrCoEtSEwZdkpow6JLUhEGXpCYGBT3J3iSPJTmZ5OAK2/cleSjJJ5PMJfnFjR9VknQuW6ftkGQLcBtwNTAPnEhyrKoeWbLb/cCxqqokVwDvBy4/HwNLklY25Az9SuBkVT1eVaeBo8C+pTtU1derqiaLlwKFJOmCGhL0bcCpJcvzk3XfJsmrknwauAd4/UoPlOTA5JLM3MLCwlrmlSSdxZCgZ4V1zzkDr6oPVtXlwPXA21d6oKo6XFWzVTU7MzOzqkElSec2JOjzwI4ly9uBJ8+2c1V9FPjJJJetczZJ0ioMCfoJYHeSXUkuAfYDx5bukOSnkmTy+UuBS4Avb/SwkqSzm/oul6o6k+Rm4F5gC3Ckqh5OctNk+yHg14HXJnkG+G/gN5e8SCpJugAyVndnZ2drbm5uTX9358F7Nniai98Tt1439giSLgJJHqyq2ZW2eaeoJDVh0CWpCYMuSU0YdElqwqBLUhMGXZKaMOiS1IRBl6QmDLokNWHQJakJgy5JTRh0SWrCoEtSEwZdkpow6JLUhEGXpCYMuiQ1YdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJgy6JDVh0CWpCYMuSU0YdElqwqBLUhMGXZKaMOiS1IRBl6QmDLokNWHQJakJgy5JTRh0SWrCoEtSEwZdkpow6JLUxKCgJ9mb5LEkJ5McXGH7q5M8NPn4WJIXb/yokqRzmRr0JFuA24BrgD3ADUn2LNvtc8AvV9UVwNuBwxs9qCTp3IacoV8JnKyqx6vqNHAU2Ld0h6r6WFX9+2TxAWD7xo4pSZpmSNC3AaeWLM9P1p3NG4APr7QhyYEkc0nmFhYWhk8pSZpqSNCzwrpaccfkV1gM+i0rba+qw1U1W1WzMzMzw6eUJE21dcA+88COJcvbgSeX75TkCuB24Jqq+vLGjCdJGmrIGfoJYHeSXUkuAfYDx5bukOTHgLuB11TVZzZ+TEnSNFPP0KvqTJKbgXuBLcCRqno4yU2T7YeAPwZ+APjLJABnqmr2/I0tSVpuyCUXquo4cHzZukNLPn8j8MaNHU2StBreKSpJTRh0SWrCoEtSEwZdkpow6JLUhEGXpCYMuiQ1YdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJgy6JDVh0CWpCYMuSU0YdElqwqBLUhMGXZKaMOiS1IRBl6QmDLokNWHQJakJgy5JTRh0SWrCoEtSEwZdkpow6JLUhEGXpCYMuiQ1YdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJgy6JDUxKOhJ9iZ5LMnJJAdX2H55ko8n+Z8kf7DxY0qSptk6bYckW4DbgKuBeeBEkmNV9ciS3b4CvBm4/nwMKUmabsgZ+pXAyap6vKpOA0eBfUt3qKqnq+oE8Mx5mFGSNMCQoG8DTi1Znp+sW7UkB5LMJZlbWFhYy0NIks5iSNCzwrpay5NV1eGqmq2q2ZmZmbU8hCTpLIYEfR7YsWR5O/Dk+RlHkrRWQ4J+AtidZFeSS4D9wLHzO5YkabWmvsulqs4kuRm4F9gCHKmqh5PcNNl+KMkPA3PA84FvJXkrsKeqvnb+RpckLTU16ABVdRw4vmzdoSWff5HFSzGSpJF4p6gkNWHQJakJgy5JTRh0SWrCoEtSEwZdkpow6JLUhEGXpCYMuiQ1YdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJgy6JDVh0CWpCYMuSU0YdElqwqBLUhMGXZKaMOiS1IRBl6QmDLokNWHQJakJgy5JTWwdewANs/PgPaM99xO3Xjfac0sazjN0SWrCoEtSEwZdkpow6JLUhC+KaqqxXpD1xVhpdTxDl6QmPEOXVuB3JdqMPEOXpCYMuiQ14SUXXbTGvDtW2ow8Q5ekJgadoSfZC7wL2ALcXlW3LtueyfZrgW8Ar6uqT2zwrFJ7vhh74XT8+UhTz9CTbAFuA64B9gA3JNmzbLdrgN2TjwPAuzd4TknSFEMuuVwJnKyqx6vqNHAU2Ldsn33AX9eiB4AXJPmRDZ5VknQOQy65bANOLVmeB35uwD7bgKeW7pTkAItn8ABfT/LYqqYd12XAl8YeYgN0OI4OxwAX0XHkHWv+qxfNMazTBT2Odfz3Bvjxs20YEvSssK7WsA9VdRg4POA5LzpJ5qpqduw51qvDcXQ4BuhxHB2OAfocx5BLLvPAjiXL24En17CPJOk8GhL0E8DuJLuSXALsB44t2+cY8Nosugr4alU9tfyBJEnnz9RLLlV1JsnNwL0svm3xSFU9nOSmyfZDwHEW37J4ksW3Lf7O+Rt5NJvyUtEKOhxHh2OAHsfR4RigyXGk6jmXuiVJm5B3ikpSEwZdkpow6FMk2ZHkH5M8muThJG8Ze6a1SrIlyT8n+dDYs6xVkhckuSvJpyf/T35+7JlWK8nvTb6WPpXkfUm+a+yZhkhyJMnTST61ZN33J7kvyWcnf37fmDNOc5Zj+NPJ19NDST6Y5AUjjrguBn26M8DvV9VPA1cBv7vCjz7YLN4CPDr2EOv0LuDvq+py4MVssuNJsg14MzBbVS9i8Y0G+8edarA7gL3L1h0E7q+q3cD9k+WL2R089xjuA15UVVcAnwHedqGH2igGfYqqeurZHzRWVf/JYkC2jTvV6iXZDlwH3D72LGuV5PnALwHvAaiq01X1H6MOtTZbge9OshV4Hpvkno2q+ijwlWWr9wHvnXz+XuD6CznTaq10DFX1kao6M1l8gMX7aDYlg74KSXYCLwH+aeRR1uKdwB8C3xp5jvX4CWAB+KvJpaPbk1w69lCrUVVfAP4M+DyLPxrjq1X1kXGnWpcfevaek8mfPzjyPOv1euDDYw+xVgZ9oCTfA3wAeGtVfW3seVYjySuAp6vqwbFnWaetwEuBd1fVS4D/4uL/Fv/bTK4x7wN2AT8KXJrkt8edSgBJ/ojFS6x3jj3LWhn0AZJ8J4sxv7Oq7h57njV4GfDKJE+w+NMyfzXJ34w70prMA/NV9ex3SHexGPjN5NeAz1XVQlU9A9wN/MLIM63Hvz37k1Unfz498jxrkuRG4BXAq2sT35xj0KeY/PKO9wCPVtWfjz3PWlTV26pqe1XtZPEFuH+oqk13VlhVXwROJXnhZNXLgUdGHGktPg9cleR5k6+tl7PJXthd5hhw4+TzG4G/G3GWNZn8Ap9bgFdW1TfGnmc9DPp0LwNew+JZ7ScnH9eOPdT/Y28C7kzyEPAzwJ+MO87qTL67uAv4BPAvLP4b3BS3nSd5H/Bx4IVJ5pO8AbgVuDrJZ4GrJ8sXrbMcw18A3wvcN/n3fWjUIdfBW/8lqQnP0CWpCYMuSU0YdElqwqBLUhMGXZKaMOiS1IRBl6Qm/hcl5fYvyI1UCwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "\"\"\"\n", " Plotting stuff\n", "\"\"\"\n", "import matplotlib.pyplot as plt\n", "\n", "plt.hist(X_NC,density=True, label=\"NC\")" ] }, { "cell_type": "code", "execution_count": 23, "id": "7d5d4b62", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(array([ 0., 3., 44., 35., 23., 23., 11., 2., 2., 1., 2., 2., 1.,\n", " 0., 3., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", " 1., 0.]),\n", " array([ 0. , 0.5, 1. , 1.5, 2. , 2.5, 3. , 3.5, 4. , 4.5, 5. ,\n", " 5.5, 6. , 6.5, 7. , 7.5, 8. , 8.5, 9. , 9.5, 10. , 10.5,\n", " 11. , 11.5, 12. , 12.5, 13. , 13.5, 14. ]),\n", " )" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAMm0lEQVR4nO3df6idh13H8ffHpLJ1U9rSmxqb4p0SupXi2hFmteAfzQLRlKZ/OOhwJWAl/2zayWSmDgT/kYAyJyhK6GoDKx2lqzS0qAvZyhBq3W3Xn2YzQ2OXGZu7jbmp4Kz7+sd9Mu5u7s05ufeee+63e78gPOd57jnn+XJz7jvPfc55SKoKSVI/PzLtASRJq2PAJakpAy5JTRlwSWrKgEtSU1s3cmdXX311zc7ObuQuJam9Z5999utVNbN0+4YGfHZ2lrm5uY3cpSS1l+Rfl9vuKRRJasqAS1JTBlySmjLgktSUAZekpgy4JDVlwCWpKQMuSU0ZcElqakOvxNyMZg89OfZ9Tx/eN8FJJOnSeAQuSU0ZcElqyoBLUlMGXJKaMuCS1JQBl6SmDLgkNWXAJakpAy5JTRlwSWrKgEtSUwZckpoy4JLUlAGXpKbGDniSLUm+mOSJYf2qJMeTnBqWV05uTEnSUpdyBH4vcHLR+iHgRFXtBE4M65KkDTJWwJPsAPYB9y/avB84Otw+Cty5rpNJki5q3CPwjwMfAb63aNs1VXUWYFhuW9/RJEkXMzLgSW4HzlXVs6vZQZKDSeaSzM3Pz6/mKSRJyxjnCPxW4I4kp4FPAbcl+STwWpLtAMPy3HIPrqojVbWrqnbNzMys09iSpJEBr6r7qmpHVc0CdwGfrar3A8eAA8PdDgCPT2xKSdIF1vI58MPAniSngD3DuiRpg2y9lDtX1VPAU8PtbwC7138kSdI4vBJTkpoy4JLUlAGXpKYMuCQ1ZcAlqSkDLklNGXBJasqAS1JTBlySmjLgktSUAZekpgy4JDVlwCWpKQMuSU0ZcElqyoBLUlMGXJKaMuCS1NQl/ZdqP+xmDz059n1PH943wUkkySNwSWrLgEtSUwZckpoy4JLUlAGXpKYMuCQ1ZcAlqSkDLklNGXBJasqAS1JTBlySmjLgktSUAZekpgy4JDVlwCWpKQMuSU0ZcElqyoBLUlMGXJKaGhnwJG9K8g9JXkjySpLfH7ZfleR4klPD8srJjytJOm+cI/D/AW6rqncCNwF7k9wCHAJOVNVO4MSwLknaICMDXgv+c1i9bPhTwH7g6LD9KHDnJAaUJC1vrHPgSbYkeR44BxyvqmeAa6rqLMCw3DaxKSVJFxgr4FX1f1V1E7ADeHeSG8fdQZKDSeaSzM3Pz69yTEnSUpf0KZSq+hbwFLAXeC3JdoBheW6Fxxypql1VtWtmZmZt00qSvm+cT6HMJLliuP1m4D3Al4BjwIHhbgeAxyc0oyRpGVvHuM924GiSLSwE/5GqeiLJ08AjSe4BXgXeO8E5JUlLjAx4Vb0I3LzM9m8AuycxlCRpNK/ElKSmDLgkNWXAJakpAy5JTRlwSWrKgEtSUwZckpoy4JLUlAGXpKYMuCQ1ZcAlqSkDLklNGXBJasqAS1JTBlySmjLgktSUAZekpgy4JDVlwCWpKQMuSU0ZcElqyoBLUlMGXJKaMuCS1JQBl6SmDLgkNWXAJakpAy5JTW2d9gBvVLOHnlz35zx9eN+6P6ekvjwCl6SmDLgkNWXAJakpAy5JTRlwSWrKgEtSUwZckpoy4JLUlAGXpKYMuCQ1NTLgSa5L8rkkJ5O8kuTeYftVSY4nOTUsr5z8uJKk88Y5An8d+HBVvQO4BfhAkhuAQ8CJqtoJnBjWJUkbZGTAq+psVT033P4OcBK4FtgPHB3udhS4c0IzSpKWcUnnwJPMAjcDzwDXVNVZWIg8sG2FxxxMMpdkbn5+fo3jSpLOGzvgSd4KfBr4UFV9e9zHVdWRqtpVVbtmZmZWM6MkaRljBTzJZSzE+6GqemzY/FqS7cPXtwPnJjOiJGk543wKJcAngJNV9bFFXzoGHBhuHwAeX//xJEkrGed/5LkVuBt4Kcnzw7bfBQ4DjyS5B3gVeO9EJpQkLWtkwKvq74Cs8OXd6zuOJGlcXokpSU0ZcElqyoBLUlMGXJKaMuCS1JQBl6SmDLgkNWXAJakpAy5JTRlwSWrKgEtSUwZckpoy4JLUlAGXpKYMuCQ1ZcAlqSkDLklNGXBJasqAS1JTBlySmjLgktSUAZekpgy4JDVlwCWpKQMuSU0ZcElqyoBLUlMGXJKaMuCS1JQBl6SmDLgkNWXAJakpAy5JTRlwSWrKgEtSUwZckpraOu0BNL7ZQ0+Ofd/Th/dNcBJJm4FH4JLU1MiAJ3kgybkkLy/adlWS40lODcsrJzumJGmpcY7AHwT2Ltl2CDhRVTuBE8O6JGkDjQx4VX0e+OaSzfuBo8Pto8Cd6zuWJGmU1Z4Dv6aqzgIMy20r3THJwSRzSebm5+dXuTtJ0lITfxOzqo5U1a6q2jUzMzPp3UnSD43VBvy1JNsBhuW59RtJkjSO1Qb8GHBguH0AeHx9xpEkjWucjxE+DDwNXJ/kTJJ7gMPAniSngD3DuiRpA428ErOq3rfCl3av8yySpEvglZiS1JQBl6SmDLgkNWXAJakpAy5JTRlwSWrKgEtSUwZckpoy4JLUlAGXpKYMuCQ1ZcAlqSkDLklNGXBJasqAS1JTBlySmjLgktSUAZekpgy4JDVlwCWpKQMuSU0ZcElqyoBLUlMGXJKaMuCS1JQBl6SmDLgkNbV12gNMwuyhJ6c9gpq5lNfM6cP7JjiJND6PwCWpKQMuSU0ZcElq6g15DlyTex/gUs7/boYZpDcyj8AlqSkDLklNGXBJasqAS1JTvompS7IZLpIad4bN8GanFwj10+nvzCNwSWpqTQFPsjfJl5N8Jcmh9RpKkjTaqgOeZAvwZ8AvATcA70tyw3oNJkm6uLUcgb8b+EpV/XNVfRf4FLB/fcaSJI2SqlrdA5NfAfZW1a8P63cDP1dVH1xyv4PAwWH1euDLq5z1auDrq3zsNHSat9Os0GveTrNCr3k7zQprm/enqmpm6ca1fAoly2y74F+DqjoCHFnDfhZ2lsxV1a61Ps9G6TRvp1mh17ydZoVe83aaFSYz71pOoZwBrlu0vgP4t7WNI0ka11oC/gVgZ5K3JflR4C7g2PqMJUkaZdWnUKrq9SQfBP4W2AI8UFWvrNtkF1rzaZgN1mneTrNCr3k7zQq95u00K0xg3lW/iSlJmi6vxJSkpgy4JDXVIuBdLtlPcl2SzyU5meSVJPdOe6ZRkmxJ8sUkT0x7llGSXJHk0SRfGr7HPz/tmS4myW8Nr4OXkzyc5E3Tnum8JA8kOZfk5UXbrkpyPMmpYXnlNGdcbIV5/3B4LbyY5K+SXDHFEb9vuVkXfe23k1SSq9djX5s+4M0u2X8d+HBVvQO4BfjAJp71vHuBk9MeYkx/AvxNVb0deCebeO4k1wK/CeyqqhtZeKP/rulO9QMeBPYu2XYIOFFVO4ETw/pm8SAXznscuLGqfhb4J+C+jR5qBQ9y4awkuQ7YA7y6Xjva9AGn0SX7VXW2qp4bbn+HhcBcO92pVpZkB7APuH/as4yS5MeBXwQ+AVBV362qb011qNG2Am9OshW4nE10nURVfR745pLN+4Gjw+2jwJ0bOdPFLDdvVX2mql4fVv+ehWtRpm6F7y3AHwMfYZkLHlerQ8CvBb66aP0MmziK5yWZBW4GnpnyKBfzcRZeUN+b8hzj+GlgHvjL4ZTP/UneMu2hVlJVXwP+iIWjrbPAf1TVZ6Y71UjXVNVZWDgYAbZNeZ5L8WvAX097iJUkuQP4WlW9sJ7P2yHgY12yv5kkeSvwaeBDVfXtac+znCS3A+eq6tlpzzKmrcC7gD+vqpuB/2Jz/Yr/A4bzx/uBtwE/CbwlyfunO9UbU5KPsnD68qFpz7KcJJcDHwV+b72fu0PAW12yn+QyFuL9UFU9Nu15LuJW4I4kp1k4LXVbkk9Od6SLOgOcqarzv9E8ykLQN6v3AP9SVfNV9b/AY8AvTHmmUV5Lsh1gWJ6b8jwjJTkA3A78am3ei1p+hoV/yF8Yft52AM8l+Ym1PnGHgLe5ZD9JWDhHe7KqPjbteS6mqu6rqh1VNcvC9/SzVbVpjxCr6t+Brya5fti0G/jHKY40yqvALUkuH14Xu9nEb7oOjgEHhtsHgMenOMtISfYCvwPcUVX/Pe15VlJVL1XVtqqaHX7ezgDvGl7Ta7LpAz68SXH+kv2TwCMTvmR/LW4F7mbhaPb54c8vT3uoN5DfAB5K8iJwE/AH0x1nZcNvCo8CzwEvsfCztmku/U7yMPA0cH2SM0nuAQ4De5KcYuHTEoenOeNiK8z7p8CPAceHn7W/mOqQgxVmncy+Nu9vHZKki9n0R+CSpOUZcElqyoBLUlMGXJKaMuCS1JQBl6SmDLgkNfX/1avzBsGF4JsAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "bins =[0,0.5,1,1.5,2,2.5, 3, 3.5, 4, 4.5, 5, 5.5, 6, 6.5,7, 7.5, 8,8.5, 9,9.5, 10, 10.5, 11, 11.5, 12,12.5, 13, 13.5, 14 ]\n", "plt.hist(X_NC, bins, label=\"NC\")" ] }, { "cell_type": "code", "execution_count": 36, "id": "92a3e7e2", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Fitting 106 distributions: 11%|█▍ | 12/106 [00:01<00:11, 8.53it/s]WARNING:root:SKIPPED kstwo distribution (taking more than 30 seconds)\n", "Fitting 106 distributions: 21%|██▋ | 22/106 [00:02<00:07, 10.95it/s]WARNING:root:SKIPPED rv_continuous distribution (taking more than 30 seconds)\n", "Fitting 106 distributions: 23%|██▉ | 24/106 [00:02<00:09, 9.10it/s]WARNING:root:SKIPPED rv_histogram distribution (taking more than 30 seconds)\n", "Fitting 106 distributions: 94%|███████████▎| 100/106 [00:22<00:06, 1.08s/it]WARNING:root:SKIPPED kappa4 distribution (taking more than 30 seconds)\n", "Fitting 106 distributions: 95%|███████████▍| 101/106 [00:31<00:16, 3.27s/it]/home/marija/.conda/envs/pytorch/lib/python3.7/site-packages/scipy/integrate/quadpack.py:880: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.\n", " If increasing the limit yields no improvement it is advised to analyze \n", " the integrand in order to determine the difficulties. If the position of a \n", " local difficulty can be determined (singularity, discontinuity) one will \n", " probably gain from splitting up the interval and calling the integrator \n", " on the subranges. Perhaps a special-purpose integrator should be used.\n", " **opt)\n", "WARNING:root:SKIPPED levy_stable distribution (taking more than 30 seconds)\n", "Fitting 106 distributions: 96%|███████████▌| 102/106 [00:31<00:09, 2.47s/it]WARNING:root:SKIPPED recipinvgauss distribution (taking more than 30 seconds)\n", "Fitting 106 distributions: 97%|███████████▋| 103/106 [00:32<00:05, 1.89s/it]WARNING:root:SKIPPED studentized_range distribution (taking more than 30 seconds)\n", "Fitting 106 distributions: 98%|███████████▊| 104/106 [00:32<00:03, 1.56s/it]WARNING:root:SKIPPED vonmises distribution (taking more than 30 seconds)\n", "Fitting 106 distributions: 99%|███████████▉| 105/106 [00:33<00:01, 1.32s/it]WARNING:root:SKIPPED vonmises_line distribution (taking more than 30 seconds)\n", "Fitting 106 distributions: 100%|████████████| 106/106 [00:33<00:00, 3.13it/s]\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
sumsquare_erroraicbickl_divks_statisticks_pvalue
laplace_asymmetric0.8340691007.300261-788.521547inf0.0942020.121887
burr120.836352888.541745-783.063512inf0.0628400.555837
skewcauchy0.857340801.601692-784.283549inf0.0837100.218161
johnsonsu0.861766884.977607-778.453770inf0.0683450.448399
exponnorm0.8747511006.719459-781.187455inf0.0958520.110524
\n", "
" ], "text/plain": [ " sumsquare_error aic bic kl_div \\\n", "laplace_asymmetric 0.834069 1007.300261 -788.521547 inf \n", "burr12 0.836352 888.541745 -783.063512 inf \n", "skewcauchy 0.857340 801.601692 -784.283549 inf \n", "johnsonsu 0.861766 884.977607 -778.453770 inf \n", "exponnorm 0.874751 1006.719459 -781.187455 inf \n", "\n", " ks_statistic ks_pvalue \n", "laplace_asymmetric 0.094202 0.121887 \n", "burr12 0.062840 0.555837 \n", "skewcauchy 0.083710 0.218161 \n", "johnsonsu 0.068345 0.448399 \n", "exponnorm 0.095852 0.110524 " ] }, "execution_count": 36, "metadata": {}, "output_type": "execute_result" }, { "name": "stderr", "output_type": "stream", "text": [ "/home/marija/.conda/envs/pytorch/lib/python3.7/site-packages/scipy/integrate/quadpack.py:880: IntegrationWarning: The integral is probably divergent, or slowly convergent.\n", " **opt)\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABDyUlEQVR4nO3deXxU1fn48c+5s2aysoZVQ1gNBAhLQBANKotLQVutaxWsRatW+1V+gntbpdoWW9TWWre6FPe9ilq3oLgBYd/3JRAgJGSZZCaznd8fd4gBAyQwYZLL8369YjL33jn3OSE+OTn3nucqrTVCCCFaPiPeAQghhIgNSehCCGERktCFEMIiJKELIYRFSEIXQgiLsMfrxG3bttUZGRnxOn2jVVVVkZiYGO8wjpkV+mGFPoA1+mGFPkDL6kdBQcFerXW7+vbFLaFnZGSwcOHCeJ2+0fLz88nLy4t3GMfMCv2wQh/AGv2wQh+gZfVDKbX1UPtkykUIISxCEroQQliEJHQhhLCIuM2hCyEgGAxSWFiI3++PdyhHJTU1ldWrV8c7jGPWHPvhdrvp0qULDoejwe+RhC5EHBUWFpKcnExGRgZKqXiH02iVlZUkJyfHO4xj1tz6obWmpKSEwsJCunXr1uD3yZSLEHHk9/tp06ZNi0zmoukopWjTpk2j/3I7YkJXSj2rlNqjlFpxiP1KKfWoUmqDUmqZUmpQoyIQ4gQnyVzU52h+LhoyQn8OGH+Y/ecAPaMfU4B/NjoKIYQQx+yICV1r/SVQephDJgIvaNN3QJpSqmOsAhRCCNEwsbgo2hnYXud1YXRb0cEHKqWmYI7iSU9PJz8/Pwanb7zlO8prv87unNqg93i93rjFG0tW6IcV+gBmP1JTU6msrIxrHB07dqSo6Ef/uzbIOeecw4wZMxg0qGXPtIbD4Zj8O7z//vv06NGDPn361Lv/mWeeISEhgcsvv7xB7fn9/kb9rMciodc30VPvY5C01k8CTwIMGTJEx2up7aTpH9R+veWKhsXQkpYGH44V+mGFPoDZD7fb3SzurjjaGJRSJCYmNos+HItY3eXy8ccf43A4GDp06I/2hUIhfvvb3zaqPbfbTU5OToOPj0VCLwS61nndBdgZg3aFOLH8rmF/LTa+3fIjH4P5F8PEiRPZt28fwWCQBx54gIkTJ7JlyxbGjx/PsGHDWLx4Mb169eKFF17A4/Ec8P5f//rXLFiwAJ/Px0UXXcTvf/97ABYsWMAtt9xCVVUVLpeLzz77DI/Hw/Tp08nPz6empoYbb7yR6667rlFxVVVV8fOf/5zCwkLC4TD33HMPbdu25e9//ztvv/02AJ988gn//Oc/eeutt0hKSuLGG2/k008/pVWrVvzxj3/k9ttvZ9u2bTz44INccsklPPfcc7zzzjuEw2FWrFjBbbfdRiAQ4MUXX8TlcjFnzhxat27Nxo0bufHGGykuLsbj8fDUU09RWlrKe++9x9y5c3nggQd48803+eUvf8mIESP4+uuvmTBhApWVlSQlJTF16lQ2bNjA9ddfT3FxMTabjddff53u3bsfwz90bG5bfA+4Knq3y3CgXGt9dH+/CSHixu128/bbb7No0SK++OILbrvtNvY/c3jt2rVMmTKFZcuWkZKSwuOPP/6j98+YMYOFCxeybNky5s6dy7JlywgEAlxyySU88sgjLF26lE8//ZSEhASeeeYZUlNTWbBgAQsWLOCpp55i8+bNjYrro48+olOnTixdupQVK1Ywfvx4zjzzTFavXk1xcTEA//73v5k8eTJgVlTMy8ujoKCA5ORk7r77bj755BPefvttZsyYUXu+FStW8NJLLzF//nzuuusuPB4Pixcv5tRTT+WFF14AYMqUKTz22GMUFBQwc+ZMbrjhBkaMGMGECRP4y1/+wpIlS2qTc1lZGXPnzuW22247oF9XXHEFN954I0uXLuWbb76hY8djv/R4xBG6UuplIA9oq5QqBO4DHABa6yeAOcC5wAagGph8zFEJcSJq4Ei6qWitufPOO/nyyy8xDIMdO3awe/duALp27crIkSMBuPLKK3n00UeZOnXqAe9/7bXXePLJJwmFQhQVFbFq1SqUUnTs2LF2CiIlJQWA//3vfyxbtow33ngDgPLyctavX1/vIppDxZWdnc3UqVOZNm0a559/PqNGjQLgF7/4Bf/5z3+YPHky3377bW0SdjqdjB9v3rCXnZ2Ny+XC4XCQnZ3Ntm3bas83evRokpOTSU5OJjU1lZ/85Ce171m2bBler5dvvvmGiy++uPY9NTU1h/y+XnLJJT/aVllZyY4dO7jwwgsB85dWLBwxoWutLzvCfg3cGJNohBBxM3v2bIqLiykoKMDhcJCRkVG7sOXge6IPfr1582ZmzpzJggULaNWqFZMmTcLv96O1rvd+aq01jz32GOPGjTvquHr16kVBQQFz5szhjjvuYOzYsdx7771MnjyZn/zkJ7jdbi6++GLsdjPNORyO2lgMw8DlctV+HQqFas+3f/uhjotEIqSlpbFkyZIjxg7UW2d9/18+sSYrRYUQgDlKbt++PQ6Hgy+++IKtW38ou71t2za+/fZbAF5++WVOO+20A95bUVFBYmIiqamp7N69mw8//BCAPn36sHPnThYsWACYI9NQKMS4ceP45z//STAYBGDdunVUVVU1Kq6dO3fi8Xi48sormTp1KosWLQKgU6dOdOrUiQceeIBJkybF7hsUlZKSQrdu3Xj99dcBMzkvXboUMC8uN+RumZSUFLp06cI777wDmCP86urqY45NEroQAjDndBcuXMiQIUOYPXv2AbfenXLKKTz//PP079+f0tJSfv3rXx/w3gEDBpCTk0Pfvn255ppraqdnnE4nr776Kr/5zW8YMGAAY8aMwe/3c+2115KVlcWgQYPo168f11133QGj5IbEtXz5cnJzcxk4cCAzZszg7rvvPuA9Xbt2JSsrK9bfJsD8q+GZZ55hwIAB9O3bl3fffReASy+9lL/85S/k5OSwcePGw7bx4osv8uijj9K/f39GjBjBrl27jjku1VRD/yMZMmSIjtcTizLq3rb40HkNeo+VbpVr6f2wQh/A7Ed6ejqnnHJKvEM5rC1btnD++eezYsWPq380t6JW+910003k5OTwy1/+skHHN9d+rF69+kc/H0qpAq31kPqOl2qLQghLGTx4MImJiTz88MPxDuW4k4QuhDisjIyMekfnsbZ8+XJ+8YtfHLDN5XLx/fffN6qdgoKCWIbVokhCF0I0C9nZ2Q2+c0TUTy6KCiGERUhCF0IIi5CELoQQFiEJXYgT3JYtW+jXr1/M2ispKWH06NEkJSVx00031W6vrq7mvPPOo0+fPvTt25fp06fH7JzCJAldCHHUDl4MFAqFcLvd3H///cycOfNHx0+dOpU1a9awePFivv7669oVpSI25C4XIQShUIirr776gPK4WVlZLFy4kLZt27Jw4UKmTp1Kfn4+v/vd79i5cydbtmwhNTWVvn371r5u27YtL730EqeddhobNmw44Bwej4fRo0cD5grSQYMGUVhYGI/uWpYkdCGaiezns5uk3eVXLz/iMWvXruWZZ55h5MiRXHPNNfWWx62roKCAefPmEQqFePjhh2tfJyQkNCimsrIy/vvf/3LLLbc06HjRMDLlIoT4UXncefPmHfb4CRMmHJC8D359OKFQiMsuu4ybb76ZzMzMow9a/IiM0IVoJhoykm4q9ZXHtdvtRCIRgNoyuvsdXBK2vhKxhzJlyhR69uzZ6MexiSOTEboQot7yuBkZGbXL6N98882YnOfuu++mvLycWbNmxaQ9cSBJ6EKIesvj3nfffdxyyy2MGjUKm83WqPYyMjK49dZbee655+jSpQurVq2isLCQGTNmsGrVKgYNGsTAgQN5+umnm6hHJyaZchHiBJeRkcGqVat+tH3UqFGsW7fuR9t/97vfHfY1mPe21yde5bpPFDJCF0IIi5CELoQQFiEJXQghLEISuhBCWIQkdCGEsAhJ6EIIYRGS0IUQP5KRkcHevXvjHUaD5OXlsXDhwniH0SxIQhdCCIs44RN6xvQPaj+EOBFVVVVx3nnnMWDAAPr168err75au8/n8zF+/HieeuopqqqquOaaaxg6dCg5OTm8++67AJx77rksW7YMgJycHP7whz8AcM8999SuBP3zn/9MdnY2AwYMqH2wxVNPPcXQoUMZMGAAP/vZz6iurgZg0qRJvPHGG7UxJCUl1X5dXzsAr7/+Orm5ufTq1YuvvvoKMBdG1X3o9MiRI2vjtCpZKSpEM7G6zylN0u4pa1Yfdv9HH31Ep06d+OADc1BTXl7OtGnT8Hq9XHrppVx11VVcddVV3HnnnZx55pk8++yzlJWVkZuby5dffsnpp5/OV199RUZGBna7na+//hqAefPmceWVV/Lhhx/yzjvv8P333+PxeCgtLQXgpz/9Kb/61a8As8bLM888w29+85tDxnmodsCs4Dh//nzmzJnD73//ez799FOuvfZannvuOWbNmsW6deuoqamhf//+x/S9bO5O+BG6ECe67OxsPv30U6ZNm8ZXX31FamoqABMnTmTy5MlcddVVAPzvf//joYceYuDAgeTl5eH3+yksLGTUqFF8+eWXzJs3j/POOw+v10t1dTVbtmyhd+/efPrpp0yePBmPxwNA69atAVixYgWjRo0iOzub2bNns3LlysPGeah2wPzlADB48ODasgMXX3wx77//PsFgkGeffZZJkybF7HvWXMkIXYhm4kgj6abSq1cvCgoKmDNnDnfccQdjx44FzCmKDz/8kMsvvxylFFpr3nzzTXr37l373srKSlwuFwsXLiQzM5MxY8awd+9ennrqKQYPHgyY9VsOLs8L5tTKO++8w4ABA3juuefIz88HOKBsr9aaQCBw2HYAXC4XADabrfaxeB6PhzFjxvDuu+/y2muvnRAXTmWELsQJbufOnXg8Hq688kqmTp3KokWLAPjDH/5AmzZtuOGGGwAYN24cjz32WG2BrcWLFwPm4+S6du3Ka6+9xvDhwxk1ahQzZ85k1KhRAIwdO5Znn322do58/1RJZWUlHTt2JBgMMnv27Np46pbtfffddwkGg4dt53CuvfZabr75ZoYOHXrAiN6qJKELcYJbvnw5ubm5DBw4kBkzZnD33XfX7ps1axZ+v5/bb7+de+65h2AwSP/+/enXrx/33HNP7XGjRo0iPT0dj8fDqFGjaqdiAMaPH8+ECRMYMmQIAwcOrH149P3338+wYcMYM2YMffr0qW3rV7/6FXPnziU3N5fvv/++9uEZh2rncAYPHkxKSgqTJ0+OyfequVMNKWeplBoPPALYgKe11g8dtD8V+A9wEuY0zkyt9b8P1+aQIUN0vP4EOtQdLVseOu+Q78nPzycvL6+JIjp+rNAPK/QBzH6kp6dzyilNczH0eKisrCQ5OTneYRzSzp07ycvLY82aNRjGocevzbUfq1ev/tHPh1KqQGs9pL7jjzhCV0rZgH8A5wBZwGVKqayDDrsRWKW1HgDkAQ8rpZyND18IIWLjhRdeYNiwYcyYMeOwydxKGnJRNBfYoLXeBKCUegWYCNStiK+BZGVesUgCSoFQjGMVQogG23+75YmkIQm9M7C9zutCYNhBx/wdeA/YCSQDl2itIwc3pJSaAkwBSE9Pr72qfbzdll3/75rDxeP1euMWbyxZoR9W6AOY/UhNTaWysjLeoRy1cDjcouPfr7n2w+/3N+pnvSEJvb77hA6eeB8HLAHOBLoDnyilvtJaVxzwJq2fBJ4Ecw49XvOgkw41h35F3iHfY6V525beDyv0Acx+uN3uZjl321DNde65sZprP9xuNzk5OQ0+viEJvRDoWud1F8yReF2TgYe0eYV1g1JqM9AHmN/gSJpA3Yufh7vgKYQQVtCQKwULgJ5KqW7RC52XYk6v1LUNOAtAKZUO9AY2xTLQo5FRXkRWyeZ4hyGEEMfFERO61joE3AR8DKwGXtNar1RKXa+Uuj562P3ACKXUcuAzYJrWOu61N+/9/t/85avH6VFWGO9QhGjWRowYcch9+fn5nH/++ccxGnG0GrT0X2s9B5hz0LYn6ny9Exgb29COXWt/BQaa65a/i9ZTDrlsWIgT3TfffBPvEEQMWPbmTK01roh5N0u/ks1UfvRRnCMSovlKSkpCa83/+3//j379+pGdnX1AGV2v18tFF11Enz59uOKKK2qX//fr14/77ruPQYMGkZ2dzZo1awCYO3cuAwcOZODAgeTk5FBZWXnI9vdf5K6v/enTp5OVlUX//v2ZOnUqcPjyuic6yxbn0tGCPvvt+ctMkkaPxnC74xSREIf3j+s/b5J2b3zizAYd99Zbb7FkyRKWLl3K3r17GTp0KKeffjpg1m1ZuXIlnTp1YuTIkXz99decdtppALRt25ZFixbx+OOPM3PmTJ5++mlmzpzJP/7xD0aOHInX68Xtdje6/aysLN5++23WrFmDUoqysrIm+f5YiXVH6H4/ANV2F5tSOhLcuZPS556Lb1BCNGPz5s3jsssuw2azkZ6ezhlnnMGCBQsAyM3NpUuXLhiGwcCBA2tL1EL9pWtHjhzJrbfeyqOPPkpZWRl2u73R7aekpOB2u7n22mt56623asvmikOz7Ag94q8BwG9z8q/sifzp6yfY++RTtLryyjhHJkT9GjqSbiqHq+u0vzwtHFiitu6+utunT5/Oeeedx5w5cxg+fDiffvppo9u32+3Mnz+fzz77jFdeeYW///3vfP7554csryusPEIPmAk9YLOzrF0PEgYMQFdX44uWBhVCHOj000/n1VdfJRwOU1xczJdffklubu5RtbVx40ays7OZNm0aQ4YMYc2aNY1u3+v1Ul5ezrnnnsusWbNqHyd3qPK6wsIj9P1TLgHDAYAndyi+pUupLlgE9IxjZEI0P0opLrzwQr799lsGDBiAUoo///nPdOjQofZCZ2PMmjWLL774ApvNRlZWFueccw5Op7NR7VdWVjJx4kT8fj9aa/72t78BZnndiRMnkpuby1lnnVVbXlc0sHxuU2jq8rm+FSvZctFFbEjtzG9G/x/Lx3ko/PUNJAwZTF6Xy+p9j5TPbRms0AdoPuVzS0pKGDRoEFu3bm30e5vrkvnGaq79iHn53JZK15gj9BpbdIQ+aBAA/mXLcYSlEKQQYNYLP/XUU2tvCRQtm2WnXCLRKZegYXbRlpqKq2dPatavp2fZdla16RbP8IRoFjp16sS6deviHYaIEQuP0M0r3/tH6AAJg81Rel+p7yKEsCALJ/ToCN32wx8hnsHmtJMkdCGEFVl3yqXGvG2xJnqXyz1f38PushX8H5BVugWlI2hl2d9nQogTkGUzmvbvvw/dAYaf9za+x7eRDQTbpZEc9HFS5e44RyiEELFl3YQenXIJ2BzYErYSiT4Rb9PJ5oo0mXYRQliNZRN6JHpR9BzH9zg8G2u3z2tTApgVGIUQLVM4HI53CM2SZRP6/pWineylJHnWApDiTGFVF3MhlTzJSIgf/Oc//yE3N5eBAwdy3XXX8f3339O/f3/8fj9VVVX07duXFStWkJ+fz+mnn86FF15IVlYWv/3tb2vrqrz88stkZ2fTr18/pk2bVtt2UlISd911FwMGDGD48OHs3m1Od06aNImbb76ZESNGkJmZWVsS93BldkePHs3ll19OdnY2+fn5nHHGGfz85z+nV69eTJ8+ndmzZ5Obm0t2djYbN27kRGPdi6J+HwBhuyaQsAeF4qacm3iwZgZel0G6r4y21WXs9aTFN1Ahoh6+pGmeCnTbq+8fdv/q1at59dVX+frrr3E4HNxwww2sXbuWCRMmcPfdd+Pz+bjyyivp168f+fn5zJ8/n1WrVnHyySczZswY3nrrLUaMGMG0adMoKCigVatWjB07lnfeeYcLLriAqqoqhg8fzowZM7j99tt56qmnuPvuuwEoKipi3rx5rFmzhgkTJnDRRRcdtszu/PnzWbFiBd26dSM/P5+lS5eyevVqWrduTWZmJtdeey3z58/nkUce4bHHHmPWrFlN8j1triyb0LWvCoBdLhtaaXLsw8mN5OF2/I1NHb303wKZFTsloYsT3meffUZBQQFDhw4FwOfz0b59e+69916GDh2K2+3m0UcfrT0+NzeXzMxMAC666CLmzZuHw+EgLy+Pdu3aAXDFFVfw5ZdfcsEFF+B0OmsfYTd48GA++eST2rYuuOACDMMgKyurduR+qDK7KSkp5Obm0q3bD4sChw4dSseOHQHo3r07Y8eaD07Lzs7miy++aKpvWbNl3YRebSb0rW4baBi66Kd89NVqxp31U7a1e4H+W6BbeRHzO2TFN1Ahoo40km4qWmuuvvpqHnzwwQO279q1C6/XSzAYxO/31xbBOvhRjkqpw5bGdTgcte85VOnd/XHU/Vyfgwtx1X2/YRi1rw3DOOA8JwrLzqFH/NUAbHTbaV3dCaPK/Ic+aUku29qZ96Z3q9gZt/iEaC7OOuss3njjDfbs2QNAaWkpW7duZcqUKdx///1cccUVB8yJz58/n82bNxOJRHjzzTc57bTTGDZsGHPnzmXv3r2Ew2FefvllzjjjjKOKJ5ZlfE801h2hR6dcNrkddC3rU7s9UAJG6jjgv3SrLIxTdEI0H1lZWTzwwAOMHTuWSCSCw+Fg4sSJ2O12Lr/8csLhMCNGjODzzz/HMAxOPfVUpk+fzvLlyzn11FO58MILMQyDBx98kNGjR6O15txzz2XixIlHFU8sy/ieaCxbPnfbVT+nav5y/nSRQbeaG2hbeQp9R3Vi5Vc7CakQI7+7H7e/hAvP/yPBaL0XKZ/bMlihD9B8yuc2Rn5+PjNnzuT9983poeZadraxmms/pHxu1P6VoiGbgzbeHoBm2MRMeg/rgF3bWZ51CYbWnFS5J76BCiFEjFg3oUdXinrCPVHaQXvHJhLcMOJnPfCrENXJWXiTusg8uhCNkJeXVzs6F82PZRP6/pWiqQFz/ryrcxHsWoYnxckGp/kMwn1pvcmoKIpbjELA4e/qECeuo/m5sGxC19Fqi2185vzTSa4lsH0+AFttNgD2tepJt8ptcYlPCAC3201JSYkkdXEArTUlJSW43e5Gvc+yd7mEawL4XWkkBzsQIUS6Yy0UzgduYJvdrANRltqD7PUyhy7ip0uXLhQWFlJcXBzvUI6K3+9vdNJpjppjP9xuN126dGnUeyyb0CPBEKWt+gOw2+HHpsKwfQEAXgPKHF7SSMLmaE9qTSXlruZ3hVtYn8PhOGDlY0uTn59PTk5OvMM4Zlbph2WnXCLBEKWtzemWlXYHuFKgohDKzXvPtznNWi/7WvWkW8WuuMUphBCxYtmETjBCWVpPADbbNXQdZm7f+o35Kfoko7K0XmSUy50uQoiWz5IJXYfD6Igi4ExBE6Hc0JAx0ty5ZR4A25QTgLLU7mR6t8crVCGEiBlrJvSaGoKOJABqbH60Ak4+zdy59WsAqg2F17aXiM1FxxOvho8QwoIsmdAjByV0gGDrPvhUCpRsoB37ANiaUAaA294RIyJPQBFCtGwNSuhKqfFKqbVKqQ1KqemHOCZPKbVEKbVSKTU3tmE2jq6pIeA071rxGwEMHeal393BM+sGUBVyMMwwi/xsilberEjtRaeqkniFK4QQMXHEhK6UsgH/AM4BsoDLlFJZBx2TBjwOTNBa9wUujn2oDadragg4zITuM8IMLF/G3m1bqAnBirIODDNWA7CVRNARKlK7kSl3ugghWriGjNBzgQ1a601a6wDwCnBwXczLgbe01tsAtNZxXa0T8dcQjI7QayJVDC37oarjsrIO5CozofvCqYR1ERHDQY8ab1xiFUKIWGnIwqLOQN3bQAqBYQcd0wtwKKXygWTgEa31Cwc3pJSaAkwBSE9PJz8//yhCPjL75i21I/SuvmUYOkRqRg98JXuoqASnr5q7skvx21PwfVdEYrAzPRy2w8bj9XqbLN7jyQr9sEIfwBr9sEIfwDr9aEhCV/VsO7jwhB0YDJwFJADfKqW+01qvO+BNWj8JPAlmPfSmqmld7VnAR87tRIKFGN5thJSNi2+7gzVfz2XeKy+wdF8HCtZu4KNILmc7fOQEIRBuc9ga21aqwd3S+2GFPoA1+mGFPoB1+tGQKZdCoGud112Ag1fiFAIfaa2rtNZ7gS+BAbEJsfEi1ZUEHckEfV8CUJCaQ2r7dPqNHoNhKDZ52zA0vBKAlUkJAAQTTiIcLeglhBAtUUMS+gKgp1Kqm1LKCVwKvHfQMe8Co5RSdqWUB3NKZnVsQ224sLeMGkciOmw+RXxR6kAAEtNa0aNvbzSKdpXmRdBduj1GyEvAlUrpsg3xClkIIY7ZERO61joE3AR8jJmkX9Nar1RKXa+Uuj56zGrgI2AZMB94Wmu9ounCPjxfZTE1Djug8RkJhKLL/AH6n/dzAIrKE0jTFehgG1RoKwCbFqyrrzkhhGgRGlRtUWs9B5hz0LYnDnr9F+AvsQvt6FVX7CVg6wiA1550wL6TBgzB7QjjDboYHVjA266zqDSK8NCXHRvlXnQhRMtlyZWi3opywsp8YlHFQQldGQbViSkA9PaZUyw73OYti+UVzaseshBCNIY1E3plDTpiJmmvPfFH+5clmGV17VXVgGZVogt0BJ9qTzgYOZ6hCiFEzFgyoVf5IuhIJQBeW9KP9n/rGojdCFMecNMnvJktrh54qnejDTt7Nu493uEKIURMWDKh+/zqh4Rezwg9pBy4Pebt9SN8iwiEOmKv2QzA5vlyYVQI0TJZMqHXBO0/TLnUM0IH2JnQCYAO1bsAG9VqBwDb1+4+LjEKIUSsWTKhB0MutN4/Qq8/oX/lHgKArxqcuoZdrgoAykst+5hVIYTFWTahE6kCoMrmqfeYTfaupLhqCGkbp/kL2OhxYwv5CepkqsplxagQouWxZkLXdiACykbYOPSIO+Ax59f7+dexydODlEpzgdHuzRXHI0whhIgpSyb0kI52S9kOe9zqBPMh0u6qSnY6e5NUuQmA7auLmjQ+IYRoCpacMA5Hi0FGjB8Sesb0D350XL5zKH2MZfiCDjoGS9G2QgC2rSwC+tfbdt12tjx03hFjaezxQghxtCw5Qg9jLg4KHWa6BcCrEklJNJ8leoZvAUbbagAqSwzCIVlgJIRoWSyX0MORMBGCAAQM1xGP3+kxb1/sXF1IYp8MPNW70NrG3kJ5gpEQomWxXELfV1GOjt7hUm1LOOLxX7hzUWhqfNC29yBSys0FRrs2lTdpnEIIEWuWS+jFpfvQ2hxdV9iTj3j8WiODNgk+NApVVU5qhZnQt62L62NRhRCi0SyX0EtLKyC67H+fPbUB71B4E9MA2LZyAfbwNgCK1pc2UYRCCNE0LJfQy/ZV1i77L3a0btB7FiSYT8vbsrUE0jW2kI9glQ3vPllgJIRoOayX0HfvBSIY2t7AETrMsw8k1enHHzIIdEonpWILIPPoQoiWxXIJvWqvOVXiiNioOMSy/4OFsGMkmo+pq9EhUivMBUa7NpU1SYxCCNEULJfQa0rLAHCGFSHV8HVTaxJ6ALCnvJLU6J0uW9cXxzw+IYRoKpZL6KEK84KoKxQBpRr8vv85R+A2gpT5FSq4A3SEskK/PMFICNFiWC6h62ofAO5wsFHvK1MpZHYwR/S7OiaQWL0LIoo92ypjHqMQQjQF6yX0gB+AhFDj71DpOTQXgBKPi9Rycx69aENZzGITQoimZLmETsgcoSdGPzdGxrhJOI0QlThIqFwDwPa1JTENTwghmoqlErqOaAibBbZStL/R77e3OZnM9ua3pEabJXR3bigjEpZ5dCFE82ephF7jD0F0UVGSoY+qjV5DhwJQnGzD7SsmEkAKdQkhWgRLJfSq0kogDNhwO4+u1HvGuEk4VJhyj5uk8uUA7FhXFrMYhRCiqVgqoRfv2Q2AgROb03lUbTjadSMz3fw6UrMOgMJ1UtdFCNH8WSqhlxSbFRINHCjX0SV0gF6DBwPgdZjPFt2xfp85Py+EEM2YpRJ6ecleAGzahnEMCb3b+GuwE6YiwY6q2UrYr/FLWRchRDNnqYReUVYGgC1ioFzuo27H0T6Tbm3NEbmzYgEA1VIeXQjRzFkqofsqzSkSewSU++gTOkDWyBEA1LADrTXePXLrohCiebNUQq/xmrcXOsIaI6FhlRYPpdt51+GMhPDbw+hwMd7iiMyjCyGatQYldKXUeKXUWqXUBqXU9MMcN1QpFVZKXRS7EBsuWG0uKnKEIyj3sSV0W2o6mUnm1+GaxRCwUVpUdawhCiFEkzliQldK2YB/AOcAWcBlSqmsQxz3J+DjWAfZUGG/udzfGQod8wgdIGuEOe0SCaxD6wiFa/Ydc5tCCNFUGjJCzwU2aK03aa0DwCvAxHqO+w3wJhC3y4eRGrMglzMYRCUkHnN7J116M4k1ATRBIqGtbF2595jbFEKIptKQ5ZSdge11XhcCw+oeoJTqDFwInAkMPVRDSqkpwBSA9PR08vPzGxnu4UUCPhTgCtawtbiE27JDjXp/ffF0CfhZ63ISDqxm+9qTuW14qPbXYEPirxtDrPt7NLxeb7OI41hYoQ9gjX5YoQ9gnX40JKHX95SIg68OzgKmaa3D6jAPldBaPwk8CTBkyBCdl5fXsCgbqOCpJwFwB2rIzMrmp8sbt/x/yxU/jmf9271Zu6uISGA9OnQ2byxxstUROeTxB5s0/YPDtn+85efnE+vv+/FmhT6ANfphhT6AdfrRkCmXQqBrndddgJ0HHTMEeEUptQW4CHhcKXVBLAJslEgAAHfAh/Ikx6TJdqPG0crrA8KEA2vpFrLUjUFCCAtpSHZaAPRUSnVTSjmBS4H36h6gte6mtc7QWmcAbwA3aK3fiXWwh6MjEYiYc+gJNVWoxNSYtOvu35+upeZTi8I1y+nWuFkcIYQ4bo6Y0LXWIeAmzLtXVgOvaa1XKqWuV0pd39QBNpSvshJzJsiJI1yDkZAQk3Zd3bvTqSaIPRRGh3fRumYfSbLGSAjRDDVokllrPQeYc9C2Jw5x7KRjD6vxKvaaxVaUcqHQ0aX/x/48UGWzkdQ3my47NrGlXRrhmmV0C57Oclf4mNsWQohYssyEcHmJWeLWnBUCw+2KWdsJgwbRtcQsKxAOrCEzEIhZ20IIESuWSej7iosBMKJ/dChX7BJ64rBckmuCJARDQICTq9ejDrrPJ2P6B7UfQggRD5ZJ6GWl5sOcDW0DYpvQE3Jy0DYbPYuiD4yuWUXHsGW+dUIIi7BMVqrYZ0652LV5H7wRw4RuJCQQ7JZBx7IqIsqGDu8mt7owZu0LIUQsWCahV1eaF0Xt0WuVx1o+92CBXr2waY1hmLe4nORdgdZSfVEI0XxYJqHXLZ0LsR2hg5nQAXrs2IVGoWs2U7h8fUzPIYQQx8IyCT3oM0vnuoLmED3WI/RgZibK4aDHbi++xFZAhG//81xMzyGEEMfCMgk95DNL57pCYQyPG2WzxfYETicJAwZgAN5E85dH4bY1tTXYhRAi3iyT0CMBPwDOUABb8rGXzq2PJzcXgDa7S8GejtYBFr749yY5lxBCNJZlEroOmgndFajBSI5NYa6DeYaZVYP7bq9hb9vWACye971ZR0YIIeLMMgmdsFmYyx3wY0tNa5JTJAwcQMCwk7FvHxvbFYORjC9Qw8ZPXmmS8wkhRGNYIqGHQyGIBAGFK+DD1qpVk5zHcLlY0+okDDSJuytR7gEAfPvm69DIWxhlZakQItYskdD93mgRLuXGGarGlta6yc61qH1vAAZujLAz3Q4qgT3lQUb4FzXZOYUQoiEsldCVcmMP+TBSYlMLvT7fdTSfj51buJ2iThuwuwcDcFrZd40epQshRCxZIqFXlZmVEFFujEgAW0pKk51ra3IHijxtSAtUcVqah1BCT1BufH6DYTVLm+y8QghxJJZI6JWlZQAo5UABttSmS+goVTtKz9vqZkP7FdhdgwA4o+ybpjuvEEIcgTUSesk+AFS0dK7RhCN0gO869AXA+e1SQj33YnPngHLi8ymG+Jc16bmFEOJQLJHQy/ZFS+dGu2Nrwjl0gJVtulHpSCCwaRMTOmezz1OKLTpKP7PsKynaJYSIC0sk9Moys3SuLWKWzm3SKRcgbNhYkH4KAEM2RNjcYQl292AMw47PZ7Dx45ea9PxCCFEfSyT06gqzdK4zumCzKS+K7rd/Hr1m7jy6D22LVg6UayQA815/iUgo1OQxCCFEXZZI6P4qs3SuM2hmdCO1aadcAAra9wGHA9+ixfw0YzSFaWuwuwbgshmUeBUrX3u0yWMQQoi6LJHQg9VVALijpXNtTVTLpa5qh5vE3FyIROi4eDsr0tahlJ0az+kAfPPhJ7UlfYUQ4nho8Qk9Y/oHlJWbC4ucoRA1djvd7v44pkvqM6Z/wPId5T9qM3n8OADK332PdZE2VLhKcNlySHQqvAEbf7z1hga3L2UAhBDHqsUndABHxCzM5QoG8Dli+6Siw0kZPx7ldlM9fz6ti05iebv5KKXYl3gmAKn79tAptPu4xSOEOLFZIqEbkSAArmANVc6E43ZeW3IyyWPGAHDW9iUstbsIGH5SGEBaoiKsDS4rfee4xSOEOLG1+IRujwQxdBiw4Qz6qLR7juv50y68AICzty2kunwga9oVALA2eQIOFSZQFWGUb/5xjUkIcWJq8QndEzYvPCojEUfIR5mzaZ5WdMjzDxuGvWNHOlaX0rd4J4vd1WgidA6cjGpl3j45ouR7bDp4XOMSQpx4WnxCT4wmdFQi9mA1FY6k43p+ZbOROnECAGO2LaTYO5CtaauxazuL3RNJcfqpDtq5uuKd4xqXEOLE0+IT+g8jdA+OkI9Kx/GbQ98v7YILABi1YymuGgcLU4oA6FfTmvWtswFI3reHvasXHPfYhBAnjpaf0EM+AJRKxBb24XUe3zl0AGdGBitbZ5AQDjBqx1K2VPdle8o6XBEnO+2n0jolSEQbfDzrD7KCVAjRZFp8Qk+unXJxYYuE8MZhhA7w0cnmA6Qv3PgVOpjKNyk7ARjsS+H51EtJtAfYVaYpePp3cYlPCGF9LT6hp0RH6AZmYa54JfS5XXLY606hW0URg/esZVt1NltT1+CMODmZBLa27QHA1/mLKF23OC4xCiGsrUEJXSk1Xim1Vim1QSk1vZ79VyillkU/vlFKDYh9qPVLCpvL/m3RkrXe43gfel1Bm513u48C4KL1X0DEwzeJZhXInOo03nWdQ9+uNsLaYM7Me+WuFyFEzB0xoSulbMA/gHOALOAypVTWQYdtBs7QWvcH7geejHWgh+KJmAndGTITpNdx/ObQD16yPydjONV2FwP3bqRHWSGFvn5sTl2DI+JkhC4nb+pfSXEG2F2umVL2SoPaFEKIhmrICD0X2KC13qS1DgCvABPrHqC1/kZrvS/68jugS2zDPDR32Jxy8dT4AeJyl8t+1Y4E5mQMB+Ci9fmgHXzlqSZChAGVHdhV7eHcqy5DobGXecnzfRe3WIUQ1mNvwDGdge11XhcCww5z/C+BD+vboZSaAkwBSE9PJz8/v2FRHoLWGmfETOhJ0cqGvxjgJJh4dHeSHCqe27JDpCeYn4/EdfIIIo99xek7l2LvMobqtB5s0UvJ3JXDq0/OJXdcBv26J7B8o5+hexcwLC+TkKd1o+I5Fl6vt0naPZ6s0AewRj+s0AewTj8aktBVPdvqfcaaUmo0ZkI/rb79WusniU7HDBkyROfl5TUsykPwe70seiICOEmqNmuiz9yQRMhoSLd+bMsV9cczafoH3JYd4uHlDWm3Dbd1zuHs7QVUvZfP3wZdQoK9Lb+0e0kqa08wnMDZv3uBHddNoLTajefr1/lr+vVoZWtwPMciPz+fY/2+x5sV+gDW6IcV+gDW6UdDplwKga51XncBdh58kFKqP/A0MFFrXRKb8A6vqiz6cGjDg7umnIDNdtTJPJZe7n02IWVw9raFZJbtwBdqy9et1wOw/L978AcNXm07kQRbgGqfwY37Zsc5YiGEFTQkoS8AeiqluimlnMClwHt1D1BKnQS8BfxCa70u9mHWr7o8Om1vJOIKlONzuI/XqQ9rZ1I73u82EgPNtSvfB61ZHOhOWUoR7pok/v38+6yxZbIyPQeDCKq8igurPo532EKIFu6ICV1rHQJuAj4GVgOvaa1XKqWuV0pdHz3sXqAN8LhSaolSamGTRVxH7QhdJeIMVFDlbB4JHWB2nzFUOhLIKV7P0N2rAQcjLskkQhi9LI2TjD184DqdULtWAJxcvIFBNSviG7QQokVr0H3oWus5WuteWuvuWusZ0W1PaK2fiH59rda6ldZ6YPRjSFMGvV/F3ujMjuHCHvJR4Ti+lRYPx+v08FJvs1b6tSvexxYJM3rwqURy9qAwOMebglP5+UfSZaSlhglrg7N2f8HJwR1xjlwI0VK16JWiZbv3AhAxbCigvBkldID3M0ewM7ENJ3n3cN7mbwC4fvJFVKYUkxJoxdm6GIC/tppMa48ff9jOpXvepHV43+GaFUKIerXohF5aHH28m2HedNOcRugAIcPOU/1+AsDkVXMIbN2Ky+nivGtzCKkgfct60cuxmSqVyL/aXUmqq4bqgIMpe/6DK+KPc/RCiJamRSf0in3mlIuNMBC/Zf+H813HfnzRJQd3OMjOO+5Eh8Nk9+rFvLRtAIwvzaS1vYhdRjteS59Ioj1Ald/Ob/c+SzggSV0I0XAtMqHvXxq/Kzrl4oo+JDpehbmO5PH+F1LqSsa3aBGlz78AwIJIJzYkbccVTuCCqkSctnLW2bqR3yEPlxGiukrx/h2/IBwIxDl6IURL0SIT+n7OsDmCTYxWXDyedVwaw+v0MCvnYgCKZ82iZsMGUPC+rRWlrlLa+NI5L1wGBPjOMYAFHYfjNEJsKPQx564riQSlkJcQ4shabEJXOoI9bI7MkwOVAFQ2wymX/RZ0yCL1Zz9FBwIU3nwLnqCPoLLxlltRY/jpUdaHUc61QJi5zqEs7piLywixbls1H9x5hUy/CCGOqMUmdHfEj0KDcpPsqwCa75TLful33ImrZw8CmzZx+8KXMHSEfcrNnOQKNBGG78llcMJSIMLnzmH87LqrcUaT+tu3XUqwsizeXRBCNGMtNqF7wj88ei7RVw403ymX/WxJiXR5/HFsqakM272aq1Z9BMAGlcynKeYtjKOLTqWvZxGg6Zh3GT//zXV47EG27gnx+tQr8ZUUxbEHQojmrOUm9NAPD4d2VZmFuaqaydL/w3F27UrnR2YRVgaXrP+cM7cVALDESGFe8m4UBuOKTqVn0nwiOkL6iAu5dNodpDiDFJXBy1OvYd/6RfHthBCiWWqxCT0l+qQibSRgrzbn0ivj8IDoo5E4fDj/yp4AwK2LX+W0HUsB+NZIoSBxDzZt5yc7TuP+l2cRjARp1T+PS3//EG09QfZV23jpvrvYlv96PLsghGiGWmxCTw2bF0JDdju2QPTxc818Dr2u/2aexku9z8amI0xbOJtTd64ABZ/bk1niKcGm7bT+Mpt7n5uJL+QjOTOHy2Y+S2Y78IdtvPnPf7P03/ejdb2VjIUQJ6CWm9AjZkLHBobWVDoSCNoc8Q2qkV7sM47Xeo7GriPcseDF2qT+icPDd4klGNjoPH8odz/xMHuq9+Bs05mJf3uDIVlpRDD49KPv+fCOywhUSqkAIUQLTuhJYXMO3anNVaKrWmfEMZqjpBT/zjqXN7ufjkOHuWv+80zYOA8UfOXw0PucVigMMleM4KGHn2PRzsUYDjdn3Psi55w7BLsKs3qzl9k3X87e5fPi3RshRJy12ISeEF1MlFhjLrpZ2aZbPMM5ekrxdL+f8GKfsdjQ/Hr5O1y37B0MHeHsiTmceuVJRIww3XcM4a1ZC3ll8etoIOvq33HltFtpkxCktNrG7Bl/pOCJO9DhcLx7JISIk/g/3ucouaOrRFO85gXRWCT0jOkfHHMbR9W+UrzUZyxFiW34v0WvccGmeXT17mHgb6socyfTwaO5ImKnY0V3tv67jDtHPMC0C35Dm5wxPNK+mBv2/YeKCoP8L5azYekFjP9/M0jN7F/vuW/LDpHXpL0UQsRLix2hO6KrRNPKqwkYdtandT3CO5q/L7oO5q6RUyh3ehi8Zx3/+OKv5OxZxy674pp7R+PqHCYpkEbn/FO596+P8dW2eewzUpnR5kZ2dcjEYw9RWKp5/q5pLHx8GpGg1IER4kTSIhO6ocMYkSCgSPRXsq5VV4K2FvvHxgGWt+3OjaNvY2nb7rSuqeSP3zzJlOXvkuAIcc2dZ9F7TBsAem8byeezNnFSq09RNi+vJ4zj6r88Rq/OdoIRG3PnruQ/109k5zf/jXOPhBDHS4tM6AnRVaIoD66aCla01PnzQyhJSOXOkdfx/CnjCSuDCzd+xcbzf0LV3HzO/tkAJt6ag0oO0ra6CxdvPo+xtnV4UgpwdcjkJw+/zYWXn0OKM0ix18bLj/yLD27/ORmhwnh3SwjRxFpkQq+yJ6HaTMGVciWuQDmrWlsroQNElMErvc/mljNuZl1aF0JFRRTecCPbr7uedpTwqwfOpvvpaQD03zOSSUXZ3Pz4fcwr/JrMiTcy6R+zye3fHpuKsGZrNRcUvsvNZc/h8MktjkJYVYtM6A5qcEWSMHBihKpZ3frkeIfUZDamdeH/zriZ9DvvxPB48M6dy+aJF1B8392ceXo6LybVUOSsJDGYStaKM/nir1u59bn7WFq9hVF3PcvkP/yO3l0chLWB3ucj+NmrfDPzRmr27Y5314QQMdYiE3qaqiLsWEmrsvVsS+mAt4Us+T9aEWXQ+qpf0P1/H9PqiivAbqf8nXfYOG4cVy95iXmBPfw3sRqVEqK1ryPdvzuDD/68ilufvY+ViSHOm/kWyzsNIs0TIBCx8e2CrTx14yS+enAK3sL18e6eECJGWuSVxLC207voTU5asZv3u50a73COG3vbtnS4525aT7qavX//O+Xvf8CZhYs4s3ARi9r15Nw7bmKt6kTBJ1tpW92FtvO78Pmyrbzc/SO+b2sn3/1r/tjqe+yblrF9n435S3ZSsPQWTslMJedn19B+8Jh4d1EIcQxa5Ai9lBSKilsDsKJNZpyjOf6cXbvS6U9/osf/PuadzNPw2xwMKl7PrltvofUff8VPO65k8BkJkBSklT+dPivz+NXOQYxzrOCz5DADH5zFZTdNomdHO2GtWLGxkhf//AgvX3ceq1/6E8Gq8nh3UQhxFFrkCF3pCFmlWwBYacELog3l6NyZf/W/gNl9xnLW9gJurl5BYONGyp55ilSeYlyvPpSf+lOWVLWGvR767xkJe+ClDd9SkbGN/pdfxER7GtvefZWV6/ays8zGzne/wvV+Pr27t6bv+IvpeOoElNEif+8LccJpkQm9U1UJiUE/uxPS2OtJi3c4ced1eni3+yhmPfggvoICyt/7LxUff0xw3Ro86/7ICGBTu36syzgDe1Im7aq60m5lV7wrIzybspFQRhZZo9uSsWY1hQtXs7vSxrJ15Sxb9zSp/3qC3qd0pfeYi2k36GxJ7kI0Yy0yoe9Iascl5/6e9tUnzi14DSlLoJTCM2QIniFD6HD3XXjnfU3lJ5/gzc8ns3gFmcUrCBt29qYPYme309jnPJnOFT1hGVQtg/xEJ5W929OpVTWdtq2lfFM55TV25i8pYv6SR0lx/ZUe3dPpPnIMnUf9FJvrh4vRh4pvy0PnHbE/dY9pyPbnxice8XtxqHYaqiHf76NptzHnbYr2Rfwcj3/bFpnQwRyVWv3ulmOhnE6SzxxN8pmj0eEw51/3OEP2rOWcyrWk71pAetF8gvYE9rbpT3G7AZS0PsUcuW80SyjsMoazu8smnAnbaFO+BbWngooaO4tWlbBo1Ss4npnNSekuMvr25aQR54KOgJLRuxDx1GITumg4ZbOxsm0mK9tm0jZ7DL8+93SqFyygeuFCkhYspOOqpwljUJbWi5LWWZS26gOJHelalgVlWQCEEmuoti9HsR63dw9BP2wsCrGxaCl8upTb7UFsCXZ2u9uz1NWHRfZT0MoW554LcWKRhH4CsrdqRcrYsaSMHQtApKoK34qV+JYuxb9yJf5Vz1O5p5KytJ6UpfagLLUH1YkdaB0eAgyBBNCuSsKBdQTD6yBUjC8EVEJiZTEjKCbP+AK3W/PVH9+mQ+9s2vcfRUr3gTIHL0QTkoQuMBITSRyWS+Kw3Npt4cpKataupWb9evzr1uFd/zl7ikOUqbZUJJ+MN7krAfdg7AxGa42OlBAJFhIJ7UCHCglEqghUw/ylu2DpLnjtE+xGmHYpig4d2jDZ62Czowur7N3j2HMhrEUSuqiXLTm59gLrfj2BcFkZgS1bCGzbRvmGdRRtLaGkNMw+XwIhdzrVCVlEbE50xEsktItIeBc6tJtIeA+hiI+iMigq20cSkM0eslnE3y55A7cBSS6DqcEI1bYkdtnTKVuQSmK3PthTU+P1bRCiRZGELhrFlpZGwsCBJAwcSCpwUnR7xvQPcIRDtPXtpquvik5hL20jYXqmtMbn70ggnEjQ5iSsy9DhEiLhEnR4Lzq8jwhBqiNQ7TOftmQLVdK5ppJnZj6E0jbs2sm0iIE9orBHNEVrI7z/xJMkeVykpiTTul07PGltsCUnYSQmYiQmMmT3anw2F367i5rNmzE8iRieBAy3G+VoWc+eFaKhJKGLmAna7BQltaUoqW3ttrq3Z2mt8Zb52bFpK394Pp924S6khDRdXAahah/hYICIrkFHKsyPcDlaBQgqH8GDp961D6qAqmIo2oxSbhRubDgwtJ2ztQ17ROEIwxu3LMIRjuAMhXCEQjh1CAcRHHaF3Q42u8LhsuFw2nElOHC4XNy+toQam4OAYSdgcxAwHARtdoKG+RGw2Sl7J4jhcqEcDpTT+cOHw/HDx8GvHQ6U3Q52+V9PxF6DfqqUUuOBRwAb8LTW+qGD9qvo/nOBamCS1npRjGMVLZxSiuRWCfQZ3IdvXt9Yu/2A+82nfYATTSIBnry0N7u3b2Ff4TZWr9qCJ1CDMxTCHg5BOICOBCDiB+1Hax8aH5FoOzUK86f1sDfa2EA5UcoJyg44UT4n+O2oDp2xYSNRGyRpsGmFoRWGxvwIRfjvi/MxIhpbJIItojF0BFskbL4OhzF0GFvE/KwiYZQOYUTCKG1us6sI3xiPogwwDFA2xb98AcJABFh2wYsYNgPDMLDZDQy7DZvN/GzYDJTDjrLZUTYb2G0omx1sxoHbDBvKbgPDZu6r81rZDLDZzc/7Xyuj9jizLRsYB31WRvQ9Bs5Vq6lyucz3GwpsNlCq9jgMZV4IN4wDtitDmdsMA6X2fx1tY/9rpersq3Nc9LNS9Rx7gjtiQldK2YB/AGOAQmCBUuo9rfWqOoedgznF2hMYBvwz+lmIxlEQQBHAxdCcvpDTF/jxwqK8vLza1yFfNSXb1rFr60ZKirZTvm8fSzbsxBEOYguH8RiKSDBEJBSGcNh8kHYkhCIM2ofWvtq2dJ1QDvm47f2/LBraIWygbJilk2woZQdc5mtl1G5HGSiM2u3rUdF90SSGgQorCBvRdlX0DAql93+lgXC0b/uPAKV17X7zWG3u0z9EaR6/f7v52Xyfrm1n/3b2t6E16z/dGG1DR88FEAGtDmjvh2PqHvdDPLXfMb3/67r76sRwwPvr7Neg6+T0Hx2//z8H5H2z3YjWvPXwi6j9O+sec4jfE0qZ39MfH6uj/176gDfPCITQGj7o3BuI38KiXGCD1noTgFLqFWAiUDehTwRe0Fpr4DulVJpSqqPWuijmEQtxEHuCh/TeA0nvPbB22/81YFVeqLqSmn27CVQUE6jcZ354K/BXVfLKN2tx6aA5ytYhjEgEpTUqos3/XyMah02hw5pIRKMjRO/20Wby2v+BBkKgQ7Xn1brecMx9x/i9iBl10Gdx7KLZtmek6RZEKn24ny5AKXURMF5rfW309S+AYVrrm+oc8z7wkNZ6XvT1Z8A0rfXCg9qaAkyJvuwNrI1VR46DtsDeeAcRA1bohxX6ANbohxX6AC2rHydrrdvVt6MhI/T6fkcf/FugIcegtX4SeLIB52x2lFILtdZDjnxk82aFflihD2CNflihD2CdfjRk2V4h0LXO6y7AzqM4RgghRBNqSEJfAPRUSnVTSjmBS4H3DjrmPeAqZRoOlMv8uRBCHF9HnHLRWoeUUjcBH2Ne239Wa71SKXV9dP8TwBzMWxY3YN62OLnpQo6bFjlVVA8r9MMKfQBr9MMKfQCL9OOIF0WFEEK0DFL6TgghLEISuhBCWIQk9CNQSnVVSn2hlFqtlFqplLol3jEdLaWUTSm1OLpuoEWKLlp7Qym1Jvpvcmq8Y2ospdT/RX+WViilXlZKueMdU0MopZ5VSu1RSq2os621UuoTpdT66OdW8YzxSA7Rh79Ef56WKaXeVkqlxTHEYyIJ/chCwG1a61OA4cCNSqmsOMd0tG4BVsc7iGP0CPCR1roPMIAW1h+lVGfgZmCI1rof5o0Gl8Y3qgZ7Dhh/0LbpwGda657AZ9HXzdlz/LgPnwD9tNb9gXXAHcc7qFiRhH4EWuui/YXGtNaVmAmkc3yjajylVBfMAhJPxzuWo6WUSgFOB54B0FoHtNZlcQ3q6NiBBGUWdfHQQtZsaK2/BEoP2jwReD769fPABcczpsaqrw9a6/9pXVub4TvMdTQtkiT0RlBKZQA5wPdxDuVozAJuh9qChC1RJlAM/Ds6dfS0Uiox3kE1htZ6BzAT2AYUYa7Z+F98ozom6fvXnEQ/t49zPMfqGuDDeAdxtCShN5BSKgl4E/it1roi3vE0hlLqfGCP1rog3rEcIzswCPin1joHsyJ6c/8T/wDROeaJQDegE5ColLoyvlEJAKXUXZhTrLPjHcvRkoTeAEopB2Yyn621five8RyFkcAEpdQW4BXgTKXUf+Ib0lEpBAq11vv/QnoDM8G3JGcDm7XWxVrrIPAWMCLOMR2L3UqpjgDRz3viHM9RUUpdDZwPXKFb8OIcSehHEH14xzPAaq31X+Mdz9HQWt+hte6itc7AvAD3uda6xY0Ktda7gO1Kqd7RTWdxYBnnlmAbMFwp5Yn+bJ1FC7uwe5D3gKujX18NvBvHWI5K9AE+04AJWuvqeMdzLCShH9lI4BeYo9ol0Y9z4x3UCew3wGyl1DJgIPDH+IbTONG/Lt4AFgHLMf8fbBHLzpVSLwPfAr2VUoVKqV8CDwFjlFLrMR+C89Dh2oi3Q/Th70Ay8En0/+8n4hrkMZCl/0IIYREyQhdCCIuQhC6EEBYhCV0IISxCEroQQliEJHQhhLAISehCCGERktCFEMIi/j+ogkBp+ErkAQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "f = Fitter(X_NC)\n", "f.fit()\n", "f.summary()" ] }, { "cell_type": "code", "execution_count": 37, "id": "b6db22d9", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Fitting 3 distributions: 100%|██████████████████| 3/3 [00:02<00:00, 1.32it/s]\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
sumsquare_erroraicbickl_divks_statisticks_pvalue
lognorm0.919659921.257538-773.477696inf0.0812092.480260e-01
gamma1.778432656.279200-671.917107inf0.2854861.331359e-11
norm1.8905411676.350736-667.539883inf0.1901642.375160e-05
\n", "
" ], "text/plain": [ " sumsquare_error aic bic kl_div ks_statistic \\\n", "lognorm 0.919659 921.257538 -773.477696 inf 0.081209 \n", "gamma 1.778432 656.279200 -671.917107 inf 0.285486 \n", "norm 1.890541 1676.350736 -667.539883 inf 0.190164 \n", "\n", " ks_pvalue \n", "lognorm 2.480260e-01 \n", "gamma 1.331359e-11 \n", "norm 2.375160e-05 " ] }, "execution_count": 37, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA6P0lEQVR4nO3dd3xUVdrA8d+Zmt4IJCGUBKQbQIgUCwRFRRGwgIqiIKus6+K6lt1F14Lruyu+ui7W5bV32UUsILiKKAoCUlyKEDoBQiAE0sskU877x0xCCAmZJBOGGZ7v53M/c+/cc899DsYnJ+fee67SWiOEECLwGfwdgBBCCN+QhC6EEEFCEroQQgQJSehCCBEkJKELIUSQMPnrxPHx8TolJcVfp2+ysrIywsPD/R1GiwVDO4KhDRAc7QiGNkBgtWP9+vVHtdZt69vnt4SekpLCunXr/HX6Jlu2bBkZGRn+DqPFgqEdwdAGCI52BEMbILDaoZTa19A+GXIRQoggIQldCCGChCR0IYQIEn4bQxdCnF3sdjvZ2dnYbDZ/h3KS6OhoMjMz/R3GCUJCQujQoQNms9nrYyShCyFOi+zsbCIjI0lJSUEp5e9wTlBSUkJkZKS/w6ihtebYsWNkZ2eTmprq9XEy5CKEOC1sNhtt2rQ545L5mUgpRZs2bZr810yjCV0p9aZS6ohS6pcG9iul1AtKqV1KqU1KqQFNikAIcdaQZO695vxbedNDfxsYdYr9VwLdPMs04J9NjkIIIUSLNZrQtdY/APmnKDIOeFe7rQZilFJJvgpQCCF8JSIiwt8htCpfXBRNBg7U2s72fHeobkGl1DTcvXgSEhJYtmyZD07fdJsPFtWspyVHe3VMaWmp3+L1pWBoRzC0AYKjHU1pQ3R0NCUlJa0bkBfqi8HpdLZqbA6HA5Op6enWZrM16WfEFwm9voGeel+DpLV+FXgVID09XfvrUdspMxbVrGfd4l0MgfRo8KkEQzuCoQ0QHO1oShsyMzPPiDtJIiMj0Vrzxz/+kS+//BKlFA888ABTpkzB5XIxffp0vv/+e1JTU3G5XEydOpXx48eTkpLC5MmTWbhwIXa7nXnz5tGzZ0/y8/OZOnUqe/bsISwsjFdffZW+ffsyc+ZMcnJyyMrKIj4+nu7du7N3714OHTrEjh07eO6551i9ejVffvklycnJLFy48KRbFENCQjjvvPO8bpsvEno20LHWdgcgxwf1CiGC1Uzv/jJuer1FjZcBPvnkEzZs2MDGjRs5evQo6enpXHHFFfz4449kZWWxefNmjhw5Qq9evZg6dWrNcfHx8fz888+88sorPPvss7z++us8/vjjnHfeeXz22Wd8++233HbbbWzYsAGA9evXs2LFCkJDQ5k5cya7d+/mu+++Y+vWrQwdOpT58+fzv//7v1x77bUsWrSIa665pkXN98VtiwuA2zx3uwwBirTWJw23CCHEmWLFihVMnDgRo9FIQkICF154IWvXrmXFihVMmDABg8FAYmIiI0aMOOG46667DoCBAweSlZVVU9ett94KwCWXXMKxY8coKnL/Yhk7diyhoaE1x1955ZWYzWbS0tJwOp2MGuW+3yQtLa2mvpZotIeulPoIyADilVLZwOOAGUBrPQdYDFwF7ALKgdtbHJUQIrh52ZNuLVrXOyrc4PfVrFYrAEajEYfD0eAx1bcc1p2St/p4g8GA2WyuKWcwGGrqawlv7nKZqLVO0lqbtdYdtNZvaK3neJI5nrtbfqu17qq1TtNaB86cuEKIs9KwYcP417/+hdPpJC8vj5UrVzJo0CAuuugi5s+fj8vlIjc316sLksOGDeODDz4A3NcU4uPjiYqKauUW1E8e/RdCnHWuvfZaVq1aRb9+/VBK8Ze//IXExESuv/56li5dyrnnnkv37t0ZPHgw0dGnHu+fOXMmt99+O3379iUsLIx33nnnNLXiZJLQhRBnjdLSUsA9JPLMM8/wzDPPAMdvZTQYDDz77LNERERw7NgxBg0aRFpaGsAJY9zp6ek1vfe4uDg+//zzk841c+bMU25Xx1LfvuaShC6EELVcffXVFBYWUlVVxaOPPkpiYqK/Q/KaJHQhhKglkB/2ktkWhRAiSEhCF0KIICEJXQghgoQkdCGECBKS0IUQIkhIQhdCiCAhCV0IcdZ48skn6dmzJ5dddhkTJ07k2Wef5bXXXmP48OH069eP66+/nvLycgCmTJnCb37zG0aMGEGXLl34/vvvmTp1Kr169WLKlCk1dUZERPCnP/2JgQMHMnLkSNasWUNGRgZdunRhwYIFgPuhpIsvvpgBAwYwYMAAVq5c2Srtk/vQhRCnXdo7aa1S7+bJmxvct27dOubPn89///tfHA4HAwYMYODAgVx33XXcdNNNREZG8sgjj/DGG29wzz33AFBQUMC3337LggULGDNmDD/++COvv/46559/Phs2bKB///6UlZWRkZHB008/zbXXXssjjzzCkiVL2Lp1K5MnT2bs2LG0a9eOJUuWEBISws6dO5k4cSLr1vl+2itJ6EKIs8KKFSsYN25czXS2Y8aMAeCXX37hoYceoqSkhNLSUq644oqaY8aMGYNSirS0NBISEmqmAejTpw9ZWVn0798fi8VywjS4Vqu1Zorc6ukC7HY706dPZ8OGDRiNRnbs2NEqbZSELoQ47U7Vk24tDU2NO2XKFD744AMuuOAC3n777ROeFK093W31evV29XS3dafBrX1MdZl//OMfJCQksHHjRlwuFyEhIT5vH8gYuhDiLHHRRRexcOFCbDYbpaWlLFrkfhVlSUkJiYmJ2O32mmlwfa2oqIikpCQMBgPvvfceTqezVc4jPXQhxFnh/PPPZ+zYsfTr14/OnTuTnp5OdHQ0Tz75JJdccgkpKSmkpaW1ysui7777bq6//nrmzZvHiBEjTnrxha9IQhdCnDUefPBBZs6cSXl5OcOGDeOBBx5gwIABTJo06aQXWL/99ts16ykpKfzyyy/17jvVNLjV+7p168amTZtqvn/qqad80JqTSUIXQpw1pk2bxtatW7HZbEyePJkBAwb4OySfkoQuhDhrfPjhh/4OoVXJRVEhhAgSktCFECJISEIXQoggIQldCCGChCR0IYQIEpLQhRDCC9WP8Z/JJKELIc4aWVlZ9OrVizvvvJM+ffpw+eWXU1FRwaZNmxgyZAh9+/bl2muvpaCgAICMjAwefvhhhg8fzvPPP09GRgb33Xcfw4YNo1evXqxdu5brrruObt268cgjj/i5dXIfOikzFtWsZ80a7cdIhDh7ZPbs1Sr19tqW2WiZnTt38tFHH/Haa69xww03MH/+fGbNmsXLL7/M8OHDeeyxx3jiiSeYPXs2AIWFhXz//fcALFy4EIvFwg8//MDzzz/PuHHjWL9+PXFxcXTt2pX77ruPNm3atErbvCE9dCHEWSU1NZX+/fsDMHDgQHbv3k1RURHDhw8HYPLkyfzwww815W+88cYTjh87dizgniq3T58+JCUlYbVa6dKlCwcOHDg9jWjAWd9DF0Kcft70pFtL7WlwjUYjhYWFpyxfdyItb6bU9RfpoQshzmrR0dHExMSwfPlyAN57772a3nqgkR66EOKsN2fOHB544AHKy8vp0qULb731lr9DahZJ6EKIs0bdaXAffPBBwP2Si9WrV59Uvvbbi+puZ2RkkJGR0WBZf/BqyEUpNUoptV0ptUspNaOe/dFKqYVKqY1KqS1Kqdt9H6oQQohTaTShK6WMwMvAlUBvYKJSqnedYr8Ftmqt+wEZwN+VUhYfxyqEEOIUvOmhDwJ2aa33aK2rgLnAuDplNBCp3G9KjQDygTP/sSohxGnV0Iuaxcma82+lGjtIKTUeGKW1vsOzfSswWGs9vVaZSGAB0BOIBG7UWi+qp65pwDSAhISEgXPnzm1ywL6w+WBRvd+nJUc3eExpaSkRERGtFdJpEwztCIY2QHC0oyltiIiIICEhgejoaNx9vzOH0+nEaDT6O4waWmuKiorIzc094RV3ACNGjFivtU6v7zhvLorW9y9f97fAFcAG4BKgK7BEKbVca11cJ8hXgVcB0tPTde0LCqfTlBkn/a4BIOuWjAaPWbZsGf6K15eCoR3B0AYIjnY0pQ12u53s7GwOHjzYukE1g81mIyQkxN9hnCAkJIR+/fphNpu9PsabhJ4NdKy13QHIqVPmdmCWdnf3dyml9uLura/xOpJWII/1C3HmMJvNpKam+juMei1btozzzjvP32G0mDdj6GuBbkqpVM+FzptwD6/Uth+4FEAplQD0APb4MlAhhBCn1mgPXWvtUEpNB74CjMCbWustSqm7PPvnAE8CbyulNuMeovmT1vpoK8YthBCiDq8eLNJaLwYW1/luTq31HOBy34YmhBCiKWQuFyGECBKS0IUQIkhIQhdCiCAhCV0IIYKEJHQhhAgSktCFECJISEIXQoggIQldCCGChCR0IYQIEpLQhRAiSEhCF0KIICEJXQghgoQkdCGECBKS0IUQIkhIQhdCiCAhCV0IIYKEJHQhhAgSktCFECJISEIXQoggIQldCCGChCR0IYQIEpLQhRAiSEhCF0KIICEJXQghgoQkdCGECBKS0IUQIkhIQhdCiCAhCV0IIYKEJHQhhAgSktCFECJISEIXQoggIQldCCGChCR0IYQIEl4ldKXUKKXUdqXULqXUjAbKZCilNiiltiilvvdtmEIIIRpjaqyAUsoIvAxcBmQDa5VSC7TWW2uViQFeAUZprfcrpdq1UrxCCCEa4E0PfRCwS2u9R2tdBcwFxtUpczPwidZ6P4DW+ohvwxRCCNEYpbU+dQGlxuPued/h2b4VGKy1nl6rzGzADPQBIoHntdbv1lPXNGAaQEJCwsC5c+f6qBn123ywqGY9LTm63u9rq12mrtLSUiIiInwXnJ8EQzuCoQ0QHO0IhjZAYLVjxIgR67XW6fXta3TIBVD1fFf3t4AJGAhcCoQCq5RSq7XWO044SOtXgVcB0tPTdUZGhhenb74pMxbVrGfdklHv97XVLlPXsmXLaO14T4dgaEcwtAGCox3B0AYInnZ4M+SSDXSstd0ByKmnzH+01mVa66PAD0A/34ToCxq+fgQ+vBFKcv0djBBCtApvEvpaoJtSKlUpZQFuAhbUKfM5cLFSyqSUCgMGA5m+DbX5Jhm/gZUvwo7/wFtX0p6j/g5JCCF8rtGErrV2ANOBr3An6X9rrbcope5SSt3lKZMJ/AfYBKwBXtda/9J6YXuvt8riUdP77o3ojpC/m39b/0Jnddi/gQkhhI95dR+61nqx1rq71rqr1vqvnu/maK3n1CrzjNa6t9b6XK317FaKt0nCsPGi+UWsyg4DJsNdK6DD+XRQR/mX5UmiKPV3iEII4TNB/aTo46Z36Wo4xDZXRxg1C0Jj4NbP2ODqSqIqYJJxqb9DFEIInwnehO6o5FrjclxaMd1+D1jC3N9bI/i7YwIAk01fYcbhxyCFEMJ3AjKh2xw2fjz4I0v2LWm4UO4WLMrJHp3ELt3hhF3LXWlsc3UkQRUyxrCylaMVQojTIyATeqm9lLu+uYu/rPpLw4UObQBgs06tZ6fiDeeVANxpWszJt9ULIUTgCciEHhcSh8VgobCykHJ7ef2FcjYA8IurvoQOnzsv5IiOoZdhPxcazogbcoQQokUCMqEblIGkiCQADpUdqr/QoY0AbG4goVdh5h3H5QDcaVzs+yCFEOI0C8iEDpAYngg0kNAdVXDEPRnkFp3SYB0fOC+lQlvIMG6ki6r78KsQQgSWgE3o7cPbA5BTWk8iPrIVnFXsdiVRRmiDdRQSyULnUADGGeXiqBAisAVsQq8ecjlcVs8Tn54Lor/Ue0H0RJ+5LgRgrOFHaGTmSSGEOJN5M9viGSkp3J3Qc8rq6aF7LojWHj9PaWCGxdWu3hzRMaQaciHnZ0geeMrz1q4na9boRuNsankhhGiugO2hVw+5HCqtZwzdc0HUmx66CwNfOIe4NzZ/7LP4hBDidAvYhF7dQz/poqjTDrlbANjiSvGqrs+dF7hXfpkPLqevQhRCiNMqYBN6YngiCsWR8iM4XLUe3z+SCc5KiOtKCWFe1bVRdyXLlQCluZC1vJUiFkKI1hWwCd1sNNM2tC1O7eRIea1XmHouiJLUlPdrKD53eXrpm+f5KkQhhDitAjahAyRG1HMvuueCKO37N6muBdXDLlsXgt3W8uCEEOI0C+iEXu+96J4LoiT1b1Jdu3UyJPaFyiLY9Y2PIhRCiNMnoBP6SY//a11zQZSkvk2vsM+17s/Mum/YE0KIM19gJ/S6d7pUFICjAqzREBrb9Ap7j3N/bv8SHJU+ilIIIU6PgE7oJ92LXuL5jExsXoVtukJCGlQWw55lLQ9QCCFOo4BO6CdN0FXimQaguQkdjvfSt37egsiEEOL0C+iE3j7C00MvO4TW2rcJfdsX7lkbhRAiQAR0Qo+0RBJpjqTCUUFhZWHLh1wA2naHtr3AVgRZP/gkTiGEOB0COqFDnXvRS3PdX0YmtaxSGXYRQgSggE/oJ1wYre6hRyS0rNLqhJ75BTgdpy4rhBBniIBP6CdMo1szht7CHnq7XtCmG1Tkw74VLYxQCCFOj8BP6LUfLiqpHnJpYQ9dqeO99C2ftawuIYQ4TQI+oR8fcsmBUk8PPaIFF0Wr9bnG/Zm5UKbUFUIEhIBP6NX3oueUZIOzCkKiweLdtLmnlHAuxHWB8qOw78eW1yeEEK0s4BN6zb3o1RN0+aJ3Dp5hl2vc6zLsIoQIAAGf0OND44kwR1BgL+GQ0diye9Drqj3somXYRQhxZgv4hG5QBs5PPB+AVaEhLb/DpbbEvhCbAmVHiC7K9F29QgjRCgI+oQMMbT8UgNWhIS2/w6W2WsMu7Y7IOLoQ4szmVUJXSo1SSm1XSu1SSs04RbnzlVJOpdR434XYuCFJQwB3Qnf5agy9mmfYJf7oKrnbRQhxRms0oSuljMDLwJVAb2CiUqp3A+WeBr7ydZCNSYlKIREzBUYjO0w+rjypP8R0xlpVAPtX+7hyIYTwHW966IOAXVrrPVrrKmAuMK6ecvcA84Ej9exrVUophrqMAKyy+fj0Sh2/OPrLfN/WLYQQPqS01qcu4B4+GaW1vsOzfSswWGs9vVaZZOBD4BLgDeALrfXH9dQ1DZgGkJCQMHDu3Lm+agdFG37NI7EWeptT+U37+wHYfLCoSXWkJUfX+3146V7OX/d77KZIXk15HpcynbJ8bbVj8KZ8aystLSUiIsLfYbRIMLQBgqMdwdAGCKx2jBgxYr3WOr2+fd4MUKh6vqv7W2A28CettVOp+op7DtL6VeBVgPT0dJ2RkeHF6b2gNcdWFEBsArtdOQy9eChWo5UpMxY1qZqsWxqIRw+nbOuzhJdns3brNpa5+p+6fC21Y/CmfGtbtmwZPvt395NgaAMERzuCoQ0QPO3wZsglG+hYa7sDkFOnTDowVymVBYwHXlFKXeOLAL1Snk8bRyU97U4qnZX898h/fVu/UhxpNwyAMcaVvq1bCCF8xJuEvhboppRKVUpZgJuABbULaK1TtdYpWusU4GPgbq31Z74OtkGeOVyGaCsAq3JW+fwU1Qn9CsM6rMibjIQQZ55GE7rW2gFMx333Sibwb631FqXUXUqpu1o7QK945kEfamkDwOpDvr8bpSIsCdoPIELZuNTws8/rF0KIlvLqJj+t9WJgcZ3v5jRQdkrLw2oizzzoAyJSsFYUknksk4OlB31/nrTxkPMzY42rWOwa4vv6hRCiBYLiSdHqhB4SlczIziPRaObvaIVbDPtch0srRhg2EEWZ7+sXQogWCKqETmQiE7pPAODTXZ8CPn6yMyqJ1a5eWJWdUcY1J+xKmbGoZhFCCH8IkoTueZdoZCID2g2gS3QXjlYcxRTp+wm1PnVdBMB44w8+r1sIIVoiOBJ6afWr55JQStX00s0xP/n8VIucQyjTVgYZtsOx3T6vXwghmis4Enr1kEuEe6bFMV3HYDVaMUXsRJmP+fRU5YSw2DnYvbHhQ5/WLYQQLRH4CV3rE8bQAaKt0VyRcgUA5pg1DR3ZbB87h7tXNn4kMzAKIc4YgZ/QbUXgsoMlEsyhNV8fH3ZZB8ru01P+pHuyz9UOig/CnmU+rVsIIZor8BN6RYH7Myz2hK/7te2HsyIZg6kMc6yvnxxVfOx0PznKhg98XLcQQjRPECT0fPdn6IkJXSlFZd7lAFjbLAODzaen/cR5MaAg8wuoKPRp3UII0RxBkNA9PfTQuJN2Ocu64yhLRZnKscT59jbDg7SF1GHgrITN83xatxBCNEfgJ/Ty6oQeW89ORWXeKAAsbVagjCW+PffAye7PtW9w8ozCQghxegV+Qq8ZQz+5hw7gquiMo6QXylCFJf4735675xgIbwd5mQxS25p0qDxZKoTwtSBI6PWPoddWmXcFWivMsT9hsPjwFXUmCwycAsBtpiW+q1cIIZohCBJ6w2Po1VyVidgL01HKSUjSfMDlu/MPnALKyBWGtbSlwHf1CiFEEwV+Qi9vvIcOUHnkKlz2SIxh+3x7G2N0MvS4ErNycpPRx0M6QgjRBIGf0KuHXBoYQ6/hCqXy8DUAWNt9hTLn+y6GQXcCcLPpW4y+nuFRCCG8FAQJ/VR3uZzIUdoHe3EaylBFSOKn+OzOlNTh7HYlkaTyudywzjd1CiFEEwV+Qq8Zcmmkh+5ReXgs2hmKKWIn5rgVvolBKd5xuh9immZa5J5fRgghTrPAT+iN3LZYl3ZGYjt0PQDWdl9iDM3ySRjznMMp0BGcZ9gF+1b6pE4hhGiKwE7oLqd7ci4UhER7fZij5Fyqjl2MUi5Ckj9EGUtbHEoFIbzr6aXz4/Mtrk8IIZoqsBO6rQjQ7mRuMDbp0Mojo3CUd8ZgLiak/Vx8cSvjO47LqdAW2PkV5G5tcX1CCNEUgZ3QvbxlsX5GbAdvxuUIxxSxC2vCAnQLx77zieLf1XOlSy9dCHGamfwdQIt4xs83HDNwTTMeodeOaGzZtxLa6XUscat5Z8s7TDl3yknlUmYs4oE0B1O8OMdrztFMMn6Da+M8hq+5AIhv9Jjaj/9nzRrdlCYIIUSNAE/o7h56oY5odhXOihRsOTcQ2uFD/r7+7ySGJzIqdVS9ZY0uJ6nFhzinMJt25QUc/OP3OPLy0HY7f99zFKWh2BLG2rBudAk9xEOR7/O3yEkcDo9Dq8D+Y0gIceYL8ITu7qEX0PyEDuAo6Ystt5CQhMU8vOJhoixRXJB8AQBV+/YxYce3DPjvNj7OPkCI8/jbj4p3HK+jd506jxHJOWTzJrMoM4Wwoe05FM63ETFsGKa2bVsUrxBC1CewE3p5y3vo1ez5F/OrjFg+yPyAhxbfwz/KriZm2SYqt21jaq1y2eHx7IjtRE54G/5wWwamhAQMVgvXv7oGhSa6soyYyhJuqfqGzkW5FOSHE26zceGhXzj050cACBsyhOGOrqxMSsNuDOz/BEKIM0dgZxNPD90XCR0U90ZfS9qK5SSv2oPV8W8qAUNEBEtiuhM/pDdPlZ1DkfX4uf56/fHx7sz5J87iuFWlstTyIEkUckPhoyQePsZDMccoW7WK8tWrmcFqiixhLE4Zymddh1FsDfdBG4QQZ7PAHtj1jKG3dMilU/FhZqx9j6xrrqXL8j1YHbAhVTF7QghH5z3Ds+k3k9fr3BOSeWP26UTmOYdjUi4mRy1hceoFdPy/OXRbsZyExx5ld3R7oqvKmbhjKe98/T/cuXkBsbbiFrVDCHF2C+weeguHXOIqirh962IuOfAzBjTKbCZmwnhiJk3iX4ffZeXO+axZfh+mqPHAuU2u/0XHtVxnXM4YwypeUWMBMEZFEXfzzUzfGEWv/H3ctOMbBuVu47rdP3Bl1mqOzikibsoUDCEhzWqTEOLsFeA9dM+QSxN76Gangwk7vuW1pf/LyAPrcSoDC1MvoOuSr0l87DFCunTh8aGPc3uf23FoB6HJc9nsavrj/DnE86HzUgxK87DpwxPneFGKzDYpPD70DqZn/J7Vib0JdVaRN/t5do+6kuLFi1t8X7wQ4uwS4Am96T30XseyeOm755i6dTFhjkpWJvbh1yP/wCv9rsOcmFhTTinF/en388DABwBY5VqMNWEBNHF63Bcd11Kswxhu3AQ7/lNvmd0xHXhiyFT+dOFdWHv1wnH4MAfvf4ADd92F/eDBJp1PCHH2CvCE7v1tiyGOSu7a9BnPLn+ZTqVHOBDRlj9fcCdPDrmdQ+ENP/wz5dwpVBy8AQNGLHErCe34NhjKvQ4xnyiec4x3b/xnBthtDZbd1PYcUj+eR+ITT2CIjKTs+x/YffUY8t99D+3y4VuWhBBByauErpQapZTarpTapZSaUc/+W5RSmzzLSqVUP9+HWo9yd0IvaqSH3jM/i1e+fY5xe1aglWJu90v57Yj7+bldD69O4ygewNXG2z3TBOwkPOUVDJZcr8N83zmS7a4OUJAFq18+ZVllNBJ74w10XbyIqKuuRFdUkPu3v7H/V7/CfuiQ1+cUQpx9Gr0oqpQyAi8DlwHZwFql1AKtde3Zp/YCw7XWBUqpK4FXgcGtEXANpx2qSnBoA8WE1VvE4HIycfs3TNyxFKN2sTuqPc8NuJE9MclNPl2iSqF873RCO76DMeQwYakvcc7fcnAUn9fosQ5MPOG4jQ8tf4Mf/g79JjZY9oRpAJ57jqirruLQY49Tvmo1e8aOI/Gxx4gec3WT4xdCBD9veuiDgF1a6z1a6ypgLjCudgGt9UqtdfUbklcDHXwbZj08wy1FhAPqpN1tywt4ZsUrTNq+BKU187pl8PuM3zUrmVfTjljKs+7GXnQeymAnNPlfWBM/AWVv9NiVrnOh1xiwl7mHXrwUOXIkXRZ8TsSIEbhKSsj5wx/IeeQRXBUVzW6HECI4qcbupFBKjQdGaa3v8GzfCgzWWk9voPyDQM/q8nX2TQOmASQkJAycO3duswMPKzvAoLXTKTAn8X7K0yfsa7NzG30+/zeWinJskVFsueZGClK6NlpnWnL9c6pvPlhEQijkenKo1prtej0/ur7AiYMY2nKJcQLxqv0p609vU8n5a+/B5LSxOHE6uyMHeR+P1oQuX0HkvHkoux1HUhKFd96Bs/2pz1lXaWkpERG+eBDLf4KhDRAc7QiGNkBgtWPEiBHrtdbp9e3zJqFPAK6ok9AHaa3vqafsCOAV4CKt9bFT1Zuenq7XrWvB+zf3rYK3RrHe1Y3rq54A3EMskzP/ww07vwNgTUJPnh04kRKLd09hNjTTYfVsi3/ffOIIlcGaQ0jyXIzWI2htoCrvMqqODQPqn5s9a9ZoWPMaLH6QPB3F5ZX/SwFRTYrHtn07B++7n6o9e1ChoSQ9+STRV3s/Q+OyZcvIyMjwuvyZKBjaAMHRjmBoAwRWO5RSDSZ0b4ZcsoGOtbY7ADn1nKQv8DowrrFk7hPVd7h4LohGV5by15WvccPO73AqA2/2voqZQ6Z6ncybw1XZnvK991CVfwFKubC2+4qwlH9isB5u+KD0X0HKxbRVxTxufrfJ5wzp0YPUj+cRPW4suqKCnAcf5PD//BVdVdWClgghgoE3CX0t0E0plaqUsgA3AQtqF1BKdQI+AW7VWu+opw7fq74HnUi6F+znhWWz6X90FwXWCB668NfM637J6ZmyVpupzB1L+f6puOzRGEOzCUt9AUv8ElCOk8sbDDD2Bcq1lWuMK7nM0PS/UgxhYSTNmkXi44+B2UzB+++zb8rtOPLyfNAgIUSgajTjaa0dwHTgKyAT+LfWeotS6i6l1F2eYo8BbYBXlFIblFItGEvxkqeHbt5byTPLX6FdRSGZsZ25J+M+Nsc3Pl7ua86y7pTtuY+q/CHu3nrbpYSnzsYYXs/vt7guPOO4AYCnza+SSNP/oFFKETtxIinvv4cpMZGKn39m7/Xjqdi0qaVNEUIEKK+6sFrrxVrr7lrrrlrrv3q+m6O1nuNZv0NrHau17u9Z6h3f8SVdcpTcn6Povu4AFpeDRSlD+NNFv+FYqPcvi/Y5VwiVuddQnvVrnJXtMFiPEtbpTUKS30OZ808o+rbzCpY5+xGnSnnB8hLGJj6BWi20Xz9SP55HaPpAHEeOsO+WSRTO/8QXrRFCBJiAfFLUUVDA/he/IX9HBE6leKHf9bzUf/wZM7e4syKV8j2/w5Z7FdplwRy1hfAuzzF7/WxKq0oB0Bi43/4bDutYBhm2c5/p42afzxQfT+c33yT25pvRdjuH/vxncp96Cu2oZ8hHCBG0AjKhOwsKsR0owmh18tnFw/kydai/Q6qHCXv+MMp2P4C9qD/K4OCNX95g9Kej+SDzA1AO8oni3qrpOLXibuMChhk2NvtsymIh8bFHSXzyL2A2k//OuxyY9mucRUU+bJMQ4kwWkAnd2iWVDtclkHp5Hpnxnf0dzilpRzS2nJso2/tbzmt3Hvm2fGatmUV4l79jil7PT7oHsx3XY1CaF80v0kWddANRk8ROmEDnt9/CGBdH2cqV7L3hBip37/ZRa4QQZ7KATOgpMxaxz5mHOdxFoY70dzhecdk68s6od5idMZuu0V0xWAoIbT+P8C7PMSe8M4udA4lW5bxhfoYYSlp0rrCBA0n9eB7WXr2w79tP1o03Ufr99z5qiRDiTBWQCR0gRrmTXqEOnFe3KaW4tPOlzB87n4qcCbiq4jBYjxKSPI8/dnLyYnhnkg25zLHMxkzLxr/N7duT8sH7RF5xBa7SUg7c9RuOvfHmiXOyCyGCSuAmdMoAKCAweui1GQ1GHEUDKdv9ABU543FVtUFZ8nm1nebyjh3IjM3mCcscaOGUuYawMJL/8Rzx90wHrTnyzDNEvf02rspKH7VECHEmCciEbqWKMFVJlTZSjtXf4bSAEUdROmW776fi4I04bYkcMxl4Li6WF1L38/S8qzlQvL9FZ1AGA21/+1uSX3geFRpK6E9r2HfrbdhzjzR+sBAioARkQo+kgjJtpZBI6ptpMfAYcRSfR/neeynffzvh5QmUGQy8bzvA6E9HM33pb1mevRyXbn6PPeryy0n56EOccXHYNm0ia/x4yv/7Xx+2QQjhbwGZ0I8STZ/Kt7ig8gV/h+JjCmdZDw7vu49eWZczrrQck9Z8n/0Ddy+9m9GfjOaNzW9wtOJos2oP6dmTYw/NIOz883Hk5bHvtskUzJvn4zYIIfwlIBN6NUfj7+cIWGsqLuF/RjzHkgOHuTe/kPaGELJLs5n982xGzhvJvd/ey7IDy7C7Gp+LvTYdGUmnN98gdtIksNs5/OhjHJo5E5dM7iVEwAvohB70eo+jzYR3uKO0ksW7d/CytRsjOgwH4NsD33LPt/cwct5Inl7zNFuPbaWxqZCrKbOZxEf+TNJTT6EsFgrn/ot9k26VV9wJEeAkoZ/peo6GSfMxWiIZtm0pL2TvZ8noedw74F5So1PJt+Xzfub73PjFjYz9bCz/3PhPsoqyvKo65tpr6Pzhh5jbt8e2aRN7r7uespUrW7c9QohWIwk9EKReDFMWQnhb2Ps9bT+4gTsSL+bzcZ/z0eiPmNhzInEhcWQVZ/HKhlcY89kYrl9wPf+38f/YW7T3lFWHntuHlPkfE37RRTgLCtj/qzvIe+lltLN5k4UJIfxHEnqgaH8e3PkdJKZB/h547VLU9i85N/5cHh78MEsnLGXOyDmM7TqWCHMEOwp28NKGlxj72VjGfjaW539+nk15m+q9U8YUG0vH/5tD/N13A3D0pZfY/6s7ZH51IQJM8F5VbIaUGYvOqPprl8+aNRpiOsLUr+Czu2HrZzB3Im85rmCWYyKVWMiaNZoLky+kylnFqpxVfL3va5YdWMbeor28vvl1Xt/8OqFEMOrHyxjeYThD2g8h3Ox+0lYZjbT93T2EpQ/k4B/+SPnq1ey55lraz3qKiIsv9uG/ghCitUhCDzSWcJjwNqx8EZY+we2mrxhiyOR39uPv7LYYLQzvOJzhHYdjd9lZn7ueyXPfxBS5lQpzEZ/u+pRPd32KSZk4L+E8Lkq+iKFJQ+kR14PwCy4g9dNPyPnDHyn/6ScO3DmN2Ftvpd2DD2CwBvJDXEIEP0nogUgpuPB3kHIRe//vJnoZ9vOF5WH4vhguvBdMlpqiZoOZIUlDqMw9RmXuWKace5D2fRXLs5ez6egm1h5ey9rDa/kH/yAuJI7BiYMZnDSY85+fSdt/f03eCy9S8N57lK9eTfunZxHSu7cfGy6EOBVJ6IEseQBXV/2NR03vcZNpGXz3P7DlExjzAnQ8v54DFG1UItP6jmNa32kUVRaxKmcVK3NWsurQKg6XHebLrC/5MutLAJLikrj8Txcx8q3NsHMne2+4kfi77iL+19NQZvNpbaoQonGS0ANcGaHMcEzjc9eFfJTwIRzZCm+MhLQbYOTjEN2hwWOjrdGMSh3FqNRRaK3ZW7yXnw79xJpDa1hzeA2Hyg7xDoeYO1Fz8zLFlesdHH3pJfYtmof1z/fRc+hVmI2S2IU4U0hCDxKrXH3g7lXwwzOw8iXY/G/IXAgXTIehv230eKUUXaK70CW6CxN7TsTpcrKzcCfrDq9jXe46Pov8mZ96HOM3i10k7M3FdecMXjr/EX65ri89k/uT1jaNtPg02oe3R6lgmF9HiMAjCT2YmEPh0sdgwG2w5HH3nTA/PAM/vcrvTZfypuNK8HJ2SqPBSM+4nvSM68mk3pPQWpNVnMXGa34i//UP6bFkJ6PWOEjP/Jl3L93Auz0VKEWsNZbe8b3p06YPveN606tNL5LCkyTJC3EaSEIPRrEpcMM7sP8n+O6vsPd7fm/6hKnGL9mVlwGFA9y3QDaBUorU6FRS+6bCCzdRsWULBx99hPit27j/MxfZ50TzxkjFltgCfjz4Iz8e/LHm2GhrND1ie9Ajrgc9YnvQPbY7XWK6YDXKXTNC+JIk9GDWaTBMXgBZP7LyjQe5wLiVAYVfwvNfQ++xkD4VUi523zXTRKF9+tB13scUzvuYvNmz6bCrkMf3GDCPvox9NwxlkyGHbfnbyDyWSUFlAWsOu8flqxmVkc5RnTkn5hz3EnsOXaO70jGqI2aDjMsL0RyS0M8GKRdys/0R0hx7eLbNF/QoWwtbPnUvcV3dQzRpEyA6uUnVKqOR2JtuJOrKUeS99DIFH32EfeFXJH/1Hefecgtt7vgbxrg4cstz2Z6/ne0F29mWv42dBTvZX7KfPUV72FO0h6/3fV1Tp0mZ6BTVidToVFKiUtyf0SmkRKUQbY329b+MEEFFEvpZZLPuwtdJd9Pj6oHw87vw83uQvxu+eRy+mQkpF0HaeOgxGiLael2vMTqaxD8/TNykW8h7/gWKFy8m/623KPjoI2JvvIE2U6fWPOhUzeawsadoD7sLd7OzcCe7Cnaxp2gPOaU5NYm+rhhrDJ2iOmEtt5K5MZOOkR3pGNmRDhEdiAuJk3F6cdaThB4gfDotQXQHGPEwDPsj7FoCGz+C7f+BrOXu5Yv7oNNQ6Hk1dLsM2pzT6LBMTXyWS+ia0Z1J275myOGt5L/zLgUffkTUmDHETZ5MSI/uAPR8ZKnnSANZs+4/Xs9Dn2Kw5GGwHuX+q2LYW7SXfcX72Jq3m8LKQgrzCgFYu2HtCecPM4WRHJlMckQyHSI68PqyArQ9Fpc9lvUzJhBjjWlSwvfm3ztr1miv62vOeVujfuE/p+O/rST0s5nRBD2udC+2Iti6ALZ+DnuWwb4f3ctXD0FMZzhnpHvWx5SLITz+lNXujunAE0Om0qXwIO9at1Ly9dcUffIJRZ98QtjQIcTeNBGjy4nTYDz5YG3BVZmMqzKZu/sf/6FPmfEFylSCwXyUOwaWEpoUyoGSAxwoOUB2STYl9hJ2FuxkZ8FOAEISjlc57F8vEmoKJTE8kcSwRJIikkgMSyQhPIGEsATahbWjXVg7oixR0ssXAU0SunALiYYBt7oXWxHsXAI7/gO7lkLhPlj3hnsBaNcbOg6GTkOg4yCITa23yj0xyXSYNY2qrCzy33ufwk8/pXzVaspXreY9awTfdErnm47pXgao0I4onI4ohkaEkzEw44S9RZVFZJdmc7DkIDmlOTy1ZCUGcwHKXEBUZCll9jL2Fu095XTCIcYQ2oW1Iz40npD2DrQjEu2IwuWMQDuiPNsRaGcYMlGpOBNJQhcnC4l2j6WnjQeXE3L+C3u+g73L4cBP7qdRj2yF9W+5y4fG8q65A5t1KltdKWTqTmTpRFyepGdJSSHx0Udoe+/vKPr0UwrmzSN2124m7FzGhJ3L2DNuIVGjRxN1xeXNDjnaGk20NZo+bfoAMPP949cAts4aTXFVMYfLDp+w5Jbnkluey5HyIxwpP0KZvYz9JfvZX7If8ymuv2qt0M5wrl/wJnEhcScssSGxxIbEEhcSR4w1hhhrDFGWKIz1/TUihI9JQhenZjBCh3T3MuwP4KiEgz+7E/uBNZC9BsryGGYsYBibaw6r0BZ26fbwyefQtifEd8PY5hzibrmJ2Ntu48pfv8Ll+9ZyUc4m2L6dvO3byXvuOeZEJrAq6VzWJvRE2+0+mzMmyhJFlCWK7rHdGyxTWlVKXkUeeeV53PzWEpSpGIOpBFW9GEvdwz6mcpSplB0FO7w6t0IRZY0i1hpLtDWaGGsM0dZoivOL2bZxG1GWKKKt0URZojCE7Ee7QsEZit1pl6kVRJNIQhdNY7JC56HuBUBrKD7Ir59+jTTDXnqp/fQy7KO9yidNZcGmrDoVKFRUMo8nRLG/XTt+dPTjjnM6U7xuL6XrttC5JJfOJbnctGMpO4a8TdjgwYQPHkS3ggp2R7fH1Yo93QhLBBGWCFKjU3EUn+rlHk6UqYzF9/XnmO0YBbYC8m35HLMdo9BW6N6uzKeosogCWwHFVcUUVRZRVFl0Uk3LNiw7YTu81ujVgPf/hxBjCJGWSCIsEURaIok0u9cjzJ7Fsx5uDifC4v4MN4cTbnJ/hpnDCDeHYzLI/+pnA/mvLFpGKYjuwFeuQXzlGlTzdRSldFMHmT8+DvK2w7Fd7qVgHxRnM9QIQ9nq/gnMhshE0FdCWZ6VowcjKT9swVVSRum331L67be8AFQZjeTEtiW2oBPF+Yexpg3AktoFZTrdP8ZGtCOKXm16eVXa4XLUJPTCykL3elUR67eup22HtjUJv7iqmOW796OMFWCwYTFXYnPasFXYyKto2dujrEYr4eZwQk2hhJnDCDOFuddNYYSaQwk1nbyEmcIIMYUQagolxBRCiNG9bjVaa76vclWhtZaLyWcISeiiVRQTwXrdAwbWuT3LaYeiA9z67L/opI6QrI5yd38LFB1AFecQYjhISuJRAOxlBspyrZTnWak4aoESSDl6GL4/zMHv3U+dKoPGEmfAGm/F0jYCS1IslsR4/uws55g5ihIVDpsrwBrpXiwRYI0AS6T7ZSHm0GY9KdsUJoOJNqFtaBPa5oTvY7JjyBiQccJ3KUuP39r2y1NXUeGooKSqxL3YSyitKqXUXkpJVQll9jJKqkootbsv+pbZyyi1l1JuL6/ZLreXU+Yoo9JZSaWzsnUa+K77F4bVaCXEGILVdHzdYrQQYnJ/VpepXrcYLVgMtdY92xajBbPRXLNuMRzfNhlM7v0Gc82n2WDGbHR/GtTZfbHaq4SulBoFPA8Ygde11rPq7Fee/VcB5cAUrfXPPo5VBAOjGeK6sNzVt+aru8cfT/rdZywkniLaqQK+uOMcYkoOEVOax7vfrCWh8hjxBUUkldlwHnVgO6ZxlJuoPKqpPGqDbTbgKLCTiwBldGEOc7Jv+eeYw5yYQpwYQ1yYrC6MIS6MFhdGKxjDQzCEh6GsYTVJ/gNzORVYsGHBhhWbNlOJhUrMVGozlZhh9X73EJTR6vm0HP+sXkyeT4PZ3Xaj2bNtwuC0gaMKDCYwnJyIlFLu3rQ5jITwhJP2e0trjc1po9xeTrmjnHJ7ORWOippkX+GooMJe4f70LOWOcmwOGxWOCmxO92elo5IKZwU2h61mvaKqAgeOml8YxRQ3O05fMCojZoMZk8F0wqfZaMakTJgM9S9F+UV8/t3nmAwmjAbjSWWNynjCZ93vjAbjCdsGZTi+r3p/+DbAgNPWtCeym6LRhK6UMgIvA5cB2cBapdQCrfXWWsWuBLp5lsHAPz2fQjSJCwNHiOWIjnXfH+/x2H8WubsT8fD2pHAyMjLA5cSZd5DKrRup2rWDqn1ZVO0/SFXuMUpy8rHYHVSVGKgq8e7cBlMFBnM5BpOL9iaNwaRRJo3B6Fk3aJRRo4zuvwzytn+BMrjXqz9RoJR7u3rd/QkoDVSvwwCgbJ57P8oABiPrneDEgFMZKX/Q6r4o7T6he93g/lQ120bAs66Ue7/ylMPzqdwBhCkjYe5Aay3q+Hp1oDXrBjB4glXhoCKOf+9pRM7hwyQmJlGFxo6LKjRV2kmVZ92u3dt27cKOi0rtwoELu3ZRpZ04tMaOZ7924aj+Hve6XTs9ny7suMs7XA7s2oWzpowLp6ecxgE4GvjZgirPUp/dezY3sMc3+niubx+pmtBq5/Cmhz4I2KW13gOglJoLjANqJ/RxwLtaaw2sVkrFKKWStNaHfB6xENUMRowJnQhL6ETYiBN3pcxYRJi9gviKIhbf3BNH7mEcR4/hzD/m/iwowFlUhLOwEGdREa6yMlwOAy4HuH9z+F9hzZrTs9j9Fcop7W/ge7NnCTuNsQSC10c09C/Wcsqdg09RQKnxwCit9R2e7VuBwVrr6bXKfAHM0lqv8GwvBf6ktV5Xp65pwDTPZg9gu68achrE4/57PtAFQzuCoQ0QHO0IhjZAYLWjs9a63smWvOmh13fFqO5vAW/KoLV+FXjVi3OecZRS67TW3j7WeMYKhnYEQxsgONoRDG2A4GmHN5eEs4Hab0PoAOQ0o4wQQohW5E1CXwt0U0qlKqUswE3AgjplFgC3KbchQJGMnwshxOnV6JCL1tqhlJoOfIX7atGbWustSqm7PPvnAItx37K4C/dti7e3Xsh+E5BDRfUIhnYEQxsgONoRDG2AIGlHoxdFhRBCBIaz+7EqIYQIIpLQhRAiSEhCb4RSqqNS6julVKZSaotS6l5/x9RcSimjUuq/nucGApLnobWPlVLbPP9Nhvo7pqZSSt3n+Vn6RSn1kVIqxN8xeUMp9aZS6ohS6pda38UppZYopXZ6PmP9GWNjGmjDM56fp01KqU+VUjF+DLFFJKE3zgE8oLXuBQwBfquU6u3nmJrrXiDT30G00PPAf7TWPYF+BFh7lFLJwO+AdK31ubhvNLjJv1F57W1gVJ3vZgBLtdbdgKWe7TPZ25zchiXAuVrrvsAO4KHTHZSvSEJvhNb6UPVEY1rrEtwJpPVm12klSqkOwGjgdX/H0lxKqShgGPAGgNa6Smtd6NegmscEhCqlTLifjA+IZza01j8A+XW+Hge841l/B7jmdMbUVPW1QWv9tda6egKY1bifowlIktCbQCmVApwH/OTnUJpjNvBH3HMUBaouQB7wlmfo6HWlVLi/g2oKrfVB4FncU6Acwv3Mxtf+japFEqqfOfF8tvNzPC01FfjS30E0lyR0LymlIoD5wO+11v6dI7SJlFJXA0e01uv9HUsLmXBPUvhPrfV5QBln/p/4J/CMMY8DUoH2QLhSapJ/oxIASqk/4x5i/cDfsTSXJHQvKKXMuJP5B1rrT/wdTzNcCIxVSmUBc4FLlFLv+zekZskGsrXW1X8hfYw7wQeSkcBerXWe1toOfAJc4OeYWiJXKZUE4Pk84ud4mkUpNRm4GrhFB/DDOZLQG+F5eccbQKbW+jl/x9McWuuHtNYdtNYpuC/Afau1Drheodb6MHBAKdXD89WlnDiNcyDYDwxRSoV5frYuJcAu7NaxAJjsWZ8MfO7HWJrF8wKfPwFjtdbl/o6nJSShN+5C4FbcvdoNnuUqfwd1FrsH+EAptQnoD/zNv+E0jeevi4+Bn4HNuP8fDIjHzpVSHwGrgB5KqWyl1K+AWcBlSqmduF+CM+tUdfhbA214CYgElnj+/57j1yBbQB79F0KIICE9dCGECBKS0IUQIkhIQhdCiCAhCV0IIYKEJHQhhAgSktCFECJISEIXQogg8f/+0uj0Yo9zkQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "f = Fitter(X_NC, distributions=['gamma',\n", " 'lognorm',\n", " \"norm\"])\n", "f.fit()\n", "f.summary()" ] }, { "cell_type": "code", "execution_count": 38, "id": "1beea73b", "metadata": {}, "outputs": [], "source": [ "\"\"\"\n", "Set up a helper function for checking p-values against an alpha level, and printing result\n", "\"\"\"\n", "def check_p_val(p_val, alpha):\n", "\n", " if p_val < alpha:\n", " print('We have evidence to reject the null hypothesis.')\n", " else:\n", " print('We do not have evidence to reject the null hypothesis.')" ] }, { "cell_type": "code", "execution_count": 41, "id": "6d9288f6", "metadata": {}, "outputs": [], "source": [ "import scipy.stats as stats\n", "from scipy.stats import normaltest\n", "from scipy.stats import kstest" ] }, { "cell_type": "code", "execution_count": 50, "id": "0a8f8714", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(154,)" ] }, "execution_count": 50, "metadata": {}, "output_type": "execute_result" } ], "source": [ "L = np.log(X_NC)\n", "L.shape" ] }, { "cell_type": "code", "execution_count": 51, "id": "83ee07de", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Statistic: \t0.50 \n", "P-Value: \t9.17e-36\n", "\n", "We have evidence to reject the null hypothesis.\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/home/marija/.conda/envs/pytorch/lib/python3.7/site-packages/scipy/stats/_continuous_distns.py:4837: IntegrationWarning: The algorithm does not converge. Roundoff error is detected\n", " in the extrapolation table. It is assumed that the requested tolerance\n", " cannot be achieved, and that the returned result (if full_output = 1) is \n", " the best which can be obtained.\n", " intg = integrate.quad(f, -xi, np.pi/2, **intg_kwargs)[0]\n" ] } ], "source": [ "stat, p_val = kstest(L, 'norm')\n", "print('Statistic: \\t{:1.2f} \\nP-Value: \\t{:1.2e}\\n'.format(stat, p_val))\n", "check_p_val(p_val, alpha=0.1)" ] }, { "cell_type": "code", "execution_count": 66, "id": "90b3ec11", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "distribution is not normal\n", "distribution is not log-normal\n" ] } ], "source": [ "import scipy as sp\n", "\n", "p=sp.stats.mstats.normaltest(X_NC, axis=0).pvalue\n", "if p<0.025:\n", " print ('distribution is not normal')\n", "p=sp.stats.mstats.normaltest(np.log(X_NC), axis=0).pvalue\n", "if p<0.025:\n", " print ('distribution is not log-normal')" ] }, { "cell_type": "code", "execution_count": 67, "id": "53177ff3", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.7529322013966726\n", "0.694634906696082\n", "1.204700345809746\n", "0.186230860341136\n", "0.7529322013966726\n" ] } ], "source": [ "shape, location, scale = sp.stats.lognorm.fit(X_NC)\n", "mu, sigma = np.log(scale), shape\n", "print(shape)\n", "print(location)\n", "print(scale)\n", "print(mu)\n", "print(sigma)" ] }, { "cell_type": "code", "execution_count": 61, "id": "adc791d9", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "KstestResult(statistic=0.08120944008227199, pvalue=0.24802604184218757)" ] }, "execution_count": 61, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sp.stats.kstest(X_NC, \"lognorm\", sp.stats.lognorm.fit(X_NC))" ] }, { "cell_type": "code", "execution_count": null, "id": "bdff3710", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "id": "60af4dca", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "hide_input": false, "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.6" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false }, "varInspector": { "cols": { "lenName": 16, "lenType": 16, "lenVar": 40 }, "kernels_config": { "python": { "delete_cmd_postfix": "", "delete_cmd_prefix": "del ", "library": "var_list.py", "varRefreshCmd": "print(var_dic_list())" }, "r": { "delete_cmd_postfix": ") ", "delete_cmd_prefix": "rm(", "library": "var_list.r", "varRefreshCmd": "cat(var_dic_list()) " } }, "types_to_exclude": [ "module", "function", "builtin_function_or_method", "instance", "_Feature" ], "window_display": false } }, "nbformat": 4, "nbformat_minor": 5 }