for Alex

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAADsxJREFUeJzt3V+MZvVdx/H3p7tLrGkV0x0V94+zJuuf1ULF6ZZoY7dp\nqgtcbJr0Atq0KSnZYEpT79h4UWJ6Q9NomlrKZoIbbGLYmEDqCovUiyoXSN2hocCCkJFiGUqyQ6ut\ntSZk5OvFPNXH6cw8Z3bOzvPMr+9XMuH8+eU5n/yyfHLmnOecSVUhSWrLG8YdQJLUP8tdkhpkuUtS\ngyx3SWqQ5S5JDbLcJalBlrskNchyl6QGWe6S1KCd4zrw7t27a3p6elyHl6Rt6fHHH3+1qqZGjRtb\nuU9PTzM3Nzeuw0vStpTkX7uM87KMJDXIcpekBlnuktQgy12SGmS5S1KDLHdJapDlLkkNGlnuSU4l\nuZDk6RHj3p5kKcn7+4snSboYXc7c7wGOrjcgyQ7g08CXe8gkSdqkkU+oVtUjSaZHDPs4cB/w9h4y\njTR94sGtOMyqXrzj+rEdW5K62vQ19yR7gPcBd20+jiSpD33cUP0scFtVvT5qYJLjSeaSzC0uLvZw\naEnSavp4cdgMcDoJwG7guiRLVfWllQOrahaYBZiZmakeji1JWsWmy72qDvxwOck9wAOrFbskaeuM\nLPck9wJHgN1JFoDbgV0AVXXykqaTJF2ULt+WubHrh1XVRzaVRpLUC59QlaQGWe6S1CDLXZIaZLlL\nUoMsd0lqkOUuSQ2y3CWpQZa7JDXIcpekBlnuktQgy12SGmS5S1KDLHdJapDlLkkNstwlqUGWuyQ1\nyHKXpAZZ7pLUIMtdkhpkuUtSg0aWe5JTSS4keXqN/R9M8mSSp5I8muSq/mNKkjaiy5n7PcDRdfZ/\nA3hXVb0V+BQw20MuSdIm7Bw1oKoeSTK9zv5Hh1YfA/ZuPpYkaTP6vub+UeChnj9TkrRBI8/cu0ry\nbpbL/Z3rjDkOHAfYv39/X4eWJK3Qy5l7kiuBu4FjVfXttcZV1WxVzVTVzNTUVB+HliStYtPlnmQ/\ncD/woap6fvORJEmbNfKyTJJ7gSPA7iQLwO3ALoCqOgl8EngL8IUkAEtVNXOpAkuSRuvybZkbR+y/\nGbi5t0SSpE3zCVVJapDlLkkNstwlqUGWuyQ1yHKXpAZZ7pLUIMtdkhpkuUtSgyx3SWqQ5S5JDbLc\nJalBlrskNchyl6QGWe6S1CDLXZIaZLlLUoMsd0lqkOUuSQ2y3CWpQZa7JDXIcpekBo0s9ySnklxI\n8vQa+5Pkc0nmkzyZ5Or+Y0qSNqLLmfs9wNF19l8LHBz8HAfu2nwsSdJmjCz3qnoE+M46Q44BX6xl\njwGXJ7mir4CSpI3r45r7HuClofWFwbYfkeR4krkkc4uLiz0cWpK0mi29oVpVs1U1U1UzU1NTW3lo\nSfqx0ke5vwzsG1rfO9gmSRqTPsr9DPDhwbdmrgG+W1Wv9PC5kqSLtHPUgCT3AkeA3UkWgNuBXQBV\ndRI4C1wHzAM/AG66VGElSd2MLPequnHE/gI+1lsiSdKm+YSqJDXIcpekBlnuktQgy12SGmS5S1KD\nLHdJapDlLkkNstwlqUGWuyQ1yHKXpAZZ7pLUIMtdkhpkuUtSgyx3SWqQ5S5JDbLcJalBlrskNchy\nl6QGWe6S1CDLXZIa1KnckxxN8lyS+SQnVtn/00n+JsnXk5xPclP/USVJXY0s9yQ7gDuBa4FDwI1J\nDq0Y9jHgmaq6CjgC/EmSy3rOKknqqMuZ+2FgvqpeqKrXgNPAsRVjCnhzkgBvAr4DLPWaVJLUWZdy\n3wO8NLS+MNg27PPArwHfAp4CPlFVr/eSUJK0YX3dUP194AngF4C3AZ9P8lMrByU5nmQuydzi4mJP\nh5YkrdSl3F8G9g2t7x1sG3YTcH8tmwe+Afzqyg+qqtmqmqmqmampqYvNLEkaoUu5nwMOJjkwuEl6\nA3BmxZhvAu8BSPJzwK8AL/QZVJLU3c5RA6pqKcmtwMPADuBUVZ1Pcstg/0ngU8A9SZ4CAtxWVa9e\nwtySpHWMLHeAqjoLnF2x7eTQ8reA3+s3miTpYvmEqiQ1yHKXpAZZ7pLUIMtdkhpkuUtSgyx3SWqQ\n5S5JDbLcJalBlrskNchyl6QGWe6S1CDLXZIaZLlLUoMsd0lqkOUuSQ2y3CWpQZ3+WIf+z/SJB8dy\n3BfvuH4sx5W0PXnmLkkNstwlqUGWuyQ1qFO5Jzma5Lkk80lOrDHmSJInkpxP8g/9xpQkbcTIG6pJ\ndgB3Au8FFoBzSc5U1TNDYy4HvgAcrapvJvnZSxVYkjRalzP3w8B8Vb1QVa8Bp4FjK8Z8ALi/qr4J\nUFUX+o0pSdqILuW+B3hpaH1hsG3YLwM/k+Tvkzye5MN9BZQkbVxf33PfCfwW8B7gjcA/Jnmsqp4f\nHpTkOHAcYP/+/T0dWpK0Upcz95eBfUPrewfbhi0AD1fVf1bVq8AjwFUrP6iqZqtqpqpmpqamLjaz\nJGmELuV+DjiY5ECSy4AbgDMrxvw18M4kO5P8JPAO4Nl+o0qSuhp5WaaqlpLcCjwM7ABOVdX5JLcM\n9p+sqmeT/C3wJPA6cHdVPX0pg0uS1tbpmntVnQXOrth2csX6Z4DP9BdNknSxfEJVkhpkuUtSgyx3\nSWqQ5S5JDbLcJalBlrskNchyl6QGWe6S1CDLXZIaZLlLUoMsd0lqkOUuSQ2y3CWpQZa7JDXIcpek\nBlnuktQgy12SGmS5S1KDLHdJapDlLkkNstwlqUGdyj3J0STPJZlPcmKdcW9PspTk/f1FlCRt1Mhy\nT7IDuBO4FjgE3Jjk0BrjPg18ue+QkqSN6XLmfhiYr6oXquo14DRwbJVxHwfuAy70mE+SdBG6lPse\n4KWh9YXBtv+VZA/wPuCu9T4oyfEkc0nmFhcXN5pVktTRzp4+57PAbVX1epI1B1XVLDALMDMzUz0d\n+8fC9IkHx3bsF++4fmzHlnRxupT7y8C+ofW9g23DZoDTg2LfDVyXZKmqvtRLSknShnQp93PAwSQH\nWC71G4APDA+oqgM/XE5yD/CAxS5J4zOy3KtqKcmtwMPADuBUVZ1Pcstg/8lLnFGStEGdrrlX1Vng\n7Iptq5Z6VX1k87EkSZvhE6qS1CDLXZIaZLlLUoMsd0lqkOUuSQ2y3CWpQZa7JDXIcpekBlnuktQg\ny12SGmS5S1KDLHdJapDlLkkNstwlqUGWuyQ1yHKXpAZZ7pLUIMtdkhpkuUtSgyx3SWpQp3JPcjTJ\nc0nmk5xYZf8HkzyZ5Kkkjya5qv+okqSuRpZ7kh3AncC1wCHgxiSHVgz7BvCuqnor8Clgtu+gkqTu\nupy5Hwbmq+qFqnoNOA0cGx5QVY9W1b8NVh8D9vYbU5K0EV3KfQ/w0tD6wmDbWj4KPLTajiTHk8wl\nmVtcXOyeUpK0Ib3eUE3ybpbL/bbV9lfVbFXNVNXM1NRUn4eWJA3Z2WHMy8C+ofW9g23/T5IrgbuB\na6vq2/3E0ySYPvHgWI774h3Xj+W4Ugu6nLmfAw4mOZDkMuAG4MzwgCT7gfuBD1XV8/3HlCRtxMgz\n96paSnIr8DCwAzhVVeeT3DLYfxL4JPAW4AtJAJaqaubSxZYkrafLZRmq6ixwdsW2k0PLNwM39xtN\nknSxfEJVkhpkuUtSgyx3SWqQ5S5JDbLcJalBlrskNchyl6QGWe6S1CDLXZIaZLlLUoMsd0lqkOUu\nSQ3q9OIwaRzG9R558F3y2v48c5ekBlnuktQgy12SGmS5S1KDLHdJapDlLkkN8quQ0irG9TVMv4Kp\nvnQ6c09yNMlzSeaTnFhlf5J8brD/ySRX9x9VktTVyDP3JDuAO4H3AgvAuSRnquqZoWHXAgcHP+8A\n7hr8V9IG+BuD+tLlzP0wMF9VL1TVa8Bp4NiKMceAL9ayx4DLk1zRc1ZJUkddrrnvAV4aWl/gR8/K\nVxuzB3hlU+kkbYlxvurhx9FW/Ka0pTdUkxwHjg9Wv5/kuYv8qN3Aq/2kuuTM2r/tkhO2T9btkhO2\nT9Y1c+bTm/rcX+wyqEu5vwzsG1rfO9i20TFU1Sww2yXYepLMVdXMZj9nK5i1f9slJ2yfrNslJ2yf\nrOPO2eWa+zngYJIDSS4DbgDOrBhzBvjw4Fsz1wDfrSovyUjSmIw8c6+qpSS3Ag8DO4BTVXU+yS2D\n/SeBs8B1wDzwA+CmSxdZkjRKp2vuVXWW5QIf3nZyaLmAj/UbbV2bvrSzhczav+2SE7ZP1u2SE7ZP\n1rHmzHIvS5Ja4rtlJKlBE13u2+m1Bx2yHkny3SRPDH4+Oaacp5JcSPL0GvsnYk475JyU+dyX5CtJ\nnklyPsknVhkzKXPaJeukzOtPJPmnJF8fZP3jVcaMfV475hzPnFbVRP6wfPP2X4BfAi4Dvg4cWjHm\nOuAhIMA1wFcnOOsR4IEJmNffBa4Gnl5j/6TM6aickzKfVwBXD5bfDDw/wf9Ou2SdlHkN8KbB8i7g\nq8A1kzavHXOOZU4n+cx9O732oEvWiVBVjwDfWWfIRMxph5wToapeqaqvDZb/A3iW5aezh03KnHbJ\nOhEGc/X9wequwc/KG4Rjn9eOOcdikst9rVcabHTMVuia47cHvz4+lOTXtybahk3KnHYxUfOZZBr4\nTZbP3oZN3JyukxUmZF6T7EjyBHAB+Luqmsh57ZATxjCnk1zurfkasL+qrgT+DPjSmPNsdxM1n0ne\nBNwH/GFVfW+cWUYZkXVi5rWq/ruq3sbyE++Hk/zGuLKsp0POsczpJJd7b6892AIjc1TV937461st\nPzewK8nurYvY2aTM6bomaT6T7GK5LP+yqu5fZcjEzOmorJM0r0OZ/h34CnB0xa6JmVdYO+e45nSS\ny307vfZgZNYkP58kg+XDLM/9t7c86WiTMqfrmpT5HGT4c+DZqvrTNYZNxJx2yTpB8zqV5PLB8htZ\n/nsS/7xi2NjntUvOcc3pxP6ZvdpGrz3omPX9wB8kWQL+C7ihBrfSt1KSe1m+e787yQJwO8s3gSZq\nTjvknIj5BH4H+BDw1OC6K8AfAfuHsk7EnNIt66TM6xXAX2T5jwW9AfirqnpgAv//75JzLHPqE6qS\n1KBJviwjSbpIlrskNchyl6QGWe6S1CDLXZIaZLlLUoMsd0lqkOUuSQ36H0fucbY4ZMCYAAAAAElF\nTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x19ee7f602b0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGeBJREFUeJzt3Xl0nfV95/H315YlW4tXCeMV28GEmJCFqkBJJ6XpZpM0\n7qRJBpqBgSTjwyR00s7MSdyeNrTDaZuGJAWCseOwmUBwSXDAocYmGG/gVd4t25LlXfImWbYsy9Zy\npe/8cS+MEFqu5Hvvc+9zP69zdHTv8/ys53N+WB8eP/dZzN0REZFwGRR0ABERSTyVu4hICKncRURC\nSOUuIhJCKncRkRBSuYuIhJDKXUQkhFTuIiIhpHIXEQmhnKA2XFxc7FOmTAlq8yIiGWnr1q117l7S\n17jAyn3KlCmUlZUFtXkRkYxkZkfjGafDMiIiIaRyFxEJIZW7iEgIqdxFREJI5S4iEkJ9lruZPW1m\nZ8xsTw/rzcweM7MqM9tlZjclPqaIiPRHPHvuzwIze1k/C5ge+5oDzL/yWCIiciX6LHd3XwvU9zJk\nNvCcR20ERprZuEQFFBGR/kvEMfcJwPFO76tjyz7AzOaYWZmZldXW1g54gycbLg/4z4qIZIOUfqDq\n7gvdvdTdS0tK+rx6tkdLttXQ3qEHe4uI9CQR5V4DTOr0fmJsWdKcamhm1f4zydyEiEhGS0S5LwXu\niZ01cyvQ4O4nE/Bze7V4y7Fkb0JEJGP1eeMwM3sRuB0oNrNq4EFgCIC7LwCWAXcAVcAl4L5khe1s\nVUUtpy80M3b40FRsTkQko/RZ7u5+Vx/rHfhmwhLFqb3D+UXZcR74zPRUb1pEJO1l9BWqL5VVE/1/\ni4iIdJbR5X6s/hLvVJ0NOoaISNrJ6HIHeHGzPlgVEekq48v9jb2nONPYHHQMEZG0kvHl3tbuvLTl\neN8DRUSySMaXO8CLm4/ToStWRUTeE4pyrzl/mbd0xaqIyHtCUe4AL2yK64HgIiJZITTlvqayluP1\nl4KOISKSFkJT7h0OP9dpkSIiQIjKHeDftxynua096BgiIoELVbnXN7Xy650ngo4hIhK4UJU7wKIN\nR4KOICISuNCV+56aC2w92tsjX0VEwi905Q7w7HqdFiki2S2U5b58z0nOXND9ZkQke4Wy3Nvanec3\nau9dRLJXKMsdoue867RIEclWoS33uoutvLqjJugYIiKBCG25Azz19uGgI4iIBCLU5V55+iJrK2uD\njiEiknKhLneAJ7X3LiJZKPTlvraylgOnG4OOISKSUqEvd9CxdxHJPllR7ku211Db2BJ0DBGRlMmK\ncm+NdPDMO9p7F5HskRXlDvD8xqM0tUSCjiEikhJZU+4XmiO8qCc1iUiWyJpyh+gHq23tHUHHEBFJ\nuqwq95MNzby6Q09qEpHwi6vczWymmVWYWZWZze1m/Qgz+7WZ7TSzcjO7L/FRE2Ph2oO4e9AxRESS\nqs9yN7PBwDxgFjADuMvMZnQZ9k1gr7t/HLgd+KGZ5SY4a0JUnr7IG3tPBx1DRCSp4tlzvxmocvdD\n7t4KLAZmdxnjQJGZGVAI1ANpe2rKvFVVQUcQEUmqeMp9AnC80/vq2LLOHgc+ApwAdgPfcve0/eRy\nV3WDbigmIqGWqA9U/wTYAYwHPgE8bmbDuw4yszlmVmZmZbW1wZbr49p7F5EQi6fca4BJnd5PjC3r\n7D5giUdVAYeB67v+IHdf6O6l7l5aUlIy0MwJsflwPVuO1AeaQUQkWeIp9y3AdDObGvuQ9E5gaZcx\nx4A/ADCzscCHgUOJDJoMP35Le+8iEk59lru7R4AHgBXAPuAldy83s/vN7P7YsIeA28xsN7AS+I67\n1yUrdKKsraxl+7FzQccQEUm4nHgGufsyYFmXZQs6vT4B/HFio6XGI28eYNFXbw46hohIQmXVFard\nWaO9dxEJoawvd4juvYuIhInKHe29i0j4qNxjtPcuImGico9ZU1nL1qM6711EwkHl3snDKyqCjiAi\nkhAq9042Hqpn3QHdc0ZEMp/KvYsfaO9dREJA5d7FzuoG3ig/FXQMEZEronLvxg/fqKSjQ09rEpHM\npXLvRsXpRn61veuNL0VEMofKvQc/+k0lLZH2oGOIiAyIyr0HNecvs2j9kaBjiIgMiMq9F/NWHaTh\nUlvQMURE+k3l3ouGy208sVoP9BCRzKNy78Oz649Qc/5y0DFERPpF5d6HlkgH31++P+gYIiL9onKP\nw9KdJ3RLYBHJKCr3OLjDQ6/tDTqGiEjcVO5x2nbsPL/eeSLoGCIicVG598P3Xt9Pc5subBKR9Kdy\n74ea85f56dpDQccQEemTyr2fnlh9kBM6NVJE0pzKvZ8ut7XzT8v2BR1DRKRXKvcB+I9dJ1l/sC7o\nGCIiPVK5D9A/Lt1LpL0j6BgiIt1SuQ9QxelGFm04GnQMEZFuqdyvwL/9ppLTF5qDjiEi8gEq9ytw\nsSXC/9WVqyKShlTuV+g/dp1kbWVt0DFERN5H5Z4A3311jx7JJyJpJa5yN7OZZlZhZlVmNreHMbeb\n2Q4zKzezNYmNmd6OnL3EvLf0UA8RSR99lruZDQbmAbOAGcBdZjajy5iRwBPA5939BuBLScia1uav\nOUjl6cagY4iIAPHtud8MVLn7IXdvBRYDs7uM+QtgibsfA3D3M4mNmf7a2p25L++io8ODjiIiEle5\nTwCOd3pfHVvW2XXAKDNbbWZbzeye7n6Qmc0xszIzK6utDd+HkNuOnef5TTr3XUSCl6gPVHOA3wI+\nC/wJ8Pdmdl3XQe6+0N1L3b20pKQkQZtOL99fXsHJBt1YTESCFU+51wCTOr2fGFvWWTWwwt2b3L0O\nWAt8PDERM8vFlgh/u2R30DFEJMvFU+5bgOlmNtXMcoE7gaVdxrwK/K6Z5ZhZPnALkLW3TlxVUcsv\nyo73PVBEJEn6LHd3jwAPACuIFvZL7l5uZveb2f2xMfuA5cAuYDPwpLvvSV7s9PfQa3t1awIRCUxO\nPIPcfRmwrMuyBV3ePww8nLhome1Cc/TwzFP3/nbQUUQkC+kK1SRauf8ML2+tDjqGiGQhlXuS/cOv\ny/VYPhFJOZV7kjU2R/g/v9iJuy5uEpHUUbmnwPqDZ3nmnSNBxxCRLKJyT5F/Xb6fqjO694yIpIbK\nPUVaIh18a/EOWiN67qqIJJ/KPYXKT1zg+8v3Bx1DRLKAyj3FnnrnMGv05CYRSTKVe4q5w/9+aSd1\nF1uCjiIiIaZyD0DdxRadHikiSaVyD8jqiloWrDkUdAwRCSmVe4B++EYFZUfqg44hIiGkcg9QpMP5\nyxe3c66pNegoIhIyKveAnWxo5q9f2qHj7yKSUCr3NLC6opZHVx4IOoaIhIjKPU08tvIAqyrOBB1D\nREJC5Z4mOhz+avEOjtdfCjqKiISAyj2NNFxu4/7nt9Lc1h50FBHJcCr3NFN+4gLf/uWuoGOISIZT\nuaehpTtPMH/1waBjiEgGU7mnqYdX7NcHrCIyYCr3NNXh8D9f3K4HfIjIgKjc01hjc4SvPltGva5g\nFZF+UrmnuWP1l5jzXBktEZ1BIyLxU7lngLKj55j78u6gY4hIBlG5Z4hfba/hR29UBB1DRDKEyj2D\nPPZWFYs3Hws6hohkAJV7hvm7V/awar9OkRSR3qncM0ykw/nmz7ex4/j5oKOISBpTuWegS63t3PfM\nZqrOXAw6ioikqbjK3cxmmlmFmVWZ2dxexv22mUXM7IuJiyjdOXepjXue2sTJhstBRxGRNNRnuZvZ\nYGAeMAuYAdxlZjN6GPevwBuJDindO9HQzN1PbdZFTiLyAfHsud8MVLn7IXdvBRYDs7sZ95fAy4A+\n7UuhqjMXuefpTVxobgs6ioikkXjKfQJwvNP76tiy95jZBOA/A/MTF03itafmAvc9s4VLrZGgo4hI\nmkjUB6qPAN9x947eBpnZHDMrM7Oy2traBG1aALYePcfXni3Tgz5EBIiv3GuASZ3eT4wt66wUWGxm\nR4AvAk+Y2Z91/UHuvtDdS929tKSkZICRpScbDp3lvz+ngheR+Mp9CzDdzKaaWS5wJ7C08wB3n+ru\nU9x9CvBL4Bvu/krC00qf1h2oY87PtupGYyJZrs9yd/cI8ACwAtgHvOTu5WZ2v5ndn+yA0n9rK2uZ\n85wKXiSb5cQzyN2XAcu6LFvQw9h7rzyWXKk1lbV8fVEZP72nlKFDBgcdR0RSTFeohti6A3Xc+8xm\nmlp0Fo1ItlG5h9zGQ/Xc/ZTOgxfJNir3LLDt2HnuWriRuostQUcRkRRRuWeJ8hMX+NKCDVSfuxR0\nFBFJAZV7Fjlc18SXFmyg6kxj0FFEJMlU7lnmZEMzX1ywgbIj9UFHEZEkUrlnofOX2vjKk5tYvudU\n0FFEJElU7lmqJdLBN17YyqL1R4KOIiJJoHLPYh0ODy4t5x+WltPR4UHHEZEEUrkLz64/wpyflemW\nwSIhonIXAN7cd4Y/n7+BmvN6bJ9IGKjc5T37Tl5g9uNv60wakRBQucv71F1s5S9+uonFm48FHUVE\nroDKXT6gtb2DuUt28zdLdtMa6fXhWiKSplTu0qMXNx/jyz/ZwKmG5qCjiEg/qdylVzuOn+dzP17H\n2wfqgo4iIv2gcpc+1V1s5Z6nN/Fvv6nU+fAiGULlLnHpcHh05QHufnoTZxp1mEYk3ancpV/eqTrL\nrEfWsariTNBRRKQXKnfpt7NNrXz12S089NpePYRbJE2p3GVA3OGptw8z+/F32H/qQtBxRKQLlbtc\nkf2nGvn84+/w07WH9GGrSBpRucsVa4108E/L9nHnwo0cPdsUdBwRQeUuCbT5SD0zH1nHovVHcNde\nvEiQVO6SUJfb2nlwaTn/5ScbOVh7Meg4IllL5S5JsflIPbMeXcfjbx2grV33pxFJNZW7JE1rpIMf\nvFHJ5x57my26jbBISqncJekqTjfy5Z9s4Nu/3Mm5ptag44hkBZW7pIQ7vFRWze//cDXPbzyq0yZF\nkkzlLil1/lIbf/fKHj4/7222HtWhGpFkUblLIPbUXODP52/ggZ9v03NbRZIgrnI3s5lmVmFmVWY2\nt5v1XzGzXWa228zWm9nHEx9Vwui1XSf5zA9W8/CK/TQ2twUdRyQ0+ix3MxsMzANmATOAu8xsRpdh\nh4Hfc/cbgYeAhYkOKuHVEulg3qqD3P7wahatP6JTJ0USIJ4995uBKnc/5O6twGJgducB7r7e3c/F\n3m4EJiY2pmSDs02tPLi0nD/80Rpe3VGjq1xFrkA85T4BON7pfXVsWU++Brze3Qozm2NmZWZWVltb\nG39KySpHz17iW4t3MOvRdby593TQcUQyUkI/UDWz3yda7t/pbr27L3T3UncvLSkpSeSmJYT2n2rk\n68+VMfvxt/VwEJF+iqfca4BJnd5PjC17HzP7GPAkMNvdzyYmngjsrG7gvme28Gfz3mHlPu3Ji8Qj\nnnLfAkw3s6lmlgvcCSztPMDMJgNLgLvdvTLxMUVgx/HzfG1RGXc8uo5lu0/qQiiRXuT0NcDdI2b2\nALACGAw87e7lZnZ/bP0C4LvAGOAJMwOIuHtp8mJLNtt78gLfeGEb04oLmPPpaXzhponk5uiSDZHO\nLKgzEkpLS72srGxAf/bvX9nDzzYeTXAiyVRXFeXx326bwldumczI/Nyg44gklZltjWfnWbs7kvHO\nNLbw8IoKbvveW3z31T0c0n3kRfo+LCOSKS61tvPchqP8bONRPj29hHs/NYXbryshdqhQJKuo3CV0\n3GFNZS1rKmu5Zkw+X7llMl8unaRDNpJVdFhGQu3o2Uv887L93PLPK/nrf9/BpkM6S1eyg/bcJSu0\nRDr41fYafrW9hg+VFPDl0kl84aaJlBTlBR1NJCm05y5Z52BtE//y+n5+519W8vVFZSzfc5LWiG5W\nJuGiPXfJWpEO5819p3lz32lG5g/hszeO4ws3TeCmyaP0IaxkPJW7CNEnRL2w6RgvbDrGxFHD+NOP\nj2f2J8Zz/dXDg44mMiAqd5Euqs9dZv7qg8xffZBrryrkjhvH8bmPjeO6sUVBRxOJm8pdpBdVZy7y\n2MoDPLbyANNKCpj10auZecM4bpw4IuhoIr1SuYvE6VBtE/NWHWTeqoOMHzGUP5wxlj+ecTW3TBvN\nkME6N0HSi8pdZABONDTz3IajPLfhKEV5Ofyn64r5g+vH8nsfLqG4UKdXSvBU7iJXqLElwrLdp1i2\n+xRmcOOEEdx+XQmfvq6ET04exeBBOvNGUk/lLpJA7rCruoFd1Q089lYVRUNzuO1DY/jda4v51LXF\nTCspDDqiZAmVu0gSNTZHWFF+mhXl0SdIjR8xlFs/NIbfmTaG264tZsLIYQEnlLBSuYuk0ImGZpZs\nq2HJtuiTKieOGsbNU0dz69QxlE4ZpT17SRiVu0iAqs9dpvrc/y/74sI8Sq8ZxW9dM4qbrhnFjRNG\n6ClTMiAqd5E0UnexheXlp1hefgqA3JxB3DB+OJ+cNIpPTB7JJyeNZNLo/IBTSiZQuYuksdZIB9uP\nnWf7sfPwTnTZ6IJcbpwwgo9NHMFHY9/HjdCxe3k/lbtIhqlvan3vYSTvKi7M5SPjhnPD+BHcMH44\nHxk3nGnFBQzSaZhZS+UuEgJ1F1tZd6COdQfq3ls2dMggPjy2iA9fXcT1Vw/n+quLmD62SPewzxIq\nd5GQam7rYGd1AzurG963fHRBLtOvKuTaqwqZflUh140t4kNXFTJ2+NCAkkoyqNxFskx9UyubDtez\n6XD9+5YX5eUwraSAaSWFTC0uYFpJAVPGFDC1uICCPFVFptF/MREBordR6G5PH6CkKI8pY/KZPLqA\nqcX5TBqdz+TR+VwzpoDRBXrweDpSuYtIn2obW6htbGHLkXMfWFeQO5iJo/KZNHoYE0flM2HkMCaO\nGsaEUcMYP3KYbqQWEJW7iFyRptZ2Kk43UnG6sdv1eTmDGD9yGONGDGXciGGMHzmUscOHvvd+7PA8\nRhfk6tGGCaZyF5Gkaol0cLiuicN1TT2OyR08iKuG53FVUR5jh0fLv6Qoj5Ki6LLiwuj3MYV5ustm\nnFTuIhK41vaO2K0YLvc6bpDBqPxcxhTmUlwYLfsxBbmUFEX3/kcX5DKmIJeR+e9+H5K1/yJQuYtI\nxuhwONvUytmmVipPX+xz/OBBxohhQxiZP4TR+dHSH5k/hFH5QxiZn/veupHDoq/f/SoampPxF4Cp\n3EUktNo7nPqmVuqbWjlEz4eFujKDwrycWNEPYfjQHIbHSn/40Oj3oqE5FOZFXxcOzaEoL4eioUMo\nyBtMYV4OhXk55AT4+MW4yt3MZgKPAoOBJ939e13WW2z9HcAl4F5335bgrCIiKeEevRd/Y3ME6P1Q\nUW/ycgZRmJdDQV4O+bnR0r9x4gge/NMbEhe2B32Wu5kNBuYBfwRUA1vMbKm77+00bBYwPfZ1CzA/\n9l1EJGu1RDpoiUQPI70r0uEp2XY8/2a4Gahy90Pu3gosBmZ3GTMbeM6jNgIjzWxcgrOKiEic4jks\nMwE43ul9NR/cK+9uzATg5BWl60FBXo6uihORjDRi2JCUbCelH6ia2RxgDsDkyZMH/HPmzrqeubOu\nT1QsEZHQieewTA0wqdP7ibFl/R2Duy9091J3Ly0pKelvVhERiVM85b4FmG5mU80sF7gTWNplzFLg\nHou6FWhw96QckhERkb71eVjG3SNm9gCwguipkE+7e7mZ3R9bvwBYRvQ0yCqip0Lel7zIIiLSl7iO\nubv7MqIF3nnZgk6vHfhmYqOJiMhABXf5lIiIJI3KXUQkhFTuIiIhpHIXEQkhlbuISAhZ9ESXADZs\nVgscHeAfLwbqEhgnUdI1F6RvNuXqH+XqnzDmusbd+7wKNLByvxJmVubupUHn6Cpdc0H6ZlOu/lGu\n/snmXDosIyISQip3EZEQytRyXxh0gB6kay5I32zK1T/K1T9Zmysjj7mLiEjvMnXPXUREepHW5W5m\nM82swsyqzGxuN+vNzB6Lrd9lZjelSa7bzazBzHbEvr6bolxPm9kZM9vTw/qg5quvXCmfLzObZGar\nzGyvmZWb2be6GZPy+YozVxDzNdTMNpvZzliuf+xmTBDzFU+uQH4fY9sebGbbzey1btYld77cPS2/\niN5e+CAwDcgFdgIzuoy5A3gdMOBWYFOa5LodeC2AOfs0cBOwp4f1KZ+vOHOlfL6AccBNsddFQGWa\n/P2KJ1cQ82VAYez1EGATcGsazFc8uQL5fYxt+38BP+9u+8mer3Tec0/XB3PHkysQ7r4WqO9lSCAP\nMo8jV8q5+0l33xZ73QjsI/rc385SPl9x5kq52BxcjL0dEvvq+oFdEPMVT65AmNlE4LPAkz0MSep8\npXO59/TQ7f6OCSIXwG2xf2q9bmY3JDlTvIKYr3gFNl9mNgX4JNG9vs4Cna9eckEA8xU7xLADOAP8\nxt3TYr7iyAXB/P16BPg20NHD+qTOVzqXeybbBkx2948BPwZeCThPugtsvsysEHgZ+Ct3v5Cq7fal\nj1yBzJe7t7v7J4g+I/lmM/toKrbblzhypXy+zOxzwBl335rsbfUkncs9YQ/mTnUud7/w7j8VPfoU\nqyFmVpzkXPEIYr76FNR8mdkQogX6grsv6WZIIPPVV66g/365+3lgFTCzy6pA/371lCug+foU8Hkz\nO0L00O1nzOz5LmOSOl/pXO7p+mDuPnOZ2dVmZrHXNxOd57NJzhWPtHyQeRDzFdveU8A+d/9RD8NS\nPl/x5ApovkrMbGTs9TDgj4D9XYYFMV995gpivtz9b9x9ortPIdoRb7n7f+0yLKnzFdczVIPgafpg\n7jhzfRH4H2YWAS4Dd3rs4/FkMrMXiZ4ZUGxm1cCDRD9gCmy+4swVxHx9Crgb2B07Xgvwt8DkTrmC\nmK94cgUxX+OARWY2mGg5vuTurwX9+xhnrkB+H7uTyvnSFaoiIiGUzodlRERkgFTuIiIhpHIXEQkh\nlbuISAip3EVEQkjlLiISQip3EZEQUrmLiITQ/wNKmis4+/2UcQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x19ee83c5780>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Task1 Perform 1000 Exp-distributed\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import stats\n",
    "import numpy as np\n",
    "import math \n",
    "def rng(m=2**32, a=1103515245, c=12345):\n",
    "    rng.current = (a * rng.current + c) % m\n",
    "    return rng.current / m\n",
    "rng.current = 1\n",
    "def Exp(l=2):\n",
    "    tmp = rng(m=2**32, a=1103515245, c=12345)\n",
    "    Exp.current = (-1/l)*math.log(1-tmp)\n",
    "    return Exp.current\n",
    "random = [Exp() for i in range(1000)]\n",
    "plt.hist(random, normed=True)\n",
    "plt.show()\n",
    "dist = stats.expon()\n",
    "x = np.linspace(0, 4, 100)\n",
    "plt.fill_between(x, 0, dist.pdf(x)) \n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-0.08984736  0.04202152  0.46272007 ...,  0.00536743 -0.21344696\n",
      "   0.00464254]\n",
      " [-0.03137377  0.14853798 -0.18128694 ..., -0.01776202 -0.00871934\n",
      "   0.04421667]]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAD8CAYAAACcjGjIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAE+pJREFUeJzt3X+s3fV93/Hny9c4IwYFCI5jfm0geU2pVFhyRViG2rBA\nBtY6h2qVQFOCqlYWE1RrNE2y1Cnrn1mqrlJWGuRuKERag1olFCt1wgBtoj9GyyUiYJIQHI8KOw52\nCIMkgM213/vjfp2e3s/9cXzP95zra54P6eh+v9/P5/M+n+/92n75+z3nfE+qCkmSBq1b7QlIkk4/\nhoMkqWE4SJIahoMkqWE4SJIahoMkqdFLOCS5N8nhJHsXaU+SzybZl+TpJO8faLspyXNd284+5iNJ\nGk1fZw6fB25aov1mYGv32AF8DiDJFHB3134lcFuSK3uakyRphXoJh6p6DPjhEl22A1+oOY8D5yXZ\nAlwD7Kuq/VV1DLi/6ytJWkXrJ/Q8FwMvDqwf6LYttP2DCxVIsoO5sw42btz4gfe9733jmelp5jtP\n7l/tKay6JOMpPDXVf82x3XFgPHXrxJjmu4buvPCPP3DFak9hYp588skfVNWmYfpOKhxGVlW7gF0A\n09PTNTMzs8ozmowb1/3Kak9heOvG8I8tsG7DWWOpm3PP7b/oW8f6rwkwpn/ET7zx5ljq1rh+D2Pw\n8MyfrPYUJibJ3w7bd1LhcBC4dGD9km7bWYtslyStokm9lXU38InuXUvXAq9W1SHgCWBrksuTbABu\n7fpKklZRL2cOSb4IfBi4MMkB4D8xd1ZAVd0D7AG2AfuA14Ff7dpmk9wFPARMAfdW1bN9zEmStHK9\nhENV3bZMewF3LtK2h7nwkCSdJvyEtCSpYThIkhqGgySpYThIkhqGgySpYThIkhqGgySpYThIkhqG\ngySpYThIkhqGgySpYThIkhqGgySpYThIkhqGgySpYThIkhqT+g5pnWZy1obea647Z2PvNQGoE+Op\nOw5j+L0CY/sdrFuXsdQ98Ub/devo0d5ranGeOUiSGr2EQ5KbkjyXZF+SnQu0/4ckT3WPvUmOJ7mg\na3shyTNd20wf85EkjWbky0pJpoC7gRuBA8ATSXZX1TdP9qmq3wF+p+v/S8Anq+qHA2Wur6ofjDoX\nSVI/+jhzuAbYV1X7q+oYcD+wfYn+twFf7OF5JUlj0kc4XAy8OLB+oNvWSPJO4CbgSwObC3gkyZNJ\ndiz2JEl2JJlJMnPkyJEepi1JWsykX5D+JeAv511Suq6qrgZuBu5M8gsLDayqXVU1XVXTmzZtmsRc\nJeltq49wOAhcOrB+SbdtIbcy75JSVR3sfh4GHmDuMpUkaRX1EQ5PAFuTXJ5kA3MBsHt+pyTvAn4R\neHBg28Yk555cBj4K7O1hTpKkEYz8bqWqmk1yF/AQMAXcW1XPJrmja7+n63oL8D+r6icDwzcDDyQ5\nOZc/qqqvjTonSdJoevmEdFXtAfbM23bPvPXPA5+ft20/cFUfc5Ak9cdPSEuSGoaDJKlhOEiSGoaD\nJKlhOEiSGoaDJKlhOEiSGoaDJKlhOEiSGoaDJKlhOEiSGr3cW0njk/XjOUTrznvXWOqOQ6bG9H+Y\nd57df83jx/uvCTA7nro1lqqw7pyNvdc8cWJcs9VCPHOQJDUMB0lSw3CQJDUMB0lSw3CQJDUMB0lS\nw3CQJDV6CYckNyV5Lsm+JDsXaP9wkleTPNU9PjXsWEnS5I38CaskU8DdwI3AAeCJJLur6pvzuv55\nVf3LFY6VJE1QH2cO1wD7qmp/VR0D7ge2T2CsJGlM+giHi4EXB9YPdNvm+1CSp5N8NcnPneJYkuxI\nMpNk5siRIz1MW5K0mEm9IP114LKq+nngvwJ/eqoFqmpXVU1X1fSmTZt6n6Ak6e/0EQ4HgUsH1i/p\ntv1UVb1WVT/ulvcAZyW5cJixkqTJ6yMcngC2Jrk8yQbgVmD3YIck702Sbvma7nlfHmasJGnyRn63\nUlXNJrkLeAiYAu6tqmeT3NG13wP8a+DfJpkF3gBuraoCFhw76pwkSaPp5csCuktFe+Ztu2dg+feB\n3x92rCRpdfkJaUlSw3CQJDUMB0lSw3CQJDXG8+316k3e8Y7xFF6X3ktm/Xj+ONW7zhlL3dnzzh5L\n3XFY/8rrY6mb18ZSljp6tPea687+B73X1OI8c5AkNQwHSVLDcJAkNQwHSVLDcJAkNQwHSVLDcJAk\nNQwHSVLDcJAkNQwHSVLDcJAkNQwHSVLDcJAkNXoJhyQ3JXkuyb4kOxdo/zdJnk7yTJK/SnLVQNsL\n3fanksz0MR9J0mhGvsdykingbuBG4ADwRJLdVfXNgW7/F/jFqnolyc3ALuCDA+3XV9UPRp2LJKkf\nfZw5XAPsq6r9VXUMuB/YPtihqv6qql7pVh8HLunheSVJY9JHOFwMvDiwfqDbtphfA746sF7AI0me\nTLJjsUFJdiSZSTJz5MiRkSYsSVraRL8JLsn1zIXDdQObr6uqg0neAzyc5NtV9dj8sVW1i7nLUUxP\nT9dEJixJb1N9nDkcBC4dWL+k2/b3JPl54L8B26vq5ZPbq+pg9/Mw8ABzl6kkSauoj3B4Atia5PIk\nG4Bbgd2DHZJcBnwZ+HhVfWdg+8Yk555cBj4K7O1hTpKkEYx8WamqZpPcBTwETAH3VtWzSe7o2u8B\nPgW8G/iDJACzVTUNbAYe6LatB/6oqr426pwkSaPp5TWHqtoD7Jm37Z6B5V8Hfn2BcfuBq+Zv19/J\n1NR4Ch892n/NDRv6rwnMnv/OsdT90X/8Ue81/89VX+q9JsC/+OVPjKXu+jePjaVuxlDTFxony09I\nS5IahoMkqWE4SJIahoMkqWE4SJIahoMkqWE4SJIahoMkqWE4SJIahoMkqWE4SJIahoMkqWE4SJIa\nhoMkqWE4SJIahoMkqWE4SJIahoMkqdFLOCS5KclzSfYl2blAe5J8tmt/Osn7hx0rSZq8kcMhyRRw\nN3AzcCVwW5Ir53W7GdjaPXYAnzuFsZKkCevjzOEaYF9V7a+qY8D9wPZ5fbYDX6g5jwPnJdky5FhJ\n0oSt76HGxcCLA+sHgA8O0efiIccCkGQHc2cdXHbZZaPNeC1Zl/HUPWtD/zXXT/VfE6j14/kdvGPq\neO81v/PWT3qvCTD7zrPGUnd91Vjq1tGjY6h5rPeaWtyaeUG6qnZV1XRVTW/atGm1pyNJZ7Q+zhwO\nApcOrF/SbRumz1lDjJUkTVgfZw5PAFuTXJ5kA3ArsHten93AJ7p3LV0LvFpVh4YcK0masJHPHKpq\nNsldwEPAFHBvVT2b5I6u/R5gD7AN2Ae8DvzqUmNHnZMkaTR9XFaiqvYwFwCD2+4ZWC7gzmHHSpJW\n15p5QVqSNDmGgySpYThIkhqGgySpYThIkhqGgySpYThIkhqGgySpYThIkhqGgySpYThIkhqGgySp\nYThIkhqGgySpYThIkhqGgySp0cuX/Wh86q3ZsdTN+jEc+jfe7L8msP6VN8ZS9/hn3t17zTuP3tV7\nTYB3vPT/xlKX2ePjqTsGdezYak/hbcUzB0lSY6RwSHJBkoeTPN/9PH+BPpcm+V9Jvpnk2ST/bqDt\nt5McTPJU99g2ynwkSf0Y9cxhJ/BoVW0FHu3W55sF/n1VXQlcC9yZ5MqB9t+rqqu7h98lLUmngVHD\nYTtwX7d8H/Cx+R2q6lBVfb1b/hHwLeDiEZ9XkjRGo4bD5qo61C1/H9i8VOck/wj4J8BfD2z+jSRP\nJ7l3octSA2N3JJlJMnPkyJERpy1JWsqy4ZDkkSR7F3hsH+xXVQXUEnXOAb4E/GZVvdZt/hxwBXA1\ncAj43cXGV9WuqpququlNmzYtv2eSpBVb9v2MVXXDYm1JXkqypaoOJdkCHF6k31nMBcP/qKovD9R+\naaDPHwJfOZXJS5LGY9TLSruB27vl24EH53dIEuC/A9+qqv8yr23LwOotwN4R5yNJ6sGo4fBp4MYk\nzwM3dOskuSjJyXce/TPg48A/X+Atq59J8kySp4HrgU+OOB9JUg9G+phsVb0MfGSB7d8DtnXLfwFk\nkfEfH+X5JUnj4SekJUkNw0GS1DAcJEkNw0GS1DAcJEkNw0GS1DAcJEkNw0GS1DAcJEkNw0GS1DAc\nJEmNke6tpPE78frrY6k7dVb/h74ynv9rrDv8yljqnv3qT/ovmgVvIza6t2bHUraOHh1L3bE4fny1\nZ/C24pmDJKlhOEiSGoaDJKlhOEiSGoaDJKlhOEiSGiOFQ5ILkjyc5Pnu5/mL9Huh+67op5LMnOp4\nSdJkjXrmsBN4tKq2Ao9264u5vqqurqrpFY6XJE3IqOGwHbivW74P+NiEx0uSxmDUcNhcVYe65e8D\nmxfpV8AjSZ5MsmMF40myI8lMkpkjR46MOG1J0lKWvYdCkkeA9y7Q9FuDK1VVSWqRMtdV1cEk7wEe\nTvLtqnrsFMZTVbuAXQDT09OL9pMkjW7ZcKiqGxZrS/JSki1VdSjJFuDwIjUOdj8PJ3kAuAZ4DBhq\nvCRpska9rLQbuL1bvh14cH6HJBuTnHtyGfgosHfY8ZKkyRs1HD4N3JjkeeCGbp0kFyXZ0/XZDPxF\nkm8AfwP8WVV9banxkqTVNdJ9m6vqZeAjC2z/HrCtW94PXHUq4yVJq8tPSEuSGoaDJKlhOEiSGoaD\nJKlhOEiSGv1/y7z6VeP5MPiJN97svea6jOf/GlVj+j/Mm/3/DnL22b3XBKi33hpLXWZnx1K23jza\nf80xzVUL88xBktQwHCRJDcNBktQwHCRJDcNBktQwHCRJDcNBktQwHCRJDcNBktQwHCRJDcNBktQw\nHCRJDcNBktQYKRySXJDk4STPdz/PX6DPzyR5auDxWpLf7Np+O8nBgbZto8xHktSPUc8cdgKPVtVW\n4NFu/e+pqueq6uqquhr4APA68MBAl9872V5Ve0acjySpB6OGw3bgvm75PuBjy/T/CPDdqvrbEZ9X\nkjRGo4bD5qo61C1/H9i8TP9bgS/O2/YbSZ5Ocu9Cl6VOSrIjyUySmSNHjowwZUnScpYNhySPJNm7\nwGP7YL+qKmDRry1LsgH4V8CfDGz+HHAFcDVwCPjdxcZX1a6qmq6q6U2bNi03bUnSCJb9mtCqumGx\ntiQvJdlSVYeSbAEOL1HqZuDrVfXSQO2fLif5Q+Arw01bkjROo15W2g3c3i3fDjy4RN/bmHdJqQuU\nk24B9o44H0lSD0YNh08DNyZ5HrihWyfJRUl++s6jJBuBG4Evzxv/mSTPJHkauB745IjzkST1YNnL\nSkupqpeZewfS/O3fA7YNrP8EePcC/T4+yvNr5ero0d5rHn9rtveaAFPnbBxLXdal95InXvtR7zUB\nWDeez6vWG2+Mp+7seP4saHL8hLQkqWE4SJIahoMkqWE4SJIahoMkqWE4SJIahoMkqWE4SJIahoMk\nqWE4SJIahoMkqWE4SJIahoMkqWE4SJIahoMkqWE4SJIahoMkqWE4SJIaI4VDkl9J8mySE0mml+h3\nU5LnkuxLsnNg+wVJHk7yfPfz/FHmI0nqx6hnDnuBXwYeW6xDkingbuBm4ErgtiRXds07gUeraivw\naLcuSVplI4VDVX2rqp5bpts1wL6q2l9Vx4D7ge1d23bgvm75PuBjo8xHktSP9RN4jouBFwfWDwAf\n7JY3V9Whbvn7wObFiiTZAezoVn+cZH4oXQj8YPTpnnbWzn4dP6Xew+/XqyuYy+paO8fs1JyR+5UE\nztB9o92vfzjswGXDIckjwHsXaPqtqnpw2CdaTlVVklqifRewa7H2JDNVtejrHmuV+7X2nKn7dqbu\nF5y5+zbKfi0bDlV1w0oKDzgIXDqwfkm3DeClJFuq6lCSLcDhEZ9LktSDSbyV9Qlga5LLk2wAbgV2\nd227gdu75duB3s5EJEkrN+pbWW9JcgD4p8CfJXmo235Rkj0AVTUL3AU8BHwL+OOqerYr8WngxiTP\nAzd06yu16CWnNc79WnvO1H07U/cLztx9W/F+pWrRy/ySpLcpPyEtSWoYDpKkxpoNh1O4dccLSZ5J\n8lSSmUnOcSVGvSXJ6WrYW6WspeO13DHInM927U8nef9qzPNUDbFfH07yaneMnkryqdWY56lKcm+S\nw0n2LtK+Vo/Xcvu1suNVVWvyAfws8DPA/waml+j3AnDhas+3z/0CpoDvAlcAG4BvAFeu9tyX2a/P\nADu75Z3Af17Lx2uYYwBsA74KBLgW+OvVnndP+/Vh4CurPdcV7NsvAO8H9i7SvuaO15D7taLjtWbP\nHGq4W3esOUPu11K3JDldnWm3ShnmGGwHvlBzHgfO6z7Pczpbi3+2hlJVjwE/XKLLWjxew+zXiqzZ\ncDgFBTyS5MnuFhxngoVuSXLxKs1lWMPeKmWtHK9hjsFaPE7DzvlD3aWXryb5uclMbezW4vEa1ikf\nr0ncW2nFerp1x3VVdTDJe4CHk3y7S9pVM6lbkkzaUvs1uFK15K1STrvjpcbXgcuq6sdJtgF/Cmxd\n5TlpcSs6Xqd1ONTot+6gqg52Pw8neYC50+ZV/cemh/1a6pYkq2ap/Uoy1K1STsfjtYhhjsFpeZyW\nseycq+q1geU9Sf4gyYVVtdZvXLcWj9eyVnq8zujLSkk2Jjn35DLwUea+g2KtW+qWJKerZW+VssaO\n1zDHYDfwie5dMNcCrw5cWjtdLbtfSd6bzN3KNMk1zP078vLEZ9q/tXi8lrXi47Xar7SP8Ar9Lcxd\nEzwKvAQ81G2/CNjTLV/B3LstvgE8y9xlm1Wf+6j71a1vA77D3DtL1sJ+vZu5L3R6HngEuGCtH6+F\njgFwB3BHtxzmvujqu8AzLPGuutPpMcR+3dUdn28AjwMfWu05D7lfXwQOAW91f8d+7Qw5Xsvt14qO\nl7fPkCQ1zujLSpKklTEcJEkNw0GS1DAcJEkNw0GS1DAcJEkNw0GS1Pj/0ZX/JukDg5gAAAAASUVO\nRK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x26fa5bf1390>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Task5 Perform Box-Muller transformation\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import stats\n",
    "import numpy as np\n",
    "import random\n",
    "def gen_normal(n):\n",
    "    cur = 0\n",
    "    res1 = []\n",
    "    res2 = []\n",
    "    while cur < n:\n",
    "        x = 2*np.random.rand()-1\n",
    "        y = 2*np.random.rand()-1\n",
    "        s = np.square(x) + np.square(y) \n",
    "        if (s <= 1):\n",
    "            m = -2*np.log(s)/s\n",
    "            l = np.sqrt(m)\n",
    "            z1 = x*s\n",
    "            z2 = y*s\n",
    "            res1.append(z1)\n",
    "            res2.append(z2)\n",
    "            cur += 1\n",
    "            Y=np.vstack((res1,res2))\n",
    "    return Y\n",
    "Z=gen_normal(10000)\n",
    "print(Z)\n",
    "X=Z[0] \n",
    "Y=Z[1]\n",
    "plt.hist2d(Y,X)\n",
    "plt.axis('equal')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAADdlJREFUeJzt3X+oX/ddx/Hny7g6mY6BiWwm0UQMlrB1ddylkwm6aSVt\nx7L6A9PNjbqNUFh0g8mWWphIEVoGOsFqCDUMsRgGW2doM9Kqhf6xVXM7265pmxFiNQkbvdvUWSbr\nYt/+8f1Gv7nem3tu8r33nPvJ8wGh95zzyff7Js195X0/53M+31QVkqS2fF/fBUiSps9wl6QGGe6S\n1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXo+/t64/Xr19eWLVv6entJWpMef/zxb1TVhqXG\n9RbuW7ZsYXZ2tq+3l6Q1Kcm/dBnntIwkNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWp\nQYa7JDWotydUpSHYsu/BBc8/f9dNq1yJNF127pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrsk\nNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ3qFO5JdiY5keRkkn0XGffmJOeS/Nr0SpQkLdeS\n4Z5kHXAPcAOwHbglyfZFxt0NPDTtIiVJy9Olc98BnKyqU1X1EnAI2LXAuN8GPgu8MMX6JEmXoMuH\ndWwETk8cnwGumxyQZCNwM/A24M1Tq04aGD/cQ2vFtG6ofgr4eFW9fLFBSfYkmU0yOzc3N6W3liTN\n16VzPwtsnjjeND43aQY4lARgPXBjknNV9fnJQVV1ADgAMDMzU5datCTp4rqE+zFgW5KtjEJ9N/Du\nyQFVtfX810k+DTwwP9glSatnyXCvqnNJ9gJHgXXAwao6nuS28fX9K1yjtOoWm1uX1oounTtVdQQ4\nMu/cgqFeVbdeflmSpMvhE6qS1KBOnbu0lq3GFItLJDU0du6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y\n3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAb5YR1ac/xg\nDGlpdu6S1CDDXZIa5LSMmrEan5UqrRV27pLUIMNdkhpkuEtSg5xzl1aQyzbVFzt3SWqQ4S5JDTLc\nJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQe4to8Fyf3bp0nXq3JPsTHIi\nyckk+xa4vivJU0meSDKb5OemX6okqaslO/ck64B7gOuBM8CxJIer6pmJYX8HHK6qSnIN8Bng6pUo\nWJK0tC6d+w7gZFWdqqqXgEPArskBVfViVdX48FVAIUnqTZdw3wicnjg+Mz53gSQ3J3kOeBB4/3TK\nkyRdiqmtlqmq+6vqauBdwJ0LjUmyZzwnPzs3Nzett5YkzdMl3M8CmyeON43PLaiqHgV+Msn6Ba4d\nqKqZqprZsGHDsouVJHXTJdyPAduSbE1yFbAbODw5IMlPJcn46zcBPwB8c9rFSpK6WXK1TFWdS7IX\nOAqsAw5W1fEkt42v7wd+FXhfku8B/wX8xsQNVknSKuv0EFNVHQGOzDu3f+Lru4G7p1uaJOlSuf2A\nJDXI7QekHiy2tcLzd920ypWoVXbuktQgO3f1zg3CpOmzc5ekBhnuktQgw12SGmS4S1KDDHdJapDh\nLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSg9w4TBoQtwLWtNi5S1KDDHdJapDhLkkNMtwl\nqUGGuyQ1yNUy0hrgKhotl527JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1\nyCdUtWoWe8pS0vTZuUtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGdQr3JDuTnEhyMsm+Ba6/J8lT\nSb6S5ItJ3jj9UiVJXS0Z7knWAfcANwDbgVuSbJ837J+Bn6+qNwB3AgemXagkqbsunfsO4GRVnaqq\nl4BDwK7JAVX1xar6t/HhY8Cm6ZYpSVqOLuG+ETg9cXxmfG4xHwC+sNCFJHuSzCaZnZub616lJGlZ\npnpDNcnbGIX7xxe6XlUHqmqmqmY2bNgwzbeWJE3osrfMWWDzxPGm8bkLJLkGuBe4oaq+OZ3yJEmX\noku4HwO2JdnKKNR3A++eHJDkx4HPAe+tqq9OvUqtKW4QJvVvyXCvqnNJ9gJHgXXAwao6nuS28fX9\nwCeAHwH+LAnAuaqaWbmyJcHi/5A+f9dNq1yJhqbTlr9VdQQ4Mu/c/omvPwh8cLqlSZIulU+oSlKD\nDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchw\nl6QGddryV1qIH8oxXO7zLjt3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoN8iEm6\ngvhw05XDzl2SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yKWQWpL7tktrj527JDXIcJekBjktI+mi\nU28+vbo22blLUoMMd0lqkOEuSQ0y3CWpQZ3CPcnOJCeSnEyyb4HrVyf5UpLvJvnd6ZcpSVqOJVfL\nJFkH3ANcD5wBjiU5XFXPTAz7FvA7wLtWpEpJ0rJ06dx3ACer6lRVvQQcAnZNDqiqF6rqGPC9FahR\nkrRMXda5bwROTxyfAa5bmXLUJ7cZkNqxqjdUk+xJMptkdm5ubjXfWpKuKF3C/SyweeJ40/jcslXV\ngaqaqaqZDRs2XMpLSJI66BLux4BtSbYmuQrYDRxe2bIkSZdjyTn3qjqXZC9wFFgHHKyq40luG1/f\nn+S1wCzwauDlJB8BtlfVt1ewdknSIjptHFZVR4Aj887tn/j664ymayRJA+ATqpLUILf8lXRRiy2R\ndSvgYbNzl6QGGe6S1CCnZa5APokqtc/OXZIaZLhLUoMMd0lqkOEuSQ3yhqqkS+L692Gzc5ekBhnu\nktQgp2Ua5Vp26cpm5y5JDTLcJalBhrskNchwl6QGeUN1jfPGqaSF2LlLUoMMd0lqkNMykqbKbQmG\nwc5dkhpk575GeONUa50d/eqyc5ekBhnuktQgw12SGmS4S1KDDHdJapCrZST1ylU0K8Nw74l/oSWt\nJKdlJKlBhrskNchpmYHxSVRJ02C4Sxok70tdHsN9hdmJS9Nl6HdjuE+JIS5pSAx3SU2wo7+Qq2Uk\nqUGdOvckO4E/AdYB91bVXfOuZ3z9RuA7wK1V9eUp1zoITr9Ia8uV2tEvGe5J1gH3ANcDZ4BjSQ5X\n1TMTw24Ato1/XQf8+fi/kjRIrYd+l859B3Cyqk4BJDkE7AImw30X8JdVVcBjSV6T5HVV9bWpV7xK\n7NAlrWVdwn0jcHri+Az/vytfaMxGYEXC3eCVtFJWI19W46eDVV0tk2QPsGd8+GKSE0v8lvXAN1a2\nqssy5PqGXBsMu74h1wbWdzkGUVvuXvRSl/p+ost7dAn3s8DmieNN43PLHUNVHQAOdCkMIMlsVc10\nHb/ahlzfkGuDYdc35NrA+i7HkGuD6dbXZSnkMWBbkq1JrgJ2A4fnjTkMvC8jbwH+Yy3Pt0vSWrdk\n515V55LsBY4yWgp5sKqOJ7ltfH0/cITRMsiTjJZC/tbKlSxJWkqnOfeqOsIowCfP7Z/4uoAPTbc0\nYBlTOD0Zcn1Drg2GXd+QawPruxxDrg2mWF9GuSxJaonbD0hSg9ZMuCf5aJJKsr7vWs5LcmeSp5I8\nkeShJD/Wd02TknwyyXPjGu9P8pq+a5qU5NeTHE/ycpJBrGBIsjPJiSQnk+zru55JSQ4meSHJ033X\nMl+SzUkeSfLM+P/ph/uuaVKSVyb5xyRPjuv7g75rmi/JuiT/lOSBabzemgj3JJuBXwb+te9a5vlk\nVV1TVdcCDwCf6LugeR4GXl9V1wBfBW7vuZ75ngZ+BXi070Lggq02bgC2A7ck2d5vVRf4NLCz7yIW\ncQ74aFVtB94CfGhgf3bfBd5eVW8ErgV2jlf2DcmHgWen9WJrItyBPwY+BgzqBkFVfXvi8FUMr76H\nqurc+PAxRs8fDEZVPVtVSz3Itpr+d6uNqnoJOL/VxiBU1aPAt/quYyFV9bXzmwVW1X8yCqmN/Vb1\nf2rkxfHhK8a/BvP9mmQTcBNw77Rec/DhnmQXcLaqnuy7loUk+cMkp4H3MLzOfdL7gS/0XcTALbaN\nhpYhyRbgZ4B/6LeSC42nPZ4AXgAerqoh1fcpRg3sy9N6wUF8WEeSvwVeu8ClO4DfYzQl04uL1VZV\nf1NVdwB3JLkd2Av8/pDqG4+5g9GPzfetZm3j916yPrUjyQ8BnwU+Mu8n295V1X8D147vPd2f5PVV\n1fv9iyTvAF6oqseT/MK0XncQ4V5Vv7TQ+SRvALYCT462jGcT8OUkO6rq633WtoD7GD0LsKrhvlR9\nSW4F3gH8YvWw7nUZf35D0GkbDS0sySsYBft9VfW5vutZTFX9e5JHGN2/6D3cgbcC70xyI/BK4NVJ\n/qqqfvNyXnTQ0zJV9ZWq+tGq2lJVWxj9mPym1Qr2pSTZNnG4C3iur1oWMv6QlY8B76yq7/RdzxrQ\nZasNLWD8gT1/ATxbVX/Udz3zJdlwfrVYkh9k9PkUg/h+rarbq2rTOON2A39/ucEOAw/3NeCuJE8n\neYrR1NGgln8Bfwr8MPDweLnm/qV+w2pKcnOSM8DPAg8mOdpnPeObz+e32ngW+ExVHe+zpklJ/hr4\nEvDTSc4k+UDfNU14K/Be4O3jv2tPjDvRoXgd8Mj4e/UYozn3qSw5HCqfUJWkBtm5S1KDDHdJapDh\nLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhr0Pz1HomVb0B/YAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x26fa85798d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH1dJREFUeJzt3Xt0XOV57/Hvg4lJAiYORkBiGzDEheOGYBKFkiZtStok\nhmbhpAmrcMilCa1LTjhJTps0pLmdlNPTk5U0t4OD42BCuBpWwGDA4Bs2thG2Lka25YtsWb5I8kWS\ndb+ONPOcPzTOGYRsjaSZ2Xtm/z5reaG9Z78zz0jopz3vfvf7mrsjIiLRcUbQBYiISG4p+EVEIkbB\nLyISMQp+EZGIUfCLiESMgl9EJGIU/CIiEaPgFxGJGAW/iEjEnBl0ASM5//zz/dJLLw26DBGRvFFR\nUdHs7kXpHBvK4L/00kspLy8PugwRkbxhZofSPVZdPSIiEaPgFxGJGAW/iEjEKPhFRCJGwS8iEjEK\nfhGRiFHwi4hEjIJfRCRiQnkDl0iYNHb28f2nd9LSHQPggnPP4u757+btZ08OuDKR8VHwi5xGVUM7\n//BgOUfb+163f3t9O0u+UMzsC6cEVJnI+KmrR+QUXqw6xs2LXn1D6AMcbunhb35Vwst7mwKoTGRi\nFPwiI6hr6eF/PF5J70D8lMd09g/y3x/dSmPnG/8wiISZgl9kBN97puq0oX9SR98g//bsrhxUJJI5\nCn6RYZ7ddoT11el34Ty3/SjrqhuzWJFIZin4RVK09w7ww3GcwX/v6Sp6Y6N/QhAJAwW/SIofr9xD\nc1f/mNvVt/by87V7s1CRSOalFfxmNs/Mqs2sxszuGuHx+Wa23cwqzazczD6UbluRsDje0cfjZXXj\nbv/Qq4do7xnIYEUi2TFq8JvZJGAhcAMwB7jVzOYMO2wtcLW7zwW+BNw3hrYiofC7koMMxH3c7Xti\ncR4pTXsRJJHApHPGfy1Q4+617h4DlgLzUw9w9y53P/kbczbg6bYVCYPeWJxHSw9P+HkeLDnEQDyR\ngYpEsied4J8OpH7+rU/uex0z+5SZ7QGeZ+isP+22yfYLkt1E5U1NuilGcuv3FXW0ZaCb5lhHH89t\nP5KBikSyJ2MXd919mbtfCXwSuHsc7Re7e7G7FxcVpbVQvEhGuDv3v3IwY8+3ZNOBjD2XSDakE/wN\nwMyU7RnJfSNy9w3AZWZ2/ljbigRhze5GDjR3Z+z5qho6eHX/iYw9n0impRP8ZcBsM5tlZpOBW4Dl\nqQeY2bvMzJJfvxc4CziRTluRoP32lcyfod+fhecUyZRRZ+d090EzuxNYCUwC7nf3nWZ2R/LxRcCn\ngc+b2QDQC/xt8mLviG2z9F5Exqy+tYdXazN/dr5uTyMnuvqZds5ZGX9ukYlKa1pmd18BrBi2b1HK\n1z8CfpRuW5GwePq1Bnz8IzhPaTDhPLvtCH/3wVmZf3KRCdKduxJpT72WvUtOy7L43CIToeCXyNpW\n10ZtU+Yu6r7h+evb2d/UlbXnFxkvBb9EVi7OyJdt1Vm/hI+CXyJpMJ7g2W3Zv9Fq2WsNeDYuIohM\ngIJfIunlvU2cSC6enk0Nbb2UHmjJ+uuIjIWCXyLpqRx2weTytUTSoeCXyOmJDbJ2z/Gcvd4LVUc1\ncZuEioJfImd9dRN9A7kL4o6+QU3hIKGi4JfIeaHqWCReU+RUFPwSKf2Dcdbtyf3C6Kt3HSOR0Oge\nCQcFv0TKpn3NdPUP5vx1m7tilB3U6B4JBwW/REqQXS7q7pGwUPBLZAzGE6zZnbvRPMOt2nlMN3NJ\nKCj4JTI217ZkZHnF8TrS3se2+vbAXl/kJAW/RMaLO48GXQIvqrtHQkDBL5Hg7qzaGVw3z0krdyr4\nJXgKfomEyro2Gjv7gy6DA83d1DR2Bl2GRJyCXyJh7e7cj90/lTUhqkWiScEvkRDkaJ7h1uwKTy0S\nTQp+KXj1rT3sORae7pWth1tpycGU0CKnouCXghe2M+yEw0sBTBshcpKCXwre2hCG7NoQdT1J9KQV\n/GY2z8yqzazGzO4a4fHbzGy7me0wsxIzuzrlsYPJ/ZVmVp7J4kVG09k3wJba8M2Rs2FvE/2D8aDL\nkIgaNfjNbBKwELgBmAPcamZzhh12APiwu18F3A0sHvb49e4+192LM1CzSNo27G0mFsJFULpjcTaH\n8A+SREM6Z/zXAjXuXuvuMWApMD/1AHcvcffW5OZmYEZmyxQZnzCN5hkubNceJDrSCf7pQF3Kdn1y\n36ncDryQsu3AGjOrMLMFp2pkZgvMrNzMypuamtIoS+T04glnfXX4+vdPUj+/BCWjF3fN7HqGgv9b\nKbs/5O5zGeoq+oqZ/flIbd19sbsXu3txUVFRJsuSiKo41EprgJOyjeZIex+7jnQEXYZEUDrB3wDM\nTNmekdz3Omb2HuA+YL67/2GBUXdvSP63EVjGUNeRSNblw5DJdSH+RCKFK53gLwNmm9ksM5sM3AIs\nTz3AzC4GngI+5+57U/afbWZTTn4NfAyoylTxIqfz0p7wd6Wou0eCcOZoB7j7oJndCawEJgH3u/tO\nM7sj+fgi4PvANOBXZgYwmBzBcyGwLLnvTOBRd38xK+9EJEVdSw97j3cFXcaoKuvaaOmOcd7Zk4Mu\nRSJk1OAHcPcVwIph+xalfP33wN+P0K4WuHr4fpFsy4duHhi6i3fdnkY+/T4NhJPc0Z27UpDCeLfu\nqeTLHykpHAp+KTg9sUE2154Y/cCQ2LCviYEQ3mQmhUvBLwVn075mYoP5E6SdfYOUHdRdvJI7Cn4p\nOPnYdfKSFmeRHFLwS0Fx97wcG5+Pf6wkfyn4paBUNXRwvCP4tXXHqra5mwPN3UGXIRGh4JeCEuZJ\n2Uajm7kkVxT8UlDyucskTAvCS2FT8EvBON7RR9WR9qDLGLeygy109IV3UjkpHAp+KRgv7WnEPegq\nxm8w4bxcrSnJJfsU/FIwCqGPvBDeg4Sfgl8KQt9AnFdq8udu3VNZv7eJeCKPP7ZIXlDwS0Eo2d9M\n70D+L17e1jNAxaHW0Q8UmQAFvxSENQU0IkbdPZJtCn4pCOvyeBjncPk0s6jkJwW/5L2qhnaOtvcF\nXUbG1DR2ceiE7uKV7FHwS95bvavwukYK8T1JeCj4Je8VYkiuKsD3JOGh4Je8Vt/aw66jHUGXkXEV\nh1pp7Y4FXYYUKAW/5LU1BXpmHE+4LvJK1ij4Ja+tLuChj6t3HQu6BClQCn7JW+29A2ypLdwlCzfu\na6avAG5Kk/BJK/jNbJ6ZVZtZjZndNcLjt5nZdjPbYWYlZnZ1um1Fxmt9dSODBTy9QU8szis1zUGX\nIQVo1OA3s0nAQuAGYA5wq5nNGXbYAeDD7n4VcDeweAxtRcZl1c7C7eY5KQrvUXIvnTP+a4Ead691\n9xiwFJifeoC7l7j7yQlGNgMz0m0rMh6xwQQv7y38KYzX7jlOooA/1Ugw0gn+6UBdynZ9ct+p3A68\nMM62Imkp2d9MV/9g0GVkXXNXjK2HNWmbZFZGL+6a2fUMBf+3xtF2gZmVm1l5U1Phn8nJxLywIzoj\nXl6ois57ldxIJ/gbgJkp2zOS+17HzN4D3AfMd/cTY2kL4O6L3b3Y3YuLiorSqV0iKp7wgh7GOdyL\nCn7JsHSCvwyYbWazzGwycAuwPPUAM7sYeAr4nLvvHUtbkbHaUnuClgjd1drQ1sv2+ragy5ACMmrw\nu/sgcCewEtgNPOHuO83sDjO7I3nY94FpwK/MrNLMyk/XNgvvQyIkil0fKyLUtSXZZx7C1amLi4u9\nvLw86DIkhBIJ57r/WEtjZ3/QpeTUpdPeyvpvXh90GRJiZlbh7sXpHKs7dyWvVBxujVzoAxw80cPu\nApyMToKh4Je8EqXRPMO9sONo0CVIgVDwS15ZuTPCwR/BaxuSHQp+yRvb6tpoaOsNuozA7Gvsoqax\nM+gypAAo+CVvPLvtSNAlBO7ZberukYlT8EtecHeeVx83z27XHz+ZOAW/5IXSAy0cbe8LuozA1TZ1\nU9XQHnQZkucU/JIXlqub5w/U5SUTpeCX0BuMJzSiJcVz248SxhsvJX8o+CX0NtU0R2puntE0tPVS\ncUhTNcv4Kfgl9NTN80b6nshEKPgl1PoG4qzW8oNvsGLHUeJamUvGScEvoba+upHOCKy0NVbNXTFK\n9mshdhkfBb+E2pNbR1y3R4Cn9L2RcVLwS2id6OpnfXVj0GWE1otVxyKx7rBknoJfQuuZyiMMxNWP\nfSq9A3FWbNfdzDJ2Cn4Jrd9X1AddQujpeyTjoeCXUNp9tINdWnhkVGWHWjh8oifoMiTPKPgllHQm\nmx53+P1Wfa9kbBT8EjqD8QTPVGrESrqe2lqvKRxkTBT8Ejov722iuUtTNKSrvrWXzbUtQZcheUTB\nL6HzWGld0CXkncfLDgddguSRtILfzOaZWbWZ1ZjZXSM8fqWZvWpm/Wb2jWGPHTSzHWZWaWblmSpc\nCtPR9l7Waez+mK2oOkarJrKTNI0a/GY2CVgI3ADMAW41sznDDmsBvgr85BRPc727z3X34okUK4Vv\naWmd5qAZh9hgQhfEJW3pnPFfC9S4e627x4ClwPzUA9y90d3LgIEs1CgREU84j5epm2e8His9rIu8\nkpZ0gn86kPrbWJ/cly4H1phZhZktGEtxEi1rdx/nWIeWVxyv2uZuXt1/IugyJA/k4uLuh9x9LkNd\nRV8xsz8f6SAzW2Bm5WZW3tTUlIOyJGwe2aILlBOl76GkI53gbwBmpmzPSO5Li7s3JP/bCCxjqOto\npOMWu3uxuxcXFRWl+/RSIOpaeti4T3/wJ2rVrmM0dfYHXYaEXDrBXwbMNrNZZjYZuAVYns6Tm9nZ\nZjbl5NfAx4Cq8RYrhevR0sPomu7EDcSdJ8p1nUROb9Tgd/dB4E5gJbAbeMLdd5rZHWZ2B4CZXWRm\n9cA/Ad81s3ozOxe4ENhkZtuAUuB5d38xW29G8lNvLM5jpeqiyJQHXz3IQDwRdBkSYmemc5C7rwBW\nDNu3KOXrYwx1AQ3XAVw9kQKl8D25tZ62Hg0Iy5TjHf08v/0on7xmLGMwJEp0564Eyt25/5UDQZdR\ncJZs0vdUTk3BL4FaX91EbVN30GUUnB0N7ZQe0Pw9MjIFvwRKZ6bZs2RTbdAlSEgp+CUwe451sKmm\nOegyCtbqXce1SIuMSMEvgVmyUWf72ZRwdP1ERqTgl0A0tPXytBZbybrHy+po7tINXfJ6Cn4JxKL1\n+xmI646tbOsdiHOfPlnJMAp+ybnGjj7dXZpDD28+RLvuk5AUCn7JucUbaukf1J2ludLVP8hvS3TW\nL/+fgl9yqrU7xqOaniHnHig5SFf/YNBlSEgo+CWnlmw6QE8sHnQZkdPWM8BDrx4KugwJCQW/5MyJ\nrn4eKDkYdBmR9ZuNtXT2qa9fFPySQ/esq1F3Q4BaumP8ZoPu5hUFv+RIfWuPVocKgSWbDmhcvyj4\nJTd+vmYfMY3kCVx3LM49L9UEXYYETMEvWbfveCfLXtNdumHx6JbD1LdqDp8oU/BL1v14ZTVxrasY\nGrF4gp+u2ht0GRIgBb9kVUlNM6t2HQ+6DBlmWWUDlXVtQZchAVHwS9bEE84Pn90VdBkyAnf4n8t3\n4q5PYlGk4JeseXjzIaqPdwZdhpxCZV0bT23VtZcoUvBLVrR2x/jpavUjh92PXtxDt+6tiBwFv2TF\nf66upr1Xd4mGXWNnP/9XwzsjJ63gN7N5ZlZtZjVmdtcIj19pZq+aWb+ZfWMsbaXwvHa4lUd1s1be\nWLKplupj6pKLklGD38wmAQuBG4A5wK1mNmfYYS3AV4GfjKOtFJCBeIK7ntyBRm/mj4G4860nt5PQ\nDy0y0jnjvxaocfdad48BS4H5qQe4e6O7lwHDP9uP2lYKy73r9+uCbh6qrGvTBHoRkk7wTwdSl0uq\nT+5Lx0TaSp6paezinnXqL85X/7mqWnf0RkRoLu6a2QIzKzez8qampqDLkTFKJJxvP7Vd8/Hkse5Y\nnO8sqwq6DMmBdIK/AZiZsj0juS8dabd198XuXuzuxUVFRWk+vYTF4o21lB1sDboMmaCX9zbx8GYt\n2FLo0gn+MmC2mc0ys8nALcDyNJ9/Im0lT1Q1tGvulwLy78/vZn9TV9BlSBaNGvzuPgjcCawEdgNP\nuPtOM7vDzO4AMLOLzKwe+Cfgu2ZWb2bnnqpttt6M5F7fQJyvP15JLK4unkLROxDn60srGdDPtGBZ\nGOfqKC4u9vLy8qDLkDT84Jkqfqe1XAvSV66/nG9+/Mqgy5A0mVmFuxenc2xoLu5K/lm585hCv4Dd\nu34/m/Y1B12GZIGCX8blQHM333hiW9BlSBYlHL669DWOtPUGXYpkmIJfxqw3FufLD1fQqcm9Cl5L\nd4z/9shWDdMtMAp+GbN/XbaDPZrbJTIq69q4+zmtq1BIFPwyJvdtrNX6uRH00OZDPFaqifcKhYJf\n0rZq5zH+94rdQZchAfne01W8UqOLvYVAwS9pqWpo5+uPV2rWzQgbTDhffriCmkZ18+U7Bb+M6khb\nL7f/royeWDzoUiRgHX2DfPGBMpq7+oMuRSZAwS+ndaKrn88u2cLxDv2iy5C6ll4+v6SUjj6tsJav\nFPxySp19A3zht6XUNnUHXYqEzK6jHXzpt2X06lNgXlLwy4j6BuLc/kA5VQ0dQZciIVV+qJV/fLhC\nY/zzkIJf3qA3FudLD5RRerAl6FIk5DbsbeLLD1fQP6gz/3yi4JfX6e4f5Av3l1Ky/0TQpUieWLun\nkX94sIK+AYV/vlDwyx909A3wuSVbdKYvY7ZhbxNfekB9/vlCwS8ANHb08be/3szWw21BlyJ5qmT/\nCW79zWZaumNBlyKjUPALNY1dfOpXJew+qgu5MjGVdW185t4S6lq0aHuYKfgjruxgC59ZVEKDpt6V\nDKlt7uZTvyphR3170KXIKSj4I2xp6WFu+80W2np0I45kVnNXPzf/uoRnKjWhXxidGXQBknuD8QT/\n9twuHtTqWZJFfQMJvra0kt1HO/mXj1/BGWdY0CVJks74I+Z4Rx//9b4tCn3JmUUv7+cLvy3lhOb3\nCQ0Ff4Ssr27kxl9spPSAhmtKbm3c18yNv9zI5lrdHxIGCv4IiA0m+I8XdvPFB8o4oaF2EpDjHf3c\ndt8WfrZ6L4NxTfMQpLSC38zmmVm1mdWY2V0jPG5m9svk49vN7L0pjx00sx1mVmlm5ZksXka360gH\nN92ziV+/XItrLn0JWDzh/GLtPj59bwk1jV1BlxNZowa/mU0CFgI3AHOAW81szrDDbgBmJ/8tAO4d\n9vj17j7X3YsnXrKkYyCe4J6X9vHJha9ofVwJnW317fz1Lzdy38Za4lrdJ+fSOeO/Fqhx91p3jwFL\ngfnDjpkPPOhDNgNTzewdGa5V0lRxqJVP/HITP1m1l5g+UktI9Q8m+F/P7+aTC1+hqkFj/nMpneCf\nDtSlbNcn96V7jANrzKzCzBaMt1AZXVtPjO8+vYObF5VQfVxn+ZIfdjS0M3/hK9z93C46tbhLTuRi\nHP+H3L3BzC4AVpvZHnffMPyg5B+FBQAXX3xxDsoqHPGE8/DmQ/xszV7djCV5KZ5wlmw6wDOVDXzz\n41dw8/tmatx/FqVzxt8AzEzZnpHcl9Yx7n7yv43AMoa6jt7A3Re7e7G7FxcVFaVXvbB293Fu/MVG\nfrB8p0Jf8l5zV4xvPbmDmxZuYtO+5qDLKVjpBH8ZMNvMZpnZZOAWYPmwY5YDn0+O7rkOaHf3o2Z2\ntplNATCzs4GPAVUZrD+yKg61cPOiEm7/Xbm6daTgVDV08NklW7jtvs1sr9eMsZk2alePuw+a2Z3A\nSmAScL+77zSzO5KPLwJWADcCNUAP8MVk8wuBZWZ28rUedfcXM/4uIqTiUAs/X7OPjTobkgh4peYE\nN93zCn/1Xy7gq385m/fMmBp0SQXBPISDu4uLi728XEP+U23c18Sil/fzSo3ufJTo+osrirjjw5dz\n3WXTgi4ldMysIt0h85qkLcQG4gme336UxRtq2aW58kVYX93E+uomrp45lQV/dhnz3n0Rk3QReMwU\n/CHU2NHHI1sO81jpYRo7NbGVyHDb6tr4yqNbeefb3sxt113CLe+fybRzzgq6rLyhrp6QSCScjTXN\nPF52mNW7jjMQD9/PRSSsJp95Bh//44u49f0z+cDl00heV4wUdfXkkYPN3Tz1WgNPVtRrFSyRcYoN\nJnh22xGe3XaES6a9lc+8dwafvGY6M897a9ClhZLO+APQ2NnHCzuO8XRlA69pcXORrDCD919yHjfN\nfSc3vPuigu8KGssZv4I/R4539LFq5zGe33GU0gMtaF4qkdyZdIbxgcumceNV7+Cjcy6kaErh/RFQ\n8IfE7qMdvLSnkVW7jrO9vk3TIouEwBkG11z8dj4250I+cuUFzL5wStAlZYSCPyAdfQOU1DSzYV8z\n6/c0cqS9L+iSRGQUM97+Fq6/4gL+bPb5/Om7zuecs/Lz0qcu7uZI30CcikOtvLr/BCX7m9lW3665\nxUXyTH1rLw9tPsRDmw/xpknGNTPfznWXT+NPL5/GNRdP5awzJwVdYsbpjH8MOvoGeO1wG6UHTlB6\noIVt9e3EBjXfvUihOuvMM5g7cyp/Mus8rp01jbkXTw3tJwKd8WdAIuHsbexkW10blXVtVBxqpaax\nSxdlRSKkfzDBlgMtbDnQAtRwhsEfXTiF913ydubOnMrcmVN51wXn5N19Awp+YDCeoKapi11HOtjR\n0E5VQzu7jnTQHYsHXZqIhEjCYc+xTvYc6+SRLYcBmHLWmcx557lcNf1tXDXjbcx5x7lcVnROqKeS\niFzwN3b2sfdYF3uOdVCd/AFWH+9Ul42IjEtn/2DKp4Ihb37TGVxx4RSuvOhcrrhoCldeNIU/umgK\n54fkXoKCDP54wmlo7WV/U1fyXzc1jZ3sa+zSYiUiknV9Awm21bezrf71awmfd/Zk3nXBOcy+4Bwu\nKzqHy4vO5vKic5g+9S05XXGsoIJ/f1MX//hQBYdP9GiRcREJnZbuGKUHWihN+XQAQxeRL5n2Vh74\n4rW8c+pbsl5HQQV/V98gNY1dQZchIjIm/YMJ9h7voidH1xXTWXpRREQKiIJfRCRiFPwiIhGj4BcR\niRgFv4hIxCj4RUQiJq3gN7N5ZlZtZjVmdtcIj5uZ/TL5+HYze2+6bUVEJLdGDX4zmwQsBG4A5gC3\nmtmcYYfdAMxO/lsA3DuGtiIikkPpnPFfC9S4e627x4ClwPxhx8wHHvQhm4GpZvaONNuKiEgOpXPn\n7nSgLmW7HviTNI6ZnmbbjDlzknHe2ZOz9fQiIll1Zo7m6wnNlA1mtoChbiIuvvjicT3HH7/zbWz9\n3kczWZaISMFJp6unAZiZsj0juS+dY9JpC4C7L3b3YncvLioqSqMsEREZj3SCvwyYbWazzGwycAuw\nfNgxy4HPJ0f3XAe0u/vRNNuKiEgOjdrV4+6DZnYnsBKYBNzv7jvN7I7k44uAFcCNQA3QA3zxdG2z\n8k5ERCQtWmxdRKQAjGWxdd25KyISMQp+EZGIUfCLiESMgl9EJGIU/CIiERPKUT1m1gQcGmfz84Hm\nDJaTKaprbFTX2KiusSnEui5x97Tufg1l8E+EmZWnO6Qpl1TX2KiusVFdYxP1utTVIyISMQp+EZGI\nKcTgXxx0AaegusZGdY2N6hqbSNdVcH38IiJyeoV4xi8iIqdR0MFvZv9sZm5m5wddC4CZ3Z1cjL7S\nzFaZ2TuDrgnAzH5sZnuStS0zs6lB1wRgZjeb2U4zS5hZoCMwzGyemVWbWY2Z3RVkLanM7H4zazSz\nqqBrSWVmM81snZntSv4MvxZ0TQBm9mYzKzWzbcm6fhh0TSeZ2SQze83Mnsv2axVs8JvZTOBjwOGg\na0nxY3d/j7vPBZ4Dvh90QUmrgXe7+3uAvcC3A67npCrgb4ANQRZhZpOAhcANwBzgVjObE2RNKR4A\n5gVdxAgGgX929znAdcBXQvI96wc+4u5XA3OBeck1RMLga8DuXLxQwQY/8DPgX4DQXMRw946UzbMJ\nSW3uvsrdB5ObmxlaKS1w7r7b3auDrgO4Fqhx91p3jwFLgfkB1wSAu28AWoKuYzh3P+ruW5NfdzIU\naNODrQp8SFdy803Jf4H/HprZDOCvgfty8XoFGfxmNh9ocPdtQdcynJn9u5nVAbcRnjP+VF8CXgi6\niJCZDtSlbNcTghDLF2Z2KXANsCXYSoYku1QqgUZgtbuHoa6fM3SimsjFi4VmsfWxMrM1wEUjPPQd\n4F8Z6ubJudPV5e7PuPt3gO+Y2beBO4EfhKGu5DHfYegj+iO5qCnduiR/mdk5wJPA14d94g2Mu8eB\nuclrWcvM7N3uHtg1EjP7BNDo7hVm9he5eM28DX53/6uR9pvZVcAsYJuZwVC3xVYzu9bdjwVV1wge\nYWjJypwE/2h1mdnfAZ8A/tJzOMZ3DN+vIDUAM1O2ZyT3yWmY2ZsYCv1H3P2poOsZzt3bzGwdQ9dI\ngrw4/kHgJjO7EXgzcK6ZPezun83WCxZcV4+773D3C9z9Une/lKGP5e/NReiPxsxmp2zOB/YEVUsq\nM5vH0MfMm9y9J+h6QqgMmG1ms8xsMnALsDzgmkLNhs66lgC73f2nQddzkpkVnRy1ZmZvAT5KwL+H\n7v5td5+RzKtbgJeyGfpQgMEfcv/HzKrMbDtDXVGhGOIG3ANMAVYnh5ouCrogADP7lJnVAx8Anjez\nlUHUkbzwfSewkqGLlE+4+84gahnOzB4DXgWuMLN6M7s96JqSPgh8DvhI8v+pyuQZbdDeAaxL/g6W\nMdTHn/Xhk2GjO3dFRCJGZ/wiIhGj4BcRiRgFv4hIxCj4RUQiRsEvIhIxCn4RkYhR8IuIRIyCX0Qk\nYv4fd5H0a8Rg6+IAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x26fa8c13978>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Task2 Perform accept-reject method\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import stats\n",
    "import numpy as np\n",
    "import random\n",
    "import math\n",
    "def gen_cauchy(n):\n",
    "    cur = 0\n",
    "    res = []\n",
    "    while cur < n:\n",
    "        x = np.random.randn()\n",
    "        y = 2.0 * np.random.randn() * math.exp(-x**2/2)/math.sqrt(2*math.pi)\n",
    "        if y < stats.cauchy.pdf(x):\n",
    "            res.append(x)\n",
    "            cur += 1\n",
    "    return res\n",
    "Y=gen_cauchy(10000)\n",
    "plt.hist(Y, 50, normed=1, alpha=1)\n",
    "plt.show()\n",
    "dist = stats.cauchy()\n",
    "x = np.linspace(-4, 4, 100)\n",
    "plt.fill_between(x, 0, dist.pdf(x)) \n",
    "plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}