{"cells": [{"cell_type": "markdown", "metadata": {}, "source": ["# p-values\n", "\n", "Compute p-values and heavy tails estimators."]}, {"cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [{"data": {"text/html": ["<div id=\"my_id_menu_nb\">run previous cell, wait for 2 seconds</div>\n", "<script>\n", "function repeat_indent_string(n){\n", "    var a = \"\" ;\n", "    for ( ; n > 0 ; --n)\n", "        a += \"    \";\n", "    return a;\n", "}\n", "// look up into all sections and builds an automated menu //\n", "var update_menu_string = function(begin, lfirst, llast, sformat, send, keep_item, begin_format, end_format) {\n", "    var anchors = document.getElementsByClassName(\"section\");\n", "    if (anchors.length == 0) {\n", "        anchors = document.getElementsByClassName(\"text_cell_render rendered_html\");\n", "    }\n", "    var i,t;\n", "    var text_menu = begin;\n", "    var text_memo = \"<pre>\\nlength:\" + anchors.length + \"\\n\";\n", "    var ind = \"\";\n", "    var memo_level = 1;\n", "    var href;\n", "    var tags = [];\n", "    var main_item = 0;\n", "    var format_open = 0;\n", "    for (i = 0; i <= llast; i++)\n", "        tags.push(\"h\" + i);\n", "\n", "    for (i = 0; i < anchors.length; i++) {\n", "        text_memo += \"**\" + anchors[i].id + \"--\\n\";\n", "\n", "        var child = null;\n", "        for(t = 0; t < tags.length; t++) {\n", "            var r = anchors[i].getElementsByTagName(tags[t]);\n", "            if (r.length > 0) {\n", "child = r[0];\n", "break;\n", "            }\n", "        }\n", "        if (child == null) {\n", "            text_memo += \"null\\n\";\n", "            continue;\n", "        }\n", "        if (anchors[i].hasAttribute(\"id\")) {\n", "            // when converted in RST\n", "            href = anchors[i].id;\n", "            text_memo += \"#1-\" + href;\n", "            // passer \u00e0 child suivant (le chercher)\n", "        }\n", "        else if (child.hasAttribute(\"id\")) {\n", "            // in a notebook\n", "            href = child.id;\n", "            text_memo += \"#2-\" + href;\n", "        }\n", "        else {\n", "            text_memo += \"#3-\" + \"*\" + \"\\n\";\n", "            continue;\n", "        }\n", "        var title = child.textContent;\n", "        var level = parseInt(child.tagName.substring(1,2));\n", "\n", "        text_memo += \"--\" + level + \"?\" + lfirst + \"--\" + title + \"\\n\";\n", "\n", "        if ((level < lfirst) || (level > llast)) {\n", "            continue ;\n", "        }\n", "        if (title.endsWith('\u00b6')) {\n", "            title = title.substring(0,title.length-1).replace(\"<\", \"&lt;\")\n", "         .replace(\">\", \"&gt;\").replace(\"&\", \"&amp;\");\n", "        }\n", "        if (title.length == 0) {\n", "            continue;\n", "        }\n", "\n", "        while (level < memo_level) {\n", "            text_menu += end_format + \"</ul>\\n\";\n", "            format_open -= 1;\n", "            memo_level -= 1;\n", "        }\n", "        if (level == lfirst) {\n", "            main_item += 1;\n", "        }\n", "        if (keep_item != -1 && main_item != keep_item + 1) {\n", "            // alert(main_item + \" - \" + level + \" - \" + keep_item);\n", "            continue;\n", "        }\n", "        while (level > memo_level) {\n", "            text_menu += \"<ul>\\n\";\n", "            memo_level += 1;\n", "        }\n", "        text_menu += repeat_indent_string(level-2);\n", "        text_menu += begin_format + sformat.replace(\"__HREF__\", href).replace(\"__TITLE__\", title);\n", "        format_open += 1;\n", "    }\n", "    while (1 < memo_level) {\n", "        text_menu += end_format + \"</ul>\\n\";\n", "        memo_level -= 1;\n", "        format_open -= 1;\n", "    }\n", "    text_menu += send;\n", "    //text_menu += \"\\n\" + text_memo;\n", "\n", "    while (format_open > 0) {\n", "        text_menu += end_format;\n", "        format_open -= 1;\n", "    }\n", "    return text_menu;\n", "};\n", "var update_menu = function() {\n", "    var sbegin = \"\";\n", "    var sformat = '<a href=\"#__HREF__\">__TITLE__</a>';\n", "    var send = \"\";\n", "    var begin_format = '<li>';\n", "    var end_format = '</li>';\n", "    var keep_item = -1;\n", "    var text_menu = update_menu_string(sbegin, 2, 4, sformat, send, keep_item,\n", "       begin_format, end_format);\n", "    var menu = document.getElementById(\"my_id_menu_nb\");\n", "    menu.innerHTML=text_menu;\n", "};\n", "window.setTimeout(update_menu,2000);\n", "            </script>"], "text/plain": ["<IPython.core.display.HTML object>"]}, "execution_count": 2, "metadata": {}, "output_type": "execute_result"}], "source": ["from jyquickhelper import add_notebook_menu\n", "add_notebook_menu()"]}, {"cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": ["%matplotlib inline"]}, {"cell_type": "markdown", "metadata": {}, "source": ["## p-value table"]}, {"cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [{"name": "stdout", "output_type": "stream", "text": ["norm.ppf(0.025) -1.9599639845400545\n"]}, {"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.200</th>\n", "      <th>-0.100</th>\n", "      <th>-0.020</th>\n", "      <th>-0.010</th>\n", "      <th>-0.002</th>\n", "      <th>-0.001</th>\n", "      <th>0.001</th>\n", "      <th>0.002</th>\n", "      <th>0.010</th>\n", "      <th>0.020</th>\n", "      <th>0.100</th>\n", "      <th>0.200</th>\n", "    </tr>\n", "  </thead>\n", "  <tbody>\n", "    <tr>\n", "      <th>0.001</th>\n", "      <td>NaN</td>\n", "      <td>NaN</td>\n", "      <td>NaN</td>\n", "      <td>NaN</td>\n", "      <td>NaN</td>\n", "      <td>7676</td>\n", "      <td>7676</td>\n", "      <td>1919</td>\n", "      <td>77</td>\n", "      <td>20</td>\n", "      <td>1.0</td>\n", "      <td>1.0</td>\n", "    </tr>\n", "    <tr>\n", "      <th>0.002</th>\n", "      <td>NaN</td>\n", "      <td>NaN</td>\n", "      <td>NaN</td>\n", "      <td>NaN</td>\n", "      <td>3834.0</td>\n", "      <td>15336</td>\n", "      <td>15336</td>\n", "      <td>3834</td>\n", "      <td>154</td>\n", "      <td>39</td>\n", "      <td>2.0</td>\n", "      <td>1.0</td>\n", "    </tr>\n", "    <tr>\n", "      <th>0.050</th>\n", "      <td>NaN</td>\n", "      <td>NaN</td>\n", "      <td>913.0</td>\n", "      <td>3650.0</td>\n", "      <td>91235.0</td>\n", "      <td>364939</td>\n", "      <td>364939</td>\n", "      <td>91235</td>\n", "      <td>3650</td>\n", "      <td>913</td>\n", "      <td>37.0</td>\n", "      <td>10.0</td>\n", "    </tr>\n", "    <tr>\n", "      <th>0.100</th>\n", "      <td>NaN</td>\n", "      <td>70.0</td>\n", "      <td>1729.0</td>\n", "      <td>6915.0</td>\n", "      <td>172866.0</td>\n", "      <td>691463</td>\n", "      <td>691463</td>\n", "      <td>172866</td>\n", "      <td>6915</td>\n", "      <td>1729</td>\n", "      <td>70.0</td>\n", "      <td>18.0</td>\n", "    </tr>\n", "    <tr>\n", "      <th>0.150</th>\n", "      <td>NaN</td>\n", "      <td>98.0</td>\n", "      <td>2449.0</td>\n", "      <td>9796.0</td>\n", "      <td>244893.0</td>\n", "      <td>979572</td>\n", "      <td>979572</td>\n", "      <td>244893</td>\n", "      <td>9796</td>\n", "      <td>2449</td>\n", "      <td>98.0</td>\n", "      <td>25.0</td>\n", "    </tr>\n", "    <tr>\n", "      <th>0.200</th>\n", "      <td>31.0</td>\n", "      <td>123.0</td>\n", "      <td>3074.0</td>\n", "      <td>12293.0</td>\n", "      <td>307317.0</td>\n", "      <td>1229267</td>\n", "      <td>1229267</td>\n", "      <td>307317</td>\n", "      <td>12293</td>\n", "      <td>3074</td>\n", "      <td>123.0</td>\n", "      <td>31.0</td>\n", "    </tr>\n", "    <tr>\n", "      <th>0.250</th>\n", "      <td>37.0</td>\n", "      <td>145.0</td>\n", "      <td>3602.0</td>\n", "      <td>14406.0</td>\n", "      <td>360137.0</td>\n", "      <td>1440548</td>\n", "      <td>1440548</td>\n", "      <td>360137</td>\n", "      <td>14406</td>\n", "      <td>3602</td>\n", "      <td>145.0</td>\n", "      <td>37.0</td>\n", "    </tr>\n", "    <tr>\n", "      <th>0.300</th>\n", "      <td>41.0</td>\n", "      <td>162.0</td>\n", "      <td>4034.0</td>\n", "      <td>16135.0</td>\n", "      <td>403354.0</td>\n", "      <td>1613413</td>\n", "      <td>1613413</td>\n", "      <td>403354</td>\n", "      <td>16135</td>\n", "      <td>4034</td>\n", "      <td>162.0</td>\n", "      <td>41.0</td>\n", "    </tr>\n", "    <tr>\n", "      <th>0.350</th>\n", "      <td>44.0</td>\n", "      <td>175.0</td>\n", "      <td>4370.0</td>\n", "      <td>17479.0</td>\n", "      <td>436966.0</td>\n", "      <td>1747864</td>\n", "      <td>1747864</td>\n", "      <td>436966</td>\n", "      <td>17479</td>\n", "      <td>4370</td>\n", "      <td>175.0</td>\n", "      <td>44.0</td>\n", "    </tr>\n", "    <tr>\n", "      <th>0.400</th>\n", "      <td>47.0</td>\n", "      <td>185.0</td>\n", "      <td>4610.0</td>\n", "      <td>18440.0</td>\n", "      <td>460976.0</td>\n", "      <td>1843901</td>\n", "      <td>1843901</td>\n", "      <td>460976</td>\n", "      <td>18440</td>\n", "      <td>4610</td>\n", "      <td>185.0</td>\n", "      <td>47.0</td>\n", "    </tr>\n", "    <tr>\n", "      <th>0.450</th>\n", "      <td>48.0</td>\n", "      <td>191.0</td>\n", "      <td>4754.0</td>\n", "      <td>19016.0</td>\n", "      <td>475381.0</td>\n", "      <td>1901523</td>\n", "      <td>1901523</td>\n", "      <td>475381</td>\n", "      <td>19016</td>\n", "      <td>4754</td>\n", "      <td>191.0</td>\n", "      <td>48.0</td>\n", "    </tr>\n", "    <tr>\n", "      <th>0.500</th>\n", "      <td>49.0</td>\n", "      <td>193.0</td>\n", "      <td>4802.0</td>\n", "      <td>19208.0</td>\n", "      <td>480183.0</td>\n", "      <td>1920730</td>\n", "      <td>1920730</td>\n", "      <td>480183</td>\n", "      <td>19208</td>\n", "      <td>4802</td>\n", "      <td>193.0</td>\n", "      <td>49.0</td>\n", "    </tr>\n", "    <tr>\n", "      <th>0.550</th>\n", "      <td>48.0</td>\n", "      <td>191.0</td>\n", "      <td>4754.0</td>\n", "      <td>19016.0</td>\n", "      <td>475381.0</td>\n", "      <td>1901523</td>\n", "      <td>1901523</td>\n", "      <td>475381</td>\n", "      <td>19016</td>\n", "      <td>4754</td>\n", "      <td>191.0</td>\n", "      <td>48.0</td>\n", "    </tr>\n", "    <tr>\n", "      <th>0.600</th>\n", "      <td>47.0</td>\n", "      <td>185.0</td>\n", "      <td>4610.0</td>\n", "      <td>18440.0</td>\n", "      <td>460976.0</td>\n", "      <td>1843901</td>\n", "      <td>1843901</td>\n", "      <td>460976</td>\n", "      <td>18440</td>\n", "      <td>4610</td>\n", "      <td>185.0</td>\n", "      <td>47.0</td>\n", "    </tr>\n", "    <tr>\n", "      <th>0.650</th>\n", "      <td>44.0</td>\n", "      <td>175.0</td>\n", "      <td>4370.0</td>\n", "      <td>17479.0</td>\n", "      <td>436966.0</td>\n", "      <td>1747864</td>\n", "      <td>1747864</td>\n", "      <td>436966</td>\n", "      <td>17479</td>\n", "      <td>4370</td>\n", "      <td>175.0</td>\n", "      <td>44.0</td>\n", "    </tr>\n", "    <tr>\n", "      <th>0.700</th>\n", "      <td>41.0</td>\n", "      <td>162.0</td>\n", "      <td>4034.0</td>\n", "      <td>16135.0</td>\n", "      <td>403354.0</td>\n", "      <td>1613413</td>\n", "      <td>1613413</td>\n", "      <td>403354</td>\n", "      <td>16135</td>\n", "      <td>4034</td>\n", "      <td>162.0</td>\n", "      <td>41.0</td>\n", "    </tr>\n", "    <tr>\n", "      <th>0.750</th>\n", "      <td>37.0</td>\n", "      <td>145.0</td>\n", "      <td>3602.0</td>\n", "      <td>14406.0</td>\n", "      <td>360137.0</td>\n", "      <td>1440548</td>\n", "      <td>1440548</td>\n", "      <td>360137</td>\n", "      <td>14406</td>\n", "      <td>3602</td>\n", "      <td>145.0</td>\n", "      <td>37.0</td>\n", "    </tr>\n", "    <tr>\n", "      <th>0.800</th>\n", "      <td>31.0</td>\n", "      <td>123.0</td>\n", "      <td>3074.0</td>\n", "      <td>12293.0</td>\n", "      <td>307317.0</td>\n", "      <td>1229267</td>\n", "      <td>1229267</td>\n", "      <td>307317</td>\n", "      <td>12293</td>\n", "      <td>3074</td>\n", "      <td>123.0</td>\n", "      <td>31.0</td>\n", "    </tr>\n", "    <tr>\n", "      <th>0.850</th>\n", "      <td>25.0</td>\n", "      <td>98.0</td>\n", "      <td>2449.0</td>\n", "      <td>9796.0</td>\n", "      <td>244893.0</td>\n", "      <td>979572</td>\n", "      <td>979572</td>\n", "      <td>244893</td>\n", "      <td>9796</td>\n", "      <td>2449</td>\n", "      <td>98.0</td>\n", "      <td>NaN</td>\n", "    </tr>\n", "    <tr>\n", "      <th>0.900</th>\n", "      <td>18.0</td>\n", "      <td>70.0</td>\n", "      <td>1729.0</td>\n", "      <td>6915.0</td>\n", "      <td>172866.0</td>\n", "      <td>691463</td>\n", "      <td>691463</td>\n", "      <td>172866</td>\n", "      <td>6915</td>\n", "      <td>1729</td>\n", "      <td>70.0</td>\n", "      <td>NaN</td>\n", "    </tr>\n", "    <tr>\n", "      <th>0.950</th>\n", "      <td>10.0</td>\n", "      <td>37.0</td>\n", "      <td>913.0</td>\n", "      <td>3650.0</td>\n", "      <td>91235.0</td>\n", "      <td>364939</td>\n", "      <td>364939</td>\n", "      <td>91235</td>\n", "      <td>3650</td>\n", "      <td>913</td>\n", "      <td>NaN</td>\n", "      <td>NaN</td>\n", "    </tr>\n", "  </tbody>\n", "</table>\n", "</div>"], "text/plain": ["       -0.200  -0.100  -0.020   -0.010    -0.002   -0.001    0.001   0.002  \\\n", "0.001     NaN     NaN     NaN      NaN       NaN     7676     7676    1919   \n", "0.002     NaN     NaN     NaN      NaN    3834.0    15336    15336    3834   \n", "0.050     NaN     NaN   913.0   3650.0   91235.0   364939   364939   91235   \n", "0.100     NaN    70.0  1729.0   6915.0  172866.0   691463   691463  172866   \n", "0.150     NaN    98.0  2449.0   9796.0  244893.0   979572   979572  244893   \n", "0.200    31.0   123.0  3074.0  12293.0  307317.0  1229267  1229267  307317   \n", "0.250    37.0   145.0  3602.0  14406.0  360137.0  1440548  1440548  360137   \n", "0.300    41.0   162.0  4034.0  16135.0  403354.0  1613413  1613413  403354   \n", "0.350    44.0   175.0  4370.0  17479.0  436966.0  1747864  1747864  436966   \n", "0.400    47.0   185.0  4610.0  18440.0  460976.0  1843901  1843901  460976   \n", "0.450    48.0   191.0  4754.0  19016.0  475381.0  1901523  1901523  475381   \n", "0.500    49.0   193.0  4802.0  19208.0  480183.0  1920730  1920730  480183   \n", "0.550    48.0   191.0  4754.0  19016.0  475381.0  1901523  1901523  475381   \n", "0.600    47.0   185.0  4610.0  18440.0  460976.0  1843901  1843901  460976   \n", "0.650    44.0   175.0  4370.0  17479.0  436966.0  1747864  1747864  436966   \n", "0.700    41.0   162.0  4034.0  16135.0  403354.0  1613413  1613413  403354   \n", "0.750    37.0   145.0  3602.0  14406.0  360137.0  1440548  1440548  360137   \n", "0.800    31.0   123.0  3074.0  12293.0  307317.0  1229267  1229267  307317   \n", "0.850    25.0    98.0  2449.0   9796.0  244893.0   979572   979572  244893   \n", "0.900    18.0    70.0  1729.0   6915.0  172866.0   691463   691463  172866   \n", "0.950    10.0    37.0   913.0   3650.0   91235.0   364939   364939   91235   \n", "\n", "        0.010   0.020   0.100   0.200  \n", "0.001      77      20     1.0     1.0  \n", "0.002     154      39     2.0     1.0  \n", "0.050    3650     913    37.0    10.0  \n", "0.100    6915    1729    70.0    18.0  \n", "0.150    9796    2449    98.0    25.0  \n", "0.200   12293    3074   123.0    31.0  \n", "0.250   14406    3602   145.0    37.0  \n", "0.300   16135    4034   162.0    41.0  \n", "0.350   17479    4370   175.0    44.0  \n", "0.400   18440    4610   185.0    47.0  \n", "0.450   19016    4754   191.0    48.0  \n", "0.500   19208    4802   193.0    49.0  \n", "0.550   19016    4754   191.0    48.0  \n", "0.600   18440    4610   185.0    47.0  \n", "0.650   17479    4370   175.0    44.0  \n", "0.700   16135    4034   162.0    41.0  \n", "0.750   14406    3602   145.0    37.0  \n", "0.800   12293    3074   123.0    31.0  \n", "0.850    9796    2449    98.0     NaN  \n", "0.900    6915    1729    70.0     NaN  \n", "0.950    3650     913     NaN     NaN  "]}, "execution_count": 4, "metadata": {}, "output_type": "execute_result"}], "source": ["from scipy.stats import norm\n", "import pandas\n", "from pandas import DataFrame\n", "import numpy\n", "\n", "\n", "def pvalue(p, q, N):\n", "    theta   = abs(p-q)\n", "    var     = p*(1-p) \n", "    bn      = (2*N)**0.5 * theta / var**0.5\n", "    ret     = (1 - norm.cdf(bn))*2\n", "    return ret\n", "\n", "def pvalue_N(p, q, alpha):\n", "    theta   = abs(p-q)\n", "    var     = p*(1-p) \n", "    rev     = abs(norm.ppf (alpha/2))\n", "    N       = 2 * (rev * var**0.5 / theta)** 2\n", "    return int(N+1)\n", "\n", "def alphatable(ps, dps, alpha):\n", "    values = []\n", "    for p in ps :\n", "        row=[]\n", "        for dp in dps :\n", "            q = p+dp\n", "            r = pvalue_N(p,q,alpha) if 1 >= q >= 0 else numpy.nan\n", "            row.append (r)\n", "        values.append (row)\n", "    return values\n", "    \n", "def dataframe(ps,dps,table):\n", "    columns = dps\n", "    df = pandas.DataFrame(data=table, index=ps)\n", "    df.columns = dps\n", "    return df\n", "        \n", "        \n", "print (\"norm.ppf(0.025)\",norm.ppf (0.025)) # -1.9599639845400545\n", "ps  = [0.001, 0.002] + [ 0.05*i for i in range (1,20) ]\n", "dps = [ -0.2, -0.1,  -0.02, -0.01, -0.002, -0.001,\n", "         0.2,  0.1,   0.02,  0.01,  0.002,  0.001, ]\n", "dps.sort()\n", "t = alphatable(ps, dps, 0.05)\n", "dataframe (ps, dps, t)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["## p-values in 2D"]}, {"cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [{"data": {"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAD4CAYAAAAZ1BptAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2dfZAc1XXof2dmdyV21zbSLLYJoFlcISnLSQo/VCROKvHHyg7RqwLnPYcgr4wAJ4t2bT9VpZIKtlLJK/KUsp16cRRsCW1sgcxMjG3qJdZL5BAkm5cUAQe5jMHgwsggCUXESAg7FjL62vP+6G5t70zfnp6Znp7pnfOr6prp27e7T3/dc+89554rqophGIZhRFHotgCGYRhG72JKwjAMw3BiSsIwDMNwYkrCMAzDcGJKwjAMw3Ay0G0BWmFsbEzHx8e7LYZhGEau+OY3v3lMVS9qZp9cKonx8XH27dvXbTEMwzByhYgcbHYf624yDMMwnJiSMAzDMJyYkjAMwzCcmJIwDMMwnKSiJERkh4i8KCLfcWwXEfkrEdkvIo+LyH8JbVsvIs/4y/o05DEMwzDSIa2WxN3ANTHbfwO4wl+mgG0AIrIc+BPgF4GrgT8RkWUpyWQYhpFbqlUYH4dCwfutVrsjRypKQlX/GTgek+U64PPq8QhwoYhcDPw68ICqHlfVl4EHiFc2hmEYuaVahdFREPGWYhFWr65XBtUqTE3BwYOg6v1OTXVHUWRlk7gEeD60fthPc6UbhmHkhmoVxsbmC3/Xsm4dvPLK/H5zc7B3b70y2LgRTp5ceI6TJ739s25VZKUkJCJNY9LrDyAyJSL7RGTf0aNHUxXOMAyjGYKuoKA1sG4dvPRSOsc+eTL+WFm3KrJSEoeBy0LrlwJHYtLrUNVZVV2lqqsuuqipUeWGYRgLcPX3z8zAwMDC2n9tzT3cFQReayBrTp6ETZuyOVdWSmIXcKPv5fRLwI9U9QXgfuA9IrLMN1i/x08zDMNomkaFfNAttG5dfRfP6tWwbRucO7fwmLU1902b6ruC0qZUguHh+DyHDnVWhgBJY/pSEfkC8A5gDPgBnsfSIICq3ikiAnwazyh9ErhZVff5+94CfMw/1GZVvavR+VatWqUWu8kw+ouZGZid9QpxERgZgRMnvO6ec+c8g/CJE9H7Dg/D+vWwc2frBXy5DAcOeK2PTs76PDzsXSd4tglX11MgTzOIyDdVdVVT++RxjmtTEoaxOFm92jPkBkxMwJ49noLYtq29YwfKpFVEvK6l8fH5rqZ2KRTgne+E/fu9lsGKFbB5M0xOetvHxtxKolKZz5eUVpRELqPAGoaRf6pVr+smKBxHRuCppxbm2bvXUxwPPtj++dpREODJCF4hPjXVeotkaAh27EhWwMcZsJtVEK1iYTkMw2iJWrfPsbHkHjdR4wBqFUTA3r3tF/DgtSRcDA97rZa47Zs3e/8nJ73uoHJ54XHLZa92rxq/nDqVrIDv1uC5OlQ1d8tVV12lhmF0jkpFtVz2irVi0fstlbxFxPstFOqLwKEhb99GBMdOugQytLoMD6tOT3u/tdtKpXmZp6frz1UuJ7umtCmV3NdTKLR2TGCfNlnedr3Ab2UxJWEYzRNV8LsWkdYL5HK5sSzNHn96Olm+gYH6Y4cL+eAeiHSv8E9Ko2tt7ZimJLpHo6pXkje1UqmvPgTVnOnpxl/WyIiXr9FXUKl4ecPVkunp+vOHzx1cU6Ewv2+QFr7e0dF6uYIqZ6P74rqfUVW/cJU2fC2jo/PXEJxvZGShDNPT7vM1U3pEPfNyef4ZhK89rkoY3Psk54+q6tYslZHf1ZGhUwpzCnNa4KxODP0/HR4603LB38winGt4reXi8wmPN6cTxb2qpZJO82ktckZhTkVUR5eeVpjz085puXBIKxOfi2+mhJ9P8P6E36Wk72XSb7oN4u5LEkUcfUxTEt2hUmlcGIC7zRukDw05vro2qnXDwwtf2DhZo86T5LraXaJkjOoXSHMJK4qo89XKFPXM05SxtuD3z79AD8lZrzDkOa2wNvI4FdZqgVMRm+Y6/hiDpcxzDa+1wlod5kTdJa9cGcjqLRP8Y/T7Eve9ZPleur5pkejKSBPEidaqDjIl0S3CNdlmC4NG6al8teV5WZvtDM6sZMlYxmKx8fniqmspyDjNHfM1Y87qKD9SOKcFzup8QRm9+zAnIhVFmee6+hiHeNWpwGqXCmu1zHMqnNNy8fn5gi/Jve3k99Lse+mSRaTl0rxSiRerVUxJdItufpVJFpF5WdtplSw2GRudLyxTLY59wgV/uLsn+F/iRZ3mDl8htFe7j6qxC+cye2RDvKoj56/Du7akCqInnn8WcrXYLxSnJ8P1m2ZpRUmYC2zWuPzw4vzz2iVw8K7930tkLWP4frvOF5F+PuyDnmOAM8xwBwBV1jLKj9jGhzjHAF7sSm+Zo3j+/0tcxDY+xAleS3R8y+Qcol6+FcTFatCmji/MAXMUOQvMUeIoJY4izFHmADu4mRO8DqWAUuAYr2eSLzR1jnnBm3z+nfxewiSRK06WFkfdxYXcmJpq6ZCt06xW6YUlty0Js0kklzEjm8T0tGqxcE4DA+g0d3jdIHJQhbkFNshoDxuv37y2fz2LJaolEWeTmOABLS/9D4Vzda2dhUu83eP80mqXj8P+kvj596JNIu5aW8Dl/joy0tLhzoN1N3WJuBfNvJu64t3k9Xcf8Pq75aBWeP95GSoTn4tpzs9pgdORZYW7TMzOKHxeJodNQon2bpoeuXuhtqtV/iMjC+/jyIjq0qUL8wTvX+AhFveMopZmPLnCeeK8BGvf17BnWdTSKe+muOtukkolWv8NDrbvNNWKkrDYTWkgMd0GOby/vU4QzuHgwfl4POXyfMyb2vg/sDBoWishFcrl9OL1tErUtRo9QoplwOjowomJAkolOHasSblqsNhNRu6pjecThEIIlEKhsDB+fxCuIQjnfNdd9QoCFsbfbyXmzqFD7QeIa0QQxTR8jaUSbNliCqFfmJmJVhAAx+MmiO4g1pJIA2tJxBJV8y+VvG3Hjy9UBrW1/MFB7/aePt2+HMFjauWRlMuwZk10JNKJCXj4YbfyqS30r78edu+Ojvpp5JSUyoBi0T2JUSuhwWuxloTRdarVhTHwR0a8Av7MGW89qImHo1sGrYALLqgvaIP90iBwTmnUbVTbWgmCuwUFeTCnQbHoyb11a3QLyAr+PqL2pQmnN0HcLHdBRSpzmjVi9MKSK8N1TomL/jAyUh/5IrAHdtrhpJ2lUmnsODMxka/4PkaPMDHhfqES0qkBdGHolncT3oxzTwP7gdsitn8KeMxfvgf8MLTtXGjbriTnMyWRHkmcppIu3R4DNTyc7Ft1hVwyZWC0jMtntVRKfIi4wA3tur4GdEVJAEXg+8CbgCHg28DKmPwfAXaE1k80e86eUxKusQStxvNNiXBroFj0Csrwu9xrtf5Sqb6WPzgYL2dw62vHM4Svu80QOobRmDYrio1aEWlVYLqlJN4G3B9a/yjw0Zj8/wq8O7SefyXR5gvSDq6ukaShlXtlCcYiRF2P1fyNnqfNMqDRcJP0xGxeSaRhuL4EeD60fhj4xaiMIlIGLge+FkpeKiL7gLPAx1X17xz7TgFTACt6LbSEyzeyw6EDgtm9AmNvYACG+TEBvUCxCBde6Bmr47ybAkNvlMHXjMDGYibOmSL4VrpFGkoiyvdLHXlvAO5T1XCJukJVj4jIm4CvicgTqvr9ugOqzgKz4LnAtit0qric5zvpVI/nTVPrDRSMB+jwqSMZHob16+FLX5r3XjI/f8NojMs5Crzvp5ukEeDvMHBZaP1S4Igj7w2wMAKYqh7xf58FHgTemoJM2dKNoH24ax/BwK9OMTIyX7sJz+87O+u5gx47Nt9QPnbMFIRhxFGtxru+dvv7SUNJPApcISKXi8gQniLYVZtJRH4WWAY8HEpbJiJL/P9jwK8AjunQe5gutCRmZtzbVqxIHinSNQaoVHJP6n7ixLwiOHvW+z1woPsvs2HkkSASQBTlcnZyuGhbSajqWeDDwP3Ad4EvqeqTInK7iFwbyroWuNc3ngS8GdgnIt8Gvo5nk8ifkugCd97p3rZ5s1ejn56er+kXi97I4HD/ZqkE99wTrQisBWAY2RAXFrxrA+hCWFiONMg4LEe1CuvWZXpKwzDiaKMMGB+P7jpOI6BfLa2E5bBJh3LIxo3ubVnNxWIYRjqsWVOvY4aHu2+wDjAlkTOq1YVxj2rJfNYqwzBaZmbG6zoONzZEPC/BXunuNSWRM+JaEaOjni3CMIzep1qtVxDgre/e3R2ZojAlkQYuF4QOuCbEtSLijNmGYXQQ14i3mJFwmza5zRVxxuysMSWRBps3e52IYYL40hnSK81TwzAaE6cIeimohCmJNJic9EaSlcteh2IwsizlUrtaTfVwhmGkhWvauJjp5FyKQKQ3XF8DTEnkiDh7RLfjuxhGX7N8eXPpRHs1icCGDb3VK2Az06VBteq5I4QnXF6/3vuf4tOOs0f0irucYRiNcXk1bdjQe84nNpguDUZHo2cvHxnxYlikhE2lbRg9SqEQ/RGK1AVmqlbhAx+Izp7GPNZx2GC6bhGlIOLSWyDOHmFdTYbRZVwGhoj0vHg1BZiSyAlxQcCsq8kwukwTHo558WoKMCWRBgXHbXSlt0DcpCS9ZOQyjL6kCQ/HvHg1BZiSSINbb20uvUlWr3Zv64VQwoZhJCcvXk0BpiTSYOtWLw53mImJ1NwU9u51b+vFmodh9B3VKtx4o9fkV/V+b7yxzphYrcLOnfnwagow76Y0qFbhppu8GXgCBgbg7rtTqRqYV5Nh9DgJPRzHxqJd2Tvt1RRg3k3dYsOGhQoCvPUNG9o+dNwMdIZh9AgJPBzjIjj3oldTQCpKQkSuEZGnRWS/iNwWsf0mETkqIo/5y++Etq0XkWf8ZX0a8mSOayxECmMktm93bxsaavvwhmFkRJyHYi96NQW0PeJaRIrAZ4B3A4eBR0VkV8Q0pF9U1Q/X7Lsc+BNgFaDAN/19X25XrsVC3ATpO3ZkJ4dhGO0R56HYy7bFNFoSVwP7VfVZVT0N3Atcl3DfXwceUNXjvmJ4ALgmBZmyxWU0iDMmJKBRQL9e9IQwDKOeuG7jQqG3v+U0lMQlwPOh9cN+Wi3/XUQeF5H7ROSyJvdFRKZEZJ+I7Dt69GgKYqeIy/bQpk0iLqDfyEhbhzYMI00azCcxO+veNa63oBdIQ0lEVZdrfW7+LzCuqr8A7AF2NrGvl6g6q6qrVHXVRRdd1LKwHWHrVpienp9gulj01tv0aYsL6BdnqzAMI2NefTU2PYj9GUWvj3VKQ0kcBi4LrV8KHAlnUNWXVPWUv/rXwFVJ980NW7d6Hk2q3m+HnZ57uXlqGH1HjHdTIw/FXrZHQDpK4lHgChG5XESGgBuAXeEMInJxaPVa4Lv+//uB94jIMhFZBrzHT+t74kZZW0A/w8gPcV1N0PsVvra9m1T1rIh8GK9wLwI7VPVJEbkd2Kequ4D/ISLXAmeB48BN/r7HReRP8RQNwO2q6p7KqY+IG2VtAf0MIz/kuasJbMR1elSrniP0oUOe0/PmzW1VEWyUtWHkiJgPVqLNrABUKtm2JFoZcW0z06VBtQpTU3DypLd+8KC3Dr3fljQMo6OIuCt2eSgeLCxHGmzaNK8gAk6ejB9i2QDXaGobZW0YPUiMC2zeW/6mJNLANZQybohlA1yjqW2UtWH0IFu2wODgwrTBQdiy5bxnfC2u9F7DlEQadOAtmJz0+ivDc5hk3X9pGEZCJifhrrsWfrB33UWVSZYujd4l6JHudUxJpIHLfSHOrcFBtQrj495Q/U2bPPv33JwXRtgUhGH0MJOT3ofqf7BVJrnxxvohFIVCKmNtM8OURBq4/Nia9G+rVuHmmxfOW3LzzY1jOBmG0QOEa3jj49yy/mxkyI0LLsiPggBTEumwZk1z6Q42boQzZxamnTkTH8PJMIweoFqFW245X8OrHvxlTp+L7m52Dc7uVWycRBq4ppsqleDYscSHsbERhpFTasqAcZ7jIOPO7N36nm1mum7hisQXF6GvButSMowcU/OtH8I9i1AhZ6VuzsRdvNxyS7clMAwjLZbjriDeemuGgqSAKYke4fRp9zYL6GcY+aHKWo6zLHLbxES+jNZgSiIXWEA/w+hxQjW5DdyJRkQ8WrIE9uzJUqh0MCXRA9g0pYaRc7ZsOR8z5wSvicxy6lRkcs9jAf56gLgQTxaryTBygF+Tq278BjHmiFxiLYk0cE04nXAi6kOH3NssVpNh5ISHHmL9S39O9KzM8S7uvYwpiTSImbowCSsc3nKlknU1GUYumJmhuu2HnMPd9N+wIUN5UiQVJSEi14jI0yKyX0Rui9j+eyLylIg8LiJ7RaQc2nZORB7zl121++aCNgP8bd4Mw8ML04aHzWBtGLlh2zY2cCeuVgTkz6spoG0lISJF4DPAbwArgbUisrIm27eAVar6C8B9wCdD236iqlf6y7XtytMV2gzwNznpzYMbDiA5O2utCMPIEy6Ddd5Jw3B9NbBfVZ8FEJF7geuAp4IMqvr1UP5HgHUpnLd3KBajFUITocInJ00pGEZeqbI2dvvEREaCdIA0upsuAZ4PrR/201x8EPhqaH2piOwTkUdE5L2unURkys+37+jRo+1JnDZttCRqAkdaeA7DyCGb+DPiupryOD4iII2WRNSdiQxfJSLrgFXA20PJK1T1iIi8CfiaiDyhqt+vO6DqLDALXoC/9sVOkVLJHeAvBpsa2zAWBwdxTwuQV6+mgDRaEoeBy0LrlwJHajOJyGpgE3Ctqp4fVqKqR/zfZ4EHgbemIFO2vPpqc+k+GzemPjW2YRg9hebWqykgDSXxKHCFiFwuIkPADcACLyUReSuwHU9BvBhKXyYiS/z/Y8CvELJl5IYWXGCrVXeQ2LhxE4Zh5Iu8ejUFtN3dpKpnReTDwP1AEdihqk+KyO3APlXdBfw5MAp8Wby21yHfk+nNwHYRmcNTWB9X1fwpiRaIay24xk0YhtGbCHMo9Y4qRc6R98AWqUivqruB3TVpfxz6v9qx378CP5+GDF2lBZvEwYPuw23enIJMhmFkwswMaGSnjDLFduBDWYuUKjbiOg1Cwb3OMzTkHA03M+M+lI2yNoz8UK3CnXdClP/OKD9mKx/OXKa0MSWRBpOTXpCl8Gi4HTucpb33UkVjo6wNIz9s2uSeivQVRr2yIOfYHNddwOayNozFQaHg/mbLcogD9/xLT3UN2BzXhmEYGbJ8eXS6MMfmd+3pKQXRKqYkMibOHjE6mp0chmG0R7UKP/5x1JY5NrCVyYc/sihCKFh3U8bEdTVVKoui4mEYfcH4eLSXYomjHOP13kq5DAcOZClWLNbd1OPYNKWGsXhwDXo9TqlxphxhSiJDbr212xIYhpEWLnvECkKKYRGMjDUlkRYJwrnGTVSX9yBghtFPzMxEj58d5FU28zF/ZXBRjIzN93jxXqFahZtvhjNnvPWDB711SNyHlPcgYIbRL8wPoKvntfyYSb7grSySmp8ZrtNgbMwdluPYMcB7sdbFTLWUw8dgGH2Jy2ANnuvrXDiGkxmuDcAdzjWUvnGje/fp6ZTlMQyjY8TZohfYIyA+SFtOMCWRAXFhwSH/oYQNo59w2aKFuXl7REATUxj3KqYk0mDp0tj0uFbEIgjtYhh9xZo19eYGCQbQBfaIgARTGPc6piTS4NSp2PS4VsQicH4wjL6hWoWdOxfaEEVgA9vYykfqd1gExmtTEmngsjonsEbbADrDyA9RUw6rwm7+a/QOi8AjJRUlISLXiMjTIrJfRG6L2L5ERL7ob/+GiIyHtn3UT39aRH49DXkyp+C4ja50wzByR+yUw+R/0JyLtksxESkCnwF+A1gJrBWRlTXZPgi8rKo/DXwK+IS/70q8ObHfAlwDbPWPly8GHMNNBgZiQ3HETFxnGEaPETvlcK1X0yIijaru1cB+VX1WVU8D9wLX1eS5Dtjp/78PmBBvsuvrgHtV9ZSqPgfs94+XL06fdqbHGa1tgiHDyA+xUw7j0CCLoDchjSu4BHg+tH7YT4vMo6pngR8BpYT7AiAiUyKyT0T2HT16NAWxO0+VtbFGa7NHGEY+aNQjMMnfRG+cm+uMQBmShpKIMt/XWmtceZLs6yWqzqrqKlVdddFFFzUpYodx9BvdynbnLub6ahj5oWGPgKvveBH0KaehJA4Dl4XWLwWOuPKIyADwOuB4wn17ny1bIhynxZvj1oG5vhpGfujnHoE0lMSjwBUicrmIDOEZonfV5NkFrPf/vw/4mnpBo3YBN/jeT5cDVwD/loJM2fLQQ/Wubg1c3xb7i2UYi4W42STPkyA0T15pOwqsqp4VkQ8D9wNFYIeqPikitwP7VHUX8DngHhHZj9eCuMHf90kR+RLwFHAW+JCq5m+I4uxsXVKVtc7si8CWZRh9w3Z3r/Fi6E1qSCqhwlV1N7C7Ju2PQ/9fBX7Lse9mIN+dLzVD76us5SZ2Em1yscmHDCNPxNme+8FD0eq0HWATf8ZZBp3bLaCfYeSDRl1N/dBtbEqiAxxcxKMvDaOfiOtqGhnJTo5uYkoiYxZB5GDD6BviupoWKBBzgTViWbLk/N8Z7sBliwCYmspAHsMw2iZuAB3UdDVdf310Jld6jrDpS9OgUDjv8jrAGc7F+APk8HYbRl/impU4YMG3nGAK417Api/tFqG35Rzu/qRF0PI0jL4hTkHUfcuLeJyEKYkUiRsbAf3hLmcY/UA/fcumJFJkI1tw2SMmJvrDXc4wFgujjqg6IyP99S2bkkgDv+35EmPOLHv2ZCWMYRjtMjMDJ07Upw8MONxiXdOU2vSlBgDXX9+wq8kwjHxQrcKdd0Zve93rHK2INqYw7nVSCcvR9+zezSYexNXVZAZrw8gPmza5y/bjxx07iUTvZC0JA4CDB2PnuO0nI5dh5J24GehWuD7zRdySMCWRBiIsJ9rVrd+MXIaRd+Iq//04D4wpiTSIqS0sXZqhHIZhtEW1Gl/578cKnymJlDhOtOHB2YdpGEbPETdNab9OOWxKIg2WLmUFhyI3OfswDcPoOeIGSPdjVxO0qSREZLmIPCAiz/i/yyLyXCkiD4vIkyLyuIj8dmjb3SLynIg85i9XtiNP1zh1is18jGFeWZA8zCt9+2IZxmKjH7uaoP2WxG3AXlW9Atjrr9dyErhRVd8CXAP8pYhcGNr+B6p6pb881qY8XWFG/4r1fJ6TDAMKzFHmALP8bt++WIaRNxpFfe1X2lUS1wE7/f87gffWZlDV76nqM/7/I8CLwEVtnrdnWL0atvEhP/KrnF/W8PdMFr/UZekMw0jKpk3ubf081qldJfEGVX0BwP99fVxmEbkaGAK+H0re7HdDfUpEljh2RUSmRGSfiOw7evRom2KnQ7UKe/dC/SA6YZYNNnmEYeSIQ9FmRaC/xzo1nE9CRPYAb4zYtAnYqaoXhvK+rKp1dgl/28XAg8B6VX0klPYfeIpjFvi+qt7eSOhemU9ifDxu4I2imv/RlobRL7Q1JcTAAJw7V59eLMLZs6nIlwatzCfRMCyHqq6OOeEPRORiVX3BL/BfdOR7LfAPwB8FCsI/9gv+31Michfw+80I323iRmYWC0rcDHWGYfQO1Sq8/HJ9+tBQwlZElIKIS88R7XY37QLW+//XA1+pzSAiQ8DfAp9X1S/XbLvY/xU8e8Z32pQnM+KNXMoUf22WMMPICbfeGj2f9eBg/3o1BbSrJD4OvFtEngHe7a8jIqtE5LN+nuuBXwNuinB1rYrIE8ATwBjwv9qUJzPcg26UCf6JrXMb4i1hhmH0BNUqvPJK9DZXej9hc1y3iDu+i6KB7hWJrp4YhtEzjI7GK4NEReQitknYiOtOsnx5tyUwDKMBcQoisevrz/5sc+k5wpSEYRh9SyOzYWLX16efbi49R5iSaIG4F6tAqHspLhCMYRhdJy6gn0gTRmvzbjLCxBmtb2VblqIYhtEGcfW4DRuaOFDBUZS60nNE/q8gY6rVuBdL2cpHshTHMIwOsXVrE5kvuKC59BxhSqJJ4rxay45w4YZh9CYjI9HpTcdqWsQ+tKYkmiQuDMdmPpalKIZhtEG1CmfO1KcXCi3EaioWm0vPEaYkmqBadY+PKHGMSb6QrUCGYbTMxo1w+nR9+rJlLYyyNsO1AV5XU9TAGmGOLURYs+NmVDcMo2vE2RZbmnLYNbfpIpjz1JREE7i6mhSJbkXkcDS7YfQDca6vLU05vGZNc+k5wpREE7gaBkXy36Q0jH4i9bms77mnufQcYUoiITMz7obBObuNhpEbGo2ybinq64kTzaXnCCvdErItZoycub4aRn6waUqbw5REAhrVPMz11TDyQ9xkYf08TakLUxIJiDNyAeb6ahg5ITbuWsEmGIrClEQC4oxcrhGbhmH0HnFdTW1N/eLyalkEbvBtKQkRWS4iD4jIM/7vMke+c6FZ6XaF0i8XkW/4+3/Rn+o0V2zfjltTmAYxjJ4irquprSENLq+WReAG325L4jZgr6peAez116P4iape6S/XhtI/AXzK3/9l4INtypM6cc3TkRG/ebqIaxGGsViI+5ZFWnR97QPaVRLXATv9/zuB9ybdUUQEeBdwXyv7Z8Wtt7q3bd/u/1nE7m+GsViIsy1u2GD2CBftKok3qOoLAP7v6x35lorIPhF5REQCRVACfqiqwQSwh4FLXCcSkSn/GPuOHj3aptjJiQviaC+VYeSHONtiU2HB+4yBRhlEZA/wxohNMSagOlao6hEReRPwNRF5AvjPiHzODjxVnQVmAVatWpVJR18j19fzjIxEaxOzSRhGf1AqRWuhRTDwomFLQlVXq+rPRSxfAX4gIhcD+L8vOo5xxP99FngQeCtwDLhQRAJFdSlwpO0rSpG45umCCaeWLo3O5Eo3DCNT4ip8qZTj11/fXHqOaLe7aRew3v+/HvhKbQYRWSYiS/z/Y8CvAE+pqgJfB94Xt3+3iJ+BrsZW4cpoc1wbRk8QV+FLZQDd5z/fXHqOEG3DRUtESsCXgBXAIeC3VPW4iKwCNqjq74jILwPbgTk8pfSXqiihtD8AABMYSURBVPo5f/83AfcCy4FvAetU9VSj865atUr37dvXstxJGB93u8uNjNTYpOO8mBaBC5xh5J2Of6I5KQNE5JuquqqZfRraJOJQ1ZeAiYj0fcDv+P//Ffh5x/7PAle3I0OnOBQTjum8V5NhGD1PYtuiEYmNuHYwPBydXiqZV5Nh5AkL6NcepiQimJmJdlZqae5bwzC6SiYB/YYcwSJc6TnClEQEs7PubdaKMIz80MirKbXv2ZREf+Gau7ytAGCGYWROx72aAhZx1AVTEhEUHHelWMxWDsMw2iPOC916BZJhSqKGatXtzTY1la0shmEY3caURA2bNkV3N42MxMR3ccUYbiv2sGEY7eLyXjKvpuSYkqjBNT7i5MmYnTZvrjdQDQ1Z7GHD6DJbtsDg4MK0wcEOeCma4bp/WLGiufTz1I6q7KFRlobRj8zMwPr1cObMfFq5DHfd1QF7xNmzzaXnCFMSNWzeXD+Qbni4QaNg06aFbyJ463GjeAzD6BgzM7BtW33X8Zo1HTJYu1wfF4FLpCmJGiYnvXES5bJnwC6XvfXYF8s1WiduFI9hGB3DFTonbgyUEU1bsZsWK5OTTdY2isVoa7f5zBpG5lSr7gq8awyU4cZaEmngevPsjTSMzInr5e1YvW0ReziakvCpVr3w4IWC92uRIw0jn8T18nZsrNOaNc2l5wjrbsJTCFNT826uBw/Ov0w2KtMw8sPMjHtb7Findtm9u7n0HNFWS0JElovIAyLyjP+7LCLPO0XksdDyqoi81992t4g8F9p2ZTvytMqmTfXjIE6eNOckw+gFZmZgYMBzJCkWYXTU+x+kjY3Ba17j/d+2LfoYIh2eB2YRO6+0OzPdJ4HjqvpxEbkNWKaqfxiTfzmwH7hUVU+KyN3A36vqfc2cN+2Z6dqeVCons1IZRi9TrXoB+YJ4SyLpfj4d/RQLhegTiPSUG2wrM9O1a5O4Dtjp/98JvLdB/vcBX1XVuPHLmeMyZplzkmE0T2DfE2m8BPa/ahVuuWVhQL40C/WO249dwi6CSmK7SuINqvoCgP/7+gb5bwC+UJO2WUQeF5FPiciSNuVpCXNOMoxkuBw8worhAx9I3ssS2P82boTTpzsjs4hFyGmHhoZrEdkDvDFiU1M99iJyMd5c1/eHkj8K/AcwBMwCfwjc7th/CpgCWNEwRkZzlMvRL3Xi2sfSpfDqq9HphpFDqlXPJnfoECxf7r3etbM1BgX8Qw/Bzp3zdr1mK88nTzaIjdYGIrBhQwYOKKOj0XNHjI52+MQZoKotL8DTwMX+/4uBp2PybgRmY7a/A88+0fC8V111laZJpaI6PKzqvd7eMjzspScivGPtYhg9wPS0arHovZIi0a9qqeS981HfQ9wSHLfXluB6MmFkJFqIkZGMBEgGsE+bLOfb7W7aBaz3/68HvhKTdy01XU1+6wIRETx7xnfalKdpghrTyZPzNohEoTgMo4cJdwuNji6MY+Sq6b/0Etx8s9f100zNPo1u2VIpecDUIF/wvZZKCyvspRJUKnDsWIbfcG0zq1F6nmhWq4QXoATsBZ7xf5f76auAz4byjQP/DhRq9v8a8ASecqgAo0nOm1ZLou0WRIC1JIwUqVRUy2Wvxl8q1VdSwy2BqNpysy2Bdpd2WxLBN1epeNcTdZ2ZtwyaJSdlAC20JNpyge0WabnAjo+7bREHDjRxIHOBNZpkZgbuvHP+9ViyxKsNx0236WJwcGH4a9d73QmGh71w3GGbBCR3Xy2XPaNy7lvtOSkDuuECm2tcH5Jr4iHDiCPK86danR/8FQwGe8tbvO6fcNlx6lRrCgLqo9K3+/6WSvXh8l35Zme9Ucy1kZPvuSdZO+LAgUWgIBY7zTY9emFJo7tpetr96pbLTR4sJ01No3lqu0DCXR613UKDgwsf/dCQ20ic9iIyL3O53PpxBgfnu35cXV493e3TLXJSBtBCd1PXC/xWljSURFw/qtkkFj/hQrBc9tbDHkDFourEhFfQRxWk09PZ9vs3WsIVmyibRLGoWijoeYUSdQwr/NsgJ2VAK0qib20SqXYh5qQ/sl9J4vM/MNDcTJOuKUS6Qa1NAhZe84oVi6Tfv5fJSRlgNomExIUBbykUR9wLsnq1F4Es3CldG48g3JG9evV8HhGvQztK4FZjm4f3Gxubly2IluYaRhtsD5aBgYUhN6OuozYqW7BtZmb+uOFrDV/zzIzzvlVX72Cs8BIiiogyJkeZWfJZxpb+eEFaVSapyvuZWvcKBw963+pLL0V7JTY7FXH7CiJpwaEMcdLP7y3C3Pn/JXmJu868n8kbiwvu1eQ64QDjzG2Y4QDjTH4gZoh0OFJe2IASpEW9p7XPZ3TUW681yNTu0+idrVajj1v7rgRyuVi9euF1DA7Oyxe897Wyhs8btbi+kbExtxyLYD6JrncdtbK0290U12c7Pd3CAScmWusjGBqq78iOWorFhf0ArfruJvWNHB5O3p8yPZ26z2WFtVrmOYVzWuSMwjkt85xWWKsV1uoQr0bsNld/e3lVS7zYke6dYuFc48c7cFaFM5GyruQxFc76cs/pEl7xZa2/5rbenXaebbAUCtH9bnHL4GDjfWrf2UqluWsaGop+55v9HoeG5vvimr2Pcdco0nP9d5hNIumNci8t0Y6lMOkS7nR2na+Rxb0ZOZM6vxeLscetsNYv/LzCsMSLOs0dWuY5lYiCsMJaHeZE9HfJiRYK/Xrl0e7+g/xEp0furitnA6UknNNy8XmtlD6iFdbqCD86f/0Fzuo0d3T+fUnj2WaxJHmvk+6f5APP+j72GKYkElCpuA13TXs1BWThwhJ2X3GdL5wnSzlDxw23AoKCsX6XhWnDnDivKLx9407XbKGfLP8ApyPyzukE/1in5CqsVRXxjN8OZaci2bk25XlJ8l4n3T+g29cUXnoMUxIJcFVW2moZ9kFLIqpFUGHtgpZEXCug4eXxnCqo0KgbpzklUeLFBjLNX8s0d/hdPXNa5Ex8jT+413HPIov3opXFWhLZ3ccew5REopvUgefZbF9qsHTBJjFduLNxQVjTb11hrQ7yk3rxeVUrE587L0/jVoB7Ec6p0rglUeLFpmwSgR0jqPGXeHG+Syiuzz9uCd/ruGfR6ntRKDSuVZtNovdtEhMTbRQqncGURAMqFffzbLmrKXzw8Kgr8D70iYmF6cHLGDjn1zrsT0wsfGFHRqI/BH+/Cu/XcvF5FeachwzSvAFR9V0q09wxXysKdgidI67gPn/fKpUErYCY+++3JBrZJCqlj2hl4nNakmO6wM4x9NdaWvKfurCl835vx/C1TU/X11gLhfnRYlEBg4J9wjc04llEbq99L4JnOzISXTCFQ7G6IovWPujwcWvfMZfs4X2D+xMVKCosT/g409MLr2tkxFuvfelq94m7j1H3Kzhu7TU2GtRRqygGBublK5WiZa39fmsX1zcStV8PKghVNSXRiNFR9/PvMSeE88SVP1GjxovF+spNo0pnsRgvQ1ylNo2RvlEVytryy1WmGIaRnFaURN8MpqtWYd069/Zeug3BQKiDB+sDpQ0Pe3FyIP56miXu+uMCxoWDIVar3iQ0cWGmSyW4/nrYvdsGehlG1rQymK7rrYJWllZaEo1akp0iXCsO19xd4zGSDDlI2ybaqCXh6iqO6hJu9noNw8gOrLsp7ua4l1Ip+XFqA5816uaMK/CjCs4khX/a3pVJCvC4QHeGYeSDVpREX3Q3NepqWrkSnnqq8XFGRrzJ2s+caZx3eBguuCA+/HOxWB8OolBo3PUVjPRPOmfA0JB3zFq5CwW49VYv1LNhGIufzGM3ichviciTIjInIs4Ti8g1IvK0iOwXkdtC6ZeLyDdE5BkR+aKIJJzAsDk2bozfnkRBgBfzJ4mCAK9fvtH8AFHxf1asiN9neNjrw9+8OXq6x0IBpqcXxvbfscMLABdOq1S885uCMAwjjoE29/8O8N+A7a4MIlIEPgO8GzgMPCoiu1T1KeATwKdU9V4RuRP4ILCtTZnqaHUyl04TFUxw8+Z6429gvI6axWvjxvnrK5Vgyxa3EdiMw4ZhNEtbLQlV/a6qPt0g29XAflV9VlVPA/cC14mIAO8C7vPz7QTe2448vUajGb6mpurTJifds3zVzuI1OelN9h5YCjKd+N0wjL6g3ZZEEi4Bng+tHwZ+ESgBP1TVs6H0S1wHEZEpYApgRaM+mRpKpexbE8PDXq0e5t1ZA4pFT0G4unomJ62wNwyjN2ioJERkD/DGiE2bVPUrCc4RNdmCxqRHoqqzwCx4husE5z3Pli1wyy2e0bldikW48EI4ftybwAbq/9f6/luBbxhGXmmoJFR1dZvnOAxcFlq/FDgCHAMuFJEBvzURpKdOUEi3qyga9fkbhmEsNrKYme5R4Arfk2kIuAHY5fvsfh14n59vPZCkZdISk5Nw6hRMTCxMHxryPH2SjCiwPn/DMPqNdl1gf1NEDgNvA/5BRO73039KRHYD+K2EDwP3A98FvqSqT/qH+EPg90RkP56N4nPtyJOEPXsWFvynTlnBbxiG4aIvBtMZhmEYXRhMZxiGYSxuTEkYhmEYTkxJGIZhGE5MSRiGYRhOcmm4FpGjwCt4Yy16nTF6X848yAj5kDMPMkI+5MyDjJAPOQMZy6p6UTM75lJJAIjIvmat9N0gD3LmQUbIh5x5kBHyIWceZIR8yNmOjNbdZBiGYTgxJWEYhmE4ybOSmO22AAnJg5x5kBHyIWceZIR8yJkHGSEfcrYsY25tEoZhGEbnyXNLwjAMw+gwpiQMwzAMJ7lREiKyXEQeEJFn/N9ljnyfFJEnReS7IvJX/jSpvSjnChH5J1/Op0RkvNdk9PO+VkT+XUQ+nZV8oXM3lFNErhSRh/1n/riI/HZGsl0jIk+LyH4RuS1i+xIR+aK//RtZPt8aORrJ+Xv++/e4iOwVkXKvyRjK9z4RURHJ3N00iYwicr1/L58Ukb/JWkZfhkbPe4WIfF1EvuU/8zUND6qquViATwK3+f9vAz4RkeeXgYeAor88DLyj1+T0tz0IvNv/PwoM95qM/vYtwN8An+7RZ/4zwBX+/58CXgAu7LBcReD7wJuAIeDbwMqaPDPAnf7/G4AvduH+JZHzncG7B0xnLWcSGf18rwH+GXgEWNVrMgJXAN8Clvnrr+/R5z0LTPv/VwIHGh03Ny0J4Dpgp/9/J/DeiDwKLMW7QUuAQeAHmUg3T0M5RWQlMKCqDwCo6glVPZmdiInuJSJyFfAG4J8ykquWhnKq6vdU9Rn//xHgRaCpEaUtcDWwX1WfVdXTwL2+rGHCst8HTGTdqiWBnKr69dC79wjeDJE9JaPPn+JVGl7NUjifJDL+LvAZVX0ZQFVfzFhGSCanAq/1/7+OBLOB5klJvEFVXwDwf19fm0FVH8ab7e4Ff7lfVb+bqZQJ5MSr/f5QRP6P3+z7cxEp9pKMIlIA/jfwBxnKVUuSe3keEbkar4Lw/Q7LdQnwfGj9sJ8WmUe9ibd+hDexVpYkkTPMB4GvdlSiehrKKCJvBS5T1b/PUrAQSe7jzwA/IyIPicgjInJNZtLNk0TO/wms8yeL2w18pNFBG85xnSUisgd4Y8SmTQn3/2ngzczXhh4QkV9T1X9OScTgPG3JiXfffxV4K3AI+CJwEynOzJeCjDPAblV9vpMV4BTkDI5zMXAPsF5V59KQLe50EWm1vuRJ8nSaxDKIyDpgFfD2jkoUceqItPMy+pWVT+F9H90iyX0cwOtyegde+fMvIvJzqvrDDssWJomca4G7VfV/i8jbgHt8OZ3fTE8pCVVd7domIj8QkYtV9QW/QIhqzv0m8IiqnvD3+SrwS3h9mb0k52HgW6r6rL/P3/lypqYkUpDxbcCvisgMns1kSEROqKrTsNglORGR1wL/APyRqj6SpnwODgOXhdYvpb7ZHuQ5LCIDeE374xnIFiVDQJSciMhqPKX8dlU9lZFsAY1kfA3wc8CDfmXljcAuEblWVbOanjLp835EVc8Az4nI03hK49FsRDwvQyM5PwhcA17Pi4gsxQv+5+wey1N30y5gvf9/PfCViDyHgLeLyICIDOLVirLubkoi56PAMhEJ+s7fBTyVgWwBDWVU1UlVXaGq48DvA59PW0EkoKGcIjIE/C2efF/OSK5HgStE5HL//Df4soYJy/4+4GvqWwszpKGcflfOduDaLvWjx8qoqj9S1TFVHfffxUd8WbOcvzjJ8/47PCcARGQMr/vp2QxlhGRyHgImAETkzXg23KOxR83aAt+G5b4E7AWe8X+X++mrgM+GrPvb8RTDU8Bf9KKc/vq7gceBJ4C7gaFekzGU/ya6492U5JmvA84Aj4WWKzOQbQ3wPTz7xyY/7Xa8Agz/4/sysB/4N+BNWd+/hHLuwXPuCO7drl6TsSbvg2Ts3ZTwPgrwF3658wRwQ48+75V4HqDf9p/3exod08JyGIZhGE7y1N1kGIZhZIwpCcMwDMOJKQnDMAzDiSkJwzAMw4kpCcMwDMOJKQnDMAzDiSkJwzAMw8n/B9Ihg6QszLbEAAAAAElFTkSuQmCC\n", "text/plain": ["<Figure size 432x288 with 1 Axes>"]}, "metadata": {"needs_background": "light"}, "output_type": "display_data"}, {"data": {"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAD4CAYAAAAKA1qZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2df5AkZ3nfv8/O3krMLhhpVmB+3SwYkfKZSky05TLGDnZOjoX+ADshji67x0kCrW42Vl2VK6kINklRpBYCKds5AnvHgk4/btoKsSsxSiwKIxkVhEIOq8L8DiCk29NFCtytgGLvgDvtPvnj7b7t7em3++2Znp6enu+nqmtmut/peWam+33e93me93lEVUEIIWQ0GRu0AIQQQgYHlQAhhIwwVAKEEDLCUAkQQsgIQyVACCEjzPigBbAxPT2tMzMzgxaDEEKGiscee+ycql7j2r60SmBmZgZra2uDFoMQQoYKEVnP0p7mIEIIGWGoBAghZIShEiCEkBGGSoAQQkYYKgFCCBlhqAQIIZfxPGBmBhgbM4+eN2iJSL+hEiBkSOh3B+15wMICsL4OqJrHhQUqgqpDJUBISfA8YHoaEDHb1NTu1/Pz/e2gl5aACxd277twwezvluh3mp6mUikbVAKEFIhtNO95wK23AhsbO23Pn9/9OkqvHXSU06ez7U/D84Bbbtn9HTY2zPekIigPVAKEFESSuWVpCbh4Mfs5u+2g49i7N9v+NJaWgEuXOvdfvOiuvOij6D+5KAEROSEi3xeRr1mOi4h8UEQeF5GviMjfz+NzCRkmkswt3Xbm3XbQcSwvA/X67n31utnfDUnfyeX7xinN+XmalvImr5nAPQBuSDj+JgDX+tsCgGM5fS4hQ0OSuaWbzryXDjqOuTlgdRVoNk0n22ya13Nz3Z0v6Tu5fN84pRlmY2O3UhgbAxYXs8s56uSiBFT1swCeTWjyFgD3qeFRAC8UkZfk8dmEFEEeDs4kc8vyMjAx4X6uycneOmgbc3PAqVPA9rZ57OX8y8vAnj2d+ycm3JRX1tmRKnDsGGcJWSnKJ/AyAE+FXp/x9+1CRBZEZE1E1s6ePVuQaIR0srgI1Gq7I3OiDs5bbsnW0SSZW+bmgBMngEZj59jk5O7XgJGp1QI2N/NXAHkzNwfcfffu79BomO/pInsvpq5glvD851MZpKKquWwAZgB8zXLsLwH8euj1wwCuSzrfddddp4QMglZL1Ywr07dmM9u5223zHhHz2G734QtUhHbb/E6u/4XLfzUKvzeANc3Qdxc1EzgD4BWh1y8H8HRBn01IB0lRJ6ur7ufJarLI09xSdebmgMOHzUwsD7j4LZ6ilMADAN7mRwn9KoAfqeozBX02IbtIWxm7teV+rjyjc0gnKyvAyZPGSQ0Ypd0LFy4YxTI+bpTL+DidyXmFiN4P4AsA/o6InBGRt4vIYRE57Dd5EMATAB4H8FEAI/6zk36TNNJPWxlbq7l9xp49+UbnkHiC2ZOqUdCBgafdNp14VjY3dxT91pZxJj/veSM8Q8hiOypyo0+AdEu7rVqv77YH1+s79mCbnVnEHHfxCTQao2FfHgZaLdWxsXz8BlX4X1FSnwAhfSM66j9yJHmkn7YydmXFROCETQ+Tk2bkGXQX587Rnl8WVlZ2ZgjtdmdEVRY2NkbPb0AlQIaaOPu+Ld9O4MR1WRkb7lhUhyMkk5j/6Ny53Saj8OI3Fydz3jmZyg6VABlq0laVhglG+nmvjCXlJRqNdfhw2jsMeeZkKjtduFUIKQ+uN2t0pD83x05/FFlZAb79beDhh5PbjVLUF2cCZKix3ayNBkf6JJ6HHkr2HeSdk6nsUAmQ0hA4eIP4bZH09ME2+/7Ro1yUReyEfQdRv8GoDRioBMhACSdmCypnATtx3GmrPGnfJ72SZRV3FesbiAkrLR+zs7O6trY2aDFIHwkqT8UVHonSbJoblJBBEUSihQMR6vXyDTpE5DFVnXVtz5kAGQieBxw65KYAgNGK1iDlxLbSfH5+uGcFVAKkcIIRFXP0kGEiaSCyvg4cPDiceYioBEjfcVnRm8SoRWuQcpI2EFEFjh8fvhkBlQDpC+FIn4MH3Vb0hglWdtLRS8pCXCRaFFVjHhqmGQEXi5HciTrQssQeNBomvJOdPikbwTW5tLQTxWbjmF9FfWWlvzLlAWcCJBfCJp9Dh7KZewAzwmq3mZiNlJsgnLTdTs9DlKU40SChEiBdYzP5uDh8uaKXDDMuVc+2toajaA2VAOmKcPZOIJvJpxQresOr1ER2qsrHxfrZVgiF909Pmy1os7iYffkzGSqCqmdJBEVrSq0IshQfKHJjUZly02xmK9YRFHLpe7HvuEru+/fvFualL00WNlyBxlahptXq3O+yBe8Nyxh9PexVTUYMlyJERRa6R8aiMgPv7G0blUC5sVXnCm+1Wh/7tbjOPq7DdhHUdseq2rVdrdbdeV1kCiuhuO/baKhOTna+rwplsYaUVsvtkoj7a/OGSoDkRlw/G5A2E+jrxW4bnTca3XfMcR21avdKpNctUEK272vbxsdVJyZ2f4/JSc4yCsJFEYT/2n5AJUByIa1Ob9Kge1dfk6RJuiWrLaqXTrgfM4EsSqgf3zc8Y+jH/zPCuJqG+vkzUwmQXLD1O9EBamz/ERwIa4Y4TdIt/R6d99Mn4Cp7+Ifux/e1fYeJCaMkqBS6xsU01M+ZMpUAyQVbvxMeoMbiYrrodT5s01CNRvc+gSTPtU3bRe304c4zcPYCOz1CsD/t94n2EP2a+bjMZvbs2TGzhb8HlUMqabdCv8xCVAIkM3F9nHUm0PhxsvnApcNK1SQOAttsVS7RQfv2DdYEEpUxLTooi0+gqK0ID2cFaLf7dxvYoBIgmchi7ahPXNL2npuTOwOXkXceQ6BRs2W7RAdNTWUzHfXq1wiPGOJmQ1X/TxxxMa3mCZUAyUSz8WOn+7vZVG037ki/mgcaNkS03d6tIERUr7gi/n/o1q8RPQ//60ScgyhygkqAJNJuqzanzqlgS5t4UoEt96mqi6Og6CueuOHq19izx10BZImHHLWZW4R+x0qEoRIgVtr779I6Nnf3zxYlEDtVdZ3XjvgNP9RElUJ4zYHrDCA6AIgbGAQO5xG7RoowDVEJkB3CN/TkpD/yj7lPI4qgjs34ezJt8QCpHuEhbDQ6yCUIoNl0DxZotQb6VYug66i7DFAJjDLtduKqWduoH75pKDARtRt3JH8GR/lENT1qKRggZAnTrXgeJc4EqAT6R7ttn7r7m20m0MSTu2/Eit14pI+4RAdlWeeQpDAqkBupiMk0lcCo4nCjtXGgwydQx6a2cWDnBhyBKTkpmLzXOYyNDfV1Gtab4ajeWi2fr5VVCbCewDBhy2sPwFt/A2bwJMawhRk8CQ8HOt4+h/uxitvQxCkIttHEKaziNszhflPZ5eTJ4aiHR4aLuTlTNSioItRo9Ha+7e0hSNJvJ6hOdvgwsLlpNBswwNoDWTRGkRtnAhEs88h263M6NaUKbNtH+OER1NhYvsMOQrqh1eo0/WTNkRSYnYbUfxDcinETnV4AZwIVYXFxp9qVCDA/31G4d/HC+zF/7A3Y3ASA3XXuLmASS3jvzo5mE7jvPjPcUAWee46jfjI4grJc4Rqjhw+bsnOuqO7UNF1fN6Xuhqhy2/Z2tv39gkqgTISL9h47lng1eDiAY/gXiHb+YU5jr6mIrTqgGo6EJBDYRYIaoysrO2ajbrhwAThyxGoyJfFQCZSFxcWdau0OmFF+QpVrAHubY+z4yXARKAZVM4DJ6j/Y2Ng9O5ifN4Oq6enSKYTJSfuxIkWlEhg0QcHzY8d2PEQOnMbexOMiwPJyr8IRMkDm5oBz58x90WoZ8yhgHqemsp1rYwO49dZSKYKPfMR+bGmpODmoBAaJ5xk75sZGelMc2BX9c7X8ILH94cOcBJAKsbJi/FiBP+v48Wz+AwC4eLHY3jWFpPvT0SCQC7koARG5QUS+JSKPi8idMcdvFpGzIvK3/vaOPD536IiGeB450uHs7XgLDmAa38c8PKxjBooxrGMGP97TuDwwitJq0edLKk407LTZdDMdnT7df9kykOT+KCxUNEsoUdwGoAbguwBeBWACwJcB7Iu0uRnAh7Kct3Ihol0smIlb3BVdQBnOElGBBZWEdE+WqnYlSX+SlFWj28X7yBgiKprBDh2HiLwewLtV9Xf81+/0lcv7Qm1uBjCrqn/get7Z2VldW1vrSbZSMTOTaY7n4Z/jkNyHLbUM92EGQEWHkxFSajzPzLDjTKwTE8CJE+b5wsLuWbiI6XubTeNMK9CWKgnxHc2m8ZNnO588pqqzru3zMAe9DMBToddn/H1R/omIfEVE/lxEXhF3IhFZEJE1EVk7e/ZsDqKVCMdpqIcDmJZzmEc7UQEAwN5k3zAho0fgTI5GFjUaRgHMzRm/QNQMGwyG19dNlF6By3aTTEJFWK/yUAJxeiw6vfgfAGZU9e8CeAjAvXEnUtVVVZ1V1dlrrrkmB9FKhK3HbjQu2zW9xh1YmLgPG9pAWvhnvc7oH0KshCOLVM3zYHSf1rOqFpq/YXnZPhsoYqCXhxI4AyA8sn85gKfDDVR1Q1V/5r/8KIDrcvjc4WJ5uTOaoV4Hjh69vGBmaeqDuHBxPPVUjYbxiTH6h5AucO1Zjx8vJKR0bs5E80UVQVEDvTyUwBcBXCsirxSRCQA3AXgg3EBEXhJ6+WYA38zhc4eLuGiGSE+eNkCp1cwsNzyoIYRkJG5AFofqzmKz8fG+zgzismgUNdDr2TEMACJyI4D/BBMpdEJVl0XkPTBe6gdE5H0wnf9zAJ4F0FLV/5N0zso5hh1I8h3X6xz9E5Ibnmd8A1kD8ocg/noQjmGo6oOq+hpV/QVVXfb3/TtVfcB//k5V/SVV/Xuq+ltpCmAoWFw0o4McRwm2AQrNP4TkTDg9Ravl/r5jxyqXl4grhrthcdFcDFtb5nVOicDjLEY0/xDSZ1ZWjCJIitUME+QlKjiKqF/kYg7qB6U2B42P7yiAMLWaWdJOCBk+AhPR6dOZ8njlaSJaXDQDwa0t050sLGQ/9UDMQSNHnAJI2k8IKT/h1NZZTUQ5ZCntk4EhFSqBbrAl7bHtJ4QMF4GJyPWe3tjouajN6mq2/XlBJdANCwuZ9ieUBiaElJVw5tJ2O91ncOFCT1lKB2VgoBLohugooVaLtQsGpQLm54e6Ch4hJFjRlUbJspS6QCUQxXXYHs1vHlEAi4um84/LY9XjgIEQMghcooiCmuAlrWYWB5VAmKDIS4/Dds8zDp0khnDAQAgJlvbaaheE0/pubJgwUsf+w5ZIrtuSy65QCYSJyy7YxbDdpTkzgBIypIQzlQaLemwOZFXg9tudTmtLL9bv/EFUAmFsw/OMw/a05swASkgFCIeUJhX2OH/eKSLEIb1YX6ASCGMbnjsO24NMEknrTMbGmAKCkMqR1kc4mpbDeuXUqWL6idFWAlEn8I03dj0fiy70iGNsDLjvPioAQiqHy9Q+xbQ8sFDyLLUoi9z6XmM4rh5pva7aanVVe7RWSy5tyvq/hFSc/fvTa4eLxL7V1h0NRY3hftH33EHT0/Hxm90U9URy1FhJf2JCSN6Ek//YiKljbEsjPyw1hocPz4tXAEDXsZvMJEEIubx+qN22F66J8Q/kFJPSFaOpBJJiOLuM3cyYSYIQUmXCoT5xRPwDPcak9MRoKoEk9dpl7KZjJglCyKgQhPrYbMWhfmhQawSAUVMCgfvdZqRvNHoK3UnJJEEIGUUch/nPe97O8yKrCY6OEginhIijXgeOHi1WJkJI9UkZ5gddU9hN+ZOfFCfe6CiBuJQQAUUtzSOEjB4pS4FzylbTNaOhBDzPPgMQcVqax5oAhJCuSVgKPMjIIAAYL+ZjBsj11wMPP2w/7uB+D6ZrgbYOIrwATh4IIb1Rr5v0QnH7i6DaM4HFxWQF4OB+9zzg0KHBTtcIIRXF83D+fHzyuaL8AtVWAsePJx9P8QN4HnDrrfbFf6wJQAjpGs/D4ts2AcSHkCYlJs2T6pqDPC85X0OzmWrLOXIEuHjRfpw1AQghXbO0hOPb34VNCRSVbaC6M4EjR5KPO6zCsGWWAFgTgBDSI6dPQxO64KKyDVRTCSTlBgKAiYlUM9DMTPJHMKKUENITKaaEohabVlMJpHlsT5ywHgr8ALaIUqDnhcWEEALvxjYE8Sbrqani5KimTyDJY9tqJfbgaX6APXu4sJgQ0hueByzc++uxKmB8PD2mJU+qOROwBdg2GqlzrCQrUrMJ3H03ZwGEkN6wJTCo1YB77im2j6meElhcjF95MTbW8xC+qJqfhJBqYzM3b29tY+5gsWkJqqcEVlftxxx68EYj235CCMmC59mzS+/FaRPa7liYPg+qpwRsK7scV14cPWrs/mHoByCE5MXSUvwSJsE2lvGunR0FpSWonhIYs3wly8oLzzPlhkXMduQI8I537E74Rz8AISQvbHErCmAO97s1zpFqRQclTZ1iVl54HnDLLcClSzv7NjaAu+4yUaTs+AkhebN3r6WoPGI6/Kuv7rs81ZoJLC3Fm30mJ2OjgpaWdiuAgIsXmRyOENIfYmvMTDyH5bF/29n4xz/uu1+gWkrANnWyFJNJmmkxORwhpB/E1pg5MY65qz7Z2biAEWm1lIBjLc+05mnHCCEkK+HCVEtLZkawq8bMs8/Gv7HPI9JclICI3CAi3xKRx0XkzpjjV4jIx/3jfyMiM3l8bgcptTzDeB6wuRl/mokJJocjhOTH4iJw8KDxBVgjQG32/z77BXpWAiJSA/BhAG8CsA/AARHZF2n2dgA/UNVXA/gTAO/v9XNjSanlGRBX2Dmg0aBTmBCSH55n0kBEw0LLUphKNCnnvssJRF4P4N2q+jv+63cCgKq+L9TmU36bL4jIOID/B+AaTfjw2dlZXVtb60m2ODwPeNvb4v3HzaaZmhFCSF5MTcUnMQDMWPVyXzQ2ZllAIJkqzIjIY6o669o+D3PQywA8FXp9xt8X20ZVnwPwIwAda3BFZEFE1kRk7ezZs91Jk1ARPsgQavs96QwmhOTJ9dfbFQAQ8T1m9GnmRR5KIG4BdFSdubSBqq6q6qyqzl5zzTXZJQnngQ4Mb7feelkRsFIYIaQoPC+5xLlIxPeYwaeZJ3kogTMAXhF6/XIAT9va+OagnwNgcYX3QFwvf/EicORIap0ZgM5gQkh+pBU3PHw44nucmwMOHdrJblCrmdd9dlDmoQS+COBaEXmliEwAuAnAA5E2DwA45D9/K4C/TvIHdI2tl9/YSHXAsFAMISRPkgadY2Mx61c9z6QrCPKfbW2Z12VfLObb+P8AwKcAfBPAf1XVr4vIe0TkzX6zuwA0RORxAH8IoCOMtN8k2fuZII4Qkidp/fbtt1t2WiwZ/aTn6KB+0VV0UIIbfqbxY6xvdNZsGxsD7ruPswBCSH7MzNhrBlxxBfDTn0Z2eh4wP28/YYZ+ehDRQeXhyitjd3s4gM1nf9axv16nAiCE5E+S5eGuu2J2DnDBQLWUQMyyaw8HsICPYkN3R6Q2GrHryAghpGuCCHXbwN3qe7RNG4I39ZFqKYGYGM8lvBcXMNmxf2qKCoAQkh+Li8aiY+vP63WL7zGp1BjQd4dltZRATIznacQH/3NhGCEkLzwPOHbMftySwcZgKzUGAK1W30er1XIMAx3O4Rk8iXXMdDRjighCSF4kOYKBFL9u0iygi/55tB3DQEftgGW8C3XsjhgqYBEeIWRE8LxkBZCKTQlYSuLmTfWUQMQvMIf7sYrb0MQpCLbRbGzSIUwIyYUgRXQSiX7dxUX7aD9YNNZnqqcElpcva1YPBzCDJ3EQbQDASczjFGaoAAghPRP4AZIsNuPjCX5dF0dCAVSr0Dxghvif/zy8Yz/EAj56OTJoHTNYwEeBjdtAHUAI6ZW0hbyNhlEA1kFn2gkKsllXzzHsMyOn4h3COIVTrQ/EFp4nhBBXkvy5ToEnSScAunIKm9OOumPYxxoair0mVosQQvpEz4P4Pi8QC1NZJbC3cSF+P04X5nAhhFQXWz/ttBB1wAvEwlRWCSwfneoMDcV5LONd5kWk6hghhGTh6FFgYmL3vokJU084kaDI+QAXiIWppBLwPLMI7wLqqOE5ANto4hRWcRvmcL9ptL5u/ggqAkKIA9HKtQBw4oSx/4uYxxMnHPrvpaWO9UwAzLqAdrtwf2XlHMOBkg3/xnWc360AwnDpMCEkhcVFM8IPd5f1epdJKHMqKG9j5B3DcUr2AiaxhPfGv2F93fzDhBASg+d1KgDA9DOZMkAvLpqFA7aB94CKnFdunYAtMZwtWgjAzoINho0SQiIk5XdzTkR5/fXJVecHmMumcjOBq6+27J+8aH5oG8eO0T9ACOkgqaN3Grx7XrICSEwx2n8qpwSsXHll+voAOooJIRFsHb2I4+A9bWXwqVMDTWZWOSWwsRG//9lnYX7opMx8mY18hJCqEY0CuvHGTiOCCHD4cErf7XnA9LS9UwIKyxSaRKWUQNL6i8vafGEh+SSsNkPIyBJEF66vGz/A+jpw773AoUO7Q0FPnkxxIXoecMstyQoASO+PCqBSIaK2wg4i5k+7rLUXF9Oz9y0vM980ISOGrQ/JHEmeVmUGMCbqn/wkw0ndGOkQUdsgXjXSn6+smEUZNkcxF5IRMpJYowuzGgjS3jAxAXzsYxlP2h8qpQRsDpzYtNxzc8ZRbMvZTf8AISOHrQ9xjgIKnAljCV1rrea4tLgYKqUEXv3qbPsxN2fmeDZHwvo6cwwRMkIsL3caCJxC+KPOBFuSyokJ42QoiQIAKqYEHnkk2/7LJKl5moYIGRnCBoLACZwawu95xnMclw8oPMBsNEo1AwiolGM4KTNr4teMSzgUhTmGCCFhPA+4/Xbg/Hl7m5zyAWVhpB3DXZPmHwDMjGB83Pyp4+PMN0TIEBKE7ouYbXq6y0m+5wG33pqsAICB5QPKApVAQOAfSFIEgZ1va8uEmFIREDI0BP12OHR/Y8OE82dSBIuLwPw8cPFicrsB5gPKQqWUgM0clFbKcxdxniEbLFNJyNCwtBTfb1+6lCEQ8Prrk9cYBdRqA80HlIVKKQGb3T+T2yPOM2SDZSoJGQo8L3ntltM6gLREcAEipYsASqJSSiA3AtPQ9rZ5tOX3KEHeD0JIMkHcRxJOpnvX6UJqUqFyQSXggu0KKkHeD0JIMrZqjgF79lhM90ERmCAYJC0NRKMxkPKQvUIl4MLKiin+HIz8azXzesj+bEJGhfDi3aS+u9EA7r47ZuAe5BcLB4Mk0WoB584N1QwgoFLrBGq1+JDcsTGa7wkZBTzPpO9PS94JpCz9GR937zT27wceeshVxL4z0usEbr89235CSHVYXAQOHnRTAKnRm0kKIBw00m6XSgF0Q6WUwBve0BkOKmL2E0Kqi60YfJSOVBDRCjLBgoGkYJBw0MgQmn+iVEoJHDnSeRGopld3KwTbxUYI6ZrgtpqfT1cAzWak746rIBPkCRuhYJDxXt4sIlcD+DiAGQCnAPy+qv4gpt0WgK/6L0+r6pt7+Vwbtmmgy/Swr0RzEwUXG1CJkQQhg8Al5VdArPknLmwoSCEfOAtWV41pqFYzH1bBYJCeHMMi8gEAz6rqfxCROwFcpar/OqbdpqpOZTl33gnkBhrMk1u5IkJIgEvxLsBEAB09CszBMx386dNmYYDtzQNI+pYnRTuG3wLgXv/5vQB+t8fz9USjYT92/HhxcnRgW44Y1CugiYiQzKSt8hUJRW4ixvSTWpB8NOhVCbxYVZ8BAP/xRZZ2V4rImog8KiJWRSEiC367tbNnz2YW5uhR+zHVAfaxtotKZPdFefCg2UeFQEgqSX11swmcPPy/sPLgjBlkxeX7V+1UBEOS9C1XVDVxA/AQgK/FbG8B8MNI2x9YzvFS//FVML6DX0j73Ouuu067QUTV/LudW6PR1Sl7p91Wrdd3C5MkKGDat9sDEpiQ8hN3W9UnLmm7cYfbPRZszaZp22xW4p4DsKYp/Wt4S50JqOr1qvramO0TAL4nIi8BAP/x+5ZzPO0/PgHgEQCv605lpXP4sP3YwBzEcUnp0nwxFy6YkAfOCgiJpeO2amxiVW/D3MZ/Ng1c/J2BX65CIZ9Z6dUc9ACAQ/7zQwA+EW0gIleJyBX+82kAbwDwjR4/10ppnffRpHRJ2UnDhM1EQR4TKgYyIqRFVs/BwynMYBtjOPXDF2Lu0j3uJx9F008cWaYN0Q1AA8DDAL7jP17t758F8DH/+a/BhId+2X98u8u5uzUHqaqOjcXP+kS6PmX+xM1ls25TU5WYvhISR6K5R8TYdycmst0ztVqlTD9xIKM5qFK5gwJKGyoaxfND1oJIhW7/i8sxcKM3lSXVxRpZjVM4hVdmP2G9PjSFXnphpHMHBSRZWo4dK5ElJTARqQInT7qbiKJsbDCyiFQOWwjoaWQI4QxGhLtyRZAwlVQCaWa+Q4eSjw+EQCG02+7lLcMEswiGmpJhJVIFfi/itYBt/2VqtZ0AjJMnzb0xok5fFyqpBObmkheObW2VuG8MhzwAGQsk+4QVwvw8FQIpP55nKr6HQviW9U7UcX5XszrOYxnvsp+nXjelHUc42icrlVQCQPLCMSBDYelBkJeZKMz6urnJpqe5QpmUj6UlU/E9xBzuxypuQxOnINhGE6ewitswh/t3Gu3ZY0Z8HelBiTNZvMhFbr1EB+14ye1bqSKFXGm3TURELxFF4W3PHnO+ikdLkAHTbqcvyHJd2DUiET69gLwXiw0z+/fbjw1lepC5OZMIRdX4DnoxGQFm5LWx0ZlGl5C8CKVr9vQmzKw/grH5A5iZ3tx9qbnekDT35E8WjVHklsdMQFV1377OgUTlMjIEI608ZgfNpmqrZUZbwair1Rr0NyTDin9dtnFA69i034fttpmZRq/H8XHOVjOCjDOBgXf2ti0vJaC6ezY6ObmzmKyS/VurlW1q7boFP5TL1J6QAP9abOJJ65jjMlFzZ6PB66sLqAQSaLWS+7fKEJ4ZBCP6blZXRu2wtkR4lfsBSQfdKn//OhRsVcc3V3KoBBII+sPoNjaW+0eVk/CNHKcU0tJY2ExOImqGj3MAAAzASURBVObcHMkNNzYzYGz+Bkebqv9ep5kAyQUqgcQfx76NZF8VN7qzacogIsP2AzYa8TbdiYkdBRGdndCcVCxRJT05uWNvn5qK/19bLbvyd+3B221tN+5I9gmQ3KASSMDWv3FEEiLJZtat87nRsM8ywj0B/Q39w/a/upgBbco/oy2Hf28xUAkkkHQf0DYZIsks0A+nc9AjRBXFxAQjQ/Kg1/8tRvm3cUCbtaf415QQKoEUJift1zlxIC76qF7vbRFb0JOktQtmDTbfQ1R57d9fnaFnL8PoXsKHYwICUsM9yUChEkihFx8X8YnrkGxx3sFoPm2k6TpStfkebEUkep1ZxPkywp8VPHc5X9x6Dpc45V4v2l5mAb5s7dbnzMgfW1rDJQ6kSgyVgAO0TfYJ2wg9qYBO0Jnltdgty5bWkWYt/JN0vrRzJSmCXh2zrr/t+PiOwggpJ9efgSbVckAl0ANcKNtH0qKD8qi01s2W1JF2o5hs50s7V61ml6NXx6ztt02YGYUHSkkBFd3oJNJfqAS6ZGQWkpWZ6DqGqNmnV99D1o60GzOK7Xwu57LR60wg+tumTH+70cc0qZYHKoEuSRrt8OIeEFl8Dy4+gWGdCRTsyHL92kzmWU6oBLqEo5whopvoINvMYhh8AsH7C3JkuUxaeE+UFyqBLkmze9LeWQG66UjLEB1UMLaZAEf+w0FWJSDmPeVjdnZW19bWCvu8xUVThN6GiElhTsiw4nmmgNfp0yZ9//JyfDr+oATAhQs7++p1Fu0aFkTkMVWddW1f6aIyWVhZAVot+/GhLEJDiE+otgtUk2sIhctcs2pj9aESCLGyYgp21eu799frZtTkeaY0L0v0krITvVaPHNk9sgfMa1ut7aDMNQt4VZ/xQQtQNoKLPTptBnZPkYORVPg9hJSBqDlnfd3e9vTpYmQi5YU+AUdmZuJvpmbTjJQIKQu2azUOXr/Vgz6BPmEbMXEkRQZBkmnS9ZoMzJxktKEScMTmGKbDmBRNmpPXdk02GnT2kk6oBBxZXrY7jAkpgmD0Pz+f7OS1XatHj9LZSzqhEnCEYXNkkIRH/zYCMxCvVZIFOoYJGQJcnL108hKAjuFSwvUFJI20ayTN2UvTJOkWKoE+k2WlJhlNXK6RpAAEmntIL1AJ9JmlpXgn3vw8ZwWjgucBU1PGPi8C1GomV1WA7RoJr+a1OXvbbTp5SW9QCfSZpGk8ZwXVx/OAQ4eA8+d39m1vm2SFgSJwWYNCZy/pF3QM9xk69EaLaKbOzU1gYyO+ba0GPPccV6OTfKFjuGTETeOj2EaCdCiXl7j/Js62b1MAALC1ZR65BoUMlCzFB4rcBlFovl/E1RJJK1iTVIiqhHVIKk20Fk2rFV/tMWv543BFyQILh5GKgyIriwH4pwC+DmAbwGxCuxsAfAvA4wDudDl3lZRAQJZSsS51XqkI+ku0imWwdVN/nv8fKYqsSqBXc9DXAPxjAJ+1NRCRGoAPA3gTgH0ADojIvh4/dyjJ4txzSQK2umoeFxeB8XFzzvHx3ZEnpJOwKWdqytjmo79dYNqJM+doRjdaowFMTu68HhszBYxWVrr+CoTkRxaNYdsAPALLTADA6wF8KvT6nQDemXbOKs4EsuAyEwhGk0lmplE3K8SVCE4bybda7r9/eGs03Gd6hPQLFDwTcOFlAJ4KvT7j7+tARBZEZE1E1s6ePVuAaOXFxaFcq+3MBuIYhRDUYFQfxN8HsfgiZpR/yy07kTeBIzZtJL+6mj4TE9n9OkjQxjBOMmykKgEReUhEvhazvcXxMyRmX+xtqKqrqjqrqrPXXHON4+mrSdh0ZGNhYadjs5FUQjDMsEQieR4wPb3T0c/P73Ty29u7254/D1y6lP0ztraSV+jW68Dhw/GdPcsykqEjy7TBtoHmoL7Tau2YM8LRQcG+pE0k+dyuDuu0CJY400tSu6yRMO226p492c00WbdazR6d1WjQvEPKDYqMDrp8kmQlMA7gCQCvBDAB4MsAfintnFQCbiT5BMK+gSRs9u/w+9IURVJIa1o7V7t5N3b6brZAwTJskwwjhSoBAL8HY+P/GYDvBSN+AC8F8GCo3Y0Avg3guwCWXM5NJeBOeJaQ1AHbsDlKwzOINEWR1kGntUtTVElydrMF55qcVB0b08szAIZtkmEnqxJg2oiKEU1bsLycbpd2SVswNma6zygixv5tO+7aLjjejZxJ1GrAC19oQj1rNWPvbzbdfhdChhGmjRhxunFMuqQtSKuxnFZrOa2dS63m5WVgzx778bHI1dxoAPfeC5w7ZxTPc8+ZRzpsCdmBSoA4LWJLUxRJIa1p7Vzz5MzNAXffbTr3gEbDpFNWNaP8sNHn3Dl29oSkksV2VORGn0D5GHR0ECEkHdAnQAghowt9AoQQQpyhEiCEkBGGSoAQQkYYKgFCCBlhqAQIIWSEKW10kIicBRBeHzoN4NyAxHGFMuYDZcwHypgPwyZjU1Wd0zCXVglEEZG1LGFPg4Ay5gNlzAfKmA9Vl5HmIEIIGWGoBAghZIQZJiWQUEixNFDGfKCM+UAZ86HSMg6NT4AQQkj+DNNMgBBCSM5QCRBCyAhTWiUgIleLyKdF5Dv+41WWdh8Qka+LyDdF5IMiIiWUca+I/JUv4zdEZKZsMvptXyAi/1dEPlSUfK4yisgvi8gX/P/6KyLyzwqS7QYR+ZaIPC4id8Ycv0JEPu4f/5si/9sMMv6hf919RUQeFpFm2WQMtXuriKiIFB6S6SKjiPy+/1t+XUT+tGwy+n3NZ0TkS/7/fWPqSbPknS5yA/ABAHf6z+8E8P6YNr8G4PMAav72BQC/WSYZ/WOPAPht//kUgHrZZPSPHwXwpwA+VML/+jUArvWfvxTAMwBe2Ge5ajB1sV8FYALAlwHsi7RZBHDcf34TgI8X/Nu5yPhbwTUHoFVGGf12zwfwWQCPApgtm4wArgXwJQBX+a9fVEIZVwG0/Of7AJxKO29pZwIA3gLgXv/5vQB+N6aNArgS5ge5AsAemIL3RZEqo4jsAzCuqp8GAFXdVNULxYno9DtCRK4D8GIAf1WQXGFSZVTVb6vqd/znTwP4PgDnVZFd8isAHlfVJ1T1IoD/4ssaJiz7nwPYX+Rs1EVGVf1M6Jp7FMDLC5TPSUaffw8zIPhpkcL5uMh4G4APq+oPAEBVv19CGRXAC/znPwfg6bSTllkJvFhVnwEA//FF0Qaq+gUAn4EZFT4D4FOq+s0yyQgzgv2hiPw3f4r2H0WkViYZRWQMwB8B+FcFyhXG5Xe8jIj8Cozi/26f5XoZgKdCr8/4+2LbqOpzAH4EoIHicJExzNsBfLKvEnWSKqOIvA7AK1T1fxYpWAiX3/E1AF4jIp8XkUdF5IbCpDO4yPhuAPMicgbAgwDuSDvpeF7SdYOIPATg52MOLTm+/9UAfhE7I5tPi8g/UNXP5iRizzLC/Ma/AeB1AE4D+DiAmwHclYd8QC4yLgJ4UFWf6tcgNgcZg/O8BMBJAIdUdTsP2ZI+LmZfNKbapU0/cf58EZkHMAvgjX2VKOajY/ZdltEfhPwJzH0xKFx+x3EYk9BvwvQ5nxOR16rqD/ssW4CLjAcA3KOqfyQirwdw0pfReq8MVAmo6vW2YyLyPRF5iao+49/4cVOv3wPwqKpu+u/5JIBfhbErlkXGMwC+pKpP+O/5C1/G3JRADjK+HsBviMgijM9iQkQ2VdXqwBuAjBCRFwD4SwD/RlUfzUu2BM4AeEXo9cvROb0O2pwRkXGYKfizBcgW/fyAOBkhItfDKNw3qurPCpItIE3G5wN4LYBH/EHIzwN4QETerKpF1Zh1/a8fVdVLAJ4UkW/BKIUvFiOik4xvB3ADYCwlInIlTHI5q+mqzOagBwAc8p8fAvCJmDanAbxRRMZFZA/MCKdIc5CLjF8EcJWIBPbrfwjgGwXIFpAqo6rOqepeVZ0B8C8B3JenAnAgVUYRmQDw333Z/qwgub4I4FoReaX/+Tf5soYJy/5WAH+tvleuLDL6ppaPAHjzAOzYqTKq6o9UdVpVZ/xr8FFf1iKLjLv8138B42SHiEzDmIeeKJmMpwHs92X8RRif6dnEsxbp3c7oCW8AeBjAd/zHq/39swA+FvKWfwSm4/8GgD8um4z+698G8BUAXwVwD4CJsskYan8zio8Ocvmv5wFcAvC3oe2XC5DtRgDfhvE/LPn73gPTScG/yf4MwOMA/jeAVxX52znK+BBMwETwuz1QNhkjbR9BwdFBjr+jAPhjv6/5KoCbSijjPpiIyS/7//U/Sjsn00YQQsgIU2ZzECGEkD5DJUAIISMMlQAhhIwwVAKEEDLCUAkQQsgIQyVACCEjDJUAIYSMMP8fNeaA1dU9cQkAAAAASUVORK5CYII=\n", "text/plain": ["<Figure size 432x288 with 1 Axes>"]}, "metadata": {"needs_background": "light"}, "output_type": "display_data"}], "source": ["import numpy, matplotlib, random, math\n", "import matplotlib.pyplot as pylab\n", "\n", "\n", "def matrix_square_root(sigma) :\n", "    eigen, vect = numpy.linalg.eig(sigma)\n", "    dim = len(sigma)\n", "    res = numpy.identity(dim)\n", "    for i in range(0,dim) :\n", "        res[i,i] = eigen[i]**0.5\n", "    return vect * res * vect.transpose()\n", "    \n", "def chi2_level (alpha = 0.95) :\n", "    N = 1000\n", "    x = [ random.gauss(0,1) for _ in range(0,N) ]\n", "    y = [ random.gauss(0,1) for _ in range(0,N) ]\n", "    r = map ( lambda c : (c[0]**2+c[1]**2)**0.5, zip(x,y))\n", "    r = list(r)\n", "    r.sort()\n", "    res = r [ int (alpha * N) ]\n", "    return res\n", "    \n", "def square_figure(mat, a) :    \n", "    x = [ ]\n", "    y = [ ]\n", "    for i in range (0,100) :\n", "        x.append( a * mat[0][0]**0.5 ) \n", "        y.append( (random.random ()-0.5) * a * mat[1][1]**0.5*2 )\n", "        x.append( -a * mat[0][0]**0.5 ) \n", "        y.append( (random.random ()-0.5) * a * mat[1][1]**0.5*2 )\n", "\n", "        y.append( a * mat[1][1]**0.5 ) \n", "        x.append( (random.random ()-0.5) * a * mat[0][0]**0.5*2 )\n", "        y.append( -a * mat[1][1]**0.5 ) \n", "        x.append( (random.random ()-0.5) * a * mat[0][0]**0.5*2 )\n", "        \n", "    pylab.plot(x,y, 'ro')\n", "    \n", "    x = [ ]\n", "    y = [ ]\n", "    for i in range (0,100) :\n", "        x.append( a ) \n", "        y.append( (random.random ()-0.5) * a*2 )\n", "        x.append( -a ) \n", "        y.append( (random.random ()-0.5) * a*2  )\n", "        \n", "        y.append( a ) \n", "        x.append( (random.random ()-0.5) * a*2  )\n", "        y.append( -a ) \n", "        x.append( (random.random ()-0.5) * a*2  )\n", "        \n", "    xs,ys = [],[]\n", "    for a,b in zip (x,y) :\n", "        ar = numpy.matrix( [ [a], [b] ] ).transpose()\n", "        we = ar * root\n", "        xs.append( we [0,0] )\n", "        ys.append( we [0,1] )\n", "        \n", "    pylab.plot(xs,ys, 'bo')\n", "    pylab.show()\n", "\n", "def circle_figure (mat, a) :\n", "    x = [ ]\n", "    y = [ ]\n", "    for i in range (0,200) :\n", "        z = random.random() * math.pi * 2\n", "        i = a * mat[0][0]**0.5 * math.cos(z)\n", "        j = a * mat[0][0]**0.5 * math.sin(z)\n", "        x.append ( i )\n", "        y.append ( j )\n", "    pylab.plot(x,y, 'ro')\n", "    \n", "    x = [ ]\n", "    y = [ ]\n", "    for i in range (0,200) :\n", "        z = random.random() * math.pi * 2\n", "        i = a * math.cos(z)\n", "        j = a * math.sin(z)\n", "        x.append ( i )\n", "        y.append ( j )\n", "        \n", "    xs,ys = [],[]\n", "    for a,b in zip (x,y) :\n", "        ar = numpy.matrix( [ [a], [b] ] ).transpose()\n", "        we = ar * root\n", "        xs.append( we [0,0] )\n", "        ys.append( we [0,1] )\n", "        \n", "    pylab.plot(xs,ys, 'bo')\n", "    pylab.show()\n", "\n", "level = chi2_level ()\n", "mat   = [ [0.1, 0.05], [0.05, 0.2] ]\n", "npmat = numpy.matrix(mat)\n", "root  = matrix_square_root (npmat)\n", "square_figure (mat, 1.96)\n", "circle_figure (mat, level)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["## p-value ratio"]}, {"cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [{"name": "stdout", "output_type": "stream", "text": ["0.9487\n"]}], "source": ["import random, math\n", "\n", "def densite_gauss (mu, sigma, x) :\n", "    e = -(x - mu)**2 / (sigma**2 * 2)\n", "    d = 1. / ((2*math.pi)**0.5 * sigma)\n", "    return d * math.exp (e)\n", "    \n", "def simulation_vector (N, mu, sigma) :\n", "    return [ random.gauss(mu,sigma) for n in range(N) ]\n", "    \n", "def ratio (vector, x, fdensite) :\n", "    under = 0\n", "    above = 0\n", "    fx    = fdensite(x)\n", "    for u in vector:\n", "        f = fdensite (u)\n", "        if f >= fx:\n", "            above += 1\n", "        else:\n", "            under += 1\n", "    return float(above) / float (above + under)\n", "        \n", "x     = 1.96\n", "N     = 10000\n", "mu    = 0\n", "sigma = 1\n", "\n", "v = simulation_vector(N, mu, sigma)\n", "g = ratio(v, x, lambda y: densite_gauss (mu, sigma, y) )\n", "print (g)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["## p-values and EM\n", "\n", "See [Applying the EM Algorithm: Binomial Mixtures](http://statisticalrecipes.blogspot.fr/2012/04/applying-em-algorithm-binomial-mixtures.html)."]}, {"cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [{"name": "stdout", "output_type": "stream", "text": ["------- sample\n", "ave 0.38\n", "mea 0.373\n", "min lk -3393.2292120130046 0.373\n"]}, {"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>average</th>\n", "      <th>pi</th>\n", "      <th>p</th>\n", "      <th>q</th>\n", "      <th>likelihood</th>\n", "    </tr>\n", "  </thead>\n", "  <tbody>\n", "    <tr>\n", "      <th>0</th>\n", "      <td>0.373</td>\n", "      <td>0.000324</td>\n", "      <td>0.341877</td>\n", "      <td>0.373010</td>\n", "      <td>-9358.705695</td>\n", "    </tr>\n", "    <tr>\n", "      <th>1</th>\n", "      <td>0.373</td>\n", "      <td>0.863747</td>\n", "      <td>0.284788</td>\n", "      <td>0.932204</td>\n", "      <td>-4531.967709</td>\n", "    </tr>\n", "    <tr>\n", "      <th>2</th>\n", "      <td>0.373</td>\n", "      <td>0.936083</td>\n", "      <td>0.346101</td>\n", "      <td>0.766941</td>\n", "      <td>-4490.512057</td>\n", "    </tr>\n", "    <tr>\n", "      <th>3</th>\n", "      <td>0.373</td>\n", "      <td>0.123023</td>\n", "      <td>0.290964</td>\n", "      <td>0.384508</td>\n", "      <td>-3563.557269</td>\n", "    </tr>\n", "    <tr>\n", "      <th>4</th>\n", "      <td>0.373</td>\n", "      <td>0.538835</td>\n", "      <td>0.053584</td>\n", "      <td>0.746213</td>\n", "      <td>-3487.438442</td>\n", "    </tr>\n", "    <tr>\n", "      <th>5</th>\n", "      <td>0.373</td>\n", "      <td>0.346351</td>\n", "      <td>0.057880</td>\n", "      <td>0.539974</td>\n", "      <td>-3302.391944</td>\n", "    </tr>\n", "    <tr>\n", "      <th>6</th>\n", "      <td>0.373</td>\n", "      <td>0.797540</td>\n", "      <td>0.376491</td>\n", "      <td>0.359248</td>\n", "      <td>-3144.938682</td>\n", "    </tr>\n", "    <tr>\n", "      <th>7</th>\n", "      <td>0.373</td>\n", "      <td>0.392520</td>\n", "      <td>0.592563</td>\n", "      <td>0.231131</td>\n", "      <td>-2902.915478</td>\n", "    </tr>\n", "    <tr>\n", "      <th>8</th>\n", "      <td>0.373</td>\n", "      <td>0.390241</td>\n", "      <td>0.459488</td>\n", "      <td>0.317648</td>\n", "      <td>-2778.903072</td>\n", "    </tr>\n", "    <tr>\n", "      <th>9</th>\n", "      <td>0.373</td>\n", "      <td>0.609127</td>\n", "      <td>0.338062</td>\n", "      <td>0.427447</td>\n", "      <td>-2764.987703</td>\n", "    </tr>\n", "  </tbody>\n", "</table>\n", "</div>"], "text/plain": ["   average        pi         p         q   likelihood\n", "0    0.373  0.000324  0.341877  0.373010 -9358.705695\n", "1    0.373  0.863747  0.284788  0.932204 -4531.967709\n", "2    0.373  0.936083  0.346101  0.766941 -4490.512057\n", "3    0.373  0.123023  0.290964  0.384508 -3563.557269\n", "4    0.373  0.538835  0.053584  0.746213 -3487.438442\n", "5    0.373  0.346351  0.057880  0.539974 -3302.391944\n", "6    0.373  0.797540  0.376491  0.359248 -3144.938682\n", "7    0.373  0.392520  0.592563  0.231131 -2902.915478\n", "8    0.373  0.390241  0.459488  0.317648 -2778.903072\n", "9    0.373  0.609127  0.338062  0.427447 -2764.987703"]}, "execution_count": 7, "metadata": {}, "output_type": "execute_result"}], "source": ["from scipy.stats import norm\n", "import random, math\n", "\n", "\n", "def average_std_deviation(sample):\n", "    mean = 0.\n", "    var  = 0.\n", "    for x in sample:\n", "        mean += x\n", "        var += x*x\n", "    mean /= len(sample)\n", "    var /= len(sample)\n", "    var -= mean*mean\n", "    return mean,var ** 0.5\n", "\n", "def bootsample(sample):\n", "    n = len(sample)-1\n", "    return [ sample[ random.randint(0,n) ] for _ in sample ]\n", "    \n", "def bootstrap_difference(sampleX, sampleY, draws=2000, confidence=0.05):\n", "    diff = [ ]\n", "    for n in range (0,draws) :\n", "        if n % 1000 == 0: \n", "            print(n)\n", "        sx = bootsample(sampleX)\n", "        sy = bootsample(sampleY)\n", "        px = sum(sx) * 1.0/ len(sx)\n", "        py = sum(sy) * 1.0/ len(sy)\n", "        diff.append (px-py) \n", "    diff.sort()\n", "    n = int(len(diff) * confidence / 2)\n", "    av = sum(diff) / len(diff)\n", "    return av, diff [n], diff [len(diff)-n]\n", "\n", "# generation of a sample\n", "\n", "def generate_obs(p):\n", "    x = random.random()\n", "    if x <= p : return 1\n", "    else : return 0\n", "\n", "def generate_n_obs(p, n):\n", "    return [ generate_obs(p) for i in range (0,n) ]\n", "        \n", "# std deviation\n", "\n", "def diff_std_deviation(px, py):\n", "    s = px*(1-px) + py*(1-py)\n", "    return px, py, s**0.5\n", "\n", "def pvalue(diff, std, N):\n", "    theta = abs(diff)\n", "    bn = (2*N)**0.5 * theta / std\n", "    pv = (1 - norm.cdf(bn))*2\n", "    return pv\n", "    \n", "def omega_i (X, pi, p, q) :\n", "    np = p * pi     if X == 1 else (1-p)*pi\n", "    nq = q * (1-pi) if X == 1 else (1-q)*(1-pi)\n", "    return np / (np + nq)\n", "    \n", "def likelihood (X, pi, p, q) :\n", "    np = p * pi     if X == 1 else (1-p)*pi\n", "    nq = q * (1-pi) if X == 1 else (1-q)*(1-pi)\n", "    return math.log(np) + math.log(nq)\n", "    \n", "def algoEM (sample):\n", "    p   = random.random()\n", "    q   = random.random()\n", "    pi  = random.random()\n", "    iter   = 0\n", "    while iter < 10 :\n", "        lk = sum ( [ likelihood (x, pi, p, q) for x in sample ] )\n", "        wi  = [ omega_i (x, pi, p, q) for x in sample ]\n", "        sw  = sum(wi)\n", "        pin = sum(wi) / len(wi)\n", "        pn  = sum([ x * w     for x,w in zip (sample,wi) ]) / sw\n", "        qn  = sum([ x * (1-w) for x,w in zip (sample,wi) ]) / (len(wi) - sw)\n", "        \n", "        pi,p,q = pin,pn,qn\n", "        iter += 1\n", "        \n", "    lk = sum ( [ likelihood (x, pi, p, q) for x in sample ] )\n", "    return pi,p,q, lk\n", "\n", "\n", "# mix\n", "p,q   = 0.20, 0.80\n", "pi    = 0.7\n", "N     = 1000\n", "na    = int(N * pi)\n", "nb    = N - na\n", "\n", "print(\"------- sample\")\n", "sampleX = generate_n_obs(p, na) + generate_n_obs (q, nb)\n", "random.shuffle(sampleX)\n", "print(\"ave\", p * pi + q*(1-pi))\n", "print(\"mea\", sum(sampleX)*1./len(sampleX))\n", "\n", "lk = sum ( [ likelihood (x, pi, p, q) for x in sampleX ] )\n", "print (\"min lk\", lk, sum (sampleX)*1. / len(sampleX))\n", "res = []\n", "for k in range (0, 10) :\n", "    r = algoEM (sampleX)\n", "    res.append ( (r[-1], r) )\n", "res.sort ()\n", "\n", "rows = []\n", "for r in res:\n", "    pi,p,q,lk = r[1]\n", "    rows.append( [p * pi + q*(1-pi)] + list(r[1]))\n", "\n", "df = pandas.DataFrame(data=rows)\n", "df.columns = [\"average\", \"pi\", \"p\", \"q\", \"likelihood\"]\n", "df"]}, {"cell_type": "markdown", "metadata": {}, "source": ["## p-value and heavy tail"]}, {"cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [{"data": {"text/plain": ["[357621, 148, 18, 1812876449, 36150]"]}, "execution_count": 8, "metadata": {}, "output_type": "execute_result"}], "source": ["from scipy.stats import norm, zipf\n", "import sys\n", "\n", "\n", "def generate_n_obs_zipf (tail_index, n) :\n", "    return list(zipf.rvs(tail_index, size=n))    \n", "                     \n", "def hill_estimator (sample) :\n", "    sample = list(sample)\n", "    sample.sort(reverse=True)\n", "    end = len(sample)/10\n", "    end = min(end,100)\n", "    s = 0.\n", "    res = []\n", "    for k in range (0,end) :\n", "        s += math.log(sample[k])\n", "        h = (s - (k+1)*math.log(sample[k+1]))/(k+1)\n", "        h = 1./h\n", "        res.append( [k, h] )\n", "    return res\n", "    \n", "# mix\n", "tail_index  = 1.05\n", "N           = 10000\n", "\n", "sample = generate_n_obs_zipf(tail_index, N)\n", "sample[:5]"]}, {"cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [{"name": "stdout", "output_type": "stream", "text": ["[9999, 55186871.0339, 233342554.46156308]\n"]}, {"data": {"image/png": "iVBORw0KGgoAAAANSUhEUgAAAWsAAAEDCAYAAADz4SVPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZwcdZ3/8denjzkzmVyTg4SQcIUEApoMEAzGCIKEU1jkUFlx3c26HivuevFz/QVW96HygH2gP1fcCIgHAnJ4RUU5gkLkcBKEhFwEQsLknEySyUzm6uP7++NbM9MzmaMnmZ7pmryfj0c/urvq29Xf6up+17e+VV1lzjlERCS/RYa6AiIi0jeFtYhICCisRURCQGEtIhICCmsRkRBQWIuIhEDOwtrM7jWz3Wa2JouyU81suZm9bGavmtnFuaqXiEgY5bJlfR9wUZZl/wP4uXPuncB1wPdyVSkRkTDKWVg75/4M7M0cZmYnmNnjZrbSzJ41s1PaigMjg8flwPZc1UtEJIxig/x+S4FPOOdeN7Oz8S3o84BbgD+a2WeAUuB9g1wvEZG8NmhhbWYjgHcBD5tZ2+DC4P564D7n3B1mdg7wEzM7zTmXHqz6iYjks8FsWUeA/c65d3Qz7uME/dvOuefNrAgYB+wexPqJiOStQTt0zzl3ANhsZh8EMO+MYPRW4Pxg+EygCKgZrLqJiOQ7y9VZ98zsAWAhvoW8C1gCPA3cBUwC4sCDzrn/NLNZwA+AEfidjV90zv0xJxUTEQmhnIW1iIgMHP2DUUQkBHKyg3HcuHFu2rRpuZi0iMiwtHLlyj3OuYqexuckrKdNm0ZVVVUuJi0iMiyZ2ZbexqsbREQkBLIKazP7nJm9ZmZrzOyB4DhoEREZJH2GtZlNBv4VqHTOnQZE8SdbEhGRQZJtn3UMKDazBFCCTrQkIv2QSCSorq6mubl5qKsy5IqKipgyZQrxeLxfr+szrJ1z28zsdvy/DJuAP3b3hxUzWwwsBpg6dWq/KiEiw1t1dTVlZWVMmzaNjHMDHXWcc9TW1lJdXc306dP79dpsukFGA1cA04FjgFIz+0g3lVjqnKt0zlVWVPR49ImIHIWam5sZO3bsUR3UAGbG2LFjD2sLI5sdjO8DNjvnapxzCeAx/NnzRESydrQHdZvD/RyyCeutwDwzKzH/LucD6w7r3fqy+VnY83pOJi0iEmZ9hrVz7kXgEWAVsDp4zdKc1OZHl8J3K3MyaRGRgXDTTTfx5z//udcyy5YtY8mSJQP6vlkdZ+2cW+KcO8U5d5pz7gbnXMuA1kJEJAT27t3LCy+8wIIFC3otd8kll/DrX/+axsbGAXtv/YNRRI4KH/jAB5g7dy6nnnoqS5cu5a677uKLX/xi+/j77ruPz3zmMwB87Wtf45RTTuGCCy7g+uuv5/bbbwfgkUce4aKLOq4DPm3aNJYsWcKcOXOYPXs269evB3y/9MKFC1m2bNmA1X+wr8EoIke5W3/zGmu3HxjQac46ZiRLLju11zL33nsvY8aMoampiTPPPJOnnnqK+fPnc9tttwHw0EMP8ZWvfIWqqioeffRRXn75ZZLJJHPmzGHu3LkArFixgquvvrrTdMeNG8eqVav43ve+x+23387dd98NQGVlJc8++yzXXHPNgMyjWtYiclT4zne+wxlnnMG8efN4++232bx5M8cffzwvvPACtbW1bNiwgfnz5/Pcc89xxRVXUFxcTFlZGZdddln7NHbs2EHXQ5OvuuoqAObOnctbb73VPnz8+PFs3z5w/x9Uy1pEBlVfLeBceOaZZ3jyySd5/vnnKSkpYeHChTQ3N3Pttdfy85//nFNOOYUrr7wSM6O3C7IUFxcfcox0YaG/7nc0GiWZTLYPb25upri4eMDmQS1rERn26urqGD16NCUlJaxfv54XXngB8K3iX/7ylzzwwANce+21AJx77rn85je/obm5mYaGBn7729+2T2fmzJls2rQpq/fcuHEjp5122oDNg8JaRIa9iy66iGQyyemnn85Xv/pV5s2bB8Do0aOZNWsWW7Zs4ayzzgLgzDPP5PLLL+eMM87gqquuorKykvLycsAf5fHMM89k9Z7Lly/nkksuGbB5yMk1GCsrK91hXXzglvLgvm5gKyQiQ2rdunXMnDlzqKuRtYaGBkaMGEFjYyMLFixg6dKlzJkzB/At72XLljFq1KgeX79r1y4+9KEP8dRTT3U7vrvPw8xWOud6/KOJ+qxFRLpYvHgxa9eupbm5mY9+9KPtQQ1wxx13sHXr1l7DeuvWrdxxxx0DWieFtYhIFz/72c96HHf22Wf3+fozzzxzIKsDqM9aRCQUFNYiIiGgsBYRCQGFtYhICCisRURCQGEtIkelO++8s8dTmN533318+tOf7lT2xz/+ca/TW716NTfeeONAVrEThbWIHJV6C+tMyWSSe++9lw996EO9lps9ezbV1dVs3bp1oKrYSZ/HWZvZDOChjEHHA//XOXdnTmokIsPb778MO1cP7DQnzoZF3+xx9MGDB7nmmmuorq4mlUrxwQ9+kO3bt/Pe976XcePGsXz5cn74wx/yjW98g0mTJnHyySe3n6Dp6aefZs6cOcRiPi4XLlzI2WefzfLly9m/fz/33HMP7373uwG47LLLePDBBzudJ3ugZHNZrw3OuXc4594BzAUagV8MeE1ERHLk8ccf55hjjuGVV15hzZo13HTTTRxzzDEsX76c5cuXs2PHDpYsWcKKFSt44oknWLt2bftrV6xY0X4+6zbJZJKXXnqJO++8k1tvvbV9eNs5rHOhv/9gPB94wzm3JReVEZGjQC8t4FyZPXs2n//85/nSl77EpZde2t4SbvPiiy+ycOHC9nNVX3vttWzcuBHw57Dueh6PwTqHdab+hvV1wAO5qIiISK6cfPLJrFy5kt/97nfcfPPNXHjhhYeUMbNuXzuU57DOlPUORjMrAC4HHu5h/GIzqzKzqpqamoGqn4jIEdu+fTslJSV85CMf4fOf/zyrVq2irKyM+vp6wJ/v45lnnqG2tpZEIsHDD3fE3FCewzpTf1rWi4BVzrld3Y10zi0FloI/ReoA1E1EZECsXr2aL3zhC0QiEeLxOHfddRfPP/88ixYtYtKkSSxfvpxbbrmFc845h0mTJjFnzhxSqRQAixYt4oYbbsjqfQb6HNaZsj6ftZk9CPzBOffDvsrqfNYikils57Pu6sorr+S2227jpJNO6rFMS0sL73nPe3juuefajxzpyeGczzqrbhAzKwEuAB7LpryIyHDyzW9+kx07dvRaZuvWrXzzm9/sM6gPV1ZTdc41AmNzUgMROSo453rciZfvZsyYwYwZM3otc9JJJ/Xa8m5zuFfn0j8YRSTnioqKqK2tPeygGi6cc9TW1lJUVNTv1+pKMSKSc1OmTKG6uhodKeZXXFOmTOn36xTWIpJz8Xic6dOnD3U1Qk3dICIiIaCwFhEJAYW1iEgIKKxFREJAYS0iEgIKaxGREFBYi4iEgMJaRCQEFNYiIiGgsBYRCQGFtYhICCisRURCQGEtIhICCmsRkRDI9rJeo8zsETNbb2brzOycXFdMREQ6ZHs+628DjzvnrjazAqAkh3USEZEu+gxrMxsJLABuBHDOtQKtua2WiIhkyqYb5HigBvihmb1sZnebWWnXQma22MyqzKxKl+4RERlY2YR1DJgD3OWceydwEPhy10LOuaXOuUrnXGVFRcUAV1NE5OiWTVhXA9XOuReD54/gw1tERAZJn2HtnNsJvG1mM4JB5wNrc1orERHpJNujQT4D3B8cCfIm8LHcVUlERLrKKqydc38DKnNcFxER6YH+wSgiEgIKaxGREFBYi4iEgMJaRCQEFNYiIiGgsBYRCQGFtYhICCisRURCQGEtIhICCmsRkRBQWIuIhIDCWkQkBBTWIiIhoLAWEQkBhbWISAgorEVEQiCriw+Y2VtAPZACks45XYhARGQQZXtZL4D3Ouf25KwmIiLSI3WDiIiEQLZh7YA/mtlKM1vcXQEzW2xmVWZWVVNTM3A1FBGRrMN6vnNuDrAI+JSZLehawDm31DlX6ZyrrKioGNBKiogc7bIKa+fc9uB+N/AL4KxcVkpERDrrM6zNrNTMytoeAxcCa3JdMRER6ZDN0SATgF+YWVv5nznnHs9prUREpJM+w9o59yZwxiDURUREeqBD90REQkBhLSISAgprEZEQUFiLiISAwlpEJAQU1iIiIaCwFhEJAYW1iEgIKKxFREIgP8PauaGugYhIXsnPsE6nhroGIiJ5JT/D2imsRUQy5WdYp5NDXQMRkbyisBYRCYE8DWt1g4iIZFJYi4iEQNZhbWZRM3vZzJblskKAukFERLroT8v6s8C6XFWkE4W1iEgnWYW1mU0BLgHuzm11AgprEZFOsm1Z3wl8EUj3VMDMFptZlZlV1dTUHFmtXI9vIyJyVOozrM3sUmC3c25lb+Wcc0udc5XOucqKioojq5Va1iIinWTTsp4PXG5mbwEPAueZ2U9zWiuFtYhIJ32GtXPuZufcFOfcNOA64Gnn3EdyWiuFtYhIJ3l6nLXCWkQkU6w/hZ1zzwDP5KQmmdLawSgikkktaxGREFBYi4iEgMJaRCQE8jOsdfEBEZFO8jOsddY9EZFO8jSs1Q0iIpJJYS0iEgL5GdapxFDXQEQkr+RpWLcOdQ1ERPJKXoV1ypl/kGwZ2oqIiOSZvArrVuL+wRNLhrYiIiJ5Js/COjhVSUvd0FZERCTP5FlYx4e6CiIieSmvwjpJdKirICKSl/IqrN1QV0BEJE/lVViLiEj3srlgbpGZvWRmr5jZa2Z262BUDKd2tohIm2xa1i3Aec65M4B3ABeZ2bzcVgv9MUZEJEM2F8x1zrmG4Gk8uOW+2dt6MOdvISISFln1WZtZ1Mz+BuwGnnDOvdhNmcVmVmVmVTU1NUdeM4W1iEi7rMLaOZdyzr0DmAKcZWandVNmqXOu0jlXWVFRceQ1SzQe+TRERIaJfh0N4pzbj7+6+UU5qU2m1oa+y4iIHCWyORqkwsxGBY+LgfcB63NVoa3poFXeqpa1iEibbFrWk4DlZvYq8Fd8n/WyXFWo/S/nNTlbH4iIhE6srwLOuVeBdw5CXQDYQzknsh0ifVZNROSokXf/YNzjyv2DlgNDWxERkTySd2Hd6AohEoemfUNdFRGRvJF3Ye0wKB6lsBYRyZB3YQ1A8Who2j/UtRARyRt5HNZqWYuItMnPsC5SN4iISKb8DGt1g4iIdJLHYb13qGshIpI38jOsyyb4c4O06PwgIiKQxT8Yh0TZJH9fvxMKTzzy6dW+Af9vjn98zqdh7AkwYgKccsmRT1tEZBDkeVjvgHEZYZ1KwOY/wQnng5kf9vDHYOyJcN5Xep7eyvs6Hj//3Y7HF30L5n1iwKotAyDZAjUbINEE42dC0cihrpHks2QLWBSiQZSlUxCJBo/TvjvVIlAwAvZshAPb/C1WBPFi2PcWNO6FZDOkk3540Sj/X4+icogWgEtDXbW/j8ahsMyfDqP2DWjYBQ27oXCEf23pOLjw6zmZ1fwPa4A3lvsfcOEI+NWn4IzrYfYH4cTz4bXHfJmz/9l/UN0pCv7CfskdsPpR2PoX//zxL/nzZr/73+BPt8GMRTBxdu7mS7qXaPIXm9j8J3jkHzqPG3M8TDwdJpwGo4+DcSf74a0NfljxKP9jqVkPf73b/3DHHO/LjTq247u0/rd+eNNeaKn3P766tyFeCqVjYeRk/8Nt2O2/dxb1P+YJp/ppjTvZNxyKR/u6xks6GgzZqH0DGmth+8swepqvy6ipECvsCJgdr/jpjpoKO9fArtU+LMomwoiJvmzNBv++sSJfrvxYiBX0/zPf9xY0H/CfX902KCiF0gp/X1jWv3kbTAd2wK7X4LVf+O9L3dt+eEEZuJT/PceKAeeDPNuLWkVi/maR/p1LP1YM5ZOD1zVB4Ui4sL8zleVb5WayR2hkl7D+yQf8/WXf8fevPOBv8z/b8ZoVd/a8Rlv9sL+fc6O/PXUrTHs3rP65f7xztQ/95f8F194PMy8d6DkKlzWP+uCYPBfGz4LXfgl/+Y5vWcSK4KQLYcIsH4Qv/xQqToEplTB5DoyeDnvf9NMZMd6X373OB05BqQ+fh2+ETU9C+RSo3dR9Hc79Nygo8fXY/jKs/WX35UZP88HTJlroW0gulf38RgsOveZnJOan8+Zyf9+VRWHkMT7YRkzwK5IRE2HTEz7wS8bAqON8K27Xa77lduhE6D5MehreHfP1iMaheEywgjrGr8xc2tejab9fBttX+VBrOdD79KMFfp5GjIfS8T7QLQqJg74V2trgQ6l8sn+vdNI3pErHB6+p8Cu1RBPsXuvrkGj2y7+u2r9HyRi/gt293u+jKhnb0aqNFXacG6hum18W21f5Fd6eDZ3reu7n/DJvrvO3gtKO0G2b/8KRfgWdbPIr/omz/fjWBj+vFTM6TzOVDKa33y835/z3adzJ/nvS0uCnN2K8n/4gyauwbl+XF5b5zZYDQViXjA1aJas6v2DFtzsev/QDOP06/yW5/2q4+HaYNt+Pq9vm7yNR32K48Gv++Ynn+y/Iyz/pmM5DH/b35/2HD4y2Tarhzjm/j6D6r4e2bjOVjoNXftZ52JvLof1Cb/0ImsKR/sdevx3eeYPfAjpuPpxy8aFlmw/4H/qejT7gIzH/A9rxiv9Rlk+Fs/4JTv0AJFth32Zf/sA23xqNFcKMS/yKpGxix4qkdJxvzTbuhcY9UDIO4kV+GPhW7/4t/n33bIQ9r8PW5+GE83wIth70jYp1y/zr24yc7AOydlPH57HoNh8mo6f7FuHezX7aNeuhYqYPttOv8UGwdzNMXwBTzvTTrd/puwHTCR9o0QL/m9j3lp9G/Q5f112vwYbHIdXigznZ5MO07Bj/OU2e40PWOT8PLfX+sykYAQdrOkKqYbffxK97G3at8a9NJXxDKtnqQ2v3Ol/mSC/JGon7+epL6Xhf55Mv8oE78XS/khho0Zjf2iode+i4eHHHlvogy6uwhozFPmqqD47WRv/Fb6zt3Pec6aq74bF/hO/Ph3mf8l/6+y6GL272f65prfctha6bdpFo0Fp3voW4+E++m2XXGnj66/529ifggv/0P47eNg1TSf/D3LYS3noOjnsX/O7z/os487KOgDv1Sjjrn2HqvPza1PzhxR3dQwCnXwtz/t7/IOt3+hAcNdV/UZ2DA9t9GB7cDSe93wfOtpU+GBt2+tZ2OuUf793sf1g163041myAG37Rv1ZJ0UgomuVb9H2JBa2lri0mgClzu39NTz/OaMzvkB57gu8m601L0JrN7Gd3bgCW88n9K+6cr0ck6vtt08nD6yrJRrLVrxAicf9daKjx9037fQDHin3LOZ30LWCL+JVFY61fQU441W+hNdf5FUDzAf+bLSj1W1bFo/00kk1DFpL5wpwb+AuVV1ZWuqqqqn6/bseS6fwpdTrXff1X8MD1sOF3MG7GoZs+7/0PwPluC4Bb6uCPX/Wb6plOvABOfr8PzbZyPWk+0PEjS6d9t8ijH+9cpnAkfPhhH7Rtki3w4Id9a6s/lyKbcBpUfgxmXzM0O9GSLfDmn3xL6dF/gpp1fvisK3xIT1/YsdNGRHLOzFY65yp7Gt/nr9HMjgV+DEwE0sBS59y3e3/VABg93d93DWrwLdXxp8CWFX7zDXzXxtRz4MHr/fNZV8DaX/k+RIAP3tf7+2UGZiQCs6+Gky7wm7e/+qQf3nIA7n2/D9qZl/tyL/+0c5/p6Om+NR6N+9bDzMv94y0r4NSrfCtr9cPw13vgt/8OTyzx/efbV/nNu7q3/Q7UE87reYfpQPj6+M7PLQqfehHGnZS79xSRw5ZN0ykJ/LtzbpWZlQErzewJ59zanNZszPTOz2df43cIQscOnxt+2XkT85SL4aofwJO3wt/dAwdrYctzftzUc/pfh6JyeOeH/W3Hq35Hy0s/8C34XWs6l73wv3zQnfz+7qc19oSOx3NvhDkfhW2roOoe+Nv9fviqH/n7N54GzB+SCL7P/rS/gxe+57uEdr7asQMmEvMriJGToPZNOFDtV2DHzfddMccv9F0QkYz/Pz15a8fji77lNzHP/Vz/Px8RGTTZXNZrB7AjeFxvZuuAycDghnU0Du/5EvzpW37PO3TfF3j6Nf4GcOMy351Rs9H3WR+JSaf7+wu/5ndkvfqQ33E59kS/khg9vX99k2a+/3TKXH8Uy+51MPE0v3fcIv5oiTWPQu3r/vb2C/51+zZ3nk7FTH+IV+3rHcNaG3y51//QMWzEBN99M/5UeO6//bAPP+K3HkQk7/WrU9LMpuGvx/hi7yUHwLiuO4cM3vt//C1bZr5FOtBGTYUFXxi46ZWM6ThyZerZ/v7YM2Hhl/1OOvBhvOs1362SaPQt6sw94a0H/eFRmTvJ6nf6fum3nvWt+G0v+64hgCuXKqhFQiTrsDazEcCjwE3OuUMukGhmi4HFAFOnTj3ympVP6fIGRz7J0DHr2Mk3fqa/Qfd79gtK/S1T2UQ441p/a3NgB2yr8oexiUhoZHUiJzOL44P6fufcY92Vcc4tdc5VOucqKyoqjrxmZvC51+Bd/+qfN+sCugNi5CS/gzaSn+fwEpHu9fmLNTMD7gHWOef+O/dVylA+xf9bDuDtlwb1rUVE8kk2zav5wA3AeWb2t+DWzV/McmRKcNhhf/4+LCIyzGRzNMhzDGWPcbzYn68j89A3EZGjTDj+ona0n1hJRI562sskIhICCmsRkRBQWIuIhIDCWkQkBBTWIiIhoLAWEQkBhbWISAgorEVEQkBhLSISAgprEZEQUFiLiISAwlpEJAQU1iIiIZCXYe2cG+oqiIjklbwM62RaYS0ikimvwtrwIZ1IpYe4JiIi+SWbazDea2a7zWzNYFTIYSy4bflgvJWISGhk07K+D7gox/XoZE9D62C+nYhI3uszrJ1zfwb2DkJdRESkBwPWZ21mi82sysyqampqBmqyIiLCAIa1c26pc67SOVdZUVFxxNNrTqQGoFYiIsNDXh0NEjFrf/xmzcEhrImISH7Jq7DO9Mn7Vw51FURE8kY2h+49ADwPzDCzajP7eO6rBSOKYoPxNiIiodBnIjrnrh+MirSZPq4UdsLOupbBfFsRkbyWd90gkaDbek+DwlpEpE3ehTXAx8+dTjRiJPW3cxERIE/DesbEMlJpx4Zd9UNdFRGRvJCXYV1SEAXgiu+u6LPs9v1NHGxJ9lkunXZsqT3Y3r2ik0WJSJjk5SEXl8yexKd5mWTasbq6jt31zazeVsen3nsidz+7mW89vh6ARadN5PdrdgLwm0+fyymTyohHI6TTjkik45jt6n2NnPutQ08OddWcyXziPSdw8oQyaupbMINxIwoHZybzWCrt2F3fzJjSAuKRCMm0w+FoTqQZURgjYmDBMfGNrX5FWRyPtg9rOx95Ku2IRSM0taYojEWIRAznHHsaWnmjpoGG5iRF8Sh1TQkONCc42JJk695Gtu1rIpl2lBXFKC+OU1YUp7w4zsjiGMXxKPFohHg0QkHMKIpFiUaM5mSarXsbaUmkKIhFKC2IUVoYo6woRsSM1lSa1qS/NSdS1DUlKIpHcTiKYlFiUcM5qG9O0Nia8uOco7QwxsjiOCOL/PuPLIpTFI9SGI/Qdtr1eNTXIxIx0mlHfUuSptYURfEIhbEokQhsrW2kviVJIpkmHotQEI1QFI9QEI2STKcxM9LOEYsYpYUxDKhrShCNGKm0Ix6NEI0YiVSaVNpRXBClJB6jpNB/Hv1V15Rg38FW0s6RSjuiEaMg5utbGPf1K4h2LLNU2tGYSNHYkiLlHBGDoliU4oIohbFI+7LvL+ccaQdp53DBPXR+7gCXBofjYGuK5oS/1dS30JxIcbAlhQuWQyziP6eWpC/T2JoilXZEzGhKpGhsTdLYmmqfdlNriuZkmmhQfTMjHrX271jaOdLOkUj6+1jwHolUmvrmJI2tSZoTaSIR/z+RkUVx/ufDcw7rs+hLXoa1mfG+meN5ct1uLvvuc+3D73zy9U7l2oIa6FQOoCge4acfP5t3HDuKi7/9bPvwWZNGsnbHAf/61Tt5bNU2Jo8qZtv+JsAv8ETKsXjB8RTFIkwsL2ZieSHjy4oYU1rAMaOKe6x3SzLFnzfuYeeBZqaNLWHlln2MKo4zeXQJG3fVUxyPMrI4zvEVpZxQMYLy4vjhf0i9cM71+ePZe7CVqrf2UhiPsn1/Eys27WFHXTONrSk27qonFZxTPGLQ9fTisYhRHGz91Df7sC6MRSguiFIQjVDT0NIeZGbgnL+PRyO0JnvfoolFjCmjiymKR9m6N01dU4L65gSJVP6f47wg5n/Eg33tjHjUKApWYtGIETUjYpBIOwqiEZxztKZcEDCJQ5Znb6IRaw/O3hTGIhTG/ErDOXB0E7iHPD/cOT58ZlBaEMPMz1vbCieV9mEcMb9CTKTStCTTRCNGxPwtHjWSaUcylSYWjVBWFKMkWGkmUo6U8593ruRlWAMsXnACT67b3eP4L7x/Bs+9vofn36zlsU++iwdf2srPq6rbxzcn0lz9/ec7vWbzNy7uFGJ7D7byyMq3+dFftgAw97jRrNyyD4Clf36z2/c9cfwIxpQU+Nc3thIx2Lir4bDmsaKskBMqSqkoK2J/YyvHjinxrbkC/yUYVVJAY2sy+ILDmu11NCdSjC4pYNv+JvY0tFAY8y280sIY+xpbaWxNse9gK0XxKGVFseALFSMWMUYWxymOR3nhzVpqDx56ZsNjxxQzpqSA980cz4wJfiulOZmiNem/nPHgh59MOxqakz4YIsbokjh1TQlakmkaW1OMKIxRWhglGomQSqcpiEZJpf2XPx71oT59XCmTyotoTaYpiEUYN6KQ0sIYhbEIpYWdv5bOOVqSPribEykSqTStSUdrKk1LIkVjIkVBNEJ5cZypY0toTaY52JKkoSVJQ3OStPNBWhCN+PtYhNKCKC3JNJFgR3bbymBkUay9RZVMO1qD9z3QnOBAU5IDzQlaEilaMlY6yXTQQkukKIxHKSv0Ld6WRJrmZIp02lFR5lf4saiRTPn5aWvtu6DFlg4m2dZaKy+Og0FBNNLeopuJk98AAAflSURBVC6IRTJaiSkaW5I0JlI0BS3IZNqRSqdJpX2It6bS4KAwHqUgapQU+u9CSUGMirJC4lHzrfpgXluSft5akun20IqaEYtGKI5H25eNmT8lRFMiRXPCb620rYjNfCvTgEjE35tZMBwMa98661oWgucZZS0oGzG/HEcU+mVUXhxnVEmcEYUxDCORTpNMOZLpNEXxKEXxKMXxaPDZ+s+upCBvI69PlotLaFVWVrqqqqp+v273LdN5a/S7OOuz95NIpTnpK7/vNH7K6GKq9zWx/msXURSPdjuNv7yxh2TKccaUUXzp0Vd5/DXf+l72mXM5bXJ5j+99sCXZ/kVsax1t3XuQbfubeXtvI6/vqqclmWbb/iZera6jrinR6fVtrfNPLjyBGRPLGFkUZ+2OA8yaNJJk0F9+9vSxFMUjbKltZFNNA2/sbmBTTQPrd9TTFGy+J1M+0Fq7aaEVRCNMHVvCgaYE5cVxjh1T0r6J39yaYkSwph83ohAzaGhOUt+cpL4lQTLlaGhJ0pRIcbAlydnTx3Lp6ZMYU1pAUTzKrEkjO3UdicjgMrOVzrnKHsfna1iDX3NHI8aV31vBdWdO5SPzjhvoqh6RbLobsp1O2vnNsrZpptKO2oYWzIzy4nh7d0RBLC/3CYvIEeorrPN6m6Ct9bzsM+8e4pp0byCCum06mTs4wAf3+JFFAzJ9EQk/NdNEREJAYS0iEgIKaxGREFBYi4iEgMJaRCQEFNYiIiGQVVib2UVmtsHMNpnZl3NdKRER6Syby3pFgf8BFgGzgOvNbFYuKpOLP+iIiAwH2bSszwI2OefedM61Ag8CVwx0Rer27maC7aMpoVOXioh0lU1YTwbeznheHQzrxMwWm1mVmVXV1NT0uyJl5eN4fcIiJiz4WL9fKyIy3GXzd/Pu/lN9SH+Fc24psBT8uUH6W5FINMJJ//Jgf18mInJUyKZlXQ0cm/F8CrA9N9UREZHuZBPWfwVOMrPpZlYAXAf8OrfVEhGRTH12gzjnkmb2aeAPQBS41zn3Ws5rJiIi7bI6Rapz7nfA73JcFxER6YH+wSgiEgIKaxGREFBYi4iEgMJaRCQEcnLBXDOrAbYc5svHAXsGsDphoHke/o62+QXNc38d55yr6GlkTsL6SJhZVW9X+B2ONM/D39E2v6B5HmjqBhERCQGFtYhICORjWC8d6goMAc3z8He0zS9ongdU3vVZi4jIofKxZS0iIl0orEVEQiBvwno4XZTXzI41s+Vmts7MXjOzzwbDx5jZE2b2enA/OhhuZvadYN5fNbM5GdP6aFD+dTP76FDNUzbMLGpmL5vZsuD5dDN7Maj7Q8EpdjGzwuD5pmD8tIxp3BwM32Bm7x+aOcmemY0ys0fMbH2wvM8ZzsvZzD4XfKfXmNkDZlY0HJezmd1rZrvNbE3GsAFbrmY218xWB6/5jpl1d5GXzpxzQ37Dn3r1DeB4oAB4BZg11PU6gvmZBMwJHpcBG/EXG74N+HIw/MvAt4LHFwO/x1+VZx7wYjB8DPBmcD86eDx6qOevl/n+N+BnwLLg+c+B64LH3wf+JXj8SeD7wePrgIeCx7OCZV8ITA++E9Ghnq8+5vlHwD8GjwuAUcN1OeMv57cZKM5YvjcOx+UMLADmAGsyhg3YcgVeAs4JXvN7YFGfdRrqDyWo+DnAHzKe3wzcPNT1GsD5+xVwAbABmBQMmwRsCB7/L3B9RvkNwfjrgf/NGN6pXD7d8FcQego4D1gWfAn3ALGuyxh/bvRzgsexoJx1Xe6Z5fLxBowMwsu6DB+Wy5mO67GOCZbbMuD9w3U5A9O6hPWALNdg3PqM4Z3K9XTLl26QrC7KG0bBpt87gReBCc65HQDB/figWE/zH6bP5U7gi0Db5enHAvudc8ngeWbd2+crGF8XlA/T/ILfEqwBfhh0/9xtZqUM0+XsnNsG3A5sBXbgl9tKhv9ybjNQy3Vy8Ljr8F7lS1hndVHesDGzEcCjwE3OuQO9Fe1mmOtleF4xs0uB3c65lZmDuynq+hgXivnNEMNvKt/lnHsncBC/edyTUM930Ed7Bb7r4higFFjUTdHhtpz70t/5PKz5z5ewHnYX5TWzOD6o73fOPRYM3mVmk4Lxk4DdwfCe5j8sn8t84HIzewt4EN8VcicwyszarkaUWff2+QrGlwN7Cc/8tqkGqp1zLwbPH8GH93Bdzu8DNjvnapxzCeAx4F0M/+XcZqCWa3XwuOvwXuVLWA+ri/IGe3bvAdY55/47Y9SvgbY9wh/F92W3Df/7YK/yPKAu2Mz6A3ChmY0OWjUXBsPyinPuZufcFOfcNPyye9o592FgOXB1UKzr/LZ9DlcH5V0w/LrgKILpwEn4HTF5yTm3E3jbzGYEg84H1jJMlzO++2OemZUE3/G2+R3WyznDgCzXYFy9mc0LPse/z5hWz4a6Ez+jk/1i/FETbwBfGer6HOG8nIvfrHkV+FtwuxjfX/cU8HpwPyYob8D/BPO+GqjMmNY/AJuC28eGet6ymPeFdBwNcjz+R7gJeBgoDIYXBc83BeOPz3j9V4LPYQNZ7CEf6hvwDqAqWNa/xO/1H7bLGbgVWA+sAX6CP6Jj2C1n4AF8v3wC3xL++EAuV6Ay+AzfAL5Ll53U3d30d3MRkRDIl24QERHphcJaRCQEFNYiIiGgsBYRCQGFtYhICCisRURCQGEtIhIC/x/+K066/q1KSgAAAABJRU5ErkJggg==\n", "text/plain": ["<Figure size 432x288 with 1 Axes>"]}, "metadata": {"needs_background": "light"}, "output_type": "display_data"}], "source": ["import pandas\n", "\n", "\n", "def graph_XY(curves, xlabel=None, ylabel=None, marker=True,\n", "             link_point=False, title=None, format_date=\"%Y-%m-%d\",\n", "             legend_loc=0, figsize=None, ax=None):\n", "    if ax is None:\n", "        import matplotlib.pyplot as plt  # pylint: disable=C0415\n", "        fig, ax = plt.subplots(1, 1, figsize=figsize)\n", "\n", "    smarker = {(True, True): 'o-', (True, False): 'o', (False, True): '-',\n", "               # (False, False) :''\n", "               }[marker, link_point]\n", "    has_date = False\n", "    for xf, yf, label in curves:\n", "        ax.plot(xf, yf, smarker, label=label)\n", "    ax.legend(loc=legend_loc)\n", "    return ax\n", "\n", "\n", "def draw_variance(sample) :\n", "    avg = 0.\n", "    std = 0.\n", "    n   = 0.\n", "    w   = 1.\n", "    add = [] \n", "    for i,x in enumerate(sample) :\n", "        x    = float (x)\n", "        avg += x * w\n", "        std += x*x * w\n", "        n   += w\n", "        val  = (std/n - (avg/n)**2)**0.5\n", "        add.append ( [ i, avg/n, val] )\n", "    \n", "    print(add[-1])\n", "    table = pandas.DataFrame(add, columns=[\"index\", \"avg(n)\", \"std(n)\"])\n", "    return graph_XY([\n", "            [table['index'], table[\"avg(n)\"], \"avg(n)\"],\n", "            [table['index'], table[\"std(n)\"], \"std(n)\"],\n", "        ], marker=False, link_point=True)  \n", "\n", "draw_variance(sample);"]}, {"cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [{"data": {"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de3hc9X3n8ff3zE0jWTfL8gULW3ZwwIQQxxGEBMJDoWmApgnZXhKebANNUod9kjZturshyfM0yTbZzW6habvbhtKEArlQEkgakiZZEpJAuy0E2Rgw2OBLfJERkmxZtmTd5vLdP+bIyEbGsmbkkc75vJ5nHs385syc7+GYz/zmd37njLk7IiISLUG1CxARkcpTuIuIRJDCXUQkghTuIiIRpHAXEYmgZLULAFi0aJG3t7dXuwwRkXll48aNB9y9darn5kS4t7e309nZWe0yRETmFTPbc7LnNCwjIhJBCncRkQhSuIuIRNCcGHMXEam0XC5HV1cXo6Oj1S6lbDU1NbS1tZFKpab9GoW7iERSV1cX9fX1tLe3Y2bVLmfG3J2DBw/S1dXFqlWrpv06DcuISCSNjo7S0tIyr4MdwMxoaWk57W8gCncRiaz5HuwTZrIdkQn3YtG59/G9HBnNVbsUEZGqi0y4/+iZF/n4/U/zzcf3VbsUEREAFixYcNzjO++8k4985CMA3Hbbbdx9990A3Hjjjdx3330AXHHFFRU5qTMSB1Tdnb97eCcAT+wdqHI1IiKndtNNN83q+0ei5/7YL/t5susw9ZkkG/ccqnY5IiKn9JnPfIZbbrll1t4/Ej33v3t4Jy11aTZcvpr/8cNtvDAwwllN2WqXJSJzxGe/9wzPvnCkou95/lkNfPo3XvOKy4yMjLBu3bpjj/v7+3nHO95R0TpOZt6H+3MvDvKz5/r4k7e+mje/ahEAG/ccUriLSNVls1k2b9587PGdd955xi6SOO/D/fZHdpFNJfiPl6xkQU2SmlTApr2H+I3XnVXt0kRkjjhVDzuK5vWYe/fhEb67eT/vvuhsmuvSpBIBr2trYpPG3UUk5uZ1uA+O5nnj6oV84LKXTsldv7KZZ144wmiuUMXKRESqy9y92jXQ0dHhlRqH+smzPXzw7k6++aE3cfGqhRV5TxGZf7Zu3cratWurXUbFTLU9ZrbR3TumWn5e99yn8voVTQBs2quhGRGJr8iFe8uCDKsW1WncXURiLXLhDqXe+6a9h5gLQ04iUj1RyYCZbMcpw93M7jCzXjPbMqntXjPbHN52m9nmsL3dzEYmPXfbaVdUAW9Y2cyBoXH29Y9UY/UiMgfU1NRw8ODBeR/wE9dzr6mpOa3XTWee+53A/wHunrSyd0/cN7NbgcOTlt/p7uuoovUrmgHYuLefFS211SxFRKqkra2Nrq4u+vr6ql1K2SZ+iel0nDLc3f0RM2uf6jkrXWT4d4ArT2uts2zFwlKg9x4Zq3IlIlItqVTqtH65KGrKHXN/C9Dj7tsnta0ysyfM7GEze8vJXmhmG8ys08w6K/3Jmk6WNms8X6zo+4qIzBflhvv1wD2THncDK9z99cDHgG+YWcNUL3T32929w907WltbyyzjeMnAMINcQeEuIvE043A3syTwH4B7J9rcfczdD4b3NwI7gVeXW+QMaiOdCBhTuItITJXTc/9VYJu7d000mFmrmSXC+6uBNcCu8kqcmXQi0LCMiMTWdKZC3gP8O3CumXWZ2QfCp97D8UMyAJcDT5nZk8B9wE3u3l/JgqcrnVS4i0h8TWe2zPUnab9xirb7gfvLL6t86WSgMXcRia1InqEK6rmLSLxFNtxTiYBx9dxFJKYiG+46oCoicRbdcE8GjBfm9zUlRERmKtrhntevMYlIPEU33DUsIyIxFt1wT+qAqojEV3TDPRGQy2vMXUTiKbrhrp67iMRYZMM9pTF3EYmxyIZ7OhkwpnAXkZiKbLhndG0ZEYmxyIa7ri0jInEW2XBPJUwHVEUktiIb7ulEgkLRKRQ1HVJE4ie64R7+SLbG3UUkjiIf7poxIyJxFN1wTxiADqqKSCxN5zdU7zCzXjPbMqntM2a238w2h7drJz33CTPbYWbPmdnbZqvwU5noueugqojE0XR67ncCV0/R/kV3XxfefgBgZudT+uHs14Sv+VszS1Sq2NNxbMxdPXcRiaFThru7PwL0T/P93gn8o7uPufsvgR3AxWXUN2PpROkzRT13EYmjcsbcP2JmT4XDNs1h23Jg36RlusK2lzGzDWbWaWadfX19ZZQxtZTG3EUkxmYa7l8CXgWsA7qBW8N2m2LZKSeau/vt7t7h7h2tra0zLOPkNFtGROJsRuHu7j3uXnD3IvD3vDT00gWcPWnRNuCF8kqcGc1zF5E4m1G4m9mySQ/fBUzMpHkAeI+ZZcxsFbAG+EV5Jc5MZmK2jHruIhJDyVMtYGb3AFcAi8ysC/g0cIWZraM05LIb+BCAuz9jZt8EngXywIfdvSq/Up1KKNxFJL5OGe7ufv0UzV95heU/D3y+nKIqQfPcRSTOInyGqnruIhJf0Q139dxFJMaiG+7quYtIjEU33DVbRkRiLPrhrmEZEYmh6IZ7QhcOE5H4imy4JwLDTD13EYmnyIa7mZFOBBpzF5FYimy4Q2ncXRcOE5E4inS4Z5KBLhwmIrEU6XBPaVhGRGIq0uGeTgY6oCoisRTtcFfPXURiKtrhrjF3EYmpSId7KqHZMiIST5EO93RSwzIiEk+RDveMDqiKSExFOtzTCY25i0g8nTLczewOM+s1sy2T2v7czLaZ2VNm9h0zawrb281sxMw2h7fbZrP4U9E8dxGJq+n03O8Erj6h7cfABe5+IfA88IlJz+1093Xh7abKlDkzGnMXkbg6Zbi7+yNA/wltD7p7Pnz4KNA2C7WVTeEuInFViTH39wM/nPR4lZk9YWYPm9lbTvYiM9tgZp1m1tnX11eBMl6udIaqz8p7i4jMZWWFu5l9CsgDXw+buoEV7v564GPAN8ysYarXuvvt7t7h7h2tra3llHFSpTNUC7Py3iIic9mMw93MbgDeDrzX3R3A3cfc/WB4fyOwE3h1JQqdCV1bRkTiakbhbmZXAx8H3uHuw5PaW80sEd5fDawBdlWi0JnQtWVEJK6Sp1rAzO4BrgAWmVkX8GlKs2MywI/NDODRcGbM5cB/M7M8UABucvf+Kd/4DEgnA4oOhaKTCKxaZYiInHGnDHd3v36K5q+cZNn7gfvLLapSUuGPZI/ni2TTiSpXIyJy5kT7DNXkS+EuIhInsQj3sYJmzIhIvEQ63DPhsExOc91FJGYiHe6pZOkgqoZlRCRuIh3u6UTpIKrCXUTiJtrhrgOqIhJT8Qh3naUqIjET6XBPJTTmLiLxFOlwz6jnLiIxFelw1wFVEYmraId7cmKeu8JdROIl0uGuMXcRiatIh7umQopIXMUi3Mc0LCMiMRPpcM/ogKqIxFSkw33i2jI6oCoicRPpcE8nNOYuIvEU6XBPJgICU7iLSPxMK9zN7A4z6zWzLZPaFprZj81se/i3OWw3M/trM9thZk+Z2frZKn460slAZ6iKSOxMt+d+J3D1CW03Aw+5+xrgofAxwDXAmvC2AfhS+WXOXCoRqOcuIrEzrXB390eA/hOa3wncFd6/C7huUvvdXvIo0GRmyypR7Exk1HMXkRgqZ8x9ibt3A4R/F4fty4F9k5brCtuOY2YbzKzTzDr7+vrKKOOVpdVzF5EYmo0DqjZF28t+xNTdb3f3DnfvaG1tnYUyStJJhbuIxE854d4zMdwS/u0N27uAsyct1wa8UMZ6ypJKBJrnLiKxU064PwDcEN6/AfjupPb3hbNmLgEOTwzfVIN67iISR8npLGRm9wBXAIvMrAv4NPAF4Jtm9gFgL/Db4eI/AK4FdgDDwO9VuObToqmQIhJH0wp3d7/+JE9dNcWyDny4nKIqKZ0IGFPPXURiJtJnqEKp564xdxGJm+iHu6ZCikgMRT/cdUBVRGIoHuGuYRkRiZnIh3sqEZBTz11EYiby4a6eu4jEUfTDXVMhRSSGIh/uGR1QFZEYiny469oyIhJHkQ/3dDKg6JBXwItIjMQi3AEdVBWRWIl+uCfCcNe4u4jESOTDPaWeu4jEUOTDPaOeu4jEUOTD/diYu8JdRGIkPuGuYRkRiZHIh3sqHJbJ5V/2G90iIpEV+XB/qedeqHIlIiJnzrR+Zm8qZnYucO+kptXAnwJNwO8DfWH7J939BzOusEwTUyF1fRkRiZMZh7u7PwesAzCzBLAf+A6lH8T+orvfUpEKy6QDqiISR5UalrkK2Onueyr0fhUz0XPPFTTmLiLxUalwfw9wz6THHzGzp8zsDjNrnuoFZrbBzDrNrLOvr2+qRSpCPXcRiaOyw93M0sA7gG+FTV8CXkVpyKYbuHWq17n77e7e4e4dra2t5ZZxUjqgKiJxVIme+zXAJnfvAXD3HncvuHsR+Hvg4gqsY8bUcxeROKpEuF/PpCEZM1s26bl3AVsqsI4ZSyUMULiLSLzMeLYMgJnVAm8FPjSp+X+Z2TrAgd0nPHfGZRIJAMZ1QFVEYqSscHf3YaDlhLbfLauiCtOwjIjEUXzOUFW4i0iMRD7cE4GRTgQM5/LVLkVE5IyJfLgDLG2soXtgtNpliIicMbEI9+VNWfYPjFS7DBGRMyYe4d6cpevQcLXLEBE5Y+IR7k1ZegfHdFBVRGIjFuHe1pzFHboPa2hGROIhFuG+vDkLwP5DCncRiYdYhHtbUy0AXTqoKiIxEYtwX9pYgxl0qecuIjERi3BPJwOWNtRoWEZEYiMW4Q4Tc901HVJE4iE+4d6sE5lEJD7iE+5NWboHRikUdelfEYm++IR7c5Z80ek5MvU1Znb0DrKzb+gMVyUiMjtiE+5tzaXpkFMNzRSKzg13PM7HvvnkmS5LRGRWxCbclzed/ESmn27rZf/ACFv2H+bomC4NLCLzX/zCfYqe+1cf3UNgpR78k/sGznRpIiIVV3a4m9luM3vazDabWWfYttDMfmxm28O/zeWXWp5sOkFLXfplV4fcfeAojzzfx/svXYUZPL77UJUqFBGpnEr13H/F3de5e0f4+GbgIXdfAzwUPq660qV/j++5f+3RPSQD4/cvX825S+rp3NNfpepERCpntoZl3gncFd6/C7hultZzWtpOmOs+Ml7gWxu7eNtrlrKkoYaL2heyac8h8gVdGlhE5rdKhLsDD5rZRjPbELYtcfdugPDv4hNfZGYbzKzTzDr7+voqUMapLW/Ksv/QCO6lue7fe+oFDo/k+N03rQSgo72Zo+MFtr04eEbqERGZLZUI90vdfT1wDfBhM7t8Oi9y99vdvcPdO1pbWytQxqktb8oyli9yYGickfECtz+yizWLF/DGVQsBuKi99Ldz9/FDM0Wd+CQi80zZ4e7uL4R/e4HvABcDPWa2DCD821vueiph+aS57p/6p6fZ2TfEJ399LWYGwFlNWc5qrOHxPS8dVP3Ud57m6r96hMHRXFVqFhGZibLC3czqzKx+4j7wa8AW4AHghnCxG4DvlrOeSmkLf7Tj1gef49ub9vPRq9bwK+ceP2LU0b6Qzt39uDsPP9/H1x/by/M9Q3z2e89Wo2QRkRkpt+e+BPhXM3sS+AXwz+7+I+ALwFvNbDvw1vBx1U38ItO/bD/Ar5zbyh9eueZly1zU3kzPkTG29w7xqe88zepFdWy4fDX3beziR1u6z3TJIiIzkiznxe6+C3jdFO0HgavKee/Z0FCToqk2RUNNir989+sJAnvZMm9YWRp3/09f20jXoRHu3XAJ61c28+87D/KJbz/N+hXNLG6oOdOli4iclticoTrh9t/t4J4Nl9BYm5ry+XOX1lOfSbKz7yjXX3w2b1zdQioR8MV3r2MkV+Dmbz99hisWETl9sQv3i1ctPHYpgqkkAuONqxfSWp/h5mvWHms/Z/ECNlz+Kn66rZf+o+NnolQRkRmLXbhPx5//1uv4/h9cRmP2+N79pa9qAeCJvbpEgYjMbQr3KTTXpVkyxbj6hW1NJANj4x6Fu4jMbQr305BNJzj/rAY2qecuInOcwv00rV/RzJP7Duv6MyIypyncT9P6lc2M5HT9GRGZ2xTup2n9iiYAjbuLyJymcD9Ny5uyLK7PaNxdROY0hftpMjPesLJZ4S4ic5rCfQbWr2hmX/8IvYOj1S5FRGRKCvcZWL+y9JOwm/box7RFZG5SuM/ABcsbSCcCDc2IyJylcJ+BTDLBBcsb2KQZMyIyRyncZ2j9imae2DfAhrs7ueNff8nzPZr3LiJzR1nXc4+zG97czpHRHI/u6ufBZ3sAWHd2E+9700qufe0yalKJKlcoInFm7tX/8eeOjg7v7Oysdhkztn9ghP+75UW+9tgedvUdpaUuzR9ceQ7vvWQlqYS+HInI7DCzje7eMeVzCvfKcXf+346D/O3Pd/BvOw+yalEdH7/6PN72miXHfoRbRKRSXincZ9ytNLOzzexnZrbVzJ4xs4+G7Z8xs/1mtjm8XTvTdcw3ZsZlaxbx9Q++kTtu7CARGDd9bSPX/vW/8t3N+49dbKxYdA4P55gLH6wiEk0z7rmb2TJgmbtvMrN6YCNwHfA7wJC73zLd94pKz/1E+UKR7zyxn797ZBc7eodY2lBDMmH0HBklV3D+9O3n8/7LVlW7TBGZp16p5z7jA6ru3g10h/cHzWwrsHym7xdFyUTAb3eczW+ub+MnW3u4f1MXtekkSxtrePyX/dz64HNc89qlLGs8+c/+iYjMREXG3M2sHXgEuAD4GHAjcAToBP7E3V82IdzMNgAbAFasWPGGPXv2lF3HfLL34DBv/eLDXLV2MX/73jdUuxwRmYdmZcx90psvAO4H/sjdjwBfAl4FrKPUs791qte5++3u3uHuHa2treWWMe+saKnlD648hx88/SI/f6632uWISMSUFe5mlqIU7F93928DuHuPuxfcvQj8PXBx+WVG0+9fvprVi+r49APPsH9ghD0Hj7K9Z1C/8iQiZZvxmLuV5vZ9Bdjq7n8xqX1ZOB4P8C5gS3klRlcmmeDPrruA9375MS79wk+PtXesbOYffu8i6mtS036vQtEZGB4nk0qwIKNz00TirpzZMpcB/wI8DUx0NT8JXE9pSMaB3cCHJoX9lKI6W2a6Hnm+j32HhsmmEvQfHecLP9zGBcsbuev9F9OYnTrg9x4c5idbe3hoWw9buwcZGB6n6FCbTvCff+1cbnhzO4lAc+tFokwnMc0zDz7zIh/+xibOXVrPzVevZWBknP6j4+w9OMz23iG29wzywuHSteTPWbyAi9qbaV2QYWFdmp8/38fPn+vjwrZGPnfdBVzY1nTS9Rwdy/Ns9xG29wyxs2+IvsEx3n/ZKtadffLXiMjcoXCfh372XC8f+upGxvMvjb9nkgHnLF7AmsULuLCtiavWLmZlS91xr3N3vv9UN5/93jMcGBpndWsdbz1/CR0rFzIwPE7PkVH29g/zVNdhnu8ZpBju/ppUQDoRMDxe4OZrzuMDl63SWbUic5zCfZ7afeAoLwyM0BL2yhfWpac91HJ4OMc/bd7PT7b28O87D5IvvrSfW+rSXLC8kXVnN3FhWyOvXlLP8qYsg6N5/st9T/Lgsz1cdd5ivvCbF9Jan5mtzRORMincY+7IaI7tPUMsWpBmSUPNK16x0t256992899/sI1MMuAPr1rDDW9uJ53UBdBE5hqFu5y2XX1D/Nn3n+Vnz/WxurWOXzt/KStbalm5sJaO9oUKe5E5QOEuM/bTbT3c+uDzPN8zSK5Q+rdy6Tkt3HHjRWSSuma9SDXNyrVlJB6uPG8JV563hELR6T48wk+e7eEz33uWP753M//7+vWabikyRyncZVoSgdHWXMuNl64iX3Q+989baa7dwueuu2Bas2pGcwX6BscoFJ2zmrKvOKwzMl6gd3CUniNj9B8dA4xEYCQDI5sunaTVUJOirTlLoA8XkSkp3OW0ffAtqzkwNM5tD+/k0V0HacymqMskj4VufU2SsXyR/QMjvDAwQvfhUQ6P5I69PjBY1piltT6Du5MvOuP5IodHchweyTGWn97lF+prknSsbKajfSFXnreY85bWa/qmSEhj7jIj7s5tD+9i875DDI8XGBrLMzSa58hojsHRPKlEwPKmLGc1ZVnWWMPi+gyLGzIEZuw7NMK+/mEODI2RCIyEGalEQGM2RWNtisZsiiUNNSxpKE0BNYyiO7lCkZFwXQPDOZ7YN0Dn7n629w4B0N5Sy9tes5SGbIr+o+McOjrOWL5IoegU3EvrSZbm8wMUww+WhbUp1q9s5g0rm1nelNUHhMwbOqAqkdY3OMaPn+3hh1u6j83pr00naK5NU5MKSARGYIY7jBeKx04MSwSl4Z6eI6MMjxcAaMymaK3P0LogQ206wWD4gQWls4HXLmvgnMULWNJQQ2t9hkUL0jqwLFWjcJfYGB7PE5i94lz+E+ULRba9OEjn7n529A1xYHCcA0NjDI8XaMgmqa9JUSw6z/UM0nVo5LjXBgavXlLP69qaWLusHoDRfJHRXIGR8QLD4S1fLFL00reFZGBkkgHpZEAyCI59yEwMURWKzliuyHCuwMh4nqJDOlFa3gzG86UPqGw6wdplDaxdVs/ZzbXkCs5YvkDRnZpUgtp0aaispS6tYxMRpXAXqZDB0Ry7DwzTNzRK3+AYXYdGeKrrME92DTAwnDtu2UwyoDZdCtlkovTtITCOHWMYyxfJF0qhny8WCcyOBX3ptUmyqQRB8FKgFx3S4dDS4FiOff0jJ6n0JamEsawxy+L6DA7kCkVyBadQLA1ZFb3Uli+Uhr4mPmAKRWdhXZrlTVmWN2epSSUohkNc2VSC5toUTbVpUgljNFdkLF8o1Vlw8oUi44UiuUJpO91L/z1qUgmy6QSN2RRN2RQN2RQ1qYBMMkFdJsmqljoaa6d/NdS401RIkQqpr0nx2rZGoPG4dnenb7B0DKEmlaAmlTgj00SPjObY1j1I9+ERMslSSAaBMTJeYCSX58hInu7Do7wwMHKsvoaaJIkgIBl+kASBkQpKxz2SCQvbAwKDg0fH2T8wwuO7+xnLF0mEH1DDuQKHR3JM1TdMJwJSCSMZfttIJwKCgNIHQG7im8zJO5Wt9RlWL6pjYV2axvADYKLWRGCkkwE1ydJ/44V1aRY3ZFhcnyGdDCgWS9+OzCAwwwxq00nq0onYHUtRuItUgJmxuKHmjK+3oSbFxasWnvH1AhSLzpHRHPmiH/tgSSXslCHq7gyPFxgYyTE4mmMsVxrGOjKaZ1ffEDt6h/jlgaNs7x3icLhMvlD6xjDTgYZUwmjMpmmtz3BWYw1nNWVpmvQNwcMhs6KD4xilD7GJb1tBeOB/8uoDK+33iQ+RieXTyeDYf4/G2hQt4XWhGrMp6tLJMzZEpnAXkRkJAqOpNn3arzMz6jJJ6jJJ4MQfh1/yiq8tFp3xQnhMI1fg4NA4fYNj9A6Okiv4sTAGKDoU3BkeyzMwkmNguLTs/oFROvcc4shojskxmwgmfTCFIT8xbFVJpaG6BJlkgkwy4Kq1i/nUr59f2ZWgcBeReSQIjJqgNCTTROl8iTNh4liDUfpwcnecsLdfLH0QePhhkguPp4zmSt9ODg6Nc3BojMHRPINjeY6O5RnNFRgLl1s6S9ugcBcROYUgMILj+vlzf/xel/YTEYmgWQt3M7vazJ4zsx1mdvNsrUdERF5uVsLdzBLA3wDXAOcD15tZ5Y8YiIjIlGar534xsMPdd7n7OPCPwDtnaV0iInKC2Qr35cC+SY+7wrZjzGyDmXWaWWdfX98slSEiEk+zFe5THUo+braou9/u7h3u3tHa2jpLZYiIxNNshXsXcPakx23AC7O0LhEROcFshfvjwBozW2VmaeA9wAOztC4RETnBrF0V0syuBf4SSAB3uPvnX2HZPmBPGatbBBwo4/XzURy3GeK53drm+Djd7V7p7lOOa8+JS/6Wy8w6T3bZy6iK4zZDPLdb2xwfldxunaEqIhJBCncRkQiKSrjfXu0CqiCO2wzx3G5tc3xUbLsjMeYuIiLHi0rPXUREJlG4i4hE0LwO9zhcVtjMzjazn5nZVjN7xsw+GrYvNLMfm9n28G9ztWudDWaWMLMnzOz74eNVZvZYuN33hifJRYaZNZnZfWa2Ldznb4rDvjazPw7/fW8xs3vMrCaK+9rM7jCzXjPbMqltyv1rJX8d5ttTZrb+dNY1b8M9RpcVzgN/4u5rgUuAD4fbeTPwkLuvAR4KH0fRR4Gtkx7/T+CL4XYfAj5Qlapmz18BP3L384DXUdr2SO9rM1sO/CHQ4e4XUDrx8T1Ec1/fCVx9QtvJ9u81wJrwtgH40umsaN6GOzG5rLC7d7v7pvD+IKX/2ZdT2ta7wsXuAq6rToWzx8zagF8Hvhw+NuBK4L5wkUhtt5k1AJcDXwFw93F3HyAG+5rST35mzSwJ1ALdRHBfu/sjQP8JzSfbv+8E7vaSR4EmM1s23XXN53A/5WWFo8bM2oHXA48BS9y9G0ofAMDi6lU2a/4S+K9AMXzcAgy4ez58HLV9vhroA/4hHIr6spnVEfF97e77gVuAvZRC/TCwkWjv68lOtn/Lyrj5HO6nvKxwlJjZAuB+4I/c/Ui165ltZvZ2oNfdN05unmLRKO3zJLAe+JK7vx44SsSGYKYSjjG/E1gFnAXUURqSOFGU9vV0lPXvfT6He2wuK2xmKUrB/nV3/3bY3DPxFS3821ut+mbJpcA7zGw3pSG3Kyn15JvCr+4QvX3eBXS5+2Ph4/sohX3U9/WvAr909z53zwHfBt5MtPf1ZCfbv2Vl3HwO91hcVjgcZ/4KsNXd/2LSUw8AN4T3bwC+e6Zrm03u/gl3b3P3dkr79qfu/l7gZ8BvhYtFarvd/UVgn5mdGzZdBTxLxPc1peGYS8ysNvz3PrHdkd3XJzjZ/n0AeF84a+YS4PDE8M20uPu8vQHXAs8DO4FPVbueWdrGyyh9FXsK2BzerqU0/vwQsD38u7Datc7if4MrgO+H91cDvwB2AN8CMtWur8Lbug7oDPf3PwHNcdjXwGeBbcAW4KtAJor7GriH0nGFHKWe+QdOtn8pDcv8TZhvT1VNYVYAAAA6SURBVFOaTTTtdenyAyIiETSfh2VEROQkFO4iIhGkcBcRiSCFu4hIBCncRUQiSOEuIhJBCncRkQj6/8vOvqF0NUv9AAAAAElFTkSuQmCC\n", "text/plain": ["<Figure size 432x288 with 1 Axes>"]}, "metadata": {"needs_background": "light"}, "output_type": "display_data"}], "source": ["def draw_hill_estimator (sample) :\n", "    res = hill_estimator(sample)\n", "    table = DataFrame(res, columns=[\"d\", \"hill\"])\n", "    return graph_XY(\n", "        [[table['d'], table['hill'], \"Hill\"],],\n", "        marker=False, link_point=True)\n", "\n", "draw_hill_estimator(sample);"]}, {"cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [{"name": "stderr", "output_type": "stream", "text": ["c:\\python372_x64\\lib\\site-packages\\pandas\\core\\series.py:679: RuntimeWarning: divide by zero encountered in log\n", "  result = getattr(ufunc, method)(*inputs, **kwargs)\n"]}, {"data": {"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD5CAYAAAA3Os7hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXxU1f3/8ddnJvtCgBAgQMK+LwkQNgEFFEVUQEVBqhUKolTcqrZqba18bX8uVVwrxQ2wVlAQKgpurFoECRQEWcMmYQ1bIITs5/fHDRhCSIZkMneWz/PxmMedmXvn3rej8/Hk3HPvEWMMSimlfJ/D7gBKKaXcQwu6Ukr5CS3oSinlJ7SgK6WUn9CCrpRSfkILulJK+YmgijYQkTBgORBavP1sY8xTpbYJBWYAXYGjwAhjzO7y9lunTh3TpEmTyqVWSqkAtWbNmiPGmLiy1lVY0IFcYIAxJktEgoHvRGShMWZliW3GAseNMS1EZCTwHDCivJ02adKE1NRUF/8RlFJKAYjInoutq7DLxViyil8GFz9KX400FJhe/Hw2cKWISCWyKqWUqiSX+tBFxCki64DDwNfGmFWlNmkI7AUwxhQAmUBsGfsZLyKpIpKakZFRteRKKaXO41JBN8YUGmOSgUZAdxHpUGqTslrjF9xTwBgz1RiTYoxJiYsrswtIKaVUJbnSh36OMeaEiCwFBgEbS6xKBxKAdBEJAmKAY+4KqZQKHPn5+aSnp5OTk2N3FFuFhYXRqFEjgoODXf6MK6Nc4oD84mIeDlyFddKzpE+BO4HvgeHAYqN3/VJKVUJ6ejrR0dE0adKEQD0VZ4zh6NGjpKen07RpU5c/50qXSzywRER+BFZj9aF/JiKTRGRI8TbvALEikgb8DnjsEvMrpRQAOTk5xMbGBmwxBxARYmNjL/mvlApb6MaYH4HOZbz/5xLPc4BbLunISil1EYFczM+qzHfge1eKnj4KCx+DvGy7kyillFfxvYK+aymsmgLvXQsn99udRinlh+bOnUtycvJ5D4fDwQcffMDw4cMr/Pyrr75K27Zt+dWvfuWBtL8Qu85dpqSkmEpfKbp1IcwZByFRcNuH0LCLe8MppWyzefNm2rZta3eM80ydOpUPPviAJUuW4HBU3A5u06YNCxcuvKQTmmUp67sQkTXGmJSytve9FjpwOL4/jP0KnCFWS33jJ3ZHUkr5qW3btjFp0iTef/99fv75Zzp0sC7DmTZtGkOHDmXQoEG0bt2ap59+GoB77rmHnTt3MmTIECZPnuzRrJc0Dt0bLNhwgN99tI7/G9qBW+5aDLNuh9ljIGMr9HsM9GSKUn7j6fk/sWn/Sbfus12DGjx1Q3uXts3Pz2fUqFH8/e9/JzExkd27d5+3/ocffmDjxo1ERETQrVs3rrvuOqZMmcIXX3zBkiVLqFOnjluzV8TnWugpTWqRnFCTR2f/yO+/2E/OqLmQNAqWPWsVdj1ZqpRykz/96U+0b9+ekSNHlrl+4MCBxMbGEh4ezk033cR3333n4YTn87kWet3oMP41tgeTv9nGG0t2sGHfSf4x6gWa1m0DXz8Fx3fDyH9DjQZ2R1VKVZGrLenqsHTpUubMmcPatWsvuk3poYV2D7f0uRY6QJDTwaPXtOG90d04kHmGG17/Lwtr3GoV8iPb4a0BsO/i/xKUUqo8x48fZ8yYMcyYMYPo6OiLbvf1119z7Ngxzpw5w7x58+jdu7cHU17IJwv6Wf3b1OXz+/vSvG4UEz5Yy9PbG5N/50JwBMF7g+GnuXZHVEr5oClTpnD48GEmTJhw3tDFWbNmnbddnz59uOOOO0hOTubmm28mJaXMwSce45vDFkvJKyjibws2M23Fbjon1uQfwxKIXzgO9q6Cfk/AFb/Xk6VK+QhvHLZYlmnTppGamsrrr79ebccIiGGLpYUEOfjLkPa8MaoL2w9lMfjtLSzr9Q4k3QZL/wZzxkL+GbtjKqVUtfKLgn7WdZ3i+XRib+rVCGP0+z/yYsSDFF35tDVO/b3BcPKA3RGVUn5i9OjR1do6rwy/KugAzeKimPvb3gzv0ojXluzg9i09yRw63Rqn/tYA2L/O7ohKKVUt/K6gA4SHOHnhliSeH96JNXuOM3BBJD8O+ggcTnh3EPw0z+6ISinldn5Z0M+6NSWBeff2JjI0iBvnnGR6h3cx9TvBx3fCsudB5+BQSvkRvy7oAG3ja/DpxN5c074eTy3KYILzKfLa3wpL/qonS5VSfsXvCzpAdFgwb4zqwlM3tGNRWiYDdozkQLfHfjlZeuqg3RGVUl7m0KFDjBo1imbNmtG1a1d69erF3LnefW1LQBR0sC7JHdO7KR/d3YsiA1esSGZJ58mYjK0wtb+eLFVKnWOMYdiwYVx++eXs3LmTNWvWMHPmTNLT0+2OVq6AKehndU6sxef39+WyFrGM+b4uzzV4hSJxWCdLN/3H7nhKKS+wePFiQkJCuOeee86917hxY+677z52795N37596dKlC126dGHFihWAde+X66+//tz2EydOZNq0aQA89thjtGvXjk6dOvHII48A8PHHH9OhQweSkpK4/PLL3ZLb527O5Q61IkN4985uvLlsBy9+tZXU2n/l/dhXCf/o19D/Sbj8Eb2yVClvsPAxOLjBvfus3xGufbbcTX766Se6dCl74py6devy9ddfExYWxvbt27ntttso76r3Y8eOMXfuXLZs2YKIcOLECQAmTZrEl19+ScOGDc+9V1UB10I/y+EQ7u3fgn+N68Hu3Ch67n+QPY2GwJJnrNmQ0lPhwI/W+PVju6zp7k4fgZyTUJCrI2SUCiD33nsvSUlJdOvWjfz8fO666y46duzILbfcwqZNm8r9bI0aNQgLC2PcuHF88sknREREANC7d29Gjx7NW2+9RWFhoVtyBmQLvaTLmtdhwf19mPjh/7gibQRTmsRzzcapyMbZFX/YGfLLIyi0eBkGweEQHAHBYdayrPecIdZNxJzB1tLhLF4GlXgdXOp1EIRGQ9121n6U8ncVtKSrS/v27ZkzZ86512+88QZHjhwhJSWFyZMnU69ePdavX09RURFhYdZvMSgoiKKionOfycnJOff+Dz/8wKJFi5g5cyavv/46ixcvZsqUKaxatYrPP/+c5ORk1q1bR2xsbJVyB3xBB6hbI4x/j+vBi19v456lwrV1OzKpX03iwrFa44X5UJgLBXnWsjDvl+cFedbrs88LcqyhkPnZkHsKsg5bz/NzipdnrG2rwhEEddtCg87WIz4Jwmv/Uvydwdb/NEIi3PL9KBVoBgwYwBNPPMGbb77JhAkTAMjOtibPyczMpFGjRjgcDqZPn36udd24cWM2bdpEbm4uOTk5LFq0iD59+pCVlUV2djaDBw+mZ8+etGjRAoAdO3bQo0cPevTowfz589m7d68WdHcJcjr4w6A2pDSuxUOz1jHsq2Bm3d2TRrWqoSgWFVr/EygqKH4UlnhexuvCEs+zj8KBdbD/f7B5PqydUfYxxAlNL4cON0Gb6yGitvv/OZTyUyLCvHnzeOihh3j++eeJi4sjMjKS5557ji5dunDzzTfz8ccf079/fyIjIwFISEjg1ltvpVOnTrRs2ZLOnTsDcOrUKYYOHUpOTg7GmHPzjD766KNs374dYwxXXnklSUlJVc/tD7fPdbeN+zIZ9dZKYiKCmTW+Fw1qhtsdqWzGwIk91kmjvNPWXxJnC39mujVq5/guq9XefAAM+BPEd7I7tVLl8pXb53rCpd4+t8IWuogkADOA+kARMNUY80qpbfoB/wF2Fb/1iTFm0iWn9xIdGsbw/tge3P72Kka9tZKZ43tRP8YL+6xFoFYT61GWq/4CB9ZbE32s+wCm9oPLJsIVj2l3jFJ+yJVRLgXAw8aYtkBP4F4RaVfGdt8aY5KLHz5bzM9KSqjJ9LHdOZKVx6i3VnL4ZI7dkS6dCDRIhoFPw8TVkDwK/vsKvHkZrJ8J27+BXcth72o4sdfq2lFK+awKW+jGmAPAgeLnp0RkM9AQKH+sjh/okliLaWO68et3f2DU26v48K6exEWH2h2rcsJrwdDXodMImP8AzL37wm3EaU2uHV4LwmtCjUbQ92Go08LzeVVAM8bYPuGy3SrTHX5Jfegi0gRYDnQwxpws8X4/YA6QDuwHHjHG/FTG58cD4wESExO77tmz55ID22HVzqOMfm81ibUj+HB8T2pHhtgdqWoKcq3x9QW51oibvGw4uQ8y90LmPjhzHHJOwKFN1vq+j0CfB62hmUpVs127dhEdHU1sbGzAFnVjDEePHuXUqVM0bdr0vHXl9aG7XNBFJApYBvzVGPNJqXU1gCJjTJaIDAZeMca0LG9/3nxStCwr0o4wZtpqmsVF8eFdPagZ4eNF3RWnDsGXj8PGOVarvfmV0Ooaa9SM9sGrapKfn096evq5cdyBKiwsjEaNGhEcHHze+1Uu6CISDHwGfGmMecmF7XcDKcaYIxfbxtcKOsC32zMYOz2VVvWi+GBsT2Iigiv+kD/YuRTWz4K0r+F0BoTGQNII6wKn0GhoeTWE1bA7pVIBoUoFXay/eaYDx4wxD15km/rAIWOMEZHuwGygsSln575Y0AGWbD3M3TPW0DY+mvfH9aBGWIAUdYCiIvj5e1jznjUksjDPej/xMhj9OTgC9k4SSnlMVQt6H+BbYAPWsEWAJ4BEAGPMFBGZCEzAGhFzBvidMWZFefv11YIOsGjzIe751xo6NIxhxm+6Ex1IRf2s/DNWX/uWz2HBI3DN36DXvXanUsrvuaUP3d18uaADfPnTQe79YC3JCTWZ/pvuRIYG6EW3xsDMUZC2CK5/yRol07gPBAXAOQalbFBeQde/kSvpmvb1efW2zvxv7wnGTFtNdl6AjuEWgetftk6a/udeeP9GmPUrHdOulA20oFfB4I7xvDwimdTdxxg3PZUzee65BabPia4H96+Fe3+Aq5+B7V/BZw/qLYaV8rAA7SdwnxuSGlBYZHjoo3WMfz+Vt36dQliw0+5YnhcSCXGtrUdOJix/AbKPWRcz6Y3BlPIIbaG7wbDODXlheBLfpR3hjndWkX482+5I9ur/R+sk6fav4LWu8J+JcGyn3amU8nta0N1keNdGvDKyM5sPnOLal79l3v/22R3JPiLWiJdxX0OLK2HjJzD7N9oFo1Q104LuRkOSGrDwgb60rh/Ng7PWcf+H/yMzO9/uWPZp0BlufhsGP2/dv33rArsTKeXXtKC7WULtCGaO78kjV7diwYYDXPvKclbsuOgFs4Gh00io3RwWP6NdL0pVIy3o1SDI6WDigJbMmXAZYcFOfvX2Kv7fgs3kFgToKBhnkHVv9sOb4dXOMP0GOJJmdyql/I4W9GqUlFCTz+7vw6juifxz+U6GvbGCbYdO2R3LHu2GWEMbr3oa9q+Ht/pbV5oqpdxGC3o1iwgJ4q83duSdO1M4fDKH61/7jvf+u4uiogA8QVi7mXUb3hEzIPekNbGGUspttKB7yJVt6/HFg5fTt0Udnp6/iXv+tYa8gqKKP+iPGnWzJtPYu8ruJEr5FS3oHhQXHcrbd6bwp+vb8dWmQzw0ax0FhQFY1EMioX5HLehKuZleKephIsLYPk0xxvDM55sJD3Hy/M2dcDgCbGaWhB7wv/ehMB+cAXi3SqWqgbbQbTKubzMevKols9ekM+mzTZWaP9CnJXSH/Gw4tNHuJEr5DW2h2+iBK1tyOreAt77dRWSok0evaWN3JM9J6GEttyywLkBSSlWZFnQbiQhPDG5LVm4hbyzZQWRoEL/t18LuWJ5RMwHaDYPvJltDGut3tDuRUj5Pu1xsJiI8M6wDQ5Mb8PwXW5m+YrfdkTzn+snWnRgXPGp3EqX8ghZ0L+B0CH+/JYmB7erx1Kc/8XHqXrsjeUZEbeg5wZqnVG8JoFSVaUH3EsFOB6+P6kzflnX4w5wf+WLjAbsjeUaH4dZyw2x7cyjlB7Sge5HQICf/vKMryQk1eWDmOtbvPWF3pOpXMwEa94Y102FZ8aQYSqlK0YLuZSJCgnjr1ynUrRHKuBmp7D9xxu5I1a/7XZB9BJY8A9Ouh6wMuxMp5ZO0oHuh2KhQ3rmzGzl5hYybnsrpXD+fcLn9jfDkIbhjHhxNgy8ftzuRUj5JC7qXalUvmtdGdWbLwZM8OGtdYNzMq3l/6yTphtlwUC84UupSaUH3Yv1a1+XP17fj602HeO7LLXbH8Yw+D0JoNKx80+4kSvkcvbDIy915WRPSMrL457KdNI+L4taUBLsjVa/wWtD4MtiXancSpXxOhS10EUkQkSUisllEfhKRB8rYRkTkVRFJE5EfRaRL9cQNPCLCUze0p2/LOvxx7gZW7jxqd6TqF58MGVshN8vuJEr5FFe6XAqAh40xbYGewL0i0q7UNtcCLYsf4wH9e9mNrDHqXUisHcGY91bz6MfrWbXzqP/e0KtBZ8DAwQ12J1HKp1RY0I0xB4wxa4ufnwI2Aw1LbTYUmGEsK4GaIhLv9rQBLCY8mBlje3B9p3gWbDjAiKkrueKFpbzyzXb2Hsu2O557NUi2lvv/Z28OpXzMJZ0UFZEmQGeg9MwEDYGS16unc2HRR0TGi0iqiKRmZOhY40vVsGY4L9ySxOonr2LyiCQSaofz8qJt9H1+CS9/s81/WuzR9SG6AWxdALkBOgerUpXgckEXkShgDvCgMeZk6dVlfOSC6mKMmWqMSTHGpMTFxV1aUnVOREgQN3ZuxAfjevLt7/szLLkBL3+znWe/2OI/Rb3bb2D3t/Bmbzi6w+40SvkElwq6iARjFfMPjDGflLFJOlBy+EUjYH/V46mKNKoVwUu3JnN7z0T+uWwnT8/3k8kyLn8UfvMl5GXBR7+2O41SPsGVUS4CvANsNsa8dJHNPgV+XTzapSeQaYwJkLtL2c/hEP5vaAfG9mnKtBW7eWLuRv+4ECmxJ3S/Gw79BHl+dp5AqWrgyjj03sAdwAYRWVf83hNAIoAxZgqwABgMpAHZwBj3R1XlERGevK4t4cFOXl+SRm5+Ic/e3ImQIB+/dqxuG8DAka06s5FSFaiwoBtjvqPsPvKS2xjgXneFUpUjIjxyTWvCgh38/attpB8/wxu/6kJcdKjd0SqvbvEI2cNbtKArVQEfb76pskwc0JJXRibz474TDHn9O9++DW+tpuAMgcOb7E6ilNfTgu6nhiY3ZM6Ey3CIcMs/v2f2mnS7I1WOMwjqtLbGpOu90pUqlxZ0P9a+QQzz7+tD18RaPPLxer7ZdMjuSJVTv6M1hPHvLa1JMJRSZdKC7udqR4YwY2x3msdF8tcFm8krKLI70qUbOAmGvwcNU2D123anUcpraUEPAMFOB09e145dR04z4/vddse5dFFx0OEmaNYPsg5BQZ7diZTySlrQA0S/1nH0bVmHVxdt5/hpHy2IMQ0BA6f0EgelyqIFPUCICH+6vh1ZuQW8/M02u+NUTkwja5npoyd4lapmWtADSKt60Yzqkci/Vv1M2mEfvOlVjeKCfnKfvTmU8lJa0APMQ1e1IiLEyR/nbuREto91vcQU38Azc2/52ykVoLSgB5jYqFCevK4tqXuOM+DFZcxa/bPv3PclJNKaoi5TW+hKlUULegAa0S2Rz+7rQ/O4SP4wZwM3T1lB6m4fuWgnphGkvmMNXywssDuNUl5FC3qAahtfg4/u7sVLtyaRfvwMw6d8z5j3fmDjvky7o5Wv3TCrL/3zh+GLP9idRimvInbdOzslJcWkpurM7t4gO6+A6Sv2MGXZDjLP5POXG9oxundTu2NdnDHw+e9gzXSYuBpim9udSCmPEZE1xpiUstZpC10RERLEhH7NWf77/lzdrh5/mb+Jf63cY3esixOBfo9by7Uz7E6jlNfQgq7OiQkP5vVRXRjQpi5PztvIR6u9eDRJVF1I7AVp39idRCmvoQVdnSckyME/ftWFy1vF8YdPfiTtcJbdkS6u5UA4tFFHvShVTAu6ukBYsJOXbk3CKcJHqV7cSm/cx1oeWFf+dkoFCC3oqkx1okK5sm1dPlmbTn6hl96hsUYDa3nqoL05lPISWtDVRd2aksCRrDwWbzlsd5SyRcYBYt2BUSmlBV1d3BWt4qgbHcrH3trt4gyCyDraQleqmBZ0dVFBTgc3d23Ekq0ZHD6ZY3ecskXV1xa6UsW0oKty3dK1EYVFhnf/u5sCb+xLj66nBV2pYkF2B1DerVlcFP1axzFl2Q5mrf6ZQR3qM7hjPL2axRLk9IL2QFQ9OLTJ7hRKeQUt6KpCU27vytKtGSzYcIBP1+3nwx/2UisimN8PasPIbgmIiH3hourB6cNQVAQOL/gfjFI20oKuKhQW7GRQh/oM6lCfnPxClm3L4L3/7uLxTzawZs9x/m9oB8JDnPaEi46HogLIOvjLMEalApQ2adQlCQt2ck37+nwwrif3X9mS2WvSuezZRTy7cAtHs3I9Hyixh7Xcsdjzx1bKy1RY0EXkXRE5LCIbL7K+n4hkisi64sef3R9TeRunQ/jdwFbMmdCLHk1jmbp8Bze9uYK9x7I9G6R+J4hJgC2fe/a4SnkhV1ro04BBFWzzrTEmufgxqeqxlK/o2rg2U+7oyuwJl3H8dB63/vN7dmZ48P4vItD2BusmXVkZnjuuUl6owoJujFkO+Mh0NsouXRJrMXN8L/IKirj97VXsO3HGcwfvOgYK82DtNM8dUykv5K4+9F4isl5EFopI+4ttJCLjRSRVRFIzMrQ15W/aNajBjLHdOZVbwO1vr2L7oVOeOXBcK0joAdu+9MzxlPJS7ijoa4HGxpgk4DVg3sU2NMZMNcakGGNS4uLi3HBo5W3aN4hh2phuHDudx7WvfMuzC7eQneeBuT8bdLHGoxcVVv+xlPJSVS7oxpiTxpis4ucLgGARqVPlZMpndW1cm8UPX8Gwzg2ZsmwHA19aztKt1XyDr/odIf80HNtVvcdRyotVuaCLSH0pvrJERLoX7/NoVferfFtsVCh/vyWJj+7uRUSIk9HvrWbS/E3kFlRTC7p+R2t58Mfq2b9SPsCVYYsfAt8DrUUkXUTGisg9InJP8SbDgY0ish54FRhp7Jp5Wnmd7k1rM/++Poy+rAnv/ncXw95YUT2jYOLaQHAk7P7W/ftWykeIXbU3JSXFpKam2nJsZY9Fmw/xyMfrqRMVyoIH+hLs7nvBfDwadn8HD28Fh01XripVzURkjTEmpax1eqWo8pgr29bjheFJbD+cxTvf7aKoyM2NibZD4HQGpGtDQQUmLejKo65qV4+r2tbl2YVbaPHHBTz/xRbc9ldiYk9rqXOMqgClN+dSHvfC8CRmrt7Lxv2Z/GPpDr7dfoQWdaN4emh7aoQFV37H0fHWtHQH1rsvrFI+RAu68rhakSFM6NccYwztG9Rg6dYMPvtxPz8fy2ZYcgMGd4wnNir00ncsAvFJWtBVwNIuF2UbEeG3/Vrw0d29mDwimQ3pmfzpPz9x53s/kJNfyeGNDTrD4c2Q66GrVJXyIlrQlVe4vlMD1j91NVNu78LGfSd54pMNletbb9IHTCHs+d79IZXyctrlorxGeIiTQR3ieeiqVkz+ZhtZuQXUjgwBoEvjWtyaklDxThJ6gDMEdi2DVldXc2KlvIsWdOV17hvQgoMnc1i8xZr8Oa+giJmr91I7IoR+rePKn8s0OFz70VXA0oKuvI7DIfy/mzoC1uX8OfmFXPvKt4ybkUqnRjHMmXBZ+RclxbaEnUs9klUpb6J96MrrhQU7+fddPZjYvwU/pmfy6MfreW3Rdt5YksbhUzkXfiC2OZzaD7kenGhDKS+gLXTlE+Jjwnn46lbsPZ7NvHX7z72/M+M0L96adP7Gsc2t5bGdEN/JgymVspe20JXPEBFeHpFM2l+vJe2v1zKmdxPmrdt34TymsS2s5dE0z4dUykZa0JVPERGCnA6CnA7GX94Mh8A/l+84f6MaDa1l1iHPB1TKRlrQlc+Kjwnn5i6N+NfKn2n6+Oc8OW+DtSIsxlqeOWFfOKVsoAVd+bRHrmnNQ1e1olezWD5Y9TM/H822bp0bGgM5WtBVYNGCrnxanahQHriqJU8PaY8xcPkLS/hpfyaEx2gLXQUcLejKL7SsF83kEdZol7lr90FYTW2hq4CjBV35jRs7N+KqtvX47McDmLCacOa43ZGU8igt6Mqv3JAUz8GTORwzkdrlogKOFnTlVwa2q0dYsIMvd+SQf/qY3XGU8igt6MqvRIQE8eR17ThJFORkgk2ToCtlBy3oyu/c3rMx9erVI9jkMXXxRv6xNI31e7X7Rfk/vZeL8kstmjSBozDjmzWkmzjaxh9g4QN97Y6lVLXSFrrySx3bdQBg0dhm3H9lS7YcPElmdr7NqZSqXlrQlX+qmQhA6Km9XNY8FmNgwItLueqlZWXfclcpP1BhQReRd0XksIhsvMh6EZFXRSRNRH4UkS7uj6nUJYppBAic+JmUxrW454rm9GweS9rhLL7fcdTudEpVC1da6NOAQeWsvxZoWfwYD7xZ9VhKVVFQKNRoACd+Jsjp4LFr2/DyiGTCgh2k7j5OZnb+uUelJqNWygtVeFLUGLNcRJqUs8lQYIaxfhUrRaSmiMQbYw64KaNSlROTACd+Pvcy2OmgY8MY3l+5h/dX7jn3/pjeTXjqhvZ2JFTKrdwxyqUhsLfE6/Ti9y4o6CIyHqsVT2JiohsOrVQ5ajSAA+vOe+uZYR35b9qRc68/St3LD7v0AiTlH9xR0KWM98r8G9YYMxWYCpCSkqJ/56rqFR0P276wLi4S6z/T1vWjaV0/+twm+0+c4f2VeygsMjgdZf2nrJTvcEdBTwcSSrxuBOy/yLZKeU50fcjPhtyTv0x6UUqretHkFhSR9PRX51omDw5sxdg+TT2XUyk3cUdB/xSYKCIzgR5ApvafK68QHW8tTx28aEG/pkN9dh09TW5+EQCfb9jP8m0ZWtCVT6qwoIvIh0A/oI6IpANPAcEAxpgpwAJgMJAGZANjqiusUpckur61PHUA4lqXuUlMeDB/GNTm3OtDJ3PYfOCkJ9Ip5XaujHK5rYL1BrjXbYmUcpdzLXTXJ4tOqB3BV5sOMrX0xNPAwHb1aVon0l3plHI7vZeL8l81igt6iaGLFemcWJOCIs/C/mEAAAyLSURBVMPfFmy5YN3mA6eYPCLZXemUcjst6Mp/hURCTCJkXFicL+aa9vXZPGkQhUXnD8K6451VHMzUWwYo76YFXfm3uNaQsfWSPhIW7LzgvfoxYWw9eMpdqZSqFlrQlX+r2wZ2LYeiQnBcWKhd3k10GIu3HOblb7ZdsC7Y6eC27onUjgypSlKlqkwLuvJvdVpDYS6c2AO1m1V6N50TazLj+928/M32MtdHhQZx52VNKr1/pdxBC7ryb2eL+LGdVSroQ5MbMiSpwQXvFxYZWvxxIcez8yq9b6XcRQu68m+xza3lsV1V3pXIhbcGCHIKNcKCOKGTZygvoBNcKP8WVQ+CI6wWejWpGRFC5hkt6Mp+2kJX/k0EajWFlf+A9jdBQje3H6JmRDDz1+9n8ZbDgHWS9M3bu9CtSW23H0up8mgLXfm/vr+zlntXVsvuHxrYitt7NubGzg0Z3DGeI1m5bNyXWS3HUqo82kJX/q/jcJj/IGSmV8vu+7euS//WdQHIKyjiwx9+JiunoFqOpVR5tIWuAkNMo2or6CWFBDkIDXKQlasFXXmeFnQVGGIaeqSggzUmXQu6soN2uajAUKMhpH0DhzZBvXbVeqiosCBW7z7Gc1/8cg+ZyBAn4/o2K/O2Akq5ixZ0FRia9IG102HFa3Djm9V6qM4JNVmw4SDvfGuNfS8yhoIiQ+fEWvRuUadaj60CmxZ0FRg63Qor34Qs1++NXlkvj+zMyyN/eb1xXybXv/Yd2XmF1X5sFdi0D10Fjqi6cPqwxw97tpvlTL4WdFW9tKCrwBEZB1kZHj9seIhV0HO0ha6qmRZ0FTii6sLpDCgq8uhhw4Ksn5m20FV104KuAkdUPTCFcOaYRw97toWuBV1VNz0pqgJHlHU1J589BCPe99hhw4Ksgv7KN9t5a3nZNwkLDXLw9p3daNeghsdyKf+jBV0Fjmb9reUlzDHqDg6H8Jcb2pGWkVXm+lM5Bfxn3X62HTqlBV1ViRZ0FTjCa0KXO2HrQo8fenTvphddt/dYNv9Zt5/8Qs/27Sv/o33oKrCE14KcE2CM3UnOCSk+aZpf6D2ZlG/Sgq4CS3hNKMyD/DN2Jzkn2Hm2oGsLXVWNSwVdRAaJyFYRSRORx8pYP1pEMkRkXfFjnPujKuUGYTWt5Znj9uYoIdhpTW2nBV1VVYV96CLiBN4ABgLpwGoR+dQYs6nUprOMMROrIaNS7hNey1rmnLDuwOgFzrbQ87SgqypypYXeHUgzxuw0xuQBM4Gh1RtLqWoSfraFfsLeHCWc63Ip0D50VTWujHJpCOwt8Tod6FHGdjeLyOXANuAhY8ze0huIyHhgPEBiYuKlp1Wqqrywy8XpEJwOYX36CWat/rmS+3Bwdft61AgLdnM65UtcKehSxnulmxLzgQ+NMbkicg8wHRhwwYeMmQpMBUhJSdHmiPK8sy30HO+a8zM+JozFWw6fm2i6Mv6S067c4ZHK/7lS0NOBhBKvGwH7S25gjDla4uVbwHNVj6ZUNQiJtpZ5p+3NUcpXD13Oiez8Sn02J7+QAS8uI6dA++ADnSsFfTXQUkSaAvuAkcCokhuISLwx5kDxyyHAZremVMpdQiKtZV7ZV23aJSIkiIiQyl3nl1N8j5jCIv2jN9BV+F+QMaZARCYCXwJO4F1jzE8iMglINcZ8CtwvIkOAAuAYMLoaMytVeUGhIE6va6FXhdNh9YoaL7pYStnDpSaBMWYBsKDUe38u8fxx4HH3RlOqGohASJTXtdCrwiFWQddRj0qvFFWBJyTSzwq6tSzUFnrA04KuAk9olF91uYgIDtEuF6UFXQWikEi/KuhgdbvoSVGlBV0FnhD/aqGDdc917XJRWtBV4AmJhNxTdqdwK6vLxe4Uym5a0FXg8cMuF6d2uSh0xiIViEIiIfsI7FwKzfrZHMY9HA5h26FTfJx6wS2UqrZfEfq3qUvtyBC37ldVDy3oKvDUTLTu5TJjGDyebo168XF1o0P5dvsRvt1+xO37vm9ACx6+urXb96vcTwu6Cjx9H7FuL7fkGWvmIj8o6PPv68PRrDy373fg5GXnbi2gvJ8WdBV4RCAqznpe6P4iaIeIkCAiarv/5xzkcOjJVh+iJ0VVYHKGWsvCXHtzeDkB9Fyr79CCrgJTUPFJvgL/aKFXFxEo0ia6z9CCrgKTttBd4nCUNb+N8lZa0FVgcha30AsrN6lEoLC6XLSF7iu0oKvAdK7LRVvo5XGIaEH3IVrQVWDSLheXiIieFPUhWtBVYNKToi4RvUeMT9GCrgLTuRa6FvTy6H3WfYsWdBWYzp0U1YJeHu1D9y1a0FVg0pOiLhG0y8WXaEFXgUlPirpET4r6Fi3oKjAFFRd0PSlaLodD+9B9iRZ0FZicwdZS+9DLJWgfui/Rgq4Ck3a5uMQh1p2GlW/Qgq4C09kWuna5lMuhfeg+xaWCLiKDRGSriKSJyGNlrA8VkVnF61eJSBN3B1XKrUSsVrq20Mund1v0KRUWdBFxAm8A1wLtgNtEpF2pzcYCx40xLYDJwHPuDqqU2zlD9OZcFXCI9rn4ElemOOkOpBljdgKIyExgKLCpxDZDgb8UP58NvC4iYvT0uPJmQSGw7t+wY7H1ukkfuO5FezN5GYfAsm0ZDHxpmd1R/MqIbgmM69vM7ft1paA3BEpOJZ4O9LjYNsaYAhHJBGKB82asFZHxwHiAxMTESkZWyk36Pgx7V/3yOqaRfVm81G96N2X59gy7Y/idOlGh1bJfVwp6WXe4L93ydmUbjDFTgakAKSkp2npX9up1r/VQFzWyeyIju2vjy1e4clI0HUgo8boRsP9i24hIEBADHHNHQKWUUq5xpaCvBlqKSFMRCQFGAp+W2uZT4M7i58OBxdp/rpRSnlVhl0txn/hE4EvACbxrjPlJRCYBqcaYT4F3gPdFJA2rZT6yOkMrpZS6kCt96BhjFgALSr335xLPc4Bb3BtNKaXUpdArRZVSyk9oQVdKKT+hBV0ppfyEFnSllPITYtfoQhHJAPZU8uN1KHUVqhfy9ozeng+8P6Pmqzpvz+iN+RobY+LKWmFbQa8KEUk1xqTYnaM83p7R2/OB92fUfFXn7Rm9PV9p2uWilFJ+Qgu6Ukr5CV8t6FPtDuACb8/o7fnA+zNqvqrz9ozenu88PtmHrpRS6kK+2kJXSilVihZ0pZTyE15d0L19cmoX8o0WkQwRWVf8GOfhfO+KyGER2XiR9SIirxbn/1FEungyn4sZ+4lIZonv8M9lbVeN+RJEZImIbBaRn0TkgTK2se17dDGf3d9hmIj8ICLrizM+XcY2tv2WXcxn62/ZZcYYr3xg3ap3B9AMCAHWA+1KbfNbYErx85HALC/LNxp43cbv8HKgC7DxIusHAwuxZpzqCazywoz9gM9s/A7jgS7Fz6OBbWX8e7bte3Qxn93foQBRxc+DgVVAz1Lb2PlbdiWfrb9lVx/e3EI/Nzm1MSYPODs5dUlDgenFz2cDV4pIWdPh2ZXPVsaY5ZQ/c9RQYIaxrARqiki8Z9JZXMhoK2PMAWPM2uLnp4DNWHPolmTb9+hiPlsVfy9ZxS+Dix+lR2PY9lt2MZ9P8OaCXtbk1KX/Qz1vcmrg7OTUnuBKPoCbi/8Mny0iCWWst5Or/wx261X85/BCEWlvV4jiboDOWC24krzieywnH9j8HYqIU0TWAYeBr40xF/0Obfgtu5IPvPu3DHh3QXfb5NTVxJVjzweaGGM6Ad/wSwvEW9j5/blqLda9K5KA14B5doQQkShgDvCgMeZk6dVlfMSj32MF+Wz/Do0xhcaYZKw5ibuLSIdSm9j6HbqQz9t/y4B3F3Rvn5y6wnzGmKPGmNzil28BXT2UzVWufMe2MsacPPvnsLFmzgoWkTqezCAiwVjF8gNjzCdlbGLr91hRPm/4DktkOQEsBQaVWuUVE81fLJ8P/JYB7y7o3j45dYX5SvWjDsHq3/QmnwK/Lh6l0RPINMYcsDtUSSJS/2xfqoh0x/pv9qgHjy9Yc+ZuNsa8dJHNbPseXcnnBd9hnIjULH4eDlwFbCm1mW2/ZVfy+cBvGXBxTlE7GC+fnNrFfPeLyBCgoDjfaE/lAxCRD7FGONQRkXTgKawTPhhjpmDNEzsYSAOygTGezOdixuHABBEpAM4AIz34P22A3sAdwIbiPlaAJ4DEEhnt/B5dyWf3dxgPTBcRJ9b/TD4yxnzmLb9lF/PZ+lt2lV76r5RSfsKbu1yUUkpdAi3oSinlJ7SgK6WUn9CCrpRSfkILulJK+Qkt6Eop5Se0oCullJ/4/wbZpoYkTLZBAAAAAElFTkSuQmCC\n", "text/plain": ["<Figure size 432x288 with 1 Axes>"]}, "metadata": {"needs_background": "light"}, "output_type": "display_data"}], "source": ["def draw_heavy_tail (sample) :\n", "    table   = DataFrame([[_] for _ in sample ], columns=[\"obs\"])\n", "    std = 1\n", "\n", "    normal = norm.rvs(size = len(sample))\n", "    normal = [ x*std  for x in normal ]\n", "    nortbl = DataFrame([ [_] for _ in normal ], columns=[\"obs\"])\n", "    nortbl[\"iobs\"] = (nortbl['obs'] * 10).astype(numpy.int64)\n", "\n", "    histon = nortbl[[\"iobs\", \"obs\"]].groupby('iobs', as_index=False).count()\n", "    histon.columns = ['iobs', 'nb']\n", "    histon = histon.sort_values(\"nb\", ascending=False).reset_index(drop=True)\n", "    \n", "    table[\"one\"] = 1\n", "    histo = table.groupby('obs', as_index=False).count()\n", "    histo.columns = ['obs', 'nb']                            \n", "    histo = histo.sort_values('nb', ascending=False).reset_index(drop=True)\n", "    histo.reset_index(drop=True, inplace=True)\n", "    histo[\"index\"] = histo.index + 1\n", "    \n", "    vec  = list(histon[\"nb\"])\n", "    vec += [0,] * len(histo)\n", "    histo['nb_normal'] = vec[:len(histo)]\n", "    \n", "    histo[\"log(index)\"] = numpy.log(histo[\"index\"]) / numpy.log(10)\n", "    histo[\"log(nb)\"] = numpy.log(histo[\"nb\"]) / numpy.log(10)\n", "    histo[\"log(nb_normal)\"] = numpy.log(histo[\"nb_normal\"]) / numpy.log(10)\n", "    return graph_XY ([\n", "            [histo[\"log(index)\"], histo[\"log(nb)\"], \"Zipf\"],\n", "            [histo[\"log(index)\"], histo[\"log(nb_normal)\"], \"Gauss\"], ],\n", "                     marker=False, link_point=True)  \n", "\n", "draw_heavy_tail(sample);"]}, {"cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": []}], "metadata": {"kernelspec": {"display_name": "Python 3", "language": "python", "name": "python3"}, "language_info": {"codemirror_mode": {"name": "ipython", "version": 3}, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.2"}}, "nbformat": 4, "nbformat_minor": 2}