{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.stats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's write a routine to compute the covariance between two variables. Remember, this is caluclated using: $\\sigma_{xy} = {\\sum (x_i-\\bar x)(y_i- \\bar y) \\over n-1}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cov(x,y) :\n",
    "    \"\"\"Compute the covariance between two sets of values\"\"\"\n",
    "    n=len(x)\n",
    "    return ((x*y).sum()-x.sum()*y.sum()/n ) / (n-1)    # add code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you were to make a 2D probability distribution function with independent variables, using a uniform distribution for each variable, what would you expect to get for the covariance? How about the full covariance matrix? What would you expect for a normal distribution?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>ANSWER HERE:</b> The variance for a uniform distribution (between 0 and 1) is, as derived in previous notebook, 1/12. If the variables are independent, then the covariances should be zero. If it is a normal distribution, the variances should be the Gaussian $\\sigma^2$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's generate a 2D probabilty distribution function with uncorrelated variables, and use your routine to caculate the covariance. Let's compare it to the [numpy.cov()](https://docs.scipy.org/doc/numpy/reference/generated/numpy.cov.html) routine, which calculates the covariance matrix. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "covariance:  0.0021789884310597827\n",
      "covariance matrix:\n",
      "[[1.06299978 0.00217899]\n",
      " [0.00217899 1.01959439]]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQ0AAAD4CAYAAAD2OrMWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2dbahm11XH/+t55j5JnplI2icB0zb3uQpSrAVtexGFoqXTlhDEqmhBb0tBy9ApYopfbB1QRAYsFWlQggyhUHMv+kVLUSOtwQpSbO2dksakTUssc6exxSa3H/oyeZnObD+cOb3nntkva7+ds8856webmXnmvOzzsv9nrbXX3puUUhAEQeAy67sCgiAMCxENQRC8ENEQBMELEQ1BELwQ0RAEwYsTfZz0zjvvVFtbW32cWhAEJhcvXnxOKXVX+/deRGNrawv7+/t9nFoQBCZEdKD7XdwTQRC8ENEQBMELEQ1BELwQ0RAEwQsRDUEQvBDREKbH3h6wtQXMZtWfe3t912hQ9NLlKgi9sbcHnDkDXLlS/fvgoPo3AOzs9FevARFtaRDRrUT0X0T0RSJ6koj+JEXFBCEL584dCUbNlSvV7wKLFO7JiwDerJT6aQA/A+BeIvq5BMcVhsRQTP7Ll/1+F24iWjRUxfdu/HPjRpGZfaZEbfIfHABKHZn8JQrH5qbf78JNJAmEEtGciB4D8C0A/6qU+pxmmzNEtE9E+88++2yK0wqlMCST//x5YLk8/ttyWf3uw1AsqxwopZIVAHcA+DSA19q2e8Mb3qCEEUGkVGVjHC9EfddMz+6uUut1Vb/1uvq37/7L5fFrXS79j1M4APaVpv2SSjxHKBH9MYDvK6X+3LTN9va2kgFrI2Jrq3JJ2qzXwKVLXdcmPxO5XiK6qJTabv+eovfkLiK648bfbwPwFgBPxR5XGBCpTP6hMPFgaoqYxt0APk1EjwP4PKqYxj8lOK4wFHZ2gAsXqi8tUfXnhQvjzXuYeDA1OrlLKfU4gNclqIswZHZ2xisSbc6fP54gBozbsmohaeRjYMqR/D6YmmXVQkRj6KTOkRAB4rGzUwU9r1+v/pyIYAAiGsMnZY5EqUlaImRFIaIxdFJG8kMEiNOgYxp9qUI2ZXTJG7mLJHclZL3WJ1at1/7H8k3S4iQ5xSZCpbw+wQsYkrvE0hg6KXMkfLsSOZaJaZv77+dZHyaL6eBA3JW+0ClJ7iKWRmJi06Kbx/GxCjiWiWmbdjGdx2RpjDyFuwRgsDRENITj+AgQx3XgNHqby6ETMnFXOsEkGuKeCMfx6UrkuEa6bUzoXJF2ToTPvkIWRDRKoc9uxdBzc5KcdNusVvrjmWInTSFbr/32FdKjMz9yF3FPWvQ51Dr03DFxlJjrnciw9BKAxDQKps9uxZBzp2i4saKTIvArWDGJRvL5NDjIfBotZrOq6bUhqkzy0s49kfkkpk62+TSEBPQ51Np0jtnMP3dCgpGTQESjBEITtFIET029G9eumdO1Jz6fxOTR+Sy5i8Q0NPj66SkDgru7Ss3n/NiGBCMnASQQWgApA3ipg6ch405SBiMluFkcIhp9k/rrnHoG8D57cELujYhMdkQ0UtF8WVerqqRKufYh9njtRnf2bH/dqL7XIu5RJ4hopMA1DiJ2cFdsXWITpM6e7Sdhy/feyHD5ThDRSAFn8JXpxc3xonf1Zeecd7UKP6ZvfYa2ONNAEdFIAWeYd8yENV1huw4XuutYLPzvh+uYtnvThaUhMRMRjSTEWBpKlfMimq6DyF0nn6HusdaLbducAlySwPdINtEAcA+q9Vu/DOBJAPe79hmsaMTENErB5kpwGjl3Up26xNwPm5Dk7PL1yVkZMTlF424Ar7/x99sBfBXAa2z7DFY0lArvPclZD+65ORPauNwJH0tjtYq7vtTBWp9zScykO/cEwCcAvNW2zaBFowRCzedY96o+N8faiLW6bC5UauuOK4RiaaQXDQBbAC4D+BHN/50BsA9gf3Nzs5urHiuhgUBXY+c2QI6FwbV8TFaDjxsU25i7EMEBkl00AJwCcBHAr7m2FUsjEleXo6kx2r6oKbpsfRqxySWoBcfHDYp1G0znms/NblApQe2MZBUNABsAPgng9znbD0I0Sn4pbJaGzXXxcWtcQUhbDIDTiG2iUMcu2ucwieV8Hh9w9XH3JtK7kjMQSgD+BsBHuPsULxqlvRQ+Kd8u14Ujhpzrd42MdZ3H5RLUx67/NF13qufj85GYSEZqTtF4IwAF4HEAj90o99n2KV40SnopfHsRONmSrgbismSavUftxC6TlbBYHO9psnX72gTBdxh/DiaSkSrJXT6U9FL4ChjH0nBZETYroL3vxsbN3c4cQdjYsGeS2q637+dji4GMyEUxiYbM3KWj65mpbDNw+U6t55oFjLOUouk65/Ob9716FTh16midFAA4PNTv397vpZfc29UcHBzdG1P9lOpm+YeQ2c7GhE5JcpfiLY0uYxquc4W4SiFdmW33RVcnkxXQ3Nc3zdy3mNyfLp5V+x737SZlBuKeeNJV70kKd8JFSIq07vo5AuabZh5SXPWpXYXcz65vNykzJtEQ98REe3lCIM8KaC73g7OKmY29vcpkPjioXulr127ehrvK/H333bw0YnvfLiYXPjiort9W52vXqus9OLjZZUi1mt1UJ1jWKUnuMghLo0lOd8X0tYwZt8E5vitxSRfwbAcuiSpXwbVv7VKYrBzfUgccOeNFclhurmsdSTAU4p5EkLMLdndX34uwsZHm5bO5CybT3ScuobsHJtfOZ2AY57w+QlO7DClFut1TxE2f74MAd1tEI4bcvmvMrFcuXI2rnQMR2hi5NF/elKLAFTfbeX3jREOxMgLrKqIRQ+5kr5yixPm6m9LPfRpjCCFB4JDStNpc429S1b0kAusqohFD7q9K7heQY0GExBtqseOYu8061OdaraoGbbuvIdaPTRh3d+3Xw2VIPSeBdRXRiCVnF2xXpi53jgqOUPjMa2GzGNrp5XXDbt9rWx3r7Vx1r+uYwh0US0NEozd0X+BcuQW6xhuSV+HK9Wg3eldaefPlNQkop6FzhKO+t+3g82KRJ6ZRwmhpiWmMiD6Cae2X2FcwXFZISByCk1G6WvFG3nKmM9zdvdktCumtcglCScFS6T0ZCSWYuDbhOHVK/7vJ0gjNw2h2d7oau22Oj7bFZrq3JqsldRZpCc83AhGN0kgVkAs9d934TMJQf325Q99jSu0a2GIXrrE2uqH4uiDr2bO8OqWwCIYULNUgolESLjM6dDAaZ1tf92G10h+j+bWezeKFYz43WwCu9VhcY1DqexrSExNjdYilMWLRSBms4hzL9uKmmGbO1ENgCyiaStst0HWTdlFs+EyYHBL4dT0Tm8tUSkwjABENE32MRQjNSuR8uVIlRNVFF4Dsupi+zD6WQ32M0OCvKV2eE5ztu/ckEBENEylNSO6xQs/J8ZFT94j4Wiapi25QnFLVb775JUrxB+PZ7nHsc+wD6T1JiO1F0WG7+dzAV6h1w3lJU85nMZvxj5dTXHRfb9/rbFtjujiNb4r5UAKdkqeRGFv3YRvXzff58nCUv72NbRbympSWRj15MbdRp7ZyTPfQ9zw+7qZPAxuKpSEZoYmxvWxtQgZYpY6PuNYy5XYpNuvcthJmsyOXgBMj8dk2tDQDsr77Neunu3c6ceaY8kMJdMrYk8RwhKB+gWwvZk2qwFfoV4zTqJoziOuWIajLqVNHprstYaopZvX9SC0aMQHZ+jmYRDim4XOfd58BUbE0EmP7WnC/nDnM0dAeltQN9sQJ3r3IOTfoYqHUyZNxx7BlgeZ+pn1bJCXGNAB8FMC3ADzB2b4o0VDK/BXw8edTE5rLYdovplHX9+TWW/MJQ2yxZbaGlJTBzBJiH6X1ngD4BQCvH6xomOAMx85ljnKyRlNkfE69dGFp2N6jgsnungDYGp1ohH4hXC6PTxq47YW3naPv/ApbsS2l0HU5fTq/62Cz/koLmjboXTQAnAGwD2B/c3Ozk4u2wu3yTJlP4RrerauT6ViuL2TO7s+xldXq5omAUmLLKymte7ZB76LRLL1ZGs1GyJ15KsTN8I0f1AKgq5Mpum86Vu2Ld7Fo0ZiKKfM09l2ocT2vAhHR4Pj6qVQ/5Ve+fjnbU+Wb3I96VGrfjXCIxeYuxPaAlBAM9cQkGtNZYU238HEb02pnvphWIlut/I9V1+n5549+OzwEvvMdYLE4vu1sVv3fwYH+WMslcPasfvFioWrGzYWwm3AWzrbhWph7SOiUxLcA+FsA3wRwFcAzAH7Htn0vlgZ3AZ5YdF+k2vQN6dlYr+0xkuaQdY7FopQ7a7QL1ybkHCnq5bpPOnfBFpD2cS8GNuIVk0/ucpnsuZdZrAXJZ3RmXafY0a3tF9t1L1LnPOgabsjsXynGt9hcu+ZzquF0fesYmEDoENEwWQD1g0/1UF0N3PXS6+rE8Ydtx2y/2DkFgVNCMjzr+TtjM0NtxXe1OVPPl26iombq/kBERERDqW7U39XAbVZGaNevyxrhDP2eemkPYONYQb7bm55foYhodEXK4fPt45oEz2Zuh8ZSYsrQunt93bf28woV44J7TpRSIhqdYmvgoUPeTcd3BfZOnkwz8a9PQ6jr3rcYhDZenzlHOdvbSsHWhohGKehyLkwJXDohKXVsiU8sxqeExDB8JlCuF0nium86QY+5zoLdFBGNEjBZGbap+0MbQ5+lOQYm1sqxzfdha9TcyYgWC15PTuhatdw6F4iIho2uuseGHoi85Rb/htCla9RcIc1HXDmD51zvxNmzR8eZzW6+btuyD4WmkptEYzoZoSb29oAzZ6osSqWqP8+cqX4PPd7WVpWdubV1/DipMk774sUX+dseHFTXf/26e9v5PLhKx7h27egZHh767RfD3h7wsY8dHef6deDEiSoDmAhYr4H3vOfmLOGazc2483eNTklyl6IsjZRjAkJ7TmJzD4bWWxHype+7NNebDXmHBjg8HgZLg6r/65bt7W21v7/f+Xm1zGbV42tDxPtKNtna0o/7WK+BS5eOrJrmGIblErjtNt6XkUhfV6EbTPee8w6ZtrEdt2eI6KJSarv9u7gnJtMwxGQ0uR/17zs7wIULlYjUZuuFC8C3v+0+9nxe7MtVLLV7kBvOO2TaZr1OX5/MiGikHH3IeXl2diqr4/r16k+g+gq5uH59kC9YbywWwAsvpBNa2whlzjvk+57ZYmN9o/NZcpeiYhpKpes90Q1Gc3XVcRdTNk3DL0VfUvbaNFe0r+Mv7dm+OMl5pjldObk4PawRC+lyTUz7oen6+l2zQXG7BU0L/vTdMKdQdIPPdCUkScs3b6cOrMZOCMRERCMlthGzpgetwyUUrpei7wY15hKy1CS3xy108KBrpHTiJDERjZT4PHBb4o7vS9OcdGfoiWKllfa9rQU6ZGV6GzEupmukdOIkMZNoSCDUl70983R6OkzB0b09XgC0ST2Vn1J+dRDcvOMdxwPUOzvV7z69aJxtOdNOrlb2oGnKHr8QdEqSuwzW0nB9JdpfAFNMQwKa5RXbbPScZ8WNKXBH0IaMlJ58TKPE6dJsLsFyWS28w4lFpHYtusyoJOp2PEmXpTl2RddT0bzXoWul2J69z3Gk96RFR0rqjWtRZm6AKmXPh3TF5im53rdS320NwxKNUteICJ3KL2RmKF/haI6ylJLuvuagRCtaw7BEo6PosDeppvLLYRmIpZG+9P2+9YxJNJL0nhDRvUT0FSJ6mog+EH3AvqPDJkxjR+pIuy5VmKhaPMl0nFRcuZJuiPlUcI1L6ft9KxWdkvgUAHMA/wPgxwEsAHwRwGts+ww2psEhJJW87y9qqpJzeYHUpXbpXPElk/vQdAfn8yoI7uNyDMBFQS73BMDPA/hk498fBPBB2z6D7T1po6ujLV5Rv2Q+s4n7lL5iGqvVMJdI8F20qoYzlaDv9IAFfhRzisavA3io8e93AfgrzXZnAOwD2N/c3OzosjNievDcF7a9bonPvrNZ3LlTltksLHtyaKUZl+KKsymQWmqgv0VO0fgNjWj8pW2fwSZ3NTE9eJ+vffMlaS9LYDP167VM2rkDfVkatQCWNumxawlGn9IMiobs08S1YFYh1rVJNFIEQp8BcE/j368C8I0Exy0b04Q7167xV2W/fPlo3oR3vav67eGHgQceqF6hNrMZcPo08Mgj1fbf+141b0Q9N2XsXJehXLkCvPOdvMmEuuTw0G+uUBtKVfN++kzq4xvQJ0o3V21OdEriUwCcAPA1AD+Go0DoT9n2GbWlwYltNL+EPkOjddvnKrXVMmaXI2fxjWmEjJLODHLmaQC4D8BXUfWinHNtPwrR4E6UYmvkQ+ltqF/osQpILfTcCZG4x3O9P9wRyz25KVlFw7eMQjSUcg8qGlpvgq2MNdvUtPJ76LOzzVpuwzWuqQfhENHoklwZn6UFGodeXDOrhQhHaON2vTM9uCkm0ZD5NHLAmTPBh9WKv8zBVPGdmwSomuMjj5j///z5KtDM5ezZo+xgX+osYRMFLbQlopEDzgNeLu0zXNesVsDzz4cLxqlTYfsNifkceNnLwvZ1PSuleMdZrYAHHwyrQ83OjnloQUEp7SIaOTA94Pn8+JiVBx5wd88eHoZbLcslcMstYfsOiWvXwkVVKfMSAefOAVevuo+xXFbPMgUpl9TIhc5nyV0mGdOwzQyVMlbRnkhmrD0eqYvu+djunWnCnlTvTwFJXpBAaMf4PnjO4CluaZ5Tgqf80g42xnaFFtL4QxHRGAK7u0otFukawcbGeKfmy1Haad8cS60erKd7lgMYlGbDJBoS0yiJz3wGeOmldMe7etV/Eesp8/KXH18KEQDe+1576vjhoT7dW9eDduVK9XsMJSzXqFOS3KV4SyPErExhio41gWoIZWPjZiuvOTO4a3/uPLAxs4F1bL1A3BMmIQ8m1cPsu+EMteR0wWoxcMU3uPPAxiRpdTykXkSDS8iDSfUwbZZGVwHNkyelx0UnBr4Zmzmsgo7nzjWJhsQ02piSfWxJQCH76HzTN71Jv+3p08Bzz1WvyMmT5mOm4OpVtx8/JZo5N7fdpt9Gl0fhmk82ti6c33OhU5LcZfKWhu4rtFjoR1mePn18vy6sAOlxUccsA9tQ9pRdqa64mMQ0CqWLmEbISuSyIFL3pR7M1kUsgfsOdZj7IaLhQ+7ek5CVyMc0zD53SWUpdblKe4HzhppEg6r/65bt7W21v7/f+XmLYWuLv+r7alXFM2az6jUSukUp8/Nar6sV5lNger5EveXaENFFpdR2+3cJhPaBblDSiRP6bb/73SpIWtAoxyiWS+AVr+i7FkcsFuZh9fXiU10MIislyMlBZ37kLsW7J13Qnn3cNdXcbBZndkvi2PF70XQjbdvqnleuQWqFpZ1DYhoF4xOvCOk9WSyGu0D0cpl2MF9dmg2/lHhCYQPcTKIhMY0S6CJesVpVrg5nbMvJk8D3v1/Vq8+xK+t15QLs7OTJG1ksgNtvr8aPEB1/BstlfF7FwJGYRsl04bceHvIHw73wQvXnPffYt8sJURVk3NmpYjo5ROOll44m71Hq6BwxiVglDCjLjc78yF3EPblByIzlRHrf99Zb87gHfaWUN12DLrubY1wSU9LealWMy+EDJI2cSVdfive9r1oljdv1WvPmN+vTkx96yNwDs1jw5iPVoVT3KeXtnokUk+pyryHmXLrh8LU1o1TZq6b5oFMSbkG1juuTAK4D2ObuV6yl0VUEOyYd3Lauhm6mrnqSGM5K59zz57JqmvVtYrM01uu0g/m465bogpbcZ1rYQs8mkKP3BMBPAng1gH8fhWjkjKI3XzJXL4br5QsRMduC1T6NbrmsxsPkEg3dvXaJuU+Kfe0qmLqvOaJhqg/3PmYalZqaLKLxw4OMRTRypQv7vNSc+Rt8RcyWi0DkP64ld9et6Rp03ZG+86DW945zftM5Tc/GZ63dAcQ3ehcNAGcA7APY39zc7OSivcllaXADefWKX74TvthwCUL9ZS0pj6O+B67gYej6q5z7q0u4q60b28fFR8QKnzM0WDQAPArgCU15e2ObcVgapiHQtqX7OHB8XaLK7Od8qZoi5koIcjWQOvGrJEujfb9MjauPQXzrtf3j4mu1rVZFJXQ16d3SaJZiRUMpffZh7BeB83JzZ8xy+fLtunKO6SsCsTGNkCCwztrrozvY5M7V9z1WyAqyPkQ0uORwUVL1XADHXyhX3CPXyvUxPRbzeXi9fK2oHKU5v4nOQkh5jp7J1XvyqwCeAfAigP8D8EnOfkWLRo5gaKqXmzvjdV3qcRsuc9nX0ggJnrbv4+6u/wC89lc4NKYRWjhWQArXrZDelayWhm8pWjRyWBopzGjdC8sRI5fFsbER5m6sVpVL5btf80sdsjCUbgLfptXje69t24dkc6YQpzFbGqGlaNHIkeAVa2mYXljO17751dI11LphxL7oHKthsXB3W3Iaue3ZpbQ8Uq5dUwuPK09nKjEN31K0aCiVfohyiCnPfXlc8YE6hqBU/1MGbmyk8f1Nz8R0fSdPhrlS7XvMeR+4x9PVx7TEY0+IaPRNs3Gbvsp1D0qIUNmEqRagPgWj3eBjXTbdbOCu/Amf49f5K6kmjda5VYV2tdaIaJQG56XxfbF2d8MCcbOZ/UvskyLNaeypLZ464Gu69voLzhWqpkWUYnmKglwOH0Q0hkZIvkiMNdGefrAWibohrlZpV7RPXTgBTdN+tgzUkN60AVgRHEQ0hoTtq9jsDWm+lKFdoKavpu54GxtHvQmrlVKnTvUvFimKjVKmAuwBEY0h4TLfY0ZYmr7C3EFZ3MlxXC5DScU1vkV3v7njY5rHGZj1IaIxJHKnR9ezcYcOyuLU0xWc7Lq4RqDaXL92g9clzLn2H2CcQ0RjSKQMFNriIjZrIsbS6HqqvtnMHm9prsmaYsoBX5dloC6OiMaQiIlPtBtLaJCP83U0jQquGzKnfq5tavfGZbE0R4vWgVyTKxA7VMB3/y6WdcyAiMbQiB1sVpvRIcPmTYOydMdr1tNnRG29vynmMZ8fr6tr7RNTA9TFEmK//GJpiGhEkTvAFSIctgFlzaxDrq9tmliG4+q4GolpBPCpU7yGbju26fpMY22486b4xigkpiGi8UO6eBlCci9sMYl2HV2i53KV6sbKnWhId30uV4Yzmtcnrdxk3fh8+UMS76T3ZIKi0X7wpm7O1GanT3dqbWW4GlqzB8Xm/7u+8LUQcC2i9jliA6XNsTVtQka8NhlgQ0+JiEYsPsHJ1AEu38CoUnGNsenCuBpeM/7BbaRNyyBGMFxWXYylMVCXIiUiGrH4NMIcAS5fKye2B4YzfV27Efks1Bw7cI3rEpgSs1yCMNDgZUpENGIJ+YrmhNslGjPQrG6YPsO42+JmO35o3VLEHlyux0C7SVMiohGLqQHknE2aE6TknDu0cTan5gu9xlzJXbljDGJpiGhE07WP6xow5tNoYlyA0Lo3E61yiEYf919iGiIa3nQZTed8obkvccjXPrSBpMpmzS1s3GuR3hMRjSLgvIw+AUUXvksoxEw718VYk3YRsmASjVnAQvNCDHt7wJkzwMFB9cofHFT/3ts7vt3mJu94ly/rz7G1BcxmwJ13Ag89xDvWfA7s7gLPPQfs7PD24dQnJ/N5t+cTRDQ659w54MqV479duVL93uT8eWC5dB+vLS5tUTo8BK5e5dXt+vVwsTDVJzfXrnV7PiFONIjow0T0FBE9TkQfJ6I7UlVstJi+xO3fd3aACxeA9RogAlYrYLE4vs1yWYlLE50ocUnR4Llip6N5ratV9ff12r1P07La2rrZahPSovNZuAXA2wCcuPH3DwH4EGe/Scc0QrrymklW7VGibbixkNj1am1xGW4iGvf6Tfesnt+TO+guZ1BzhEFT5A6EolqicY+z7WRFgzNaVLdPiin02/v7TlcXUydXjwpnwmTdvB31NbhEKHf36Ui7Z7sQjX8E8E7L/58BsA9gf3Nzs5OL7gTuF8Y3s7LG1zLRnSdkeUEbtjEdtuvnTpLj2r+5DydzM3ei1kgTwYJFA8CjAJ7QlLc3tjkH4OMAyHU8NSZLw+cLE/pilTiFvs0F6voLy7mvuVPCR5pyns3SAPBuAP8JYMndZzSi4SMEoS9WiV+xkAlxXIQKHUe4xdIIIotoALgXwJcA3OWz32hEw9RodEIQ+mKV6C+7YhTceSmaAd6YwGzIREIS03CSSzSeBvB1AI/dKH/N2W8UomEb1q0TgpgXq8TI/O5u3LwUul6PnF9q6T3xJnsg1KeMQjRs3YCcYOAYXiybELrm4uAsojTwmMDQMYmGZISGYkrSUsqcVbmzA1y6VGVeXroUn30ZS2xSVDsBbb2u/g0cZaWa4GRydp1dKrA40XcFBsvmpr5RuDIYS6FON6+zR+sxMICfmO3s3Lz91pY7K3U+twuHLttVKAKxNELRpUsP6UXnjoEJwTVobbmsBKp9/4iqP2uLpW9LTNAiohGKyTTv6kWPdS24Y2BCsLkV9X168MGb79/DD1fuXQmum2BGF+jIXUYRCO2TFF18OXMLRtoFOTUggdARkcK1yOle9W2FCVmhSlC6ZXt7W+3v73d+3tEwm1Xf7zZEVc8Ml729SmguX65civPnpWELP4SILiqlttu/i6UxREwxA98uytK6gHMgc20kR0RjiAy956YruFMrCl6IaAyRsccMUlkHObuVJ4zENISyaCedAZUVFSKKqWI/E0ViGsIwSGkdpIr9CMcQ0RDKImXSmcR+siCiIZRFSutg7LGfnhDREMoitXUwhW7ljhHREMpCrIPikaHxQnnohtsLxSCWhiAIXohoCILghYiGIAheiGhMARm0JSRERGPsyKCt7piIOEeJBhH9KRE9TkSPEdGniOgVqSomJEIGbXXDhMQ5asAaEf2IUuo7N/7+ewBeo5R6r2s/GbDWITJoqxu2tsyz01+61HVtkpBlwFotGDc4CaD7IbOCHRm01Q05J2oujOiYBhGdJ6KvA9gB8EfxVRKSIoO2umFC4uwUDSJ6lIie0JS3A4BS6pxS6h4AewB+13KcM0S0T0T7zz77bLorEOxIWnY3TInwaQYAAALJSURBVEick03CQ0RrAP+slHqta1uJaQijZGQTNWeJaRDRTzT++csAnoo5njACJtLtqGUiI2pjB6z9GRG9GsB1AAcAnD0nwohJtT6sUDQyR6iQjhF2O04ZmSNUyM+Euh2njIiGkI4JdTtOGRENIR0T6nacMiIaQjokJ2QSyHR/Qlpkqr7RI5aGIAheiGgIguCFiIYgCF6IaAiC4IWIhiAIXvSSRk5Ez6IaqzJm7gTwXN+V6JEpX/9Yrn2tlLqr/WMvojEFiGhfl7c/FaZ8/WO/dnFPBEHwQkRDEAQvRDTycaHvCvTMlK9/1NcuMQ1BELwQS0MQBC9ENARB8EJEIyNE9GEieurG0pUfJ6I7+q5TbojoXiL6ChE9TUQf6Ls+XUJE9xDRp4noy0T0JBHd33edciAxjYwQ0dsA/JtS6gdE9CEAUEr9Qc/VygYRzQF8FcBbATwD4PMAflMp9aVeK9YRRHQ3gLuVUl8gotsBXATwK2O7frE0MqKU+pRS6gc3/vlZAK/qsz4d8LMAnlZKfU0p9RKAvwPw9p7r1BlKqW8qpb5w4+/fBfBlAK/st1bpEdHojt8G8C99VyIzrwTw9ca/n8EIGw0HItoC8DoAn+u3JumRmbsiIaJHAfyo5r/OKaU+cWObcwB+gGrpyjFDmt8m5/8S0SkAfw/g/a1F0keBiEYkSqm32P6fiN4N4JcAnFbjDyA9A+Cexr9fBeAbPdWlF4hoA5Vg7Cml/qHv+uRAAqEZIaJ7AfwFgF9USo1+1WsiOoEqEHoawP+iCoT+llLqyV4r1hFERAA+BuDbSqn3912fXIhoZISIngZwC4DDGz99Vik16qUrieg+AB8BMAfwUaXUZNYvIKI3AvgPAP+NaqlSAPhDpdQj/dUqPSIagiB4Ib0ngiB4IaIhCIIXIhqCIHghoiEIghciGoIgeCGiIQiCFyIagiB48f9LJdqUGeZpcQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n=1000\n",
    "# Generate a set of n randomly distributed points (uniform or gaussian)\n",
    "x=np.random.normal(size=n)\n",
    "# generate an independent set of randomly distributed points\n",
    "y=np.random.normal(size=n)\n",
    "# plot them in a square plot\n",
    "plt.figure()\n",
    "axes=plt.gca()\n",
    "axes.set_aspect('equal')\n",
    "plt.plot(x,y,'ro')\n",
    "# calculate and print the covariance. Compare with numpy.cov\n",
    "c=cov(x,y)\n",
    "print('covariance: ',c)\n",
    "print('covariance matrix:')\n",
    "c_matrix=np.cov(x,y)\n",
    "print(c_matrix)\n",
    "if not np.isclose(c,c_matrix[0,1]) : print('Problem with covariance calculation?')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How do the values compare with your expectation? How does the numpy routine calculate covariance by default (as compared to how it calculates variance)? <Br>\n",
    " \n",
    " <B> ANSWER HERE: </B> We get the expected variances as the sample size tends towards infinity. At smaller sample sizes, there is noise around the expected values.  There does not seem to be a difference between my covariance routine and numpy.cov(); the default for np.cov() is to use dof=1 (unlike np.variance()!)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's generate a perfectly correlated (or anticorrelated) data set. What do you predict for the covariance matrix in this case, again, both for a uniform distribution and a normal distribution in the variables?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b> ANSWER HERE: </b> The covariance should be equal to the variances for perfectly correlated data, and negative of the variances for perfectly uncorrelated data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "covariance:  -0.09679730934183998\n",
      "covariance matrix:\n",
      "[[ 0.09679731 -0.09679731]\n",
      " [-0.09679731  0.09679731]]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQcAAAD4CAYAAADhGCPfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAWIUlEQVR4nO3dfYxcV33G8e/jOAZZBMiLE1JgvaAYCSeqCgwRaQW02KmCK+FQQRs0NhtIaxyXKhLij0juHxXIEgW1FNSsrS0FNrFbCKjUbgkK2YWUChFgrQCJHSUOqdekiRJjXiRn1Sbe/PrHncHryZ177+6dlzszz0dazczO8dyz6/jJveec+zuKCMzMWq3qdwfMrJocDmaWyuFgZqkcDmaWyuFgZqlW97sD7VxyySUxPj7e726YDbXDhw//PCLWpb1X2XAYHx9nbm6u390wG2qS5tu958sKM0vlcDCzVA4HM0vlcDCzVB0JB0nXSXpY0qOSbk15/0WSvtx4//uSxjtxXDPrntLhIOk84DbgncBG4H2SNrY0uwn4ZURcAXwa+Juyx+XAARgfh1WrkscDB0p/pJmd1Ykzh6uBRyPisYh4FvgSsLWlzVZguvH8q8AmSVrxEQ8cgB07YH4eIpLHbdtg8+YVf6SZnasT4fBK4GdLXj/e+F5qm4g4A/wauHjFR9y9GxYWXvj92VkHhFmHdCIc0s4AWotEFGmDpB2S5iTNnTx5sv0RT5xo/97sLOza1f59MyukE+HwOPDqJa9fBTzRro2k1cDLgF+0flBETEVELSJq69alruhMjI1l92jfPo9BmJXUiXD4IbBB0mskrQFuAA61tDkETDSevwf4VpQpQbVnT/b7Ecmlh5mtWOlwaIwhfBi4G3gIuDMijkj6mKR3NZr9E3CxpEeBjwAvmO5clnodNm3KbjM/71kMsxJU1RqStVotcm+82rw5GWPIsnYtTE0lgWJm55B0OCJqae8N9grJmRm4+WbImhVdWPAlhtkKDHY4AExOwh13wPr17dv4EsNs2QY/HCC5ZDh+PD8gvFDKrLDhCIemPXuSMYYss7Pwkpf4LMIsx3CFQ72eDD5mnUEAPPMMfOADDgizDMMVDlDsEgPguefgllt60iWzQTR84dCUt1AK4NQpnz2YtTG84VBkoRTA9u2+F8MsxfCGA5xdB5ElAvbu9SyGWYvhDgdI1kHs3w9r1mS3892cZucY/nCA5BLj85+H887Lbjc11Zv+mA2A0QgHSAJiejp7qfXioldSmjWMTjhAEhA7d2a3mZ/3IKUZoxYOkIxB5M1iRLhgjI280QsHODuLkTUG4YIxNuJGMxwgOYM4c8Z3c5q1Mbrh0LRnT/Yg5fx8UgbfAWEjxuHQHKR0wRizczgcwAVjzFI4HJpcMMbsHA6HVkULxjggbMg5HFoVLRgzO+tLDBtqDoc0RQvGTEw4IGxoORyy5BWMWVz0NKcNLYdDliIFYxYWfAZhQ8nhkGdmJj8gfAZhQ8jhUMTMTFIwJuteDC+UsiHjcCiqWQ8ia5rzxIne9cesyxwOy9Gc5mx3BhGRLMOW4IILfJlhA83hsFxFziAATp+GG290QNjAcjisxNKFUlk3bJ05k195yqyiSoWDpIsk3SPpWOPxwpQ2vyPpe5KOSPqJpD8tc8zKaC6Uev757HanT3uptQ2ksmcOtwKzEbEBmG28brUAvD8irgSuA/5e0stLHrda8qpae6m1DaCy4bAVmG48nwaub20QEY9ExLHG8yeAp4F1JY9bLTt25LfxNKcNmLLhcFlEPAnQeLw0q7Gkq4E1wE/bvL9D0pykuZMnT5bsWg9NTubvrOVpThswq/MaSJoBXpHy1rL+VyjpcuAOYCIiUi/UI2IKmAKo1WqxnM/vu8nJ5HHv3vT3x8Z61xezDsgNh4hoO5om6SlJl0fEk41//E+3afdS4OvAX0XEfSvubdVNTsIjjyRjDEutXVts12+zCil7WXEImGg8nwAOtjaQtAb4GnB7RHyl5PGqr7nUujnNuX59Mu1ZryeDkuPjsGqVS85Z5Sli5Wfvki4G7gTGgBPAeyPiF5JqwM6I+DNJ24AvAEeW/NEbI+JHWZ9dq9Vibm5uxX2rnAMHkoHLhYWz31u79mxwmPWBpMMRUUt9r0w4dNPQhcP4eFKDstX69cl6CbM+yAoHr5DslXazFa5qbRXlcOiVrNkKb95rFeRw6JW8qtYRyTSoA8IqwuHQK0WrWnt3b6sIh0MvFalq7d29rSIcDv2Qt3mvl1pbBTgc+qG5eW87XmptFZC7fNq6pHkvxr59yaVEk5daW0X4zKGflu7u3brU2qzPHA79trSi1PHj5waD78WwPnI4VFXzXoz5+eSyY34etm2DSy5xSFhPOByqavfuc2/Sajp1Cj74QQeEdZ3DoaqypjOffRZuuaV3fbGR5HCoqrzpzFOnetMPG1kOh6rKuxcDkhkOj0FYl3idQ1U1Zy3e//7svTGaYxBL/4xZB/jMocrqdbj9djj//Ox2HoOwLnA4VF29Dl/4Qv7dnKdO+fLCOsrhMAiK3M0JMDHhgLCOcTgMkj17si8xFheThVMOCOsAh8MgaV5iZN3uvbDgehDWEQ6HQVOvJzdrZU1zumitdYCnMgdRc8pyYiK5lEgzP392g19PcdoK+MxhUNXrMD2dfQbhSwwrweEwyIoUrXXJOVshh8Ogy5vmXLXK9SBsRRwOw6LdvRiLi2frQXia05bB4TAsll5iSHDeeS9ss7DghVJWmMNhmCwtOdfuZi0vlLKCHA7DKqsehGcxrACHw7DKqwfhWQzLUSocJF0k6R5JxxqPF2a0famk/5H0D2WOaQU1xyDSxh7AG+dYrrJnDrcCsxGxAZhtvG7n48B/ljyeLUe7hVLeOMcKKBsOW4HpxvNp4Pq0RpLeBFwGfLPk8Wy5Wmcxlm6c430xLINi6VZsy/3D0q8i4uVLXv8yIi5sabMK+BawHdgE1CLiw20+bwewA2BsbOxN8/PzK+6b5Wjui7G0/P3atd5xa8RIOhwRtbT3cs8cJM1IejDla2vB4+8C7oqIn+U1jIipiKhFRG3dunUFP95WJG1fDM9i2BK5d2VGxOZ270l6StLlEfGkpMuBp1OaXQO8VdIu4CXAGkmnIyJrfMK6rd1shWcxrKHsmMMhYKLxfAI42NogIuoRMRYR48BHgdsdDBXQbrbC92JYQ9lw+ARwraRjwLWN10iqSfpc2c5ZFxW5F2P7dti1q/d9s0ooNSDZTbVaLebm5vrdjeF24EAyxnDiRHK2kFY4RkoqT3mQciiVGpC0IVbkXowI36w1ohwOlshaMbm46EuMEeRwsMSePdlVrSNg3z6fQYwQh4Ml6nXYuTM/ILwOYmQ4HOysyclk8LHdzVrgdRAjxKXp7VzNWYnt25MzhVa+m3Nk+MzBXqjdJcbatbBli2/WGhEOB0vXvMRYejfnxERyC/j8vIvWjgAvgrLixseTQGi1fn2yXsIGjhdBWWf4Zq2R4nCw4toNRnqQcig5HKy4tJu1XHJuaDkcrLisknM2dLzOwZanXncYjAifOVhnuWjt0PCZg3VOa9Ha5joI8NnGAPKZg3WOi9YOFYeDdU679Q7z87B6tetBDBiHg3VOXsGYvXthc9ti5lYxDgfrnLzNewFmZz1IOSAcDtY5S9dBZPEYxEBwOFhnNYvWZhWMmZ/3NOcAcDhYdzSnMNvx7d6V53Cw7pichE2bstt4mrPSHA7WPTMzsH9/9hiEb/euLIeDdVdzDKJdQPh278pyOFhv+HbvgeNwsN7Iut17165kBaXklZQV4huvrHfSbvfetStZOdnUXEkJyaCm9Y3PHKy/pqaW933rGYeD9dfiYvvve6FUX5UKB0kXSbpH0rHG44Vt2o1J+qakhyQdlTRe5rg2RPJWUnqhVN+UPXO4FZiNiA3AbON1mtuBT0XE64GrgadLHteGRd5KyoUF2LbNg5R9UDYctgLTjefTwPWtDSRtBFZHxD0AEXE6IhZa29mImpyEm2/OPoOAZJDSAdFTZcPhsoh4EqDxeGlKm9cBv5L0r5Lul/QpSan/JUjaIWlO0tzJkydLds0GxuQknDmTfzfn3r2+xOih3HCQNCPpwZSvrQWPsRp4K/BR4M3Aa4Eb0xpGxFRE1CKitm7duoIfb0OjSD0Ij0H0TG44RMTmiLgq5esg8JSkywEaj2ljCY8D90fEYxFxBvg34I2d/CFsSDQXSmVpjkF4JqPryl5WHAImGs8ngIMpbX4IXCipeSrwDuBoyePasKrXkzGIPJ7J6Lqy4fAJ4FpJx4BrG6+RVJP0OYCIWCS5pJiV9AAg4B9LHteGWXOQMs/CAtxyS/f7M6IUEf3uQ6parRZzc3P97ob1U+s+GO3s3+99MVZI0uGIqKW95xWSVl2uSdlXDgertmY9iP3727eZn/f2e13gcLDBUK/DxRe3fz/Cg5Qd5nCwwfGZz+Svg1hYgA99qDf9GXIOBxscrQVj2nnmGS+17gCHgw2W5hjE889nD1S6HkRpDgcbXFn1J9vVibDCHA42uOr1ZJaiHc9glOJwsMGWNfjoGYxSHA422IrUg/DOWivicLDB16wHEdF+FsM7ay2bw8GGS7sdtLyz1rI5HGy4tNtZ64orvHHOMjkcbLik7ax1zTUwO3t2erO5cc6VV/a3rxXncLDhs3Sh1PHjcO+96e2OHoXNm3vYscHicLDhl7UganbW05xtOBxs+OWVvZ+YcECkcDjY8MvbOGdx0QulUjgcbPhNTsLGjdltvFDqBRwONhqOHIFNm7LbeKHUORwONjpmZpJyc+3GIJorLCVPc+JwsFFTr8P0dH5FqaNHRz4gHA42eopWlDp6dKRXUjocbDQtXSiVZe9euOCCkZzJcDiY5Tl9eiSnOh0OZnnTnDCSU50OB7MjR4oFxIhNdToczCAJiLzNeyNGqialw8GsqVlyLmsGY34etm8fiVkMh4PZUpOTyQzG/v3t98WISGYxhvx2b4eDWZrmVGfWWcTs7FCfQZQKB0kXSbpH0rHG44Vt2n1S0hFJD0n6rJT1GzerkLzak0O8s1bZM4dbgdmI2ADMNl6fQ9LvAr8H/DZwFfBm4O0lj2vWG3v2ZJ89LC4O7SBl2XDYCkw3nk8D16e0CeDFwBrgRcD5wFMlj2vWG/U67NyZ3WZIN84pGw6XRcSTAI3HS1sbRMT3gG8DTza+7o6Ih9I+TNIOSXOS5k6ePFmya2YdMjmZf7v3EC6Syg0HSTOSHkz52lrkAJKuAF4PvAp4JfAOSW9LaxsRUxFRi4jaunXrlvNzmHXXzEz+zlpDtkgqNxwiYnNEXJXydRB4StLlAI3Hp1M+4t3AfRFxOiJOA98A3tLJH8KsJ5o7a7Wb4ly1aqg27y17WXEImGg8nwAOprQ5Abxd0mpJ55MMRqZeVpgNhLSNcyAZnByizXvLhsMngGslHQOubbxGUk3S5xptvgr8FHgA+DHw44j495LHNeuf1noQaZcaQzAGoYjodx9S1Wq1mJub63c3zPKtWpWcMbSS8utF9JmkwxFRS3vPKyTNymq3UKpZk3JA9+Z0OJiV1W4Moqm5N+eA3YvhcDArq2hNygHbes/hYNYJRWtSDtAgpcPBrNOGZKGUw8Gs07L25sy7y7NCHA5mndbuXoy1a5PBywHhcDDrhubWe81ByvXrk0HLer3fPStsdb87YDa06vWBCoNWPnMws1QOBzNL5XAw67cDB5LbvCt2u7fHHMz66cCBZOpzYSF53bzdG/o+XuEzB7N+2r37bDA0VeR2b4eDWT+1WzFZgZWUDgezfmq3YrICKykdDmb9lHa799q1sGVL3wcpPSBp1k/NQcfdu5NLibGxJBimp/s+SOkycWZVMz6eBEKr9euT28I7yGXizAZJu8HI+fmelpxzOJhVTdZgZLPkXA8CwuFgVjV5NSkhCYguD1I6HMyqZmlNyixd3jjH4WBWRc2alFkl57q8ktLhYFZlWSXnIBmklLpS9t7hYFZlk5PJ7t55Zmfhyis7emiHg1nVTU4mJefyBimPHu3oLIbDwWwQFB2k3LevY4OUDgezQdEcpMwS0bFBSoeD2aBJK3u/VIdu9y4VDpLeK+mIpOclpa7PbrS7TtLDkh6VdGuZY5qNvJkZ2Lix/furVnXkbs6yZw4PAn8MfKddA0nnAbcB7wQ2Au+TlPGTmVmuI0eSWYy0jXsXF5PLi+bdnCsMiFLhEBEPRcTDOc2uBh6NiMci4lngS8DWMsc1M5JZjDvuOLtxTtqCqRILpXox5vBK4GdLXj/e+J6ZlbV0d+92O3yvcAwit9iLpBngFSlv7Y6IgwWOkXLeQ2oRCUk7gB0AYxUok2U2UMbG0utArPDfUm44RETZdZmPA69e8vpVwBNtjjUFTEFS7KXkcc1Gy54955a5h1Kb9/bisuKHwAZJr5G0BrgBONSD45qNlqULpTqweW/Zqcx3S3ocuAb4uqS7G9//LUl3AUTEGeDDwN3AQ8CdEXGkzHHNrI2lYxDHj5eqOVmqwGxEfA34Wsr3nwC2LHl9F3BXmWOZWW95haSZpXI4mFkqh4OZpXI4mFmqym5qI+kkkLKi4wUuAX7e5e6UUeX+Vblv4P6VVaR/6yNiXdoblQ2HoiTNtduxpwqq3L8q9w3cv7LK9s+XFWaWyuFgZqmGIRym+t2BHFXuX5X7Bu5fWaX6N/BjDmbWHcNw5mBmXeBwMLNUAxEOeQVqJb1I0pcb739f0njF+vcRSUcl/UTSrKSczQd6278l7d4jKbKKBferf5L+pPE7PCLpn6vUP0ljkr4t6f7G3/GWtM/pUt8+L+lpSQ+2eV+SPtvo+08kvbHwh0dEpb+A84CfAq8F1gA/Bja2tNkF7Gs8vwH4csX69wfA2sbzm6vWv0a7C0gKBd8H1KrUP2ADcD9wYeP1pRXr3xRwc+P5RuB4D/v3NuCNwINt3t8CfIOkIttbgO8X/exBOHMoUqB2KzDdeP5VYJOUVpa3P/2LiG9HRLM8z30k1bB6pWiB348DnwT+t4d9g2L9+3Pgtoj4JUBEPF2x/gXw0sbzl9Gm0lk3RMR3gF9kNNkK3B6J+4CXS7q8yGcPQjgUKVD7mzaRFJf5NXBxT3q3/AK6N5Ekea/k9k/SG4BXR8R/9LBfTUV+f68DXifpu5Luk3Rdz3pXrH9/DWxrFD66C/jL3nStkBUXeC5V7KVHihSoLVzEtguWU0B3G1AD3t7VHrUcNuV7v+mfpFXAp4Ebe9WhFkV+f6tJLi1+n+Ss678kXRURv+py36BY/94HfDEi/lbSNcAdjf61KQfdUyv+tzEIZw5FCtT+po2k1SSndlmnWp1UqICupM3AbuBdEfF/Peob5PfvAuAq4F5Jx0muSw/1cFCy6N/vwYh4LiL+G3iYJCyq0r+bgDsBIuJ7wItJbnqqgsIFnl+gVwMnJQZcVgOPAa/h7IDQlS1t/oJzByTvrFj/3kAyqLWhir+/lvb30tsBySK/v+uA6cbzS0hOky+uUP++AdzYeP76xj8+9fB3OE77Ack/4twByR8U/txe/QAlf/gtwCONf2C7G9/7GMn/hSFJ6q8AjwI/AF5bsf7NAE8BP2p8HapS/1ra9jQcCv7+BPwdcBR4ALihYv3bCHy3ERw/Av6wh337F+BJ4DmSs4SbgJ3AziW/u9safX9gOX+3Xj5tZqkGYczBzPrA4WBmqRwOZpbK4WBmqRwOZpbK4WBmqRwOZpbq/wFTBOKQKJbMLgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n=100\n",
    "# generate a randomly distributed data set in one variable\n",
    "x=np.random.uniform(size=n)\n",
    "# set the other variable to be THE SAME (i.e., not an independent random sample), or THE SAME with opposite sign\n",
    "y=-x            # add code here\n",
    "plt.figure()\n",
    "axes=plt.gca()\n",
    "axes.set_aspect('equal')\n",
    "plt.plot(x,y,'ro')\n",
    "# calculate and print the covariance matrix\n",
    "c=cov(x,y)\n",
    "print('covariance: ',c)\n",
    "c_matrix=np.cov(x,y)\n",
    "print('covariance matrix:')\n",
    "print(c_matrix)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's generate a correlated data set, but not with a perfect correlation. Do this by injecting some random noise on top of the underlying correlation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "covariance:  0.08522014496508686\n",
      "covariance matrix:\n",
      "[[0.08269301 0.08522014]\n",
      " [0.08522014 1.10109394]]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAFlCAYAAAD292MqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2df4xm11nfv+ednfF6ZuwEv95WgDPvQhOimqhSyYqmatWWTETDFsWC0ipoHFISdeOxaC2hqgKNRH+g+YNWbXEliLsqCYln+JG2oqRoKypvbCWKEmBMQnASgly647hBybIuYLIx6+ye/nHneu/cOb/vc3+ce78f6Whm3rnvveece+73POc5zzlXaa1BCCEkX2Z9Z4AQQkgzKOSEEJI5FHJCCMkcCjkhhGQOhZwQQjKHQk4IIZlzqo+L3nvvvfrs2bN9XJoQQrLl6aef/iOt9Zn6570I+dmzZ3FwcNDHpQkhJFuUUoemz+laIYSQzKGQE0JI5lDICSEkcyjkhBCSORRyQgjJHAo5IYRkDoWcEEIyR0zIlVJLSqlPKaV+TeqchBBC/Eha5I8A+Lzg+QghhAQgIuRKqfsA/D0A/1nifIQQMlj294GzZ4HZrPi5v993jsSW6P80gH8O4C6h8xFCyPDY3wcuXACuXy/+Pjws/gaAra3estXYIldKfS+Ar2itn/Ycd0EpdaCUOrh69WrTyxJCTAzQWhwVOzu3Rbzk+vXi8x6RcK38DQBvU0pdAfBLAN6slNqrH6S1vqi1Pqe1PnfmzInNuwjpnrGJXmktHh4CWt+2FnMv15B47rm4zzuisZBrrX9ca32f1vosgLcD+IjW+sHGOSOkTcYoegO1FkfFxkbc5x3BOHIyTcYoegO1FjujixHW7i6wunr8s9XV4vMeERVyrfVTWuvvlTwnIa0wRtEbqLXYCV2NsLa2gIsXgcUCUKr4efFirxOdAC1yMlXGKHoDtRajSbGsuxxhbW0BV64At24VP3sWcYBCTsZAyoM/FtGrMlBrMYpUy9o3whrbxHYdrXXn6Y1vfKMmRIS9Pa1XV7UuHvsira4Wn4d8d7HQWqniZ8h3SLssFsfvZZkWi/TvNWkjAwPAgTZoqir+1y3nzp3TfGcnEeHs2cJqq7NYFMNekhezWSG1dZQqXBk26gt1gGKEdfFi4V4ZSRtRSj2ttT5X/5yuFZI3Y5y0nDKpcxcut9IE2giFnOTNGCctp0yTuQvbJOQE2giFnOTNGCctp0wbE7YTaCMUcpI3Y4jUIMeRCu8rI1Xe8Q7gzjuB+by9NtJzVAwnOwkh48M1+SndyXd4LU52EkLCyT3uussFQgPY7oFCTsgQ6VNIx7ChWJeRKgOIiqGQEzI0+hbSAViYjbnnHvPnbUSquKJiOuqQKeSEDI2+hXQAFmYj9veBP/3Tk5+vrLQTqWKLijl/vrMOmUI+RnL3b/ZN3/XXt5DmHne9swO8/PLJz++6q51oJlvk1KVL3XXIpnX7bSfutdIiI9pXoheGUH+p+41IMYQ6iKG+Z46p7oDi/12ilHg+YNlrhUI+NvoWgdwZQv0NQUhz2VDMVFe2NJ93m7cW2pJNyOlaGRt9D8ub0LdLA5Cvv5QyDWGR0wD33D5GWa8PPnjSfTEUulxRalL3thMt8hYZgkWZwhCsUK1l6y+0TLlYv31SraP5XOvl5TArvO7S6Lquha8HulYmwlAEMZahdECS9RdSplzvVxeUIliKcKxw19PaWlzHCmi9tHT7ng3gnlDIc0GiB8/RwmthYigZqfoLKdNQOrAQumxXMb7velpe1vrUqbjvlOVxXbcq+j09YxTyHJiydZaToIUSUqYhdWAuum6bruiTEEGez+O/u7rq/17PbxyikOdAzmLW1EIZYycWUqZc7nnX+Yx1pdTrVcIVY+tge7xnFPJY+hg65WKd1ZES4RxdQj7qk3Tz+fHySdRdF/XWddsMscjLPJnKbPt+U4Ev67mn55RCHkNf1mEu1lmdNvM9FnF3tamQMtqO6aKt7u1pPZuZ73Fbsdmx8eGmCUtTvWxvuzuJ+dzvI6dFnomQ93WjcnUvtGWhdF0fMZ1GbAfje8u761yueohpqymdok9Q19baq7P68Wtr9nz4wjrro6HtbX/HCpijVmx1YupQhKGQx9CniyNHC7Stjq/LDjWm04jtYPb27AJUfrfezuoib6uH0LYaKj719hc7aShVZyZ817a1C591nvKs2SZUl5dPus8EmY6Qhwqh6zi6CuJoy3LuskONueexVrDLoi0tvlCRr9dDaF5c7oSqJZoa8idZZzZ817W1i64NDennoMI0hDxUUHzHtSVMubpOQmijg2r6ANbz5LLAYjqNmGN9AhoijDaxjwmF803yLRbpIX/SdWYjJDTQhKRBUG1TTTq2RKYh5E2tk+pxQxSmqdGk4wuxLqvnassidz3wLreJT/RDF6eEXkOpOHEqfc5t1JnrntqW5rvaRZNr1/3sKVsDlPdagGkIeWjP25cPvC3LoA0XzVBcQKn5CBXI8mFuy0fuE5Ht7TBruY2Jyvo1bHmdz5tHzEiNRn0TkaH1UL+2qX6lXE1AsdJU4BmahpBLWuR95s9H2y6aMbiAQq3Laie6vX1bHJaWir9txMzFuKIjYkYNscRY+7b8+K7fZqSPJL5Ri6ncKatDXWl9vXExpiHkUj7ytpC6btsdUa4uoOrD6ptErJepzTZhExGf0Crl7kx8hHRmS0t+y7QN2rxO7Lmbzg24Ys9N97RBeVsTcgCnAfwmgN8B8FkA/8r3ncFHrbRJ7HVNx6e4aGKum+MK05RhcKqPXIoQoW1y/VDfeNeYfN3Ly3KdZmyH3GS1Z/ld15yBr+1F0KaQKwDrR78vA/gNAG9yfWfwceRDIXbIZ3voYxt316Im0ana8mxbkViP1e5j/kJCaGNdBl12VDZs7VdilWjMs7G35xbfGEsb0HplJW4yNKHuO3GtAFgF8NsA/prrOAp5IK7JpzaFuUvXk9S1XELsizxI6RwlytNUaEMn8WxuprJuYsqW0uHWv+cqbxO2t+3nNS2QWlmxH1+dNwh108WmBCOhVSEHsATg0wD+DMBPWY65AOAAwMHGxkZ0ASZJqDi14SrpyvUkZf03neg2dY4pq/RSOs3yO/X75OvQXGWp4nIdmCI1mkao1MsXY9VW6yQkX+XxsTHmrg6lnOgur9eGiKe0cd2ykL9yMuDVAJ4E8AbXcbTIA7E1tqWldsVFGtdDKOXSCBWaGMu9brGFCFeT8sR2nr4Y9RLb/Q8RqaYbRcVMJM5mhYDWR0yzmf1ehHYUoe0gpm6apKH5yE+cEPgXAP6Z6xgKeSAhjbTu6w09T1fhhL5rS3YyIULoul5I1EuqcIWWJ0bMXSJZvZ7pHsSuTPR1ULZ8t23RhnYU9bpsGqliSy4feawBZqDNyc4zAF599PudAD4G4Htd35m0kDeJWgnZl0MiSkfSreITtq47Gdv1TLvhuYQr9vwh5Yn97t5eWD7rrofY+Gjfni6uOrW12fm8ue85djVqNa2txb8Orrym6/+2BV5CbbpNIf8rAD4F4DMAngHwE77vTFbIXQ9qiHi2Ha7my2MKoULTZSio6Xoxll3K+UNIseZ9k7Sm+7myEm+Ru9pFiusmZEVrSL6aWNazmdZ33BFfD7Z826JcBLe3ncaCoKHTNAqli7hgaVeHrdH35Z+3bZwVIiqlH72tjid1UrruP15Z8buuQlPIni4pgtxUhMtr2sJL20queyQV+eSAQj4EYht8vQF0ERcsGU/tstT68s/bhCpUVKRGLDGjAld8tSnkrlxc43K9hIhV08l03zVMdbm8nObyGELqYCEdhXwIxDZ4UwNwPfBVS6x6XOjmQq481ifPQixSV8fVBbETYSGdpMSIxeVTNsU221Y9ukY8pkiPkPaWMsJImUytun7qbcnlWw81hmIX5zRNrraxtiY2iqOQh9C2r9b2AKes1PStSEvZ7tOVx5R44hTRk7wHoQ99PfLCJgwu4fRFb4TWS4jPO3TyO1aIXNjKZTIY5nP3a9mAuAir6gghtCzliCS2M0tNsSOgoYcfhqRBCnlX0ROmB8Ln5/TlU+qhDbHiY8Q5pk5tnVMXu//VQw9NERWlVZwavVG957Z8uKIwbK4IiZTauW9umvPr62BcriJf+/Ld03qIn3Rdue5N7PUS3KDTEfIuIwak2NsL30RIIv7V5rJpuqDGVjbf/fBFMNSFNvS+hvrIQ0MPSxEyLU6RcAe4UvmSYwnhWVmJW7Eqdd3y2tXr1e+rr83GdGZtdHqua6WEdUYyHiF3PdCpVnVoPG5bxHQiEqJgOm+oeye1w3MNzVPKFBOb7YtaiRWq+tL9mI4gNc1mcotr1tfjhFxyUU/VGo/xrZvmaIDj+8e3Wf9tpMla5D6hTvXJuh7C0Mpu4tuNsXIlw8mqebcdX8+DbWOi7W23WMfGILf0IETVf+i1JS3WlFSG4NXLUXY4pv/52kQVyfJV25Mrqqn+t29f9raX1EunxMitcQi5T6hTwn9cjbStFXl1Yrb1bOInjZ18q9at71hXLLzrvjV5AKVGS6lC1SSOWjKV7cTWkYaUz/Y6t6ZtztWeXPVm+p9rYU3KPezLig/pmCyMQ8h9Qp1ikYeESUmE7LmI3Z851ZKtf9/nlwTCfeSu+nN9p8nD5Ju0jdkGIUWoUkcV0qFx9Q7N5HpILVv1nDHnWl/3n9N0TEi+TO6x2HvY5DVuTTruhpFY4xByn2CmWMYhD6HvHE0XAqR+P7ZBhQyz68en1Fe9DKFD6NgHOjSqxyRKJiEo8xm7fD1lufnaWpqQ2VJp4UlHtpRvvokpn2+EUP4vNU/1vNi2ubCFQMaWRyo1sMRLxiHkqQ9p7DltD6yNphZ56vfb9M3afHixQmELywt5kJaXzXHAKUPsmM7eNEHqymeqcK6s3LYuJe6Z5LlSU+jzKJ1P07OyuTm8emq4qnkcQq51s0nFkHO6boDr+yahCu19m0TbmL4n1ehc1w0ZZrv26fA19rql3DTCwud+c/mIXd9r6mf1ze+0nRYLmbfFx0aRSZfX5F6yHbu+Hm6QtHFf6i+8jmA8Qq51O2Jekmodm4bXMROeqWUyfU/C2vCV19fAfR1Z01FMbNljBdNkWZoWbUk82FJimnqPJdwxXXZ+Ie0mZO5HYgfG1DT5lZ22B0pKzFOt45Rl9hKdkc3f2+TBDBlNhHQWLsujaaSP7fummO6QEFXfvdvbMy/akhJg0wRom/uF1P3K5f2SFKeuRDLFRbW05L53XeQ9IXx2PEJumyCSeAN3SYqf3dUgTMdL7aAXsr+5a+8VW6N1TSbGTAr6XowQWs/1Y317zZgWvOzt+fcAsd07m0j49rN2Lbs35VtqEtaVH6kO35aq77yM+Y50PoaeJruy0zfp1DYpEzfV4Wv53dTXiNWJcU+kuFzq50mdtAwtl61+UwWnOimaco6qcZC6aCg2OsPVqYXM5cTUSdsTfpub4e6nLlbHSiUpa32SFnnIBEmKmyLUKkyduIm1fGLynxK2GNtgq7gmCl3lC7E8XPXbRHCanqN0MaV8f3Mz3dKsum1MG5iZXIyxdZJap6FpNgt3DS0vH48PH7qF3lTMJ+sjj3mQ2liR6bJ+XQIXm/eY/PuuaypvbIx0FVfH4YpkCXnprKt+mz40Tc5Rli1F+Nrws9b3nI9xFdXbiHTemqb6nEQuFnpsms0mHLWSsqLQR4xrImQ0UP27KsZN9/KwYZqAA+yTv007w5QYbdN5TSvzXJ2E7bqzWZgghYQ8+u5FGyFzqanaUfcVD10fLUilevx++XfXr3Tr6v5Fkr+QxzbYkOF8jGsi5Prl+eqWZ2oYVoibJeZFBL78+yxn14ZZJSFzAaZOz7UKz9VBlB2Dy81QlqdJBEib1mGKdexrw22n6jbLtnYRO1qQNoaGnBLJX8hjh1pNLHJb9EDqClBXjLmvg7BZsCU+d0dMnfk6v9jYb4mHr/r+SVvH4BKMahSPy6pzLdt2dbjr6+brxy5pj+1kmvj9U/bOttV7tY1Xt5Qtd8Ns2gGWz2LTvA4tTXrTrNDG28RHbloWbgrps13btMLMteozpLG7LJUU370tpS4CsnUAUg9gGY8e2zFU49hTF4iErv7b3DzpFggVsTJENKZsvrqov8C4Olrc3pbzkftCRofoix9CWlpyP2sWpiPkMVEf9YZWxh/bzlsl1DoNOa5JaKBrgVSM8JXHrq3dtlxLy6rE9VCa6l1ywirVigx1K1Vj7usx6KF+eFP5Q2LA29rE6fTp27+vrcVvmhZTx9WyT0XAm/rtE8hfyENEISYu0xaxYjt3iKVtGglI++Hr37etOAzt9MpzuY4rxTz04ZzP29mYKMUFkZKq99G3diGm/dnOFTppO+RUjjK7jjZRSuv77+++vGXn1aSNM/zQckNjKiZ2AtK1yMY1QRi7aCcm3tx3flfIXPn9EAEph4Fjm3AKeVBToqVc7dDmz19bG0f9dlmGWBeWVKqGDzbpuCa5ICikgcTViD1JLJ8viYlVL48PifrwxTa7YrtL10tMbLTW6THLuaYmk4m2+9t3mWLTkO+51v1OhJauldR95Se5RD/ETRCDa/FK6GrPUOrnq0ah2PYF0dpd3pCd22xlLEP6+n4Q+0gp8wYpyRbK6fqO601RfdTT9naca2nI93KIKSGePH8h98USx4qtq4JD85Mi9r6hWDVCxtZQJSbHpLZgHXNK7Sht99T2ooMy2dwrfQhW2Z7bWipfnVDv+z73lSYp5FqbrdBqKF+JtO+6fl7Tw+XbLTBk175qHmz5SwlVq6cpPzx9pZwmMyUm8kJSuf1s3+XtK03StaK1e6l2+fvamj0OvEqK7zo2aqbJREhbVthY968ISezAwtLm5rTbSVdpshZ50+07y46gfKBd/ukqvlWBZaqGBA5xNVpMXHS5wKXNRtx1+U2dfHnf+r43sUkp2Zc3514fOaYhCTmA1wB4EsDnAXwWwCO+74hb5E2Sz79uWnBjS6VQ9mXNuPJZnXwLiUSwvW1H6iG3vVhZKp06ZS7nGDZf8oWWMuWRhuRaAfCNAL7j6Pe7APw+gPtd32m0srMNa8HlFw/tPMqQvr4tcen6WV8/OWqROnd9T5uQSJxv+qZ+67fvVO7bnZPPncmcBOPIZ2iI1voPtda/ffT7i0eW+Tc3Pa/jgvLnfO65tP9VuesuYGvLf7xS4flKQbp+/uzPgBdfBB5/HLhypfhMqgzXrhX1tbEB7O4Cly758/+lLwErKzLXz5GXXwYee6yoO5I3u7ty5zKpe2oCcBbAcwDuNvzvAoADAAcbGxvRPVGrLgtbz+jaca+eymGSyyLvayWaZD21ZQnmWidMTClpfT1eAwshbXeyE8A6gKcBfL/v2E72Iw9NrrDBlG1zQ6JhYvYIn1JiVAnTVNIQ3xAEYBnArwP40ZDjxZfop75Q1hWpEnMel1D7Fgu5riMl9Dl1GKGbYtFHPPzkemHI1FMirQk5AAXggwB+OvQ7rVjk1RWRIZa0L1rF98IGk1CnrPZ0vXVnb08mxCynh+nUKX/Hs7kZF00knZTKq3PsK/ne3DTlNDSLHMDfBKABfAbAp4/Sedd3Wt/GNjSyIiVaxeVPT9kcy5W/2aw7d0Nbe2JLp2p9VPdO7yqVo78c6qrv+9R3BNeQU0LEitZat+4jj0mtvViiGpcZ2ohc17IJs8nyjn2BxJCEoNymtu985JCm8tKEpsm3rwxTkgSOQ8hLXHuJx+wT7HvdkkmwU19I0edCodBGNSaBqr/qTLKepM85n9MFMcWUgE3IT8kFMnbIa18LHB6e/PzmTeDCheL3ra3i5yOP2GNub94sfu7vAzs7t2Oaz58vYprLvx9/vDhuZ8d83evXgaWl2+ercs89wNmz5u+lolTRFKRYLIo6ePFF/7HzeXHcjRty15dmaQl41avkY61ns3bWAVy7Bpw+LX9eMmz292/rVFNM6t52amSRh/gnq+4Mn4sl5PVhKyth0RR1i7tN/63kuV3vKq2Xz7Ud65DcRUPKCxOTKQnutdJ4ZWfn7OwU1eCiurrSZwlfu1aslnNx44b/mMUCuHix+KlUYbneuuX+ThMkz33tWpj1+tJLwIMPApcvm//vuy9doRSwutp3LghxIzhizE/IQ5bM33PP7d+XltrLS8nqarHcdmurWMb++OPAH/9x+9ftmjY7Jkm0LtxdhEyE/IR8YyPueJPfWpLSEi99Xfv7hZ++7esSN0MZHRBiYyYnv/kJ+e6uf9j8wgu3f18s2svL2loxQtjZKQQcKH6nNUgI8SE4ws1PyLe2bvuibcxmRTp7tohAWV5uJy9f/Wph+R0eFlb4/n74bomEkGkjaGTmJ+RAIebnz9tDwW7evC2wjz1WTFS2vX3s9euFNR7r+mmL5WVgfb3vXBBCbLz2tWKnylPIH34YeO97w/yg5TFd+EwPD8NcP13w8svAHXcAp/JcKkDI6HnySbFT5SnkFy/2nQM7Dz4I3Hln4T/vm2vX2h+JEELSEPSR52muDT0iZEgC6ot/J4RkT34W+bd/e985CIPhb4QQF4Kj9ryE/C1vAT73ub5zQQghzXnTm8ROlZeQ25aGE0JIbjz1lNip8hJyQggZC4JzfRRyQgjpA8F9oPIS8s3NvnNwku3tvnNACMmR179e7FR5CfkTTwD33993Lm6zsgJ86EPy5x1K6CLJi/m87xyQGD7/ebFT5SXk+/vAF77Q7ByuHcdiw4G+/nX5t9AADF0k8cznwNe+1ncuSAyCz3k+Qt50e1ilgL094Lu+y37MSy/FnTOX/bnJuFGqaLvcdXOy5CPkUtvDfuQj9v8NfcUoISa0LnbiJHkhuKldPkLedHvYjY2w18QRQkgXvOMdYqfKR8ibbg+7tib7JntCCGnCpUtip8pHyM+fb/Z9Lu0nhAwJwZfQ5CPkbYT5EUJIXwi+hCYPId/fbyfMjxDSLlwTYWd3V+xUeQj5zk7fOTCzuclFGIS4YHCBna0tsVPl8WKJob7Q+PJlWhyEkHiE3yCWh0U+lBcam6DFQQiJ5fRp0dPlIeS7u+6l9YQQkhMvvCB6OhF1VEq9Tyn1FaXUMxLnO4GgL4kQQnpH2MsgZeb+PIC3Cp3LDPc1IYSMhabrYmqICLnW+qMAZMcKhBAyVn7u54qwaiE6czwrpS4opQ6UUgdXr17t6rKEEDI8btwAHnlE7HSdCbnW+qLW+pzW+tyZM2fivry/zzA/Qsi4EFzkmEcoCHctJIQQK3kI+VAXBBFCSCorK2Knkgo//EUAnwDweqXU80qpd0uc9xWGvCCIEEJSuHFD7FRSUSs/qLX+Rq31stb6Pq31z0mc9xUEN5chhJCxkYdrZWtLfG8CQggZC3kIOSC+N0EQjJQhhLTFJN/ZKbw3gZPTp4somTZXk7KTIGS6zGbAY4/JnU7sTG3T5YTnSy8B994LPPxwe9dgOCUhRIh8hLzrCc9r14D3vrfbaxJCpsGtW3mu7GzM1hZw//1954IQQmSY3MrOki9/ue8cEELI4MhLyPkCZkIIOUFeQk4IIeQEeQm54N4EhBAyFvIR8re8RXRvAkIIGQt5CPn+PnD5ct+5IISQQZKHkO/s9J0DQgiRZSYnv3kI+eFh3zkghBBZ3vMesVPlIeRLS33ngBBC5NjcBH72Z8VOl4eQ37zZdw4IIUSOT3yimPsTIg8hXyz6zgEhhMhx/bro3F8eQn7+fN85IIQQWQTfRZyHkF+61HcOCCFEFsGtufMQcsGeixBCBoHg1tx5CPk99/SdA0IIkWVrS+xUeQg5IYQQK3kIeZfv6ySEkLaZz0VPl4eQd/m+TkIIaZtHHxU9XR5CvrvLt84TQoiFPIR8awt485v7zgUhhMhw4cIEV3bu7wNPPtl3LgghRIZJrux85BHg1q2+c0GGCDdUI7kyuZWdfOkyscEN1UiuTG5lJyGEjAmlJriyUzjmkhBCeuXNbx7eyk6l1FuVUl9QSj2rlPoxiXMeQzjmkhBCemVo+5ErpZYA/AyA7wFwP4AfVErd3/S8x9jaolVOCBkPA4xa+U4Az2qt/0BrfQPALwF4QOC8x+EyfULImBhY1Mo3A/hi5e/njz47hlLqglLqQCl1cPXq1fircAdEQsiYGFjUimntvD7xgdYXtdbntNbnzpw5I3BZQgjJlNXVwUWtPA/gNZW/7wPwJYHzHoeuFULIWLh4cXBRK78F4HVKqW9RSq0AeDuADwuc9zh0rRBCxsBiISriAHCq6Qm01l9XSv0IgF8HsATgfVrrzzbOGSGEjBFBl0qJSBy51vqS1vrbtNZ/SWstn0uAy/QJIePg4x8XP2UeKzsBbo5ECBkH730vsL4+rAVBncHNkQghY+GrXwV++IfFxDwfIV8s+s4BIYTI8fLLYqs78xHy3V1glk92CSHEi9DqznyU8eMf58slCCHjQmh1Zz5CfvFi3zkghEjhG12vrnaTj74RCkXMR8g52UnIeHCNrpUC3vnOfOfF5vOiDIuFe9fW+VxsYVA+Qk4ImQZaAx/4AHD+fH6W+dJSsZ3IxkZhbT/6KLC8fPK4lRXR9yxQyAkhw+P6deBDHypcqjlZ5jdvFh3R4SHwjncUc3vvf/9xy3w+B973vsHttdINa2t954AQ0iXXrhUiuLtbuCqa0EfEm9bF4p/3vKf4u3S3PPro8PZa6YT9/SKAnhAyLS5fBj72sUIUm9BnxNtXv3pbvw4PgQsXit8nZ5ELvhKJEJIZN270nQNZhF/zBuQi5IKvRCKEkN45PBQ9XR5Czr3ICSE2lpYK/3NOL2hXaqKbZhFCuAtoHaWKSJF77gFefLHv3ISjtah7JQ8h52veyJRZXga2t4uIB25TcZxyEvTatfx86YIu4zyEXPBt04Rkx8svA489VvhVm0ZvkOEgqGt5CLlEHCkhOUMBHxerq6KvfMtDyLe2gIce6jsXhJDcGYJBuFgUK1YF48iV7qGnP3funD44OIj/4hBuAiEkb1ZWignSvjbia6C5Sqmntdbn6p/nYZETIsXKynhfULK5WSTi5saNcBHf3JRtL0oV5zt7dqLhh4KFJhPmxg3gG76h71zIsrRURLU8+2yxpJ3IcfmybKSQ1rc31bpwYSY8QFcAABw0SURBVILv7OQyfSLFtWt950CWmzdvR7U0ZT4H9vaanyeEMiY+xGU6xlGU4FL9fGpHeEkrIaNCaq7r2jXgh36oG+G8dauY+AvJey7x88vLhfsulMm9s5Mr2gjphlu3uhHO2Wx8BtrddwPvfnf4dgGTe2cnX/VGxgANktvk+EyvrJjf+FNy7VrxdqNHH/WLuWAseT5CntNbQggxUUYskDyZzYo3+7z//YUeKWXumK9fL9456pqLEY4lz6dV7e66e0JCho7WxXL7KeIzxGI7uBg/tAQrK8AHP1gI79YWcOWK2wXlGm0sFsX3J/diCaAo9N13950LQkgKV664xTzGJ18X8bZHOYuF/R2bKT7uFt6vkI+QA9wFkZxkyqt9F4t8XI77+3J7Jt24cXynwzYnZpUCzp+3W8+7u4WvO4YWNgHMS8in8oKJMfhRu3hZ9t4e8Pjj05xALCfKUoSkD971ruJnbnsmlS9Qfvhh8/+3tgpft8tnXkV4s6xKPnVyAvAPAHwWwC0A50K/98Y3vlFHs7en9fJyuS6KaahpdbW4V4tFN9dbWtJ6c1Nrpfor73zezbVms+LnYlHUcfXZkL5WG/W5WBT57aq+pOujWucunVpdNddl/b4lAOBA65Oa2tT0ewbA9wP4aMPz+NnZme5EUU7ceWfxs6v3rN68WSyj1rqb69W58872VorWrbvTp4tRiPBEmZGY+tzbK473rQg9PCz2GBnqylqXNa11+CrM8hkAihDExx8vvt/mfTOpe2wC8BTatsilLYTSumGST11aqUNI0m1zaclef6VVW7K3199oBCjyWWVpSaa+Tp/Wen29mzKUo0hfXSpl16e9Pa1XVuznFgItWeTdIT1BkMuS31T69LNfvz6tiWmtZc9365bdaj08vL3R0v5+sfGS7/rlplrSLydeXS0WvlQJWeQTUl8vvQTccUdavkpCJlbX1m7Hc/vee1DXoP39YoShFPDgg+ZXzQnup+LEpO7VBOAJFC6UenqgcsxT8FjkAC4AOABwsLGxEd8VteEHZGIaYlos3Jbtykr8PMT6etzxthHrfF5YplV/b5kXpcIs8uq52qzHkHzVRzhaa729fdIyr1vWJl+4Lbks+UhgschPfJCSQoS8mpJcK1q3e9OZmIaQSsHwHVcKVJd5KzuQkr29dDEuBbTppLitDuoCbTvOJrLVTsA0SRmTb1NnkQiFnGlcaYxzHFXB8Fm2pcD0kUet/Rapz9dcnmt7O9yyDb2WyS9tq6tUkY3JXwc+8hMfxCQA3wfgeQB/DuDLAH495HvJQt7npA7TcFJ1gm1trf/8AM3bZt0y9B1fin5TEUxJoR1I2dlsbto7ptXVQsybullMLp8qpsnI+ggjhlAX0uZm2vkt2IS80YyY1vpXtNb3aa3v0Fr/Ra31321yPif7+9NexUcKlComAstXZZ0+3XeOij2AHnro9qKQxSL+FWH1iTTXis2VlWJRSX0xynwOrK+HXzNlIZVS4VvPal0ce/myfRL0+nXg0qW4fJv42teKMD9XiJ/W7r9j8E3qzmbFBPMTT6RfIwaTuredkizyECuAFvu0Uh/WqCnNZsctwhR3wXweNpm2vm63OE3Hz+fm/FRD7vquR6Vknl2Xm0TatSJ9vkDQpo88Ng0ijpyJqa2U2lZLN0M5yTafh7kMXEZO1addHle6Bcpzbm/3W18uN0XMXIgrOsR1T1yTmjZMHaArZtw3eRpI/kLex8QOU7O0WExrYVAbyScOIda0a6QgtXirFCfJ5zR2pFDvtKqiactXXeBjltOHinOs6DvIX8iHMARkCk+lddl3PsaQlpbMYhErmm2Mak3x1b49kcp8uCzxkGOqqZy4tImmqRPz1YfUqkxBN0z+Qq51s7jVMaZy6N13PoAieqTuBpjyKKotV2BVXIbgbjT59k1L1U3HS+ff9TxU22RM6GZTn7drTUDCQqHxCHnMyrGxpyE90MBJC6btfK2vD6fs9XpoczRSWuhSz8J83mxn0erEqS9PVWHssqM3iWZoAEUqPi/CJC3ysbpWmgjR0tLwLN+qn9L1UM/n8gtBpNLKyklhc1mZ9bS9PVyjw+YTns+Px+SX0S6hZQi9n1VhjH2mY+5BPZlEM+T6TSxy13M5WR/5kMRKKtlCw2KSzf/XVyr3bXY9dMvL4Rac7yFLaRfb2/b8lUNw00SWaQ8O230dyv2o13sZFVPeq3pbqouLtAFl2r0xpA0sLR3Pe0wKiSYJrY8YXG1lslErQxxCSzRqCb9/eZ42t/0sw+F8x/n89rNZ8UA2va9N4qCbTD4NbQQUk6orYn11YBO4JlsjlB24qU5D7mH1nrtcQfO5P5rE1FELhQgG13EC+Qt5rg9PF6k6XI0ZDpvSyop5+F26C0K+H5LfskGndGL1hyzWBRC7gVLJGCbbS1x14BLW5eWT9zg0hNG1JL4qoq57VE6W2tpZiBUtGA7oFH/J6xyRv5DzVW/2VO/hm060VRelmERTMu8xFl7Tob9v+1eXpeRzF1Xz2OWGXiZhDWkvrugOn9FksnhD70OoNeoaqdnyXs4Z+ZCylEOEWtjKH4eQN5nsGFqSchVVG47UsN9mNfTpVigXtFQf4pjwy+oqyfncbFW6HjJXuetx3m3Ww8rKyTDPlM7V1NmUdeBrmyFbv8Z+N6a+Q89dHamVPnat/aORUOHtYZl+/kKek2vF9yBIDc+rDS20o1MqbMdA06RUnxN4TSZGTZOPy8v25e+mh9lXpyFtNaXz9vl7pe5L6T4LedZChCpV5JoYC9Vzb26aj1lftz9/pnbi6uBTXXQNyF/Ic5rsVMq/KCJ1+F1GutT3zYg536lTYWWoClqfoXRNXBVK2SeBTSsmUydP6z7een1VI0bKa/o69JDl4pIGTumacNVBqA/aVDZfeWK2HHCNqHwdb6yf39b50CIfuUUOtLNPdvkANJkrCBXFtbVhhtC1lUohSA1n9NWVaaIv5OUL9fzVcZ0jpZ1UI0Nsm2y5iI1ASX3GXSOqkPOYRjopbxESnsz0kb+Qt+17jEl9RC6UoVtDjJpoa7TU9SjM5981pbW18NFK3VKzCU7M+yVdVqErRjomnzHEiHH9OrF1b8tnyHliVnr6JsElQxY95C/kWg9HxIbyVpq+k8l6yzmlvD4tZgK+Lh42iy70+7ZzmFwYddePK99NfLwxYhw6txCbz5DzhK70jNkNsQPGIeRdTbitrISv4ptqqk6ODXWkUE3VobTL4o1pY7G+e5t41C26WMvQFUdvEyLX5HGqRR47IR0yoe7bR9yWD1+H2MdKTwHGIeRan2z4bQhI6cboW3yapC4s5Fjhqz4QXdaFKbbX9f7GvT3/qCulzKEiEON7jal/Uz1ILoyx5cO0d40vxNUXo+7Lp2vCtbqy2EYPE5khjEfIq7RpCS4Ww9z0KCTdcUdRPyFivrqa7ipKcUUA7btilpZOToTVXQt1a9q0fLxq6c5mRT2lbNPrGpbbfKyhvtfYejRZwRI+Xpe/v17/MddJnXStfz/2uk1W/7boMx+fkLftZunDrSJ5zfLhqVuep06dFLnU66ZMDpZC0tb9q++PHXOt0JWB1fP6Ik9c1p+ERRxb/23FOLcRU91DVMgrpE58tpzf8Ql5266DxSLvCTzbqsfqi4KbvJiitGBjQ8ZcC29Sd7cryyvRVmIfPN9cismd4avzmOF7U4tcijZcEVLn9FnJTdw5IWstBOt8fEIeaolUF32Y4n1d/rsm1mpKyqnjKIUzxrIO8U3GnrNMtp31YtpKU6Hwnc80QjKlGCs2po22uQ9IG9aohJXvy5fr/yEdQOwe7A0Zh5CH9H6m5Gu8Lj9lV8Loe6tMGXXRVX7W18Nm/ss8VX3OtgcwRiBTrPXYFXjSD55PeGJ86jGERFiZXskmLbzS/mEJi9x3jibXaOt+OshfyJv4VGMjBlIWUcQm0+ZHrobR9eigFGpbh5ny8gSXQMaIgMQKvJUVe/hg6oPnun+hLqym26n6BKU890CjMo7R5hxC2VaaWP0hzyN95DWaWqMhDbRJZxET+VFaR3Xxcn0ntA6kXi7hmpQM3X869B7EPrCpE1ESYW0umk7gVndoTLVofeISssOhdLRFE0u9qZXvayuuDbRSz23aw0eI/IW8qTUaYg02OX+shbq5Gb74ISbSQ8Jqry+YqIfupYi4pDDbXFChPnjTvZd68FLaUrmhlkSnEnLtkIl8KUuyz8iTkOs3EXLb82iKnBIifyEPaaApM8dSYXB1n3GT87gafUjUgy2FTLa5Qub29uI2Ygq1TGKHtzm4BUI71PKhDylTSKcT2tmHHCdRn0O4V656azqhanseW+qs8hdyX8PzuQJirUFbci3LrjaSJm/SCbEQYydiq/VjmkQMEdyYDqT6gmXf+WMf9jZilqUJtYxLfGWKadc+g6IaceQ6TqI+h3yvXHNAkhOqguQv5FqH9X57e8f91dWX/ZpEJcUV4epQQoZtIQ94SG9ua0Bt7RERU4bS0gyNxQ15ZZZk/LWLJp2dq0yu8qVGV9RdAHU3mGnv+XqoZptCNASL3ITr/khPqAoyDiEvcYULmob+5SZYJrGIXZ4eEmFStX5T9w537bdcrQdTmVKEJ4TYMsRYOy7LPaQe2/Tppl7LNLcQE5NcvY7L4Ki2/1A3Yd1l05Yfu28fuQ3XRKUrb6Z2mrtFDuDfAvg9AJ8B8CsAXh3yPbG9Vqq4hkmu1OTlv76Z//IBlnjXaMwmQ20h9VLhqqXSxPVSrWupcvcQG/wKrroIMRxi3IR1a7HNdtRlGw0lxYp2GU4ddVZtCfl3Azh19PtPAfipkO+J7n5oCyOTTqYGaBvqm1wb1QUtTfY26ROpunStCjU9AL76knxgQu9N1z5e15yILyZ6iG2pL0JH0yZclndHnVXrrhUA3wdgP+RY0f3IU2OayxRixdv8ujGWtiv6IDT/fU8Qufy0MdZguWVs6JC0q1C5kGv1KYS+d0p2PWdSZ4iWd5XYeYs6A5i47ULI/weAB0OOTRbypqF9phsX8r7F6sOSGmIYO2QbinhU8+ibaI6pl/KBD6mrEB+55ERnyL1YW2u+eCcWUyRUdeKy6zmTKkP1hVfxWeK+vA5g4jZZyAE8AeAZQ3qgcszOkY9cOc5zAcABgIONjY20UkguUa9PkrqEpRSXJi4c382uT4y53hIeg4SVZCt3yP4drvqMeTB8USumjjK17KaoFd+oyXV/pO6BqTM7derkPejDKh6AyHmRiBnvubNqzSIH8E4AnwCwGvodcYs8dlVlqh8sxOI05aUeHhnyoNWFK2W1mFTDi5lsDJ10Ll/U4BPEegcXek9NZTftbyNVB64OKDb+O3Wys08G4HbwItHZ9Ow+amuy860APgfgTMz3xH3k5YRntYI3N+2NK3Wz/9B9LFzhkTEPdVMRlrKSYjcHSh05pVr49S0FXGVPrc+QMplEK/QeNAk/HIJY5mCRD8CibkpbQv4sgC8C+PRReizke2ILgnxWakjjskXBmIS4Sx+bxIMhZSXFTgCmzmXETnJWU2jZU+sz1SIPvQe++x1jkfdhNeYikkOfkPWQ/4KglIYiuew5NQ8h+SkfRsDvlogRYSkrKcYyjjneV7YmYhzaCTTxj4a0g9B7ENJWbYvdfKOYrgQ1c5HMgfyFPEWUUq0cl8XeJAIgNDwsVrRsSD7UIb7v6iZZLtdWaNmauEfaiASq++tD/O2h9yC0LfpGpDm4OEgy+Qu55Eosn98x1WL3YTpfjNilXLuJlWTqxGKiUmJEPHQPlnpaWwuLFpnP7a/0a5uQeyDV1nKYdCTJ5C/krogVFymRADEWeyz1/ISKYtdDVZuwVFeolhZ67NYIsRatr55S9j4Z4tBfIn+0yEdN/kK+t2deSel66W7IOZtY7BKEhPb1QYwgdLE03HcNClVBLpOOKbg6uqF30kLkL+Ra+5co20htAF1YNy73QX0iK/a8TRp2TCcW6vtvIii+Do+ug9uMUdR8ocdj7bxqjEPI2/CTu4iJW3eFLYZcx7fYJda33WUMus8NIyEoPn85LfJx42qPE3InjUPI24hc8WET7LqoLC83W1Yv5QeWKHNZ7tjQzLatQFuHN1LrqxVytdZdRtyEJnjHIeSmWFqfj7yNmxyzUCVUPEP8zE3PVa6ADH2QJSNeJAUjVzHqm5xdELTItdZjEvK61evzI7dxk2Mm96odRoo/PqXzcUX4NHmQQwV0SIJB0b9NzoJHH7nWeixCntIQ27jJKRa5Lx8hMdMx7qCYfdtDzhtTj0MRjAk94EHk7oJg1MpIhDy1IUrf5BQfeejKvfK4elljd+4zlbnJgywRjti1YAylQxkKrI/sGYeQD6khxkatxIpbG6sSm9Sfa6QgeR1JhtKhDIW2RihDt4aHnr8IxiHkOQ+Vm4iblDA2qb8YURzKfRpKh9IWKQLVxeh0SM/k0PMXyTiEXOt8e9euRDQkHyn1F2ORN7mOJCN7iI8xlLINvbMcev4iGY+QpxIrLG0IUeo5h9AYY4V8KAyhQ2mDIbQJrYfvvhp6/iKZtpCnLG4ZgrUzpPzYIl58m5aRdhiKQA2lQ7Ex9PxFMm0hj72ZQ7z5fVuWMTH8fed1CgyljQ7ByHAx9PxFMm0hj7VehmLtDI0QgR7ZgzNYhlTPQ++4h56/CMYn5L6bY9uXI1eLPBdYd90xIoEiYdiEfIYc2d8HLlwADg8LmTg8LP7e37/9/3e9C7h2zX2e8+fNn+/uAqurxz9bXS0+J26eey7uc5LO1hZw5Qpw61bxc2sr/hz7+8DZs8BsVvwsnyGSFyZ1bzs1tsh9Vl/oEnrf0n5aO/FM2SLPrc0MyT1DgsCoXCs+H3boplZdxmF3wRDytrcXv0PlGMhRFKfc6WaKTcjzdK1sbLg/t/0/9Dw2fC6dPhlS3pRy/z1GdnaA69ePf3b9evH5UKEbbDTkKeQ+H/buLnDqlPscq6uFjzzGPzjkh3UoedvZAW7cOP7ZjRvDqKM2yVEUfQYRyYY8hXxrC7h4EVgsCmtvsSj+Lid7traAV73K/v3FAnjnO4EPfCDOgh3ywzqUvA0lH12ToyhyUn805CnkgH/G/oUXzN9Tqjj+0qV4C3bID+tQ8jaUfHRNjqLoM4imSo6RPCbHedupkwVBvokcqRc5A0W8et+TWkOZbBtKPvpgCJPNpBkDb78YVdRKCL4bkjpjb1toNISbPRQhGUo+CIll4JE8NiFXxf+65dy5c/rg4KD9C+3vF66S554rhva7u7eHjWWUR9W9sroaNrQ8e7bwqddZLAq3DSEkT2azQrrrKFW4cXtGKfW01vpc/fN8feQhuPzoTfyDU53QI2TsZDrHM24h95G6xDnTm00I8ZDjpDWmLuSpZHqzCSEeMo3kaSTkSqmfVEp9Rin1aaXU/1JKfZNUxgZN7M3OMZxpKLDu0mHdpSGxGVnXmGZAQxOAuyu//1MAj4V8r5dXvfXFwMOZBg3rLh3W3ShB21ErSqkfB7Chtd72HdtZ1MoQYIRLOqy7dFh3o6S1qBWl1K5S6osAtgD8hOO4C0qpA6XUwdWrV5teVp62hqGMcEmHdZcO625SeIVcKfWEUuoZQ3oAALTWO1rr1wDYB/AjtvNorS9qrc9prc+dOXNGrgQStLlzICNc0mHdpcO6mxReIddav0Vr/QZD+tXaob8A4O+3k82WaXPnQEa4pMO6S4d1NymaRq28rvLn2wD8XrPs9ESbw9BMw5kGAesuHdbdpGg02amU+m8AXg/gFoBDAA9prf+v73uDm+zkxBAhJANsk52ety+40Vrn6Uqps7tr3neFw1BCSAZwZSfAYSghJGsaWeSjYmuLwk0IyRJa5IQQkjkUchfcqyIc1hUhvUHXio36iyfKRUIAXTB1WFeE9Mq43xDUBIYkhsO6IqQTpvmGoCZwr4pwWFeE9Eq+Qt62T5Z7VYTDuiKkV/IUcslNrmwdAveqCId1RUi/mDYpbzs1frHEYnF8w/wyLRZx5/Ftvr+3V5xTqeInN+W3w7oipHXQ9oslYmg82TmbFbJbR6ni9UyhcJKOEJIR45rslPLJcpKOEDIC8hRyKZ8sJ+kIISMgTyGX2uSKk3SEkBGQ78pOiU2uyu/v7BTulI2NQsS5GpEQkhH5CrkU3PWQEJI5ebpWCCGEvAKFnBBCModCTgghmUMhJ4SQzKGQE0JI5lDICSEkcyjkhBCSORRyQgjJHAo5IYRkDoWcEEIyp5f9yJVSVwEYNgIP4l4AfySYnRxgmacByzwNmpR5obU+U/+wFyFvglLqwLSx+phhmacByzwN2igzXSuEEJI5FHJCCMmcHIX8Yt8Z6AGWeRqwzNNAvMzZ+cgJIYQcJ0eLnBBCSIXBCrlS6q1KqS8opZ5VSv2Y4f93KKV++ej/v6GUOtt9LmUJKPOPKqU+p5T6jFLqslJq0Uc+JfGVuXLcDyiltFIq6wiHkPIqpf7h0X3+rFLqF7rOozQB7XpDKfWkUupTR237fB/5lEQp9T6l1FeUUs9Y/q+UUv/xqE4+o5T6jkYX1FoPLgFYAvC/AXwrgBUAvwPg/toxDwN47Oj3twP45b7z3UGZvwvA6tHv21Mo89FxdwH4KIBPAjjXd75bvsevA/ApAN9w9Pdf6DvfHZT5IoDto9/vB3Cl73wLlPtvAfgOAM9Y/n8ewP8EoAC8CcBvNLneUC3y7wTwrNb6D7TWNwD8EoAHasc8AOADR7//VwCbSinVYR6l8ZZZa/2k1vr60Z+fBHBfx3mUJuQ+A8BPAvg3AF7qMnMtEFLefwzgZ7TW/w8AtNZf6TiP0oSUWQO4++j3VwH4Uof5awWt9UcBvOA45AEAH9QFnwTwaqXUN6Zeb6hC/s0Avlj5+/mjz4zHaK2/DuBPAMw7yV07hJS5yrtR9Og54y2zUuqvAniN1vrXusxYS4Tc428D8G1KqY8rpT6plHprZ7lrh5Ay/0sADyqlngdwCcA/6SZrvRL7vDs51Tg77WCyrOvhNSHH5ERweZRSDwI4B+Bvt5qj9nGWWSk1A/AfAPyjrjLUMiH3+BQK98rfQTHi+phS6g1a6z9uOW9tEVLmHwTw81rrf6eU+usAHj8q8632s9cbovo1VIv8eQCvqfx9H04Ot145Ril1CsWQzDWUGTohZYZS6i0AdgC8TWv95x3lrS18Zb4LwBsAPKWUuoLCl/jhjCc8Q9v1r2qtX9Za/x8AX0Ah7LkSUuZ3A/gQAGitPwHgNIr9SMZM0PMeylCF/LcAvE4p9S1KqRUUk5kfrh3zYQDvPPr9BwB8RB/NImSKt8xHbob/hELEc/edAp4ya63/RGt9r9b6rNb6LIp5gbdprQ/6yW5jQtr1f0cxqQ2l1L0oXC1/0GkuZQkp83MANgFAKfWXUQj51U5z2T0fBvBDR9ErbwLwJ1rrP0w+W9+zu45Z3/MAfh/FjPfO0Wf/GsWDDBQ3+78AeBbAbwL41r7z3EGZnwDwZQCfPkof7jvPbZe5duxTyDhqJfAeKwD/HsDnAPwugLf3necOynw/gI+jiGj5NIDv7jvPAmX+RQB/COBlFNb3uwE8BOChyn3+maM6+d2m7ZorOwkhJHOG6lohhBASCIWcEEIyh0JOCCGZQyEnhJDMoZATQkjmUMgJISRzKOSEEJI5FHJCCMmc/w8y9p8+kis0FwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def get_correlated_sample(n,uniform=True,xmin=0,xmax=1,mean=0,sigma=1,noise=1) :\n",
    "    \"\"\"Generated correlated sample data set\"\"\"\n",
    "    if uniform :\n",
    "        x=np.random.uniform(xmin,xmax,size=n)\n",
    "    else :\n",
    "        x=np.random.normal(mean,sigma,size=n)\n",
    "    y=x+np.random.normal(0,noise,size=len(x))   # add some random noise here, Gaussian with input noise standard deviation\n",
    "    return x,y\n",
    "\n",
    "x,y=get_correlated_sample(10000,noise=1)\n",
    "plt.figure(figsize=(6,6))\n",
    "plt.plot(x,y,'ro')\n",
    "\n",
    "# calculate and print the covariance matrix\n",
    "c=cov(x,y)\n",
    "print('covariance: ',c)\n",
    "c_matrix=np.cov(x,y)\n",
    "print('covariance matrix:')\n",
    "print(c_matrix)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How do you expect your results for the covariance matrices change if you changed the width of the underlying distributions, both for uniform and normally distributed variables?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>ANSWER HERE:</b>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Play around with your data sets, changing the widths and types (uniform and normal) of the underlying distributions to build up some intuition about the covariance matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "# code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OK, now let's work the other direction, i.e., given a covariance matrix, generate a multivariate distribution. For normal distributions, we can use the [numpy.random.multivariate_normal()](https://docs.scipy.org/doc/numpy-1.15.1/reference/generated/numpy.random.multivariate_normal.html) routine. Play around with this for different values of the input means and covariance matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/holtz/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:13: RuntimeWarning: covariance is not positive-semidefinite.\n",
      "  del sys.path[0]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x12bc42890>]"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAFlCAYAAAD292MqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAbHElEQVR4nO3db4xcV3nH8d+za6ftEiIae+mfJDsbJEQbRVQhKwSNhBChUjBRwgsqQSfIhUpWUUvTCkRDV2pfVJYqtULkBYq0CqRRPAJVgQpUuZQ/pa36goh1oCVgaKPUdgyhWdtt88dIJvbTF3cmnh3fO3Nn7rl/zr3fj7Ta7PX63jMh/ObMc59zrrm7AADxWqp7AACAYghyAIgcQQ4AkSPIASByBDkARI4gB4DI7anjovv37/f19fU6Lg0A0Tp27NgZd1+dPF5LkK+vr2t7e7uOSwNAtMzsZNpxSisAEDmCHAAiR5ADQOQIcgCIHEEOAJELEuRm9ioze9TMvm9mx83szSHOCwCYLVT74f2SvuTu7zazqyStBDovAGCGwkFuZtdIeouk35Ykd78g6ULR8wIA8glRWnmNpB1JD5nZt8zsQTN7RYDzAgByCBHkeyS9QdID7n6LpBcl3Tf5S2Z2yMy2zWx7Z2cnwGWHBgNpfV1aWkq+Dwbhzg0AEQgR5KclnXb3x4Y/P6ok2Hdx9y1333D3jdXVK7YKWMxgIB06JJ08Kbkn3w8dIswBdErhIHf3H0t62sxeNzx0u6TvFT1vLpub0vnzu4+dP58crwKfBgA0QKiulQ9JGgw7Vp6S9P5A553u1Kn5joc0+jQweiMZfRqQpH6//OsDwFCQPnJ3//awbPJ6d3+Xu/9PiPPOtLY23/GQ6v40AABDca/sPHxYWploWV9ZSY6Xrc5PAwAwJu4g7/elrS2p15PMku9bW9WUNur8NAAAY+IOcikJ7RMnpEuXku9V1afr/DQAAGPiD/K61PlpAADG1PKot9bo9wluALVjRg4AkSPIASByBDkARI4gB4DIEeQAEDmCHAAiR5ADQOQIcgCIHEEOAJEjyAEgcgQ5AESOIAeAyBHkABA5ghwAIkeQA0DkCHIAiBxBDgCRI8gBIHIEOQBEjiAHgMgR5AAQOYIcACIXT5APBtL6urS0lHwfDOoeEQA0wp66B5DLYCAdOiSdP5/8fPJk8rMk9fv1jQsAGiCOGfnm5uUQHzl/PjkOAB0XR5CfOjXfcQDokDiCfG1tvuMA0CFxBPnhw9LKyu5jKyvJcQDouDiCvN+XtrakXk8yS75vbXGjEwAUS9eKlIQ2wQ0AV4hjRt5m9McDKCieGXkb0R8PIABm5HWiPx5AAAR5neiPBxAAQV4n+uMBBECQ14n+eAABEOR1oj8eQAB0rdSN/ngABTEjjxk96AAUcEZuZsuStiX90N3vDHVeZKAHHcBQyBn5vZKOBzwfpqEHHcBQkCA3s+slvVPSgyHOhxzoQQcwFGpG/glJH5V0KesXzOyQmW2b2fbOzk6gy3ZYVq/5tddSNwc6pnCQm9mdkp5192PTfs/dt9x9w903VldXi14WaT3oV10lPfdcUi93v1w3J8yBVgsxI79N0l1mdkLSZyW9zcyOBDgvpknrQX/lK6Wf/nT371E3B1rP3D3cyczeKukjs7pWNjY2fHt7O9h1MbS0lMzEJ5lJlzKrXgAiYWbH3H1j8jh95G3C3i1AJwUNcnf/J3rIa8TeLUAnMSNvE/ZuATqJIG+iIkvv+33pxImkJn7iBCEOdABB3jSjpffjLYQf+IC0fz+94QBSsfth06Qtvb9wQTp7Nvln9lQBMIEZeRXmKZXkWWJPbziAMQR52dJKJdNWW+ZtFWRPFQBDBHnZ5t2lMK2FMA294QCGCPKyzbtL4WQL4b590t69u3+H3nAAYwjysi2y2nK8hfDMGemhh+gNB5CJIC9biNWW9IYDmIIgLxurLQGUjD7yKvT7BDeA0jAjB4DIEeQAEDmCHAAiR5ADQOQIcgCIHEEOAJEjyAEgcgQ5ij2RqAliHz9QEEHedGWH1Lzb7DZN7OMHAjB3r/yiGxsbvr29Xfl1ozMKqfFtcFdWwi7xX19Pwm9Sr5fs69J0sY8fmIOZHXP3jSuOE+QNVkVILS0lM9lJZskmXU0X+/iBOWQFOaWVJpt3L/NFLLLNbhGhS0VVjx9oIIK8iUZhl/VpKWRIhdhmN68y6tlVjh9oKIK8acbDLk3okKpym92sx97dc8/is3O2CQaokTdOVl1cSkLq8OF4Qyqrnj0S+kYu0DLUyGORVf82S386UEw91LNKQtMeSp1HTP8ugIAI8jItEizz3Lyrooc6ZDim1bMnLXojl35ydJm7V/516623eusdOeK+suKexErytbKSHA/193q93b83+ur1rjxnr+dulnyfNYair2HWObPGnTb2vPL+uwAiJmnbUzKVIC9LnmDJCti8wWuWfg2z3ddYNIzLDMfQbxJ5/l0AkSPIqzYrWEIEWZ6gLRLGZYfjop8U0jAjRwdkBTk18rLMqnVnteLNc7MvTw91kUVFZS+26feTG7iXLqXfyJ0H/eToMIK8LLOCJcSqzTw91EXCOKZwpJ8cXZY2TS/7qxOlFffppYOqSgFFSzghyx8AClFGaYUFQXWpYmfD8Wttbiaz/bW1uBcVAR2WtSBoTx2DgS4HaRUB2+8T3ECLEeR1ImABBMDNzqqxjBxAYMzIqzRZFx8tI5eYmQNYGDPyKoXoHQeACQR5lap44k+XUKYCJBHk1eKxZOGw2yHwMoK8SjGtlGw6ylTAywoHuZndYGZfN7PjZvZdM7s3xMBaiWXk4VCmAl4WomvlJUkfdvfHzeyVko6Z2Vfc/XsBzt0+9I6HsbaW/kg8ylTooMIzcnd/xt0fH/7z85KOS7qu6HmBqShTAS8LWiM3s3VJt0h6LOXPDpnZtplt7+zshLwsuogyFfCyYJtmmdnVkv5Z0mF3//y032XTLACYX9amWUFm5Ga2V9LnJA1mhTgAIKwQXSsm6VOSjrv7x4sPCQAwjxAz8tskvU/S28zs28OvAwHOC8SDVaaoUeH2Q3f/V0kWYCxAnNgMDTVjZSdQFKtMUTOCHCgqllWmlH9aiyBHs8QYNtdem368SatM2WSs1QhyNEeMYTMYSM89d+Xxq65q1ipTyj+tFmxB0DxYEIRU6+vp+6f0etKJE1WPJp+sMe/bJ505U/lwMi0tJW+Ok8ykS5eqHw8WUuqCICCIWGrN47LGdu5cteOYhb3wW40gR3Nk1ZqzjjdBLAHJJmOtRpADUr6brGm/E0tAsslYq1EjR3PUVcedXNAjJWE8HnTTfkdKbhqeOpXMxA8fJiBRiqwaOUGO5qjrZmee68Z4Ixatw81ONF9dZYo8N1ljvBGLziDI0Rx11XHz3LCM5aYmOokgR7P0+0mp4tKl5HsVteY8nwRiuamJTiLIgTyfBBb9tBDjlgOIDjc7gbLk6YYB5sDNTqBq7G+CihDkiFtTShdp46DTBRUhyBGvpuyWmDWOrK0FlpaolSMoghzxakrpImsc0pWdLpJ08WLzt+dFVAhyxCurRJG2ArOOcZw7l9zYXF6+8s+olSMgghzxylqMY1btbHfaYqF+P3ufGGrlCIQgR7wOH05Ce5J7sdnurBuok39+4MD0xUKsCkXJCHLEq99P3y1RWny2O+sGatqfP/ywdPBg9mIhVoWiZCwIQtxC70o463yLXm8wYKtbFMY2tmin0KsnZ+2JzrMvUSNWdqKdQu+YmFW3dk9m41m94dS7USOCHPGbZ8fEWTcy0+rZIydPSs8/L+3du/s49W7UjCBHd+RZCTo+w09z4YJ0zTU8+xKNQpAjDiH2VMm7EnQ0w89y7lz1e6YDUxDkaL5pM+l5An6eTawGg/QedYl6OBqHrhU0X1bL37590k9+snuWbZaEfa93ZYvfPK2DWb9rJj3yCLNw1IKuFcQrayZ99uyVpZLRxCSt/j3Pwpysa7oT4mgcghzNt2gpY7L+PU+rYtY1s26CAjUiyNF8WTPpfftm/93JmXXeVkWW1SMiBDmaL2smff/92T3fI4vO5kMvNAJKRJAjDmkz6cme78kuk8kZ9LwtjPMsNAJqRJAjbqOwdU+6SbJm0E15LBxQAtoP0Q2hd0kEakD7IbqNJ9qjxQhydANP6UGLEeToBtoJ0WIEObqBdkK0GEGO5gmx02Ea2gnRUkGC3MzuMLMfmNmTZnZfiHOio2JsEyzrjQfIqXCQm9mypE9KeoekmyS918xuKnpedFTePcObIsY3HrROiBn5GyU96e5PufsFSZ+VdHeA86KLYmsTjO2NB60UIsivk/T02M+nh8eA+cXWJhjbGw9aKUSQpz1G5YrlomZ2yMy2zWx7Z2cnwGXRSrG1Ccb2xoNWChHkpyXdMPbz9ZJ+NPlL7r7l7hvuvrG6uhrgsmilMtsEy7gpGdsbD1ppT4BzfFPSa83sRkk/lPQeSb8V4LzoqtHOhiGNbkqO6tmjm5Kj6y1q9Hc3N5NyytralY+YA0oWZNMsMzsg6ROSliV92t2nTkfYNAuVY9MstEDWplkhZuRy96OSjoY4F1AKbkqixVjZiW7gpiRajCBHN3BTEi1GkKP9BoPLC3eWl5NjbJqFFglSIwcaa7Jb5eLFyzNxQhwtwYwc9alis6kmLKGve1Otuq+P8rl75V+33nqro+OOHHFfWXFPtppKvlZWkuMhme2+xujLbLEx93rJ3+318o21qtfZ1OsjKEnbnpKpPHwZ9Sirr3tUDx8tznnhBens2eLXmSzRSEmJZladve7+9bqvj6B4+DKapYy+7rQtZZ9/Xtq7d/fvLdKtsmiJZp7XWUYJhP75TiDIUY8y+rrTwvbCBemaa4rv3bJoIOZ9nWXta07/fCcQ5KhHGX3dWaF67ly+R7xNmxEvGoh5X2eeGf8iM3b657shrXBe9hc3O+Hui908nKbXS7+x2evlG8u0m4JFbhrmeZ2zbsqWfX1EQRk3OwlytEeRsMvzJlBmIM66fpE3KbRGVpBTWkF7FNnLPE8NvN/PV6JZRFoJxEw6cCD/+NBZBDnaZdGwrfumYL8vHTyYhPeIu/Tww0ktvO7xodEIckAKc1OwaPvg0aNJeI8b3fAczcwnZR1HpxDkgFT8EXMh2genlU+OZmz3/8ADLLsHKzuBQkYrSdNWT0rzraCctgrz1KkrZ+vj8qwyRfRY2QnMI0+ZZHwWnmWem5HTyjuzauHnz0v33pv/WmgVghyYlLdMkraIZ9I8NyOnlXfSQn7S2bOLl1jYITFqlFaASXk3mlpaqrbcMauMkzbGvOddZEMwVI7SCpBX3p7tabPtMp5ANGqtPHIk+3cW6Stvwp7tKIQgBybl7dnOqmkfORJ+wdC4fl/aty/9zxbpK2exUfQIcmBS3p7yoi2LRdx/f7jNsFhsFD2CHJg0T0CXuWw/1BhnYYfE6BHkQJq6AnoeocYY6k0hb+cLHTLBEeRAKDEHVNE3hbwtm2U9QKPjaD8EQuh6C1/elk2eIVoI7YdAmWJs4Qv5CSJv5wsdMqUgyIEQYguo0CWOvJ0vdMiUgiAHQogtoEJ/gsjb+UKHTCkIciCE2AIq9CeIvJ0vdfbetxhBDoTQ5IBKq4WX8Qkib+dLDK2dkSHIgVCaGFBZtfADB+L6BIGpCHKgzbJq4UePNvcTBOZGkANlq3Oh0LRaeBM/QWAhBDlQpiJtfiHeAK69dr7jiBJBDpRp0TY/lrJjDizRB8qU9RQhs6SkkSXUUvZFr49GYok+UIdF2/xC9XnHtlAJCyHIgTItulAoVADHtlAJCyHIgTItulBo0QCevEEq0WbYAdTIgaYaDJKboqdOJTPxw4enB3DXt9LtAGrkQGzm7fNu6la6MT9wIxIEORCz8ZBM63KR6t1KlzbKShQKcjP7SzP7vpn9u5n9rZm9KtTAAMwwGZJZ6uxQaeqnhJYpOiP/iqSb3f31kv5D0seKDwlokTLLCmkhOWnvXumFF6ora0y+3iZ+SmihPUX+srt/eezHb0h6d7HhAC0yefNxVFaQwtx8nBaGZsmNzhdflM6eLef6k9Jer1n6pwX62IMKWSP/gKS/z/pDMztkZttmtr2zsxPwskBDlV1WyArDXk965JH02XqZZY201+uehPk4+tiDmxnkZvZVM3si5evusd/ZlPSSpMzPbe6+5e4b7r6xuroaZvRAk5X9HM9pveabm9l187LKGlnndaePvWQzg9zd3+7uN6d8fUGSzOygpDsl9b2OpnSgqcpeHj9tsVFWbXrW9YvU9Kd9QmC73FIV7Vq5Q9IfS7rL3WfcdQE6porl8Vm95svL08eVpmir4IED8x1HMIVWdprZk5J+RtLwboq+4e6/O+vvsbITnTHv6sxQJuvS47L+P190x8VQOzYiU9bKTpboA220SKgW3fKWLXNLxxJ9oEvSyjpS0lOeVSopWtNny9zaEORAG41uhO7bt/v42bPZde+iNX22zK0NQQ60Vb8vXX31lcezeskX3XI31N/HwqiRA21G3bpVqJEDXUTduhMIcqDNqFt3AkEOtBl1604otPshgAj0+wR3yzEjB4DIEeQALuP5mlGitAIgUfaDMFAaZuRAGy0ys+b5mtFiRg60zaIz67IfhIHSMCMH2mbRmXVTFg9Rp58bQQ60zaIz67oXDw0G0v790j33LP5wi44iyIG2WXRmXefioVE56OzZK/+MOv1MBDnQNkVm1lmPjitbWjloHHX6qQhyoG1iXJY/K6jZ5GsqulaANoptWf7aWvqj6SQ2+cqBGTnQBrF3emQ9mm7fvuZ/mmgAghyI3ehGYcydHv2+dPCgtLyc/Ly8LH3wg9KZM4R4DgQ5ELs2rMgcDKSHH5YuXkx+vngx+TmmN6MaEeRA7NqwIrMNb0Y1IsiB2DVlRWYRbXgzqhFBDsSu7hWZIbThzahGBDkQuxj7xie14c2oRgQ50AaLrshsSttiG96MasSCIKCrmvYgidgWMTUIM3Kgq+gUaQ2CHOgqOkVagyAHuopOkdYgyIGuolOkNQhyoKvoFGkNulaALqNTpBWYkQOIU1N64BuAGTmA+DStB75mzMgBxIce+F0IcgDxoQd+F4IcQHzogd+FIAcQH3rgdyHIAcSHHvhd6FoBECd64F/GjBwAIkeQA4jHaBGQmbRnT/K944uBpEBBbmYfMTM3s/0hzgeg4epYVTlaBHTyZPLzxYvJ99FioA6HeeEgN7MbJP2GpG42cAJdMx6o7vMFaZE3gLRFQCMxLAYq8c3P3L3YCcwelfTnkr4gacPdz8z6OxsbG769vV3ougBqsr5+eVY8rtdLnheaZXJZvZS0DObtNllaSt44spglzyxtoqKvfcjMjrn7xhXHiwS5md0l6XZ3v9fMTmhKkJvZIUmHJGltbe3Wk2n/IQBovqxAnRWki74BzPr7856nDkVf+1BWkM8srZjZV83siZSvuyVtSvrTPANw9y1333D3jdXV1dwDB9Awi66qLLqsPm0R0EjTFwOVvKXAzCB397e7+82TX5KeknSjpH8bzsavl/S4mf1ikJEBaKZFV1UWXVY/vghIkpaXk+8xLAYqeUuBhW92uvt33P3V7r7u7uuSTkt6g7v/OMjIADTToqsqQyyr7/eTUoS79NJLyfcTJ5od4lLpWwrQRw5gfqNAvXQpf5B2eVl9ya+9cNfKIuhaAYD5LXyzEwDQbAQ5AESOIAeAyBHkAMrF0+5LR5ADKM+8+7IQ+gshyAGUZ56n3RfZjKvjCHIAYY3PqrP2Rklbmj5P6GMXHvUGIJy0Xf7SpC1NL3k/kjZjRg4gnGl7ho9kLU0veT+SNiPIAYQzbfY8a2l6yfuRtBmlFQDhrK0tvu/2KNw3N5M3hLW1JMS7sBdLQczIAYRTdFa9yGZcIMgBBNTlHQ5rRGkFQFj9PsFdMWbkABA5ghwAIkeQA0DkCHIAiBxBDgCRI8gBIHIEOQBEjiAHgMgR5AAQOYIcACJHkANA5AhyAIgcQQ4AkSPIASByBDmA9hgMpPV1aWkp+T4Y1D2iSrAfOYB2GAykQ4cuP/z55MnkZ6n1+6MzIwfQDpubl0N85Pz55HjLEeQA2uHUqfmOtwhBDqAd1tbmO94iBDmAdjh8WFpZ2X1sZSU53nIEOYB26PelrS2p15PMku9bW62/0SnRtQKgTfr9TgT3JGbkABA5ghwAIkeQA0DkCHIAiBxBDgCRI8gBIHIEOQBEjiAHgMgR5AAQOYIcACJn7l79Rc12JJ0McKr9ks4EOE9MeM3d0cXXzWuerufuq5MHawnyUMxs29036h5HlXjN3dHF181rXgylFQCIHEEOAJGLPci36h5ADXjN3dHF181rXkDUNXIAQPwzcgDovCiD3MzuMLMfmNmTZnZf3eMpm5ndYGZfN7PjZvZdM7u37jFVycyWzexbZvZ3dY+lCmb2KjN71My+P/zf/M11j6lsZvZHw/+2nzCzz5jZz9Y9pjKY2afN7Fkze2Ls2LVm9hUz+8/h95+f97zRBbmZLUv6pKR3SLpJ0nvN7KZ6R1W6lyR92N1/VdKbJP1eB17zuHslHa97EBW6X9KX3P1XJP2aWv7azew6SX8gacPdb5a0LOk99Y6qNH8t6Y6JY/dJ+pq7v1bS14Y/zyW6IJf0RklPuvtT7n5B0mcl3V3zmErl7s+4++PDf35eyf+xr6t3VNUws+slvVPSg3WPpQpmdo2kt0j6lCS5+wV3/996R1WJPZJ+zsz2SFqR9KOax1MKd/8XSecmDt8t6eHhPz8s6V3znjfGIL9O0tNjP59WR0JNksxsXdItkh6rdySV+YSkj0q6VPdAKvIaSTuSHhqWkx40s1fUPagyufsPJf2VpFOSnpH0f+7+5XpHValfcPdnpGTSJunV854gxiC3lGOdaL0xs6slfU7SH7r7c3WPp2xmdqekZ939WN1jqdAeSW+Q9IC73yLpRS3wUTsmw5rw3ZJulPTLkl5hZvfUO6q4xBjkpyXdMPbz9Wrpx7BxZrZXSYgP3P3zdY+nIrdJusvMTigpob3NzI7UO6TSnZZ02t1Hn7geVRLsbfZ2Sf/l7jvu/lNJn5f06zWPqUr/bWa/JEnD78/Oe4IYg/ybkl5rZjea2VVKbop8seYxlcrMTEnN9Li7f7zu8VTF3T/m7te7+7qS/53/0d1bPVNz9x9LetrMXjc8dLuk79U4pCqckvQmM1sZ/rd+u1p+g3fCFyUdHP7zQUlfmPcEe4IOpwLu/pKZ/b6kf1Byd/vT7v7dmodVttskvU/Sd8zs28Njf+LuR2scE8rzIUmD4UTlKUnvr3k8pXL3x8zsUUmPK+nQ+pZausLTzD4j6a2S9pvZaUl/JukvJP2Nmf2Okje135z7vKzsBIC4xVhaAQCMIcgBIHIEOQBEjiAHgMgR5AAQOYIcACJHkANA5AhyAIjc/wMsdFZ7TiLU2QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n= 100\n",
    "\n",
    "# The desired mean values of the sample.\n",
    "mu = np.array([5.0, 0.0])\n",
    "\n",
    "# The desired covariance matrix.\n",
    "r = np.array([\n",
    "        [  3.40, -05.05 ],\n",
    "        [ -05.05,  5.50]\n",
    "            ])\n",
    "\n",
    "# Generate the random samples.\n",
    "y = np.random.multivariate_normal(mu, r, size=n)\n",
    "plt.figure(figsize=(6,6))\n",
    "plt.plot(y[:,0],y[:,1],'ro')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "OK, now let's work with quantifying correlation for distributions with arbitrary width, using the correlation coefficient, which is essentially a normalized covariance.\n",
    "\n",
    "[scipy.stats.pearsonr](https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.stats.pearsonr.html) computes the (Pearson) correlation coeffient, and returns it plus the probability of getting the value given the number of points, given the null hypothesis. [scipy.stats.spearmanr](https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.stats.spearmanr.html) does the same thing for Spearman's correlation coefficient. Let's write a quick routine that calculates both for an input data set:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pearsons:  (-0.9088473503370985, 5.3988638206682304e-39)\n",
      "Spearmans:  SpearmanrResult(correlation=-0.9007260726072607, pvalue=2.894703262889244e-37)\n"
     ]
    }
   ],
   "source": [
    "def corrcoef(x,y) :\n",
    "    \"\"\" Compute Pearsons and Spearmans correlation coefficients and probabilities\"\"\"\n",
    "    print('Pearsons: ',scipy.stats.pearsonr(x,y))\n",
    "    print('Spearmans: ',scipy.stats.spearmanr(x,y))\n",
    "    \n",
    "corrcoef(y[:,0],y[:,1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Experiment with generating sets of data with different covariance matrices and measuring the correlation coefficients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What if we want to make a nonlinear relation? One way to accomplish this is to generate a random sample from a distribution of each x point, with some nonlinear relation. Play around with the relation and the scatter around it, and see the effect on the measured correlation coefficients."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.18"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}