wind.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import datetime as dt\n",
    "import matplotlib.pyplot as plt\n",
    "from statsmodels.distributions.empirical_distribution import ECDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load wind data\n",
    "wind_df = pd.read_csv(\"German_Wind_Power/50Hertz.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Date</th>\n",
       "      <th>00:00:00</th>\n",
       "      <th>00:15:00</th>\n",
       "      <th>00:30:00</th>\n",
       "      <th>00:45:00</th>\n",
       "      <th>01:00:00</th>\n",
       "      <th>01:15:00</th>\n",
       "      <th>01:30:00</th>\n",
       "      <th>01:45:00</th>\n",
       "      <th>02:00:00</th>\n",
       "      <th>...</th>\n",
       "      <th>22:00:00</th>\n",
       "      <th>22:15:00</th>\n",
       "      <th>22:30:00</th>\n",
       "      <th>22:45:00</th>\n",
       "      <th>23:00:00</th>\n",
       "      <th>23:15:00</th>\n",
       "      <th>23:30:00</th>\n",
       "      <th>23:45:00</th>\n",
       "      <th>year</th>\n",
       "      <th>month</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2019-08-23</td>\n",
       "      <td>74</td>\n",
       "      <td>73</td>\n",
       "      <td>63</td>\n",
       "      <td>55</td>\n",
       "      <td>54</td>\n",
       "      <td>53</td>\n",
       "      <td>51</td>\n",
       "      <td>52</td>\n",
       "      <td>45</td>\n",
       "      <td>...</td>\n",
       "      <td>63</td>\n",
       "      <td>61</td>\n",
       "      <td>57</td>\n",
       "      <td>52</td>\n",
       "      <td>51</td>\n",
       "      <td>49</td>\n",
       "      <td>45</td>\n",
       "      <td>40</td>\n",
       "      <td>2019</td>\n",
       "      <td>8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2019-08-24</td>\n",
       "      <td>37</td>\n",
       "      <td>38</td>\n",
       "      <td>35</td>\n",
       "      <td>35</td>\n",
       "      <td>31</td>\n",
       "      <td>32</td>\n",
       "      <td>30</td>\n",
       "      <td>29</td>\n",
       "      <td>28</td>\n",
       "      <td>...</td>\n",
       "      <td>112</td>\n",
       "      <td>107</td>\n",
       "      <td>104</td>\n",
       "      <td>104</td>\n",
       "      <td>99</td>\n",
       "      <td>95</td>\n",
       "      <td>95</td>\n",
       "      <td>96</td>\n",
       "      <td>2019</td>\n",
       "      <td>8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2019-08-25</td>\n",
       "      <td>93</td>\n",
       "      <td>88</td>\n",
       "      <td>85</td>\n",
       "      <td>81</td>\n",
       "      <td>76</td>\n",
       "      <td>80</td>\n",
       "      <td>81</td>\n",
       "      <td>76</td>\n",
       "      <td>73</td>\n",
       "      <td>...</td>\n",
       "      <td>103</td>\n",
       "      <td>96</td>\n",
       "      <td>93</td>\n",
       "      <td>88</td>\n",
       "      <td>86</td>\n",
       "      <td>82</td>\n",
       "      <td>80</td>\n",
       "      <td>77</td>\n",
       "      <td>2019</td>\n",
       "      <td>8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2019-08-26</td>\n",
       "      <td>72</td>\n",
       "      <td>66</td>\n",
       "      <td>59</td>\n",
       "      <td>54</td>\n",
       "      <td>48</td>\n",
       "      <td>46</td>\n",
       "      <td>45</td>\n",
       "      <td>43</td>\n",
       "      <td>42</td>\n",
       "      <td>...</td>\n",
       "      <td>47</td>\n",
       "      <td>46</td>\n",
       "      <td>47</td>\n",
       "      <td>49</td>\n",
       "      <td>50</td>\n",
       "      <td>51</td>\n",
       "      <td>48</td>\n",
       "      <td>42</td>\n",
       "      <td>2019</td>\n",
       "      <td>8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2019-08-27</td>\n",
       "      <td>39</td>\n",
       "      <td>38</td>\n",
       "      <td>36</td>\n",
       "      <td>36</td>\n",
       "      <td>34</td>\n",
       "      <td>31</td>\n",
       "      <td>31</td>\n",
       "      <td>30</td>\n",
       "      <td>31</td>\n",
       "      <td>...</td>\n",
       "      <td>51</td>\n",
       "      <td>52</td>\n",
       "      <td>51</td>\n",
       "      <td>52</td>\n",
       "      <td>54</td>\n",
       "      <td>54</td>\n",
       "      <td>56</td>\n",
       "      <td>59</td>\n",
       "      <td>2019</td>\n",
       "      <td>8</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 99 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "        Date  00:00:00  00:15:00  00:30:00  00:45:00  01:00:00  01:15:00  \\\n",
       "0 2019-08-23        74        73        63        55        54        53   \n",
       "1 2019-08-24        37        38        35        35        31        32   \n",
       "2 2019-08-25        93        88        85        81        76        80   \n",
       "3 2019-08-26        72        66        59        54        48        46   \n",
       "4 2019-08-27        39        38        36        36        34        31   \n",
       "\n",
       "   01:30:00  01:45:00  02:00:00  ...  22:00:00  22:15:00  22:30:00  22:45:00  \\\n",
       "0        51        52        45  ...        63        61        57        52   \n",
       "1        30        29        28  ...       112       107       104       104   \n",
       "2        81        76        73  ...       103        96        93        88   \n",
       "3        45        43        42  ...        47        46        47        49   \n",
       "4        31        30        31  ...        51        52        51        52   \n",
       "\n",
       "   23:00:00  23:15:00  23:30:00  23:45:00  year  month  \n",
       "0        51        49        45        40  2019      8  \n",
       "1        99        95        95        96  2019      8  \n",
       "2        86        82        80        77  2019      8  \n",
       "3        50        51        48        42  2019      8  \n",
       "4        54        54        56        59  2019      8  \n",
       "\n",
       "[5 rows x 99 columns]"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "wind_df.Date = pd.to_datetime(wind_df.Date)\n",
    "wind_df[\"year\"] = pd.DatetimeIndex(wind_df[\"Date\"]).year\n",
    "wind_df[\"month\"] = pd.DatetimeIndex(wind_df[\"Date\"]).month\n",
    "wind_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Subset to seasons\n",
    "spring = wind_df.loc[wind_df[\"month\"].isin([3,4,5])]\n",
    "summer = wind_df.loc[wind_df[\"month\"].isin([6,7,8])]\n",
    "autumn = wind_df.loc[wind_df[\"month\"].isin([9,10,11])]\n",
    "winter = wind_df.loc[wind_df[\"month\"].isin([1,2,12])]\n",
    "\n",
    "seasons = [spring, summer, autumn, winter]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get ECDF and samples for each hour in each season\n",
    "n_sample = 5001\n",
    "quantile = np.linspace(0.01, 1, n_sample)\n",
    "samples = np.zeros((4, 24, n_sample-1))\n",
    "\n",
    "\n",
    "for s in range(4):\n",
    "    df = seasons[s]\n",
    "    for h in range(24):\n",
    "        i = np.arange(4*h+1, 4*h+5)\n",
    "        test = df.iloc[:, i]\n",
    "        test = test.to_numpy().flatten()\n",
    "        # Estimate ECDF\n",
    "        ecdf = ECDF(test)\n",
    "        # Inverse sampling\n",
    "        sample = np.interp(quantile, ecdf.y, ecdf.x)\n",
    "        samples[s, h, :] = sample[1:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [],
   "source": [
    "name = [\"spring\", \"summer\", \"autumn\", \"winter\"]\n",
    "# Save sample data\n",
    "for s in range(4):\n",
    "    data = samples[s, :, :]\n",
    "    filename = name[s] + \"_wind.csv\"\n",
    "    np.savetxt(filename, data, delimiter=\",\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "166.10560000000012\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.LineCollection at 0x20fbacdaa90>"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAQx0lEQVR4nO3dX4wdZ3nH8e+vzh9QCErc2K6xrdpILpWDShKtXFAqRDEihiCcm1RGonIrV75JW1ArUbtIrbiwFHqB4KJBsgLUFQTXCtBYQQJc0whVamPWJIE4jhtD3HhlEy9QxJ8L04SnF2fSnti73uPdc9g9b74faTUz73nn7PPI9m/HszNzUlVIktrya4tdgCRp+Ax3SWqQ4S5JDTLcJalBhrskNeiqxS4A4Kabbqr169cvdhmSNFaOHTv2g6paMdNrSyLc169fz+Tk5GKXIUljJcl/zfaap2UkqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBS+IO1cW0fveX573v6XvvHGIlkjQ8HrlLUoMGCvckNyR5MMnTSU4keUuS5UkOJ3mmW97YN39PklNJTia5Y3TlS5JmMuiR+yeAr1TVbwNvAk4Au4EjVbURONJtk2QTsB24GdgK3Jdk2bALlyTNbs5wT/Ja4K3ApwCq6hdV9WNgG7C/m7YfuKtb3wYcqKoLVfUscArYPNyyJUmXM8iR++uBaeAzSR5Lcn+S64BVVXUOoFuu7OavAc707T/Vjb1Mkl1JJpNMTk9PL6gJSdLLDRLuVwG3AZ+sqluBn9OdgplFZhirSwaq9lXVRFVNrFgx47PmJUnzNEi4TwFTVfVot/0gvbB/PslqgG55vm/+ur791wJnh1OuJGkQc4Z7VX0fOJPkDd3QFuAp4BCwoxvbATzUrR8Ctie5NskGYCNwdKhVS5Iua9CbmP4M+FySa4DvAX9M7wfDwSQ7geeAuwGq6niSg/R+ALwA3FNVLw698iXAG6AkLVUDhXtVPQ5MzPDSllnm7wX2zr8sSdJCeIeqJDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGDfoZqkvaQj7LVJJa5JG7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJatBA4Z7kdJLvJHk8yWQ3tjzJ4STPdMsb++bvSXIqyckkd4yqeEnSzK7kyP33q+qWqprotncDR6pqI3Ck2ybJJmA7cDOwFbgvybIh1ixJmsNCTstsA/Z36/uBu/rGD1TVhap6FjgFbF7A95EkXaFBw72AryU5lmRXN7aqqs4BdMuV3fga4EzfvlPd2Msk2ZVkMsnk9PT0/KqXJM1o0McP3F5VZ5OsBA4nefoyczPDWF0yULUP2AcwMTFxyeuSpPkb6Mi9qs52y/PAl+idZnk+yWqAbnm+mz4FrOvbfS1wdlgFS5LmNme4J7kuyfUvrQPvBJ4EDgE7umk7gIe69UPA9iTXJtkAbASODrtwSdLsBjktswr4UpKX5j9QVV9J8k3gYJKdwHPA3QBVdTzJQeAp4AXgnqp6cSTVS5JmNGe4V9X3gDfNMP5DYMss++wF9i64OknSvHiHqiQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGDfrIXwmA9bu/PO99T9975xArkXQ5HrlLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1aOBwT7IsyWNJHu62lyc5nOSZbnlj39w9SU4lOZnkjlEULkma3ZUcuX8AONG3vRs4UlUbgSPdNkk2AduBm4GtwH1Jlg2nXEnSIAYK9yRrgTuB+/uGtwH7u/X9wF194weq6kJVPQucAjYPpVpJ0kAGPXL/OPAh4Jd9Y6uq6hxAt1zZja8BzvTNm+rGXibJriSTSSanp6evtG5J0mXM+UlMSd4DnK+qY0neNsB7ZoaxumSgah+wD2BiYuKS11vnJxpJGqVBPmbvduC9Sd4NvAp4bZLPAs8nWV1V55KsBs5386eAdX37rwXODrNoSdLlzXlapqr2VNXaqlpP7xelX6+q9wOHgB3dtB3AQ936IWB7kmuTbAA2AkeHXrkkaVYL+YDse4GDSXYCzwF3A1TV8SQHgaeAF4B7qurFBVcqSRrYFYV7VT0CPNKt/xDYMsu8vcDeBdYmSZon71CVpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAbNGe5JXpXkaJInkhxP8pFufHmSw0me6ZY39u2zJ8mpJCeT3DHKBiRJlxrkyP0C8PaqehNwC7A1yZuB3cCRqtoIHOm2SbIJ2A7cDGwF7kuybAS1S5JmMWe4V8/Pus2ru68CtgH7u/H9wF3d+jbgQFVdqKpngVPA5mEWLUm6vIHOuSdZluRx4DxwuKoeBVZV1TmAbrmym74GONO3+1Q3dvF77koymWRyenp6AS1Iki521SCTqupF4JYkNwBfSvLGy0zPTG8xw3vuA/YBTExMXPK6Zrd+95cXtP/pe+8cUiWSlqorulqmqn4MPELvXPrzSVYDdMvz3bQpYF3fbmuBswstVJI0uEGullnRHbGT5NXAO4CngUPAjm7aDuChbv0QsD3JtUk2ABuBo0OuW5J0GYOcllkN7O+uePk14GBVPZzk34GDSXYCzwF3A1TV8SQHgaeAF4B7utM6kqRfkTnDvaq+Ddw6w/gPgS2z7LMX2Lvg6iRJ8+IdqpLUIMNdkho00KWQkn71FnLJq5e7yiN3SWqQR+6vQAu9CUrS0ueRuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktSgOT9DNck64B+B3wB+Ceyrqk8kWQ78E7AeOA38QVX9d7fPHmAn8CLw51X11ZFUr7GykM9uPX3vnUOsRGrfIB+Q/QLwl1X1rSTXA8eSHAb+CDhSVfcm2Q3sBv4qySZgO3Az8DrgX5L8VlW9OJoW9Eqw0A/19oeDXmnmPC1TVeeq6lvd+k+BE8AaYBuwv5u2H7irW98GHKiqC1X1LHAK2DzkuiVJlzHIkfv/SbIeuBV4FFhVVeeg9wMgycpu2hrgP/p2m+rGJGlkPO33cgP/QjXJa4AvAB+sqp9cbuoMYzXD++1KMplkcnp6etAyJEkDGCjck1xNL9g/V1Vf7IafT7K6e301cL4bnwLW9e2+Fjh78XtW1b6qmqiqiRUrVsy3fknSDOYM9yQBPgWcqKqP9b10CNjRre8AHuob357k2iQbgI3A0eGVLEmayyDn3G8H/hD4TpLHu7G/Bu4FDibZCTwH3A1QVceTHASeonelzT1eKSNJv1pzhntV/Rszn0cH2DLLPnuBvQuoS5K0AN6hKkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhp0Rc+WkTQefIqmDHdphBYastJ8eVpGkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIa5B2qekVYyJ2i3oqvceSRuyQ1yHCXpAYZ7pLUIM+5S3rFa/ERyR65S1KDDHdJapDhLkkNmjPck3w6yfkkT/aNLU9yOMkz3fLGvtf2JDmV5GSSO0ZVuCRpdoMcuf8DsPWisd3AkaraCBzptkmyCdgO3Nztc1+SZUOrVpI0kDnDvaq+AfzoouFtwP5ufT9wV9/4gaq6UFXPAqeAzcMpVZI0qPmec19VVecAuuXKbnwNcKZv3lQ3dokku5JMJpmcnp6eZxmarwMP7ObAA7sXuwxJIzLsX6hmhrGaaWJV7auqiaqaWLFixZDLkKRXtvmG+/NJVgN0y/Pd+BSwrm/eWuDs/MuTJM3HfO9QPQTsAO7tlg/1jT+Q5GPA64CNwNGFFiktpoXevSgthjnDPcnngbcBNyWZAv6WXqgfTLITeA64G6Cqjic5CDwFvADcU1Uvjqh2SUuQj1deGuYM96p63ywvbZll/l5g70KKkiQtjHeoSlKDDHdJapDhLkkNMtwlqUF+WIekJcPLTofHI3dJapDhLkkNMtwlqUGGuyQ1yHCXpAZ5tYykS3jVypVZis/T8chdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDVoZOGeZGuSk0lOJdk9qu8jSbrUSMI9yTLg74F3AZuA9yXZNIrvJUm61KiO3DcDp6rqe1X1C+AAsG1E30uSdJFRfRLTGuBM3/YU8Lv9E5LsAnZ1mz9LcnKW97oJ+MHQK1w8S6Kft7y08tH3DOPtlkRPQ2Q/S18zPeWjwPz7+c3ZXhhVuGeGsXrZRtU+YN+cb5RMVtXEsApbbK31A+31ZD9LX2s9jaKfUZ2WmQLW9W2vBc6O6HtJki4yqnD/JrAxyYYk1wDbgUMj+l6SpIuM5LRMVb2Q5E+BrwLLgE9X1fF5vt2cp27GTGv9QHs92c/S11pPQ+8nVTX3LEnSWPEOVUlqkOEuSQ1a0uE+jo8wSPLpJOeTPNk3tjzJ4STPdMsb+17b0/V3Mskdi1P17JKsS/KvSU4kOZ7kA934WPaU5FVJjiZ5ouvnI934WPbzkiTLkjyW5OFue9z7OZ3kO0keTzLZjY1tT0luSPJgkqe7f0tvGXk/VbUkv+j9Iva7wOuBa4AngE2LXdcAdb8VuA14sm/s74Dd3fpu4KPd+qaur2uBDV2/yxa7h4v6WQ3c1q1fD/xnV/dY9kTvHozXdOtXA48Cbx7Xfvr6+gvgAeDhcf8719V5GrjporGx7QnYD/xJt34NcMOo+1nKR+5j+QiDqvoG8KOLhrfR+8OlW97VN36gqi5U1bPAKXp9LxlVda6qvtWt/xQ4Qe8O5LHsqXp+1m1e3X0VY9oPQJK1wJ3A/X3DY9vPZYxlT0leS++g71MAVfWLqvoxI+5nKYf7TI8wWLNItSzUqqo6B72wBFZ242PVY5L1wK30jnbHtqfuFMbjwHngcFWNdT/Ax4EPAb/sGxvnfqD3A/drSY51jyqB8e3p9cA08Jnu1Nn9Sa5jxP0s5XCf8xEGDRibHpO8BvgC8MGq+snlps4wtqR6qqoXq+oWendOb07yxstMX9L9JHkPcL6qjg26ywxjS6afPrdX1W30nix7T5K3XmbuUu/pKnqnaj9ZVbcCP6d3GmY2Q+lnKYd7S48weD7JaoBueb4bH4sek1xNL9g/V1Vf7IbHuieA7r/GjwBbGd9+bgfem+Q0vVOXb0/yWca3HwCq6my3PA98id5piXHtaQqY6v6HCPAgvbAfaT9LOdxbeoTBIWBHt74DeKhvfHuSa5NsADYCRxehvlklCb1zhSeq6mN9L41lT0lWJLmhW3818A7gaca0n6raU1Vrq2o9vX8jX6+q9zOm/QAkuS7J9S+tA+8EnmRMe6qq7wNnkryhG9oCPMWo+1ns3yLP8Rvmd9O7OuO7wIcXu54Ba/48cA74H3o/gXcCvw4cAZ7plsv75n+46+8k8K7Frn+Gfn6P3n8Jvw083n29e1x7An4HeKzr50ngb7rxseznot7exv9fLTO2/dA7R/1E93X8pX/7Y97TLcBk9/fun4EbR92Pjx+QpAYt5dMykqR5MtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSg/4XhhXc3oMoj20AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sample = sample[sample > 0]\n",
    "md = np.median(sample)\n",
    "print(md)\n",
    "plt.hist(sample, bins=20)\n",
    "plt.vlines(x=md, ymin=0, ymax=100, colors='r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAhFUlEQVR4nO3deXxV9Z3/8dcnO2SDkBDWkKAghB0iS12qiPuCa+tSba3K2JHazrSPqdNOt+m009rpYuvCjypWayvWtlqqKO4iIpV9CUgIa0IgG5B9vfn+/sitjTHABW84d3k/H4887r3nnNy8v23y9su5ZzHnHCIiEv5ivA4gIiLBoUIXEYkQKnQRkQihQhcRiRAqdBGRCBHn1Q/OzMx0ubm5Xv14EZGwtHbt2irnXFZP6zwr9NzcXNasWePVjxcRCUtmtvdo67TLRUQkQqjQRUQihApdRCRCqNBFRCKECl1EJEIct9DNbJGZVZjZlqOsNzP7lZkVm9kmM5sa/JgiInI8gczQfwtccoz1lwKj/F/zgEc+eSwRETlRxz0O3Tm33Mxyj7HJXOBJ13kd3lVm1s/MBjvnDgQrpIhIOOvocJQcbqSovJ6i8jomDkvnnFE9nhv0iQTjxKKhQEmX16X+ZR8rdDObR+csnpycnCD8aBGR0FXT1Mbj7+7m8Xf3UNPU9uHyL513WsgWuvWwrMe7ZjjnFgILAQoKCnRnDREJOyWHGnnjgwrqW9qpb2mnsaWd+hYfDS3tNLR2LmtoaaehxUdVfQst7R1cmJ/NnLEDGZWdyqiBKaQmxfdKtmAUeikwvMvrYUBZEN5XRCRkbDtQy7LCg/zytR0fLouLMZIT40hOiO18TIwjJTGOgamJJCfGkd4nnuumDmP80PRTkjEYhb4EmG9mi4EZQI32n4tIuGtt7+DVreVsKj3C+pIjvL/7EGYwfmgaV00awm2zckmMi8Gsp50U3jhuoZvZ08B5QKaZlQLfBeIBnHMLgKXAZUAx0Ajc3lthRUR6Q1Orjw8O1lJ6uIkDNU0cqGnmj6tLaGj1ATBpWDr/Nmc0t80aQf/kBI/THl0gR7ncdJz1DrgnaIlERE4B5xyL3t3Ds2tKKCqvo6PLp3rJCbEM7teH4f378MjnppEUH+td0BPg2eVzRUS8UNPUxrLCgyx+fx/r9h1h8vB+zD//dMYNTSd3QDKD0pNIS4oLqV0pgVKhi0jUKDvSxNyH3qWyroWcjL788Jrx3HRmDjEx4VfePVGhi0jE23+kiVU7q/nFa0VU1rWw8NZpXJifHZaz8GNRoYtIRPvZK9v59RvFAGSlJvK9K/MjssxBhS4iEezn/jI3gxe/fA5jBqVGzO6VnqjQRSSidHQ4nl1bwvPry3hvVzX9+8az4huzSU6M/LqL/BGKSNTYfrCOB14vYunmgyTExfD1i0Zz26dyo6LMQYUuImHMOUdhWS0riqt4t7iKd3ZUER9rzB4zkPuvn0hmSqLXEU8pFbqIhJ3mNh9PrdrL71btZW91IwCjBqZw7wWjuP1TuSF9NmdvUqGLSNhobe9g24Fa/vuFrazde5gZeRncc97pnD9mIFmp0TUb74kKXURCQkNLO5V1LVTVt3zksbK+lf1HmthT1UDp4UY6HKQmxvHAjZOZO3mo17FDigpdRE4J5xwHapr5++5q1u87wsGa5s7Srm+hqq6Vpjbfx77HDAYkJzA4vQ+Thvdj7uQh5GUmc8GYbNL79s41xcOZCl1EgsY5R3ltC5v311BYVkPJoX9evfBATRPNbR1A5wx7SL8+ZKYmMC2nP5kpiWSmJpL1kccEMvomEBcbyK2PBVToIhIkLe0+7n16PcsKywGIMRiUlsSg9CTyh6QxZ+xAhvbrQ0FuBmMHpxEbwSf4eEWFLiInpc3XQcmhRvZWN1JYVsOSjWUUldcz//zTOX9MFvmD0+mTEB6XnY0UKnQROarW9g4ON7ZSVd/CoYZWKmpbWLfvMH/ffYjdVQ34ulxEfGpOPx68eQpXTBziYeLopkIXiVK+DseRxlaqG1qprm/lYG0TReX17CivY1dVA1V1LdQ2t3/s+1IS4yjI7c8l4waRm5lMXmZf8jJTyIjSY79DiQpdJErsP9LEg28Us3bvIarrWznc2PqRu/RA502P8zKTOSM7lXNOz2RASiIZyQlkpiR8+HxERl99UBmiVOgiUaCovI6bFq6irqWdc0dlUZCbQWZyAhnJCWSkJJKZnMDAtERyMpJJiFNZhysVukiEO1jTzE0LVxEbY7z0lXM4LSvF60jSS1ToIhGqpqmNVwoP8vNXi6huaOXpu2aqzCOcCl0kQs19cAV7qhsZntGHn1w3gZkjM7yOJL1MhS4SIZpafby5vYIVxVVsLq1hT3UjN8/I4YdXj4/I263Jx6nQRSLAW9sr+MriDdQ0tZGaFMfEYen8y6dHcuvMESrzKKJCFwlTLe0+nly5l79tKmNTaQ0js5J5+JapzBw5QKfVRykVukgY8nU47v7dWt7cXsnEYel8/6pxfPbM4STF61T7aKZCFwkjLe0+HluxmydW7qG8toUL87P5zW0FXseSEKFCFwkjf91Qxv0vb+ecUZl86/J8Lhyb7XUkCSEqdJEw4Zzjrxv2E2Pw0C1TSUvSDR7ko1ToImHge0sKeaXwIGU1zVw7ZajKXHqkQhcJUR0djh8u3carW8vZd6iRtKQ4fjB3HFdN0n00pWcBFbqZXQI8AMQCjzrnftxtfTrwFJDjf8//c849HuSsIlGl5HAjj63YzaTh/bht1ghuKBhOeh/NzOXojlvoZhYLPARcCJQCq81siXNua5fN7gG2OueuNLMsYLuZ/d4519orqUUimHOO17ZV8LNXtgPwzUvHMGPkAI9TSTgIZIY+HSh2zu0CMLPFwFyga6E7INU6T0lLAQ4BH78yvogc06pd1fzwxW1s3l/DiAF9efiWqSpzCVgghT4UKOnyuhSY0W2bB4ElQBmQCnzWOdfR/Y3MbB4wDyAnJ+dk8opEpOY2Hw+8voPfLN/FoPQk7r9+ItdOGaobScgJCeS3padziLvd54SLgQ3AEGAy8KCZpX3sm5xb6JwrcM4VZGVlnWBUkcjk63A8/NZOHnlrJ7NOG8CL957DZwqGq8zlhAUyQy8Fhnd5PYzOmXhXtwM/ds45oNjMdgNjgPeDklIkwtQ1t7GptIZlhQdZuvkgVfUtTByWzuNfOFNFLictkEJfDYwyszxgP3AjcHO3bfYBFwDvmFk2cAawK5hBRcJZdX0Lb26vZNWuatbtPcyuqgYAkuJjmD1mIJdPGMIFYweqzOUTOW6hO+fazWw+sIzOwxYXOecKzexu//oFwA+A35rZZjp30XzDOVfVi7lFwkJ5bTOPrdjNk+/tobmtg/594ynIzeCaKUMZPyyd6bkZJCfqdBAJjoB+k5xzS4Gl3ZYt6PK8DLgouNFEwlNFXTOPvbObNz6oYEdFPTEGcycP5Y6z88gfnEaMLm0rvURTA5Egcc5RWFbLFb9eAcC5o7O4duowLh6XzUjdy1NOARW6yElobvOx71Aju6sa2FvdwOb9tazaVU1lXQtm8O3L8/ni2Xlex5Qoo0IXOY7G1nYqalsor22mvK6FrWW1PP7ublra/3mqRXZaIrNGDuDs0zO5MD+b/skJHiaWaKVCF+lBYVkNT7+/j1cKy6moa/nY+um5GdwyM4fcAcnkZibrGisSElToIl3srKzn/pc/YFlh+YeHFI4fmk52ahLZaUlkpyUyMDWJtD5xuvmyhBwVukS15jYfhxpaqa5vpeRwI998bjO+DsdX54zi9rPyNPOWsKJCl6hyuKGVbz63ma0Haqmub6W+5aPXkBuSnsTiebPIGdDXo4QiJ0+FLlHj9W3l/NszG2ho9XHxuGyy05LITEkkIzmBAckJDEhJ4IxBaaToRB8JU/rNlYjX0eF4bv1+vvbsRgBevPdsxg1J9ziVSPCp0CXi7K1uYEPJEbbsr2FTaQ2FZbXUt7QzJacf379qnMpcIpYKXSJCeW0z2w/W8dSqvbyytRyAhLgYxg5O45opQ5k4LJ0rJw0hKT7W46QivUeFLmGnvqW98ySf2mYq61pY9O4eNpYc+XB9YlwMv7tjBlNy+hGvqxdKFFGhS9hobe/gzifXsLyo8mPrYgyeunMGo7NTyUxJ9CCdiPdU6BIWmlp9fO3ZDSwvquTqyUM474yBDPSf5JOdlkhKok70EVGhS8hravVxxxOreW9XNXMnD+FH106gb4J+dUW601+FhLTmNh93PbmG93ZV87MbJnHt1GFeRxIJWfrESELassKDrCiu4n+vmaAyFzkOzdAlJLX5Oli8uoSH3ywG4KJxgzxOJBL6VOgSkr701Fpe21bB1Jx+/OyGSWTo+uIix6VCl5DR7uvgrxvK+P3f97Ju3xFump7Dj64Zr6NXRAKkQhfP1TW3sWV/Ld96bjO7qhoYNTCFr104mttm5arMRU6ACl08taO8jst/vYJW/+3cpudl8My8mSpykZOgQhfPNLa28/VnN9La3sF3r8znqklDGKCzPEVOmgpdTgnnHB8crOPd4iq2H6yjqKKe7QdraW7r4MpJQ7j9rDyvI4qEPRW69LoDNU3c9eQatuyvBSAzJZHR2SncPH0Ec/IHMiNvgMcJRSKDCl16VXObj6sfepeGFh8/mDuOi8cNYmBaktexRCKSCl16RXV9C4+8tZM/vL+PxlYf37kin1tn5XodSySiqdAl6EoPN3LjwlWUHWnisgmDuWrSEOaMzfY6lkjEU6FLUNU2t3HTb1ZR09TGX/71LCYP7+d1JJGooUKXoNpcWkPJoSYWfG6qylzkFAvoaotmdomZbTezYjO77yjbnGdmG8ys0MzeDm5MCReb99cA6INPEQ8cd4ZuZrHAQ8CFQCmw2syWOOe2dtmmH/AwcIlzbp+ZDeylvBLC1u87zI9f+oBzR2cxcWi613FEok4gM/TpQLFzbpdzrhVYDMztts3NwF+cc/sAnHMVwY0p4eCXr+0gLsb48bUTiNPNmUVOuUD+6oYCJV1el/qXdTUa6G9mb5nZWjO7LVgBJTxsLDnC20WV3FAwnCH9+ngdRyQqBfKhaE9XSXI9vM804AKgD/Cema1yzhV95I3M5gHzAHJyck48rYSkplYf311SSIzB/Nmnex1HJGoFMkMvBYZ3eT0MKOthm5edcw3OuSpgOTCp+xs55xY65wqccwVZWVknm1lCzEtbDrCh5Ag/vX4SQzU7F/FMIIW+GhhlZnlmlgDcCCzpts1fgXPMLM7M+gIzgG3BjSqhqvRwEwBnj8r0OIlIdDvuLhfnXLuZzQeWAbHAIudcoZnd7V+/wDm3zcxeBjYBHcCjzrktvRlcvNfS7uPbz2/h2bWljBmUygDdJk7EUwGdWOScWwos7bZsQbfXPwV+GrxoEuoefWc3f1xTyh1n53Hv7FE6skXEYzpTVE5Km6+Dh94s5rSsZP7r8rG6w5BICNCUSk7Y1rJarn14JY2tPr5+0Rkqc5EQoRm6nJAt+2u45dG/U9PUxi8+O4lLJwz2OpKI+KnQJSC+DseCt3fy0JvFJCfG8fPPTOKaKcO8jiUiXajQ5bg2lBzhF68W8XZRJWMGpfLEF6eTrYtviYQcFboc1Zb9Nfx02XbeLqokKT6GKyYO5n+uHk+/vjo8USQUqdDlY7aW1fL8hv38duUekhNi+cYlY7h11ghSEvXrIhLK9BcqQOfNnJcVHuQv6/bzdlElcTHG+WMG8pPrJpKhE4ZEwoIKPco1trbz5geVPPxWMYVltQxKS+LrF43mczNHaNeKSJhRoUepTaVHWLh8F69vq6CpzUdmSiIP3DiZKycOISZGx5WLhCMVehTasr+Gqx58l/Q+8Vw7dSiXTxzMjLwBxKrIRcKaCj3KHGpo5YHXdwDw2OcLKMjN8DiRiASLCj3KfO2PG1i+o4q7P30ak4b38zqOiASRCj1KLC+qZOnmA7y3q5o5Ywdy36VjvI4kIkGmQo8CbxdV8vlF75OSGMfZp2dx5zkjvY4kIr1AhR4FfvbKdhJiY1j1zQt0cpBIBNPlcyPcXzfsZ1NpDXedm6cyF4lwKvQItryokn97ZgMzR2Zwz/mnex1HRHqZCj1COef4ycsfEB8bw6IvnEnfBM3ORSKdCj1CvVVUSWFZLfPOHakyF4kS+kuPMM451pcc4X+XbqNPfCx3nasjWkSihQo9gryzo5IfvLCVovJ6EmJj+J9rxpOWFO91LBE5RVToEaChpZ2fvVLEE+/tISejL/dfP5GLxw0ivY/KXCSaqNDD3Hs7q/n+3wopKq/j+mnD+PYV+aRqVi4SlVToYaqjw1FUUce8J9dQ39rOo7cVcMHYbK9jiYiHVOhhZO3ew7yzo5J1+46wft9h6prbiY0xHvt8AbPHqMxFop0KPUzsqWrgukdWAjBmUCpXTBzCtBH9mZGXwfCMvh6nE5FQoEIPA/Ut7XzlmQ0A/PlLs5g2QtcwF5GPU6GHuBc2lfFfz2/hSGMbN8/IUZmLyFGp0EPYn9eW8rVnN5IQG8Pz95zFZN2QQkSOQYUeYtp8Hby3s5qn39/HS1sOMm1Ef7535TgmDEv3OpqIhLiACt3MLgEeAGKBR51zPz7KdmcCq4DPOuf+FLSUUeSmhatYs/cwqYlx3HvBKOaffzoJcbrkjogc33EL3cxigYeAC4FSYLWZLXHObe1hu58Ay3ojaKSraWrjpc0H2Ly/hovHZfPAjVNIio/1OpaIhJFAZujTgWLn3C4AM1sMzAW2dtvuy8CfgTODmjAKvLjpAPcuXo+vwzFuSBr3nH+6ylxETlgghT4UKOnyuhSY0XUDMxsKXAPM5hiFbmbzgHkAOTk5J5o1IlXUNvOjpdvwdTie+9dPMXl4P8zM61giEoYC2TnbU7u4bq9/CXzDOec71hs55xY65wqccwVZWVkBRoxsP122nf1HmvjT3bOYktNfZS4iJy2QGXopMLzL62FAWbdtCoDF/jLKBC4zs3bn3PPBCBmp3tpewbNrS5kzdiAFuTq+XEQ+mUAKfTUwyszygP3AjcDNXTdwzuX947mZ/RZ4QWV+bJV1Ldz+29VkpSby1TmjvY4jIhHguIXunGs3s/l0Hr0SCyxyzhWa2d3+9Qt6OWNEam7z4Rx845IxjB+qY8xF5JML6Dh059xSYGm3ZT0WuXPuC588VuRbt+8wAPGx2mcuIsGhM1Y88ErhQb72x41MHt6PObqGuYgEiQr9FHttazn3/GEd44em8+Qd00lO1NUXRCQ41CanSHObj1+8VsSiFbvJH9JZ5rqBs4gEk2bop8ji9/fx/97excyRA3jyiypzEQk+zdBPgZZ2Hwve3sXQfn147PNn6mJbItIr1CynwB/XlHKwtpl7dOVEEelFapde1u7r4KE3ihmcnsT104Z5HUdEIpgKvZf9Y3Y+f7Zm5yLSu9Qwvaijw/Gr13eQnZbIZwqGH/8bREQ+ARV6L1qysYyDtc18efYo4mP1P7WI9C61TC96u6gSgGunDvU4iYhEAxV6L3ppywHmjB1I3wQdHSoivU9N0wte3nKAVwrLaW7rYHR2qtdxRCRKqNCDyDnHM6tLuO8vmwG4bMIgrpmi3S0icmqo0IPoqb/v49vPb2Hs4DQevHkKp2WleB1JRKKICj1Imtt8PPBaEcMz+vDil88mJkbXOReRU0sfigbJw2/tpKq+lXtnj1KZi4gnVOhBsPj9ffzq9R1cPnGw9pmLiGdU6EHwt01ljMxK5oHPTiZOJxCJiEfUPp/Qc+tLebe4mhl5A1TmIuIpNdAn0NTq40dLPyAvM5nvXpnvdRwRiXIq9E/g5cIDVNa18KXzTiMpPtbrOCIS5XTY4kma/4d1vLDpAAmxMcwaOcDrOCIiKvST0dTq472d1czIy+DXN01hYFqS15FERLTL5UQ1tfq444nVHG5s5c5zRqrMRSRkaIZ+gpZs3M/KndXcf/1ELszP9jqOiMiHNEM/ASuLq3jgtR2YweUTBnsdR0TkI1ToAWpt7+CLT6zGzHjqjhkkJ+ofNyISWlToAWjzdfDlp9fR3NbBd6/M56zTM72OJCLyMSr0ADz6zm6WFZbz3SvzuWjcIK/jiIj0SIUegNLDjaT3ief2s/K8jiIiclQBFbqZXWJm282s2Mzu62H9LWa2yf+10swmBT+qN5rbfLy4+QCD03V4ooiEtuMWupnFAg8BlwL5wE1m1v3CJbuBTzvnJgI/ABYGO6gXnHM8/f4+jjS2cUPBcK/jiIgcUyCHakwHip1zuwDMbDEwF9j6jw2ccyu7bL8KGBbMkF5wzvHfL2zl8Xf3MDWnHxfpmHMRCXGBFPpQoKTL61JgxjG2vwN4qacVZjYPmAeQk5MTYMRTr6K2mV+9sYOnVu3juqnD+L8bJmKmuxCJSGgLpNB7ajLX44Zm59NZ6Gf3tN45txD/7piCgoIe38NrzjnuenING0trmDaiP1+/eLTKXETCQiCFXgp03YE8DCjrvpGZTQQeBS51zlUHJ96pt3xHFRtLa7h15gh+cPV4r+OIiAQskKNcVgOjzCzPzBKAG4ElXTcwsxzgL8Ctzrmi4Mc8dR54rYik+Bi+OmeU11FERE7IcWfozrl2M5sPLANigUXOuUIzu9u/fgHwHWAA8LB/90S7c66g92L3jhU7qli37wj/cu5IBqQkeh1HROSEBHRBEufcUmBpt2ULujy/E7gzuNFOrfLaZr71/GZSEuOYd+5Ir+OIiJwwnSnq96Wn1lJe28wTXzxTs3MRCUsqdOCdHZWs23eE28/KY9qIDK/jiIiclKgv9JrGNv77b1tJiIvh7k+f5nUcEZGTFvWF/vBbxeyoqOeRW6aS3ife6zgiIictqgv9+fX7+c07uzh3dBazxwz0Oo6IyCcStYVe09TGv/9xA+OGpLPgc1N1NqiIhL2oLfQnVu6hw8F/XT6Wvgm6nZyIhL+oLPTmNh/Liyrp3zeeGSMHeB1HRCQooq7QnXP86+/XsXbfYb552Viv44iIBE3UFXpZTTNvfFDBl2eP0k0rRCSiRF2h76yoByAno6/HSUREgisqPg3cdqCWZYUH2VFez4ubD5AYF8P5Z2R5HUtEJKgivtC37K/h2kdW0trewYgBfZkzNpsrJg7W9VpEJOJEbKEXldfxzOoSfrdqLymJcXznunyunjLU61giIr0mIgt9a1ktl/3qHQAuHpfNj66ZoBm5iES8iCr0plYf3/jzJpZsLCM1KY7/vXYCV0wc4nUsEZFTIqIK/fGVu1mysYz555/OXeeO1MW2RCSqREyhd3Q4dlY0kBAXw9cvPsPrOCIip1zEHIf+i9eK+PO6UsYPSfM6ioiIJyKi0Btb21m0Yjc5GX159u5PeR1HRMQTEVHom0praGj1cde5I4mN0WVwRSQ6hX2h765q4N6n1xMfa1ycn+11HBERz4R9of/81SJqmtr4zhX5DExL8jqOiIhnwrrQnXNU17cwYkBfbp2V63UcERFPhXWhP7ZiNyt3VjN7jHa1iIiEdaEv31HF6OwU/kPHnYuIhG+hN7S0s2V/Df37JhCjI1tERMK30NfuPcyhhlaumzrM6ygiIiEhbAv9mTUlxBjMOk03eRYRgTAu9IraZgpyMxiuW8mJiAABFrqZXWJm282s2Mzu62G9mdmv/Os3mdnU4Ef9J+ccbT5HnPadi4h86LiFbmaxwEPApUA+cJOZ5Xfb7FJglP9rHvBIkHN+qKXdx11PrmFDyREyddMKEZEPBTJDnw4UO+d2OedagcXA3G7bzAWedJ1WAf3MbHCQswLwTlEVr22r4EvnncZPb5jYGz9CRCQsBVLoQ4GSLq9L/ctOdBvMbJ6ZrTGzNZWVlSeaFYD+yQlcNmEQt84cQWJc7Em9h4hIJArkBhc97ah2J7ENzrmFwEKAgoKCj60PxLQR/Zk2YtrJfKuISEQLZIZeCgzv8noYUHYS24iISC8KpNBXA6PMLM/MEoAbgSXdtlkC3OY/2mUmUOOcOxDkrCIicgzH3eXinGs3s/nAMiAWWOScKzSzu/3rFwBLgcuAYqARuL33IouISE8Cukm0c24pnaXdddmCLs8dcE9wo4mIyIkI2zNFRUTko1ToIiIRQoUuIhIhVOgiIhHCOj/P9OAHm1UCe4+yOhOoOoVxelukjQcib0waT+iLtDGd7HhGOOeyelrhWaEfi5mtcc4VeJ0jWCJtPBB5Y9J4Ql+kjak3xqNdLiIiEUKFLiISIUK10Bd6HSDIIm08EHlj0nhCX6SNKejjCcl96CIicuJCdYYuIiInSIUuIhIhQqrQj3cz6lBlZovMrMLMtnRZlmFmr5rZDv9j/y7r/tM/xu1mdrE3qY/OzIab2Ztmts3MCs3sK/7lYTkmM0sys/fNbKN/PN/3Lw/L8fyDmcWa2Xoze8H/OtzHs8fMNpvZBjNb418WtmMys35m9icz+8D/tzSr18fjnAuJLzovzbsTGAkkABuBfK9zBZj9XGAqsKXLsvuB+/zP7wN+4n+e7x9bIpDnH3Os12PoNp7BwFT/81SgyJ87LMdE5x21UvzP44G/AzPDdTxdxvXvwB+AF8L9d86fcw+Q2W1Z2I4JeAK40/88AejX2+MJpRl6IDejDknOueXAoW6L59L5fyj+x6u7LF/snGtxzu2m8xry009FzkA55w4459b5n9cB2+i8R2xYjsl1qve/jPd/OcJ0PABmNgy4HHi0y+KwHc8xhOWYzCyNzoneYwDOuVbn3BF6eTyhVOgB3Wg6jGQ7/12b/I8D/cvDapxmlgtMoXNWG7Zj8u+e2ABUAK8658J6PMAvgf8AOrosC+fxQOd/ZF8xs7VmNs+/LFzHNBKoBB737xZ71MyS6eXxhFKhB3Sj6QgQNuM0sxTgz8BXnXO1x9q0h2UhNSbnnM85N5nO+91ON7Pxx9g8pMdjZlcAFc65tYF+Sw/LQmY8XZzlnJsKXArcY2bnHmPbUB9THJ27YR9xzk0BGujcxXI0QRlPKBV6pN1outzMBgP4Hyv8y8NinGYWT2eZ/9459xf/4rAeE4D/n71vAZcQvuM5C7jKzPbQuWtytpk9RfiOBwDnXJn/sQJ4js5dDuE6plKg1P8vQYA/0VnwvTqeUCr0QG5GHU6WAJ/3P/888Ncuy280s0QzywNGAe97kO+ozMzo3Pe3zTn38y6rwnJMZpZlZv38z/sAc4APCNPxOOf+0zk3zDmXS+ffyRvOuc8RpuMBMLNkM0v9x3PgImALYTom59xBoMTMzvAvugDYSm+Px+tPgrt9KnwZnUdU7AS+5XWeE8j9NHAAaKPzv7R3AAOA14Ed/seMLtt/yz/G7cClXufvYTxn0/nPvU3ABv/XZeE6JmAisN4/ni3Ad/zLw3I83cZ2Hv88yiVsx0PnPueN/q/Cf/z9h/mYJgNr/L93zwP9e3s8OvVfRCRChNIuFxER+QRU6CIiEUKFLiISIVToIiIRQoUuIhIhVOgiIhFChS4iEiH+P2Lu8Y3gEaVXAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(ecdf.x, ecdf.y)\n",
    "# plt.scatter(sample, quantile, c='r')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "a48edb53158dd8a459f3e7f13fbcd78db9c99f511fd1a201dbc67d4baa3f7a7a"
  },
  "kernelspec": {
   "display_name": "Python 3.9.7 64-bit ('power': conda)",
   "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.9.7"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}