ex1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Klimadynamics 1\n",
    "\n",
    "#### Nils, Nastaran and Rovina"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Write a program, that evaluates a bivariate Gaussian PDF $ f(\\vec{x}; \\mu, \\Sigma) $ with given parameters $ \\vec{\\mu} $ and $ \\Sigma $ on a regular 2-dimensional grid\n",
    "   * Define the covariance matrix, $ \\Sigma $ by its variance $ \\sigma_{1} $ and $ \\sigma_ {2} $, and correlation, $ \\rho $  \n",
    "   * Apply this $ f(\\vec{x}) $ to a mesh-grid with range: \n",
    "       $ x_{1} = (-3\\sigma_{1},\\cdots,3\\sigma_{1}) $ and $ x_{2} = (-3\\sigma_{2},\\cdots,3\\sigma_{2}) $\n",
    "   * Plot a contour plot with axes $ x_{1}, x_{2} $\n",
    "   \n",
    "Note: $ \\Sigma = \\begin{pmatrix} \\sigma^{2}_{1} & \\sigma_{12} \\\\ \\sigma_{21} & \\sigma_{2}^2 \\end{pmatrix} $ and $ \\rho $ = $ \\frac{\\sigma_{12}}{\\sigma_{1}*\\sigma_{2}} $\n",
    "     where, $ \\rho \\in (-1,1) $  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from numpy.random import multivariate_normal as mv_norm\n",
    "from numpy.linalg import det, inv, eig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_multivariate_normal(x, mean, cov):\n",
    "    \"\"\"\n",
    "    mean : 1-D array_like, of length N\n",
    "\n",
    "    Mean of the N-dimensional distribution.\n",
    "\n",
    "    cov : 2-D array_like, of shape (N, N)\n",
    "\n",
    "    Covariance matrix of the distribution. Ref lecture notes eq 55\n",
    "    \"\"\"\n",
    "\n",
    "    q = mean.shape[0]\n",
    "\n",
    "    Z = np.sqrt((2*np.pi)**q*det(cov))\n",
    "\n",
    "    dist = 1/Z * np.exp(-0.5*(x-mean).T @ inv(cov)@(x-mean))\n",
    "\n",
    "    return dist\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma_1 = 5\n",
    "sigma_2 = 1\n",
    "rho = -0.8 \n",
    "\n",
    "mean = np.array([0, 0])\n",
    "\n",
    "cov = np.array([[sigma_1**2, rho*sigma_1*sigma_2], \n",
    "                [rho*sigma_1*sigma_2, sigma_2**2]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get values from multivariate normal distribution\n",
    "x_norm = get_multivariate_normal(x=np.array([0, 1]), mean=mean, cov=cov)\n",
    "\n",
    "x1 = np.linspace(-3*sigma_1, 3*sigma_1, 100)\n",
    "x2 = np.linspace(-3*sigma_2, 3*sigma_2, 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "x1_mesh, x2_mesh = np.meshgrid(x1, x2)\n",
    "    \n",
    "prob = np.zeros(shape=x1_mesh.shape)\n",
    "for i in range(prob.shape[0]):\n",
    "    for j in range(prob.shape[1]):\n",
    "\n",
    "        prob[i, j] = get_multivariate_normal(x=np.array([x1_mesh[i, j], x2_mesh[i, j]]), mean=mean, cov=cov)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# prob = np.apply_along_axis(func, 0, np.array[mesh], mu,sigma) ---> vectorize it intead of the for loop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWgAAAEICAYAAAByEW6PAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO2df7Ck1Vnnv0933zsDAxNEMCGAAc3oLqZWogi4WTVZowE2ZtTVKtA1hrVqZBdKs7qlSbD8kV1rUcto3LCwI2LCmgRTRjazcVYC7mbZbIUEiJFASMoJScyESQghwJBh5t7ufvaP95x+v6ffc/q83f3ee9++/Xymuu7b55z3fU/37Tn32895foiqwjAMw2gfna2egGEYhhHHFmjDMIyWYgu0YRhGS7EF2jAMo6XYAm0YhtFSbIE2DMNoKbZALxki8psi8mdbPY8YIvLTIvKBBq6jIvLiJuZkGFuJLdDbDBF5lh5DEXmOnv/0Vs/PIyLnuYW059tU9Z2q+sMbfN8Pishx9348ISJ/KSJnub63i8iaiBx1j4dE5D+JyPPo/NeJyGDsfX7bRs7ZWF5sgd5mqOop/gHgHwD8CLW9c7PmwQtvC7nOvT/fBuA0AH9Afb+rqqcCOBPA1QAuBfD/RGQXjfkwv8+qet2mzdxYKmyBXk5WReQ2pxIfFpGLfIeIvFBE3isiXxGRz4rIL1DfDhH5QxF5zD3+UER2uL6Xi8hhEflVEfkSgD8VkY6IvEFEPiMiXxWR94jI6e5y97ifTzkV+r1OnX6I7vcdInKXiDwpIl8WkTe59otF5MMi8pSIHBGRt4nI6rRvgqo+CeC9AF4S6TuuqvcBeA2Ab0SxWBvGpmIL9HLyGgC3o1CPBwC8DQBEpAPgfwD4OwBnA/hBAK8XkVe5865HoSgvBPCdAC4G8Gt03RcAOB3AiwDsA/ALAH4UwA8AeCGArwG40Y39fvfzNKdCP8wTFJFTAdwN4K/duS8G8DeuewDg3wE4A8D3unn+22nfBBE5A8C/BPC3qTGqehTAXQC+b9rrG8a82AK9nHxIVQ+q6gDAf0Ox2ALA9wA4U1XfrKprqvoogD8GcKXr/2kAb1bVx1X1KwB+C8DP0HWHAH5DVU+o6nMAfh7A9ap6WFVPAPhNAD9R0/zxagBfUtXfd2r2qKp+BABU9QFVvVdV+6r6OQD/FcUfgbr8kYg8heIP0REAv5QZ/xiKPzyeS516949Lp7i3YdSmzXZCY+P4Eh0fA7DTLZovAvBCt3h5ugD+rzt+IYDPU9/nXZvnK6p6nJ6/CMAdIjKktgGA59eY47kAPhPrEJFvA/AWABcBOBnF5/iBGtf0/IKq3jLF+LMBPEnP71XVfzbF+YYxE6agDeYLAD6rqqfR41RVvcL1P4Zi0fV8s2vzjKdG/AKAy8eut1NVvxgZG5vLtyb6bgLwKQB7VHU3gDcBkPzLmx4ROQXAK1H+kTKMTcMWaIP5KIBn3EbfSSLSFZGXiMj3uP53A/g1ETnT2W9/HcAkn+qbAfy2iLwIANx5e13fV1CYRL4lce77AbxARF7vNidPFZFLXN+pAJ4B8KyI/CMA/2bWF5zC3fO7Afx3FLbzP236HoaRwxZoY4SzSf8Iik3AzwJ4AsAtALwf8H8EcD+ABwF8AsDHXFuKt6LYhPyAiBwFcC+AS9y9jgH4bRQubBU7rtuc+yE3ny8B+HsAr3Dd/x7ATwE4isJG/uczv+gqv+Lm+iSA21CYTv6pqn69wXsYRi3EEvYbhmG0E1PQhmEYLWXuBVpEdorIR0Xk71zQw281MTHDMIxlZ24Th4gIgF2q+qyIrAD4EIBfVNV7m5igYRjGsjK3H7QWK/yz7umKe5hh2zAMY04aCVQRkS6K3e4XA7jRR3yNjdmHIvwXXXS/+2TsnueGs587LbaJahhbylF87QlVPXOea7zqFbv0q08OsuMeePDEnap62Tz3apJGFmjnnnWhiJyGInLsJar60NiY/QD2A8BuOV0v6fwQpDPjQiubuLepw/yYyilzLuoz3NMwtit36198Pj9qMl99coCP3vnN2XHds/7+jHnv1SSNrnSq+hSADwJozV8gwzAMBTCs8a9tzK2gReRMAOuq+pSInIQiLPZ3sufVUc8NKOW6Kj2peqeZg1O+uXtmFXbqnqasDWMmFIp1zZs42kYTJo6zALzD2aE7AN6jqu9v4LqGYRiN0UaFnKMJL44HAby0gbmU1FStM9uwG7xWoIZj846o3tS9plLWpqYNozYKxWABN/wt3ahhGEvBcAG9f9uzQGdUc1bhbqZnR3DbyUpWh5l5kRKOvcapbOOmqg0jigIY2AJtGIbRTkxBG4ZhtBAFsG426GZpxKzR4EZilIwJYyYTSMbsUQyJfNjM7GEYURTamIlDRC5Dkeu8C+AWVb1hrF9c/xUoSsq9TlU/5vo+hyKP+QBAX1UvmnSvVi/QhmEYjaDAoIH12bkT34iimMRhAPeJyAFV/SQNuxzAHve4BEWJtkuo/xWq+kSd+23tAj3rxl7svJRS3ui8HV26fvQrVLc8jKjemMJObizOsqFornmG4SIJG+FiAIdcxXuIyO0A9gLgBXovgNtcIrl7ReQ0ETlLVY9MezNL2G8YxhIgGNR4ADhDRO6nx76xC52NoqCx57BrqztGUZSAeyBy7QqtNHFE7a5zqmbpZP4WTaO0M5sNOkz8rY5NIaKWU3braezVtW3UY9cwjO1IsUlY6//4Exm7cOwi4//ZJo15mao+JiLfBOAuEfmUqt6TupkpaMMwtj2FH3QtBZ3jMIBz6fk5AB6rO0ZV/c/HAdyBwmSSpJUKekRK8UUVdtkWVcsphZxT1jMgqXtFlLfKMNIft1tPY6+Oj80EvZiSNrYxw3oKOsd9APaIyPkAvgjgShQV5pkDAK5z9ulLADytqkdEZBeAjqoedcc/DODNk27W7gXaMAyjAbyCnvs6qn0RuQ7AnSiU1K2q+rCIXOP6bwZwEIWL3SEUbnZXu9OfjyJfPlCsve9S1b+edL/FXqBTqjmmYHP9wdg5f5E5BU02aul2K/2BDTsQyFVlnbVXT+NTbTZqY5uiEAwasuiq6kEUizC33UzHCuDayHmPAvjOae612Au0YRhGTRoycWwqtkAbhrHtUQjWtJsf2DJas0Bnw7q5fxqXOG/a4HMi95JpNiRzSHwzTr2poEv34o07Z+IINhlpYzFu+pi8oThN0IttIhrblSJQZfGc1lqzQBuGYWwkTWwSbjYLvUAnNwZj7aSEA7Xs21OqfJZQcf5DTQpYfAe725HaHinshIKOKevkhqJTzsLf6kgBx5R1esPRNhGNxUZVMFBT0IZhGK1kaAq6JUTszVHVzGNzCjx1/RyxsHBuY5c7jSjsYURhU7sMqFIx26t9AExwrfJ1ZYNezEZtbCOKTcLFW+4Wb8aGYRhTYpuEW00uZDvlBeKPOWAk1p+6FhMbG1PQrER53n5sTmEDwMB5afBc2IvDK+uYqi6eVM5jVT23jdrUtNEyBuYHbRiG0T6ajCTcTFqzQGvgw1vzL10NW/DI9pxSxV45J/v9+Qkb9jTzyalpjdiNB7FkSgCkUMgxVQ2Qsqa2pH81nNpu0kZtHh9GyxiaF4dhGEb7KJIl2QJtGIbROhSCdQv1bgnZsPGIG13MrAHkTSCZEPIskVDvwJTBpgY2D/jXwIEq3cgmIW0MasIlz8965k1EN69sRZdiUPW1GMYGo4qFDFSZe8Yicq6I/G8ReUREHhaRX2xiYoZhGM0hGNZ4tI0mFHQfwC+r6sdE5FQAD4jIXWNlyOOwipq1wneOWCh31M0uopq5n1U3jdVAQUeunyNwqXOnc1svsWHYdcesirl/9M2g7JduYkPRjc0GvVC3ueQZi4RiMRX03Au0KyV+xB0fFZFHUFSwzS/QhmEYm8TSbxKKyHkAXgrgI01ed1OQiF0aKJUzq+bu5LDwqKqugxeSGrFLA2HQysC1D7rR/pFC7vfj/aysvZsc29Aj7nmBVZls1OITPqWSMc3ikmeq2mgQhSx3wn4ROQXAewG8XlWfifTvA7APAHbi5KZuaxiGkUUBrC9rLg4RWUGxOL9TVf8yNkZV9wPYDwC75fTEFv+UxAI/6hDzvEiFgjvlnFLNGtimpdrG4rBmIEtgmk8o6FFm0sAGTZ4ZfmzC7hxT1kEgS8T7I5j9FIEuM3l8mI3aaBRZznzQUvyv/hMAj6jqW+afkmEYRrMoljeS8GUAfgbAJ0Tk467tTa7ybRIdarbKdOAJMKwq1Q2DlKRGvDgChdyr2qZDBU3XquvlwT7KLB45HN6P6XP5LFbIw8q8JfCZpvO8mg6+OZDduF/TRh14gZCNmt0/XIkujZTnGm+nAXRdU9PGbCylglbVDwEL+MoNw1gaVGVpFbRhGEarKTYJLdR7dqYJWpl1czBGtKJKJviETCwaHFc3FAOzgHcxm8YNL7lh6MxAK2TC6JPZYMVtxvXpAutkamDTx8iVkPp5E9EHsgQuheXY/CYiv57JYePxLHoJlzwzdxi1sZqEhmEYraTYJFw8S+yWLtAz5YAOzk/kOg6O3T1m/eMZCT7RxIbhaJOQg1p6VTXNqjr4zOTeAhaSXoly/qR+dVO1s07vUY/VNitr194vvwIGarnj1DKralbguU1Efgl1w8YjwS1Fc0RNm5I2arD0kYSGYRhtZOkjCRvFKaK0bXIGY38qfLquPbuTOo4oZFLNw4i9OrRR01Sm+BYRU9DDAbm2uaAVZRs1qekO26udsha2UdMcxStsds1bX6fJ+MrpiUAX4XkV90iGjTtpnQsVL5ojiZeC65qyNkqsaKxhGEYLUQXWI1kX285iLtDDiF05oYrVhxbn/npO4xkicQWMmI2ZlOhwJaagq9cKrpkS1epThNJQtun7OJXAs6O8MLd31oprdUghC6VcFefdIWulDTq0Ubv2TlxhB6rXt/FL4QAXd99RAiYgUNixABcrDmDkKEwctkAbhmG0kqWMJGyKeT06AmIhzckUnnP6VEfCwlMK2h8PA1VN0xopcG7L3D7w7KjaoDv9so1Vs9Lx0NnMu+zxwfbqNecRwiH267w/ECl6wL/DPqU29e8RzStQ4yMbNanqTBImCxU3cjTpZicilwF4K4rNsFtU9YaxfnH9VwA4BuB1qvox6u8CuB/AF1X11ZPutXia3zAMY2oKE0fukb1KsbjeCOByABcAuEpELhgbdjmAPe6xD8BNY/2/COCROrO2BdowjKWgoZqEFwM4pKqPquoagNsB7B0bsxfAbVpwL4DTROQsABCRcwD8CwC31LnZ1pk4dFgjpDsTDkxtXJGaN7hGJgy2BbAJZFQdO2H2mGbz0FtT6Pccc6Njs8aQA1nctIcJN7zg8+OPeXpcDLxa1Dtp7uisOZe8depf4wActxnXi/d7E0WwiRi43EVc8mjagUueb+KXxS8iEuBioeJGjsKLo5Z77hkicj893+9y2XvOBvAFen4YwCVj14iNORtFacA/BPArAE6tM5nW2KANwzA2iikCVZ5Q1Ysm9McuMq7iomNE5NUAHlfVB0Tk5XUm08oFuswHPYVRPxOIohzsEMjaaj29WKVtTlqUVdWBgo650VVVM1Cq6SH9Vri/lvudY+A2DDusODlXEkdtO/e/4Tqp6pXyuOv6vTseAHRpQ7DTiWwicqALbx6ecGqaFTa75LmfseAWYOx/gv+WMGv1FnPDWypqmjByHAZwLj0/B8BjNcf8BIDXiMgVAHYC2C0if6aq/yp1M7NBG4ax7fFeHLlHDe4DsEdEzheRVQBXAjgwNuYAgNdKwaUAnlbVI6r6RlU9R1XPc+f9r0mLM7DVCnqqFKOT3as4/DuaRIkVW0xtJxS4RPp1GjXNQyOBKIHLnfttBK53CQWtEXs3403uQSHvoKg32ZOdgB2QDbpLCnrojruBvZzG9qoh7p0OK+yIbZrc9ML0rsVYIXs5v8OBbXp0EK/eUjvxktmll4ImAlVUtS8i1wG4E8Wic6uqPiwi17j+mwEcROFidwiFm93Vs96vlSYOwzCMJlEV9BuKJHTl/A6Otd1Mxwrg2sw1Pgjgg7l7tXqBnip4hW2LsRDwVHCKl5hc8ZqVea96PifOD2sGTp6il3+pQBSvlkO7dPx4NDZjl5akgi6PvYJmu/RwlZT9CfdzhdV+NYS9y6o6SLlKatqNEbZRx1LFSrw/CHDxbZFQ8aK/XuIl8/JYDiybnWEYRguxhP3zMoM9uk46Um+PDlRYzKMjsEGz3633uea8nnwcs1dPnn5AxGc6sFGzgmbbdK86NqqmWRzyccS7o0PuynzsFfJgrWzr0lx6xyOlvlhNR8p+dYLyWRG1LIn+srX0mY6Eigf99MKVC4zX9fKgscbiYgu0YRhGC7GE/YZhGC2mIT/oTaWdC7T/OklfMefeMExtEnoThsQrXsOFOfvscDw9APFK20Fe5nhe43JA9Ti2cQiUZg2AXfI0OnZ03ahfWnzDsENudsIbhs600TlBbZENw2GvvEGPzR296jFvKHakelwnSKk0h8QDWUorD2XTY/MXmzt8fyozngW1LDSqQN8S9huGYbQTM3FsNrlkSgD8hmEQvMIbSU6pBeHfvAnoxw7IfYvduvp0X1exJOfaxpt1sQ3FdLIlmuKKrzlI/V2tHqdEA9/XBa2w2mc1PfQVV1bp/qym/YYlu94Fm5vlJLyy5nD3YBPRu+FxoAu74fFYqf6Hi24iBm0RNV3jc2RBLYuN2aANwzBajNoC3TAJ17upkinl6hf6fg5U6UZs1IGqjqculb6vYlK2BZW2RzZqlsiZ+afU9CiohVTzKik+Zw8Wfi2pt8uHhdNcOSBkuOarfleDU4pj95M+Tb1EUIufb849MFTN1J0IYCnHkgL210+66Tk1nanYwu3mhre4LOImYSNWcxG5VUQeF5GHmrieYRhGk6g2lixpU2lKQb8dwNsA3NbQ9WoRq/BctE+R3N/XwCM1HtiYvfoihS2spkkty0r1Xuz9IcNqG6cD9SlCJaWqYx4frKDJi0J6xc06K+X8Ol1OIlVeqoyGJ9VLtnXd4QJROIk/2ZV9WHiglLkoAXuijOoydiptfMxeIGx3jinrQEmvUTSNI5d4qZH6h+bl0XIEg2X14lDVe0TkvCauZRiGsRGYDXoCIrIPRQFF7MTJ019g3tSk7CMcSy06yLhe9MkxOKheTYrLqekOVcRmpdkZ2ajZLk2X8ibRxFSi8o/ts6SmvXLurZTqsNcrj7sdUtsRyd4fdOi4ePMG6+WbONhBxy7Ue7hKqpo9N2JJngLVzP1OrUdKhRXzLtu7scRKlVcSMo2XR0xNm5fHYmK5ODK4ul77AWC3nD5NtgrDMIz5UEyTur01tNuLwzAMoyEW0YtjMRfoWULB+etmZMMwyHZHm4SjenqBm11iw3C9OE8o/7Gs8yaiq93HFbXZ3BGpxB0cZwJc+PPnNwTZrLFK5o7Vbnm84o6DrG/0dXDdmThO9EtbxIm1chdwfbX4GA2Ol/39HoVU84ZitO4ivfc+R3On2gaMRcu731kQKh50V00gTCThX9bcYW54i4ku6CZhU2527wbwYQDfLiKHReTnmriuYRhGU6jmH22jKS+Oq5q4zoaQyx3tJVmqKngsmVJmw1DI9S5Uyy5II9JWHPtbVVV15djfIvGh8qKRNwNZNZ+8Urqj7egVN+6RXO+QmvalgtYo3P34jlJBf32tiAF/brVsW1spj/sr5cfMV1fRwM2unPeojGCgoHknFHTsEi8lvrmOEi/Fu0fnJ/JJjVVy6YcTxIxueEVHakbGBmJeHIZhGC2kUMi2QG8umVDworleatIwiQ7dIpb+MghUoYrVXulRxWq2R3fWqgEdnTV2jXM21URlk6iaJrWNyJcAtiuvkIL2qhkATnFqeme3vFknIs15k+X4oFTIz66uup87Rm1HV8rjYytllqW+U9bDLqvqaqAKEm52yr9n/+VHEh/jTuVSM7nhARillZ3bDa/o8I2Z2RhNYm52hmEYLaWNNuYc22eBziZWajAUPFaFGgC8QiVvB1kvFVfH2V27HMjCCYjWi7kOqI0rbYfKupoiNEh2FFELrKbZ3uyV865uaZc+iY47kUCWE+SGsdup4mdXdo7anqLjp1dOGh0/0yvaT1DYeb9XqvFRoErg2QE6jtimg98HTTLivTGvmp7Xy6Not6CWzUYhGC6gF8f2WaANwzAmsIACekkX6FgoeMyzA8iHgner9uigyjQpQekWypp9ozuUgKjjkvB3qYyVT0QEjNmj+74tkSLUJTvikG1W1ayKvb2ZVfMp3TIj/w4n41dIzg9J1h4fegVd2p1PIdV8CnmM7FwpXsRT1P91SvLUd+XCtUN+1BJX0yOFzL+uIKZ/Mv5SdSyTtX2mA5UWqamFcj/DvDw2EdskNAzDaDELKKFtgTYMYykwBd0WIhuGs1YFz4aC9yMbhpTkWdjlzm0isllDyUTSdW52XKl7mKqkveZ/cls1h3OfgkvW6biv1Q0TNnvsIHPG83rHAAA7Zb1yDgAMnLHgxHCFzilNGLt7x0fHJ/eKie/oltf/Kr0HR7vFJuJ6l1zzyNwRmAJGbnap32c9c0etjUNn6kq64Y3aalQQd59P/hyaG97GogCGQ1ugDcMw2odiLJHLYrD9F+hsYqVMKLhTPsnNIXa/85uDiVBwWSvaO+xCFgStuNBlro5NqnlAxUI6TllzpW2uGeirn3AOZ052xGHbsSxfvCHolfOp3VIJr0jZ30Xxfq1r+XF6XvcYHT83OvabjydRUMxOUtNfdu/9UxSivkYvcp2DUvy3o2D6scRK028chleKM40bXlRNk0KOqmnbOGwU84M2DMNoK7ZALzAR1zuA1HTM9Q4Yc79z0ojd7NZLpRgLBWd7dDcSpMFh4V2u8+ft1QkbtU/xyZVPOEUoJzvyYdscfDKM2KhZNe/uPEftVXey0ygD0qmkoE920v9kculbJbXuEzpxWPpXZVc5V/o1eDWtEbt0cRzTwHOq6YQM8wFDSv1JNa3ezS6ukM0NbyOQhdwkXLzQGsMwjFnQGo8aiMhlIvJpETkkIm+I9IuI/JHrf1BEvsu17xSRj4rI34nIwyLyW7l7LY+CzoaCS3TsKJCFk9mzGop4d0Q9O4BR3LbE7NYAOu44ZYPuHZ+c+H7ICtsFuAzoHJ9YHyhThAJlsqPdlCL0OF3Me2l4WzMQquZdslZpY07tlLZrr7zZnr2DInC89wiHondpT+AJOj4G5/Eh9MJjmoOUUxi1HlHTrJC5XqT34sgYMnNeHtwu9DlKqen4TczLY2oU0Aa8OESkC+BGAD8E4DCA+0TkgKp+koZdDmCPe1wC4Cb38wSAf66qz4rICoAPicj/VNV7U/czBW0YxpIgNR5ZLgZwSFUfVdU1ALcD2Ds2Zi+A27TgXgCnichZ7vmzbsyKe0z8i788CpqZxU86VTJrUE1Tmkys5NVyxLMDKL072LOjG1S/Ztu0VsYGatorbKq0zSWpOLm+TxPKyY44bNv7N7OXBuOV8y6yUa+QVB2gbN/lPEJOFgr/Do7Xg2sCQI/8yjkN6uPuZ+kvAqwjoqaD/wKRJFJBWfFKd0DSy8N/fnrDaH9wWZ+utktJtSJqOpmu1JiNem/hGSJyPz3f7wpee84G8AV6fhiFOkZmzNkAjjgF/gCAFwO4UVU/Mmkyy7lAG4axfNRboJ9Q1Ysm9Mdk9viVk2NUdQDgQhE5DcAdIvISVX0odTNboOf1k54msZL36EgpbKeQO52qZwcQKuhet2qD7gVeHMXPIRVs5UKuXJLKJ9fnFKGc7MhHBbJvM3tpeFg170wUct3lvDNO7pTXWiW17JXzSqQtxeN0zGq6794EIY8UtkHLKOM/nRSM7dGx/xkvjZazTQdeGsPiveNvX/z5q5383zw76tNcoMphAOfS83MAPDbtGFV9SkQ+COAyAMkF2mzQhmEsBQ0Vjb0PwB4ROV9EVgFcCeDA2JgDAF7rvDkuBfC0qh4RkTOdcoaInATglQA+NelmpqANw1gOGvDiUNW+iFwH4E4UbkC3qurDInKN678ZwEEAVwA4hOJL3dXu9LMAvMPZoTsA3qOq7590P1ug65IKZOHNnVEbnUZf9Uebj33+akuBLM6c0aHADaWahl0OB/abgCtV1zsAcN+i0aPfsLK5gypt+5qBTydyOPtkRxyyzcEn3o2ONwP5q9lOytHcce/ODjL9rIBDyIv3piuzfWUPzR3FvYLNzaA8i/9Jv6Nh3BwyKj1OMqsbccmrYwKJh4UT7uORq85igSzTESkMNBOqehDFIsxtN9OxArg2ct6DAF46zb1sgTYMY/szRSBKm7AF2jNjIEsQ7ODDfSM1DQFSTqmahifWK/2cWAkRBd0L3PBIpXmXPZ4eV8/usZouNtN8vUCgrHwClClCucqKD9kGyuCTXfRtYFeXXePK+644Nb1C79sKbYquOle9Lp4etXVX6ytCTvz0uFOdx+k/Zj+ipiXY0+WglqraDjYZc4EsCde4kTsmb0ZT/6g12EMkNT1S2Jl0pUVHdA7LhzS1Sbip2AJtGMZyYAp6mzBFIEvMJY+VTzSpO4WCh5ky3bMgyT+p6cD9zinkoMp1ealhNPESzYvc94bdooMrbXPNQJ9cn1OEcrIjH7bNwSfsRhfam52CJrs0H3dHKoeLAxwdHQ1Wq45H6/Q1oT8sjwdOQX+FlNMaKeS+O5bABg06rtqmuYq6DPnrifdrLJs6KbeAiNoOPidaTbwU1Dr0Mj6XrrToqIxdWhbwLbAF2jCM7c8yJ+wXkcsAvBWFQfYWVb2hieu2gpqBLEX75Arho2RKnI50mkAWtk27425CQcds0BrYq6nd2ab7vTJ4hStt+5JUnFifU4T6ZEccsh0En7CXhnvtK4FnRznxHS7x0fMCoVyq6QGeKY+dVwor6EBNuzeES319jVTx2sApaFLCMqD3O6KmAy+PINNsRBUPSeGyqvWfKfby4DJr0fJa9dOVRjG7dGNeHJvJ3IEqlN3pcgAXALhKRC6Y97qGYRiN0lC60c2kCQU9yu4EACLiszt9cuJZ24RsWDgbpH1ynIhnRzHWQeqUVbNEFDTYy4MdTUYKmpQ/qe1hVGjQwUoAAB1cSURBVG2TZwdVB/CFXL9M3wxWyUvDpwjlorKpUO2RlwZ9s9hB6UK9mua2U0j6nE6pSde7hW16faVUwMe1PM8XIODyXn164U8PXFEDUs0yoN8dtfuXHtilg7HeTxolg/JJJ+YTHdiYq0mWogmWgNHnKJWu1JIsbR+aCPVOZW4KEJF9InK/iNy/jhPj3YZhGBuKaP7RNppQ0HWyO8Gl7NsPALvl9Ba+FYZhbFsUjYR6bzZNLNB1sjstPolAltzY4Kun+xMdVnum05zZIlmRJch45sZyRRYa2xu54dH1gw1FDoBxbbyJ2KGNt24RCs6VtrlmoK9+kjJrcNh2GYBSmip4Q3BHUB3FX4vqG9LXd/9NbE1LN7zjK2WlGJ/Heo38C48PyuO+M1Ecpfe73y/P59+DfzkyqLYVx7Ggl0S2u1xQy8gEwhuSXBV8FKlS9ltYeJ4FlIVNmDjqZHcyDMPYUpbSxJHK7jT3zNrMrGHhfmOMK7JEwsKDjcNUWLijI1Ul7DqKayVUs3aqm4eB2qbr9p2aXuuU6pIrbfuagVz5JJXDuQzbPkqtpZr2G4KsmpkVepG73Hu3m8LOv7FbXveoC10/tlrO+7mTaBPRqem1dVbYVIGG1LZXyMH+LW8iRhU0u+/RhqHb8OO2IId4zM0uEhYeS7AEzBgWvgxKuoULcI5G/KBj2Z0MwzBaxbIu0EvNDGHhQSBLJCw8UE6psHBPIqilDGSha/EFWOT7eSeCXvzrWhdSmnS+r7TdmeJ/AIdsc/CJd6NjWzOr5iHdw7eeTMr9VC2DYk7vFfU5jw3L+orPDsrjYzsLZX18vVTV/T5VQeeQe/d76Efs0sWxBD+BUlUX51fDwllBhzbqSCBLLF0pt/H5/OUlopCjanqb26XbasLIYQu0YRjLwZJ6cRieeesbOuUTDV4Zb/dqORHI4uVlN2XDpuNRab4gBD0SAEOvi9X0MRS2Xk6WPw0DKiQwCj4hX/ldZNhlYe9bu/TO7KTK4j4NKtdS/AZKAvXMajHvZ3eWNurjZI8+uk6Jl9adQqY2tkd7Ec8KuxMo8GqSJemX9+pwoMrAKfpAFZPnjwtUYdWc/Jy4D9VUYeHbFFPQhmEYbcUWaANADS+PyT6rsQRLRT8r3MmVrmMKmz0+uoENulcZG9qoZXxaI7s1AKw732WuqJ1T06kERz5sm32b2UuD7c1eOa8nspR5T5JddD6X7TptpTh+drUsVHCM1PQJUtMnnBoesM80K2Sf52qdfarLuQT2aJ9YiULQhXPBRhIv5fykg1DwiO/80vtJmw3aMAyjxdgCbRiG0U5mrEO8pdgCvQVk80nHMuAB0Sx4wZd7NoGsrVX6s2GjgQkjFigSN4H4K6+jdFdLmTt8zUCufMImDp+NjkO2OfiE3eh4QzB2Lc+KVHNXA2WNxVNWymueslq64X19RzmH/s7iv8pgjdwDTyLz0cgNr7wvh4p3Avc8F/TCJg7aROx4EwebLShoxpszwtzTiQ1DP4Y3CVE1j1l9w3ZiC7RhGMuBmTiMCjMkWUrlk47WN4y0pYjmkwZKaR0JJS9wbl1BoEtsLOWTVlbT5VhfaXsw5LHlsc/h7BMdAWXINlAGnwClG10qrNyr6QHNq0tlUnY6Nc21Fk9ZKavCPLtaHh9fLeYz2Fmq3iGr6bXiNZJAD9zwYu53nXVSvX3eMPSfAwo1Z1U79P6YkyuyAOSSF6nIUtyjZiBL0VEZuzDYJqFhGEaLsQXamEg2kEXCcYinKwVKRRRLVwqE9k+6QHlIzfVTGsYTGJWRLnQrUsXrSiHibgxX2uaagb76CacI5WRHHLbtA1DYjY7tzV458/2Zjts12kFSdyep6Z298njHanG8vlpea7CDFPS6V9CklEsBjs4OobGurV91vSuOfSALKVU6llEFcXrDe5xsiY69HTuZtrb6rW0hFXIdbIE2DMNoHwLz4jDqMmO6UuXkPIjYFvkWkbYUfsx0ycHJZjpSJhJpQxBQ0Xdqdo3yInClbV8zkBPrc4pQTnbkw7Y5+IS9NNje7OGahevD6sefEz4FdRd7xfFzK2XbYJWS6Ltw9eEq1X0sp4oh2aYH7gsBq2pW00MXTi6rVc8OgLw7Ip4dQT/oM1WOjKYpzQZP0bUW0rPDbNCGYRgtxhZoY2qmSFcaU9NCSjamppNeHhEvjNns0kDMNi1BXDh1ODXdJ9W8Rik6faXtPoVGnyClyKHYPtmRD9kGSt9moPTS6NB3W1bN3mNkmAgV79F5vsRXj1T1Otl91RUiYAU9oGOqdTBSzoHHB9uunZeHJNKVqrdRc5J/8iuPeXdoKnWpP56cOaA4zxd2WFQ/aVugDcMw2omZOAzDMNqKLdDGXNR1w4udA4S1DmNueHScDWqh47rmDuEw68CsQZuH/iUGbZR32Zk7uNJ2UDOQqp/4PM6cjY7Dtn0ACrvRMd608dygtD8M6ZWz6cO7OPa6lJuaTBz9XtE/XCHzQY+uRcXK/TGbQLpk4hh6lz3aOFTKme1NHIhtHAJhfUNfUieX+a5GhfDWmzAmoc15cYjIZQDeisK2d4uq3jDWL67/ChRZD16nqh8TkXMB3AbgBSjSme9X1bdOulcTVb0NwzDaj9Z4ZBCRLoAbAVwO4AIAV4nIBWPDLgewxz32AbjJtfcB/LKq/mMAlwK4NnJugCnoNpJxwyuaM9VZXHvODW+8PUbWDW+Un5jbaF48R7c5GFS/pk1CvzHW71OVE9os45qBvvoJbxxysiMfts3BJ7G6iayajw+oPiG9Bo1sJHLgkHTdZhwLVVLTwxVSyD3fVo4dUH/HHQup5uA4GshC31561Q3DZGIlrW4iIulmF8stvTih4A3ZoC8GcEhVHwUAEbkdwF4An6QxewHcpsWbeq+InCYiZ6nqEQBHAEBVj4rIIwDOHjs3wBS0YRjLQT0FfYaI3E+PfWNXORvAF+j5Ydc21RgROQ/ASwF8ZNKUTUG3nXmrs0Ts0kCYOKluUEvSLh0RTELh1axcfAi4cAHPiJrmUPUBudlxpW1fM5Arn3CKUJ/siEO2OfjEu9GxrTlI3EQ1A304ugZ26XLeozenS0qVytaweX4UqR2xSwOlq56yXZqCYjouNamupELBuTqLa6dUtEIpbL1yDr9dkc9dxv1uYSqE1zRhAHhCVS+a0B/7bzJ+5YljROQUAO8F8HpVfSYydoQt0IZhbHsEjZk4DgM4l56fA+CxumNEZAXF4vxOVf3L3M1sgV4kskEtk+2FgZqmzPIjOyTdKpAAWv1kB2o61h/Yo6veHTG7NAD46XIG0TBgg5MOFccnSFX7xPpAmSLUJzoCypBtoAw+CZJQZZI4DWiukZcdvDFsjw7egm61bdiL2KipTYJ+9y2jR+9Lj43fXOTBX4yTKdHEvf0/Wf19+yRWamiBvg/AHhE5H8AXAVwJ4KfGxhwAcJ2zT18C4GlVPeK8O/4EwCOq+pY6N7MF2jCM5aCBBVpV+yJyHYA7UbjZ3aqqD4vINa7/ZgAHUbjYHULhZne1O/1lAH4GwCdE5OOu7U2qejB1P1ugF5Wsz3Q9Lw9gRp/pSOgw26WDitTBWPeT5xrz8iDPjlTie3E2aK60HZSkcsn1OUUoJzvqjcK34yqQ1bRXzn3ylhgOWElGLsBqOqKsUwo7pqA7dOx9opW+OegKXYDt0d4POqgIT++hT+ifCwUH6DMTN0y3PhS8oUAVt6AeHGu7mY4VwLWR8z6E/FZPwFxeHCLykyLysIgMRWSSYd0wDGPrcNnsco+2Ma+b3UMAfhzAPQ3MxTAMY+NoIFBls5nLxKGqjwChy5axyczrhgfUDmqZxg2PzRrdILua/9mttAHlhmC4SciBLGW7N3306as+V9r2NQODyifkruaz0XHIdhB8wmYJ18xmjSGbGLx7IH+LT9iJYiaOIBmda2dvuWATsRsJZOGN1NiGYSJ3dBlklAkFLwYVTYFr3RQmjC0OZLGE/RNwDt/7AGAnTt6s2xqGYQBopwkjR3aBFpG7UST3GOd6VX1f3Rup6n4A+wFgt5y+gG/VAjCDG17QngpqGbXR5YeRX2EkhLhoLw9HLnkUmRGtx0eKsMPqkDcPfU7soJI2bxgWx4P16gYbgFEOZ5/oCMAoZLt4giq898nKvl/d3AxD3yOXSkT+eIUcuOZ1IgEwrGRpE1BZQbvQ+NFmIRBuGHaruaWVAlmiG4azhoJvJS01YeTILtCq+srNmIhhGMaGsh0XaGNBqemGV7RnglpibnhsV/YHvbiRT2IuXMN4fxnIQv2JxEre5Y6rkVD8zahiSScIma4eBylCeds8o6Zj80rZy5O26RhSnYuyWnYKehgEqtAN1lkhe7XNbRF7czdhS44FsPB7tCCh4A1GEm4q87rZ/ZiIHAbwvQD+SkTubGZahmEYzSJDzT7axrxeHHcAuKOhuRgbwYypSxHYoOesIB7Yo11QS5D+kpPNO0+BIJCFEgWxbXrdpyYtT+dQcO/l0Vmj21OlbZ8wP0ysn0h25G+biubx4pC9TEi5x74RJBWdV9Bszo54eSQVdmCDdjZmDsbpcii4nwwngaLPScQerQP+arAgoeDb1QZtGIaxHVhEE4ct0MvEFD7TwWkzVBDne0k0qxBNhdSyV9aB3Zq9PCjUWlzazUBV9yM26h1cTqq8r6+0Hab9jKcIHSU4iqUYBfK2c1b5g2p/zuMj9KOW4CeA0KOD1bBX211ui9ijE1XBo/ZotkvPYI/eslBwW6ANwzDaiSlowzCMtmILtLEwZINapDJWeW+IvuuPzqLNPt7Yi0Y8J75Gj0wb1N/hTUQOcBn6IAs2h1Q3EQdk1qAC3iPTR2Di4Ihoah+ZOIKNOVThb++BGx7dtx/pn8YNb3QSnSJxc8fIvS7lZudrPCbyQQcbhu4aQfBKZMNw5qrgG2nuUAv1NgzDaCWL6gdtC7RRI6hFwnEYCxX3u0Nc626QcMkbVjcBETlObhIGG4Y918ZueHS87jYRqRI4bxj6QJYBB68k6gTmFHSk6HdYizHmZhckhIr0T6Oqk/NyATSpiimjTUIOXuGw8MwmYWzDcJbglc0gs1ndRmyBNgxjKTAFbSw2WTe8ql26OHS2R8QDGGI26GTiJW9v5rZB3GXP26a5irV3vSuOqzbqQE27QJIuBZQMVhJquldNYMSeadmY3IiajtmlAaAzCtapnlOLiKrVpI05V5MwYtsmiR61R0eCVwD61lXHvhzYo/PDs1igimEYRnuxTUJj+zCLl0eqIAApxVFYeKLuXdzLI+HxMXASlxUy26i9wia79HCdbdQd10Z2aVbQbJvuhT8BjCXc90q0bMuGhXMa1sDO7vurhQ6C9jqKMJbgqBPpZ8Wa8ejQwLgeG1tfqm6mPdoWaMMwjDaisE1CY5uS8fLw1Cqv5dRhVDWDtBfbOYO5ZDw+hjwH5+XBVa5Xqz7TbJfmMlLKaUp7EvwsjmkqMS+PTI2wZFi4Tyi1Xm3jscGm18w26mry/7A/4tHBLyzYdnAKO1oJHFueTMk2CQ3DMNqKLdCGYRjtwwJVjO1Pwg2v7J5csSVoT9Q/xMh1joJeUrUOvVtX4JJHLne+nUwcnUhmPK5TyCaOIVUAF2fa6PSqZg8Ao5qBYeWTclq5+obRQBZ+3ZFgHUm9L/PaWtncMUyYPkZj6XMwyESoZNjQDUNtZ0L+HLZAG4axHCze+mwLtDEjs6ppL7KE1R+HHjtVmwoVp2PxoeXBphQnVnK1FIMESzyXbvATGAsbp6CXUaAKK2yuuDKqxB2pwoKx3M0ZvNILky1NVtBIuOTVVtMxdZwcSt9+Yi53M+aL3mjMxGEYhtFGFKEpbEGwBdqYn5ir1LyJl1Kh4rEAl0To8Ug9pmzUPoQ8sEvTvFbIPc/V+dP1uL3a22250nZSQefEqg9kCWzUWjkOwtoTanrkfjerG17qeBbmTUc6L4u3PtsCbRjGcmAmDsPwzJt4KRUqHgtwCYJTIqqyx3ZpPnbJlLiNah6y98dIQa+wvTqSzJ7ula0ZmEuwxG9hpDJ60ssjKJwwrJ7fpJdHUHo8k9B/CoLPR0PC2rw4DMMw2ohlszOMBJnES0wsHDiXeClI6MO3Hdmoq4rSDXDXpLYeKeSg/pVrj6hqAFB/vB5pA6IlqYKSWTn7bhDO7ueUUM0Dtrn79yBil+ZrNKmwGUn54MTGbpw9ughUaeZ1ichlAN4KoAvgFlW9YaxfXP8VAI4BeJ2qfsz13Qrg1QAeV9WX5O6V+5JlGIaxPRjWeGQQkS6AGwFcDuACAFeJyAVjwy4HsMc99gG4ifreDuCyulO2BdowjKVAVLOPGlwM4JCqPqqqawBuB7B3bMxeALdpwb0AThORswBAVe8B8GTdOZuJw9hcIpnxwm6/iRgJbhk7L+qSx6YAH8iS20TkKtdBKDmbVlxmPB7LJVV8RjwKVOENx+Dl+qxvqYCOuuaOIDyc5xoxfbAJZKPMGTn8xl/CfLWh+aDr26DPEJH76fl+Vd1Pz88G8AV6fhjAJWPXiI05G8CRutP1zLVAi8jvAfgRAGsAPgPgalV9ap5rGoZhNE/tXBxPqOpFE/ozWVVqj6nFvAr6LgBvVNW+iPwOgDcC+NU5r2ksA1OFik/hkhdLwpTdREwEunAiIK9Ag+RAFPTilTX1S5fmysrau+RJxPUOqG94TLnhRULfw43DyEZpKtBlFoU9b0DLRtHMt4XDAM6l5+cAeGyGMbWYywatqh9QVV/u8l43EcMwjHahhfdL7lGD+wDsEZHzRWQVwJUADoyNOQDgtVJwKYCnVXVq8wbQrA36XwP48wavZywLs4SK03lBuDD54Y2SMOXCxlM26i7ZmP1YVods9/XKmW3UXbZBc3uneq2Ygp5GiSbc6CTmShhzNdzEyiZbRgMK2lkLrgNwJwo3u1tV9WERucb13wzgIAoXu0Mo3Oyu9ueLyLsBvByFrfswgN9Q1T9J3S+7QIvI3QBeEOm6XlXf58ZcD6AP4J0TrrMPhcsJduLk3G0NwzCapaE9SFU9iGIR5rab6VgBXJs496pp7pVdoFX1lZP6ReRnUThe/6Bq+k+U2wndDwC75fQFjOkxNpVMqDgT2KCDS9QMG2cbdaLaOAbFmKiqBkq1HHhT0DEr635VbUtMTU+RojQgZkNO2ah9KHbMLp063oyQ6YSXz1yXHC7et4R5vTguQ7Ep+AOqeqyZKRmGYTSMorGcHpvJvDbotwHYAeAupwDuVdVr5p6VYTBzenxkw8ZTftR8E6dqA1XN3h9egbJSDuzV3Wr7MDHWzyeloKcICy/vz/OOqPyMF0iSBUlAJKgdiNIq5lqgVfXFTU3EMAxjQ1m2BdowDGNhsAXaMDaYBqu3lOMygS7AyKwgia//OqqPSCHVQXa2SO2+pIkjVtsvYtZImTpyJo7Y2NjG4Xi7O9acS96MC+HGVvXGUtqgDcMwFoKl8+IwjFYwp0tendzTozzWMSUMUtYcvh245EWUNd8r6mbHoeZTKOgYKVUbCdZJbxhmVG2rF0A1E4dhGEYrUdgCbRhbziwueRkbdTHWp+3k/khQSyLoJVDWnUioN9lfswp7Utu0xBathGoe2Z7nTaa0VbRZ4CewBdowjKVg6fygDaPVTOHxUXZPDhsPkjEFlae96qVzBlUvECBur+agFxXvWRFXyBILg24iLNzfP+WlEatfGLFda1tt0bZAG4ZhtBDVMN/IgmALtLFcNFAooGyKeH8kynNFlXXOXs3X6kQUdjCZOe3RdfyOYx4fdc+pcY8NLXk1PpcFwRZowzCWA1ugDcMwWohiYRI7MbZAG8vLTJuIk+sjhmMnB73kNhSDKiyJDcfyAs3nTwYQN2cklOhoc7CVSlUXsmqMLdCGYWx/FLZJaBgLzxSbiEwuhDw6NrWhCJ+UKKFEIwo66+PbRFi4705tEtZN0rRVSraVyn4ytkAbhrEc2AJtGNuIjI06HKque7JrXtGcsVd7EgockeZUUMyGM1Vq09i3jM1aNC1ZkmEYRjtRtDzbXhxboA1jGma2UdcPeinPSSwokVqKY1crD2cNAa9LTgHPam/eCDu1KWjDMIw2YqHehrFczGCjDodOr6rD82ssOJlrbAgZ9bt5dme+aY1SXS3EFmjDMJYDiyQ0DMNoKWaDNowlJ/U1epYNxaypIG++qGUG2UCmMmdspAlC1bw4DMMwWospaMMwomTc88KhmQ3F2DXDwROvlSN6rwRzb/ht2sadhsE8C4It0IZhbH8s3ahhGLWYwj2vPGXy4lInxDxx4tT3mpmtdnPb6vvPwFxOkiLyH0TkQRH5uIh8QERe2NTEDMMwmkJR/OHJPeogIpeJyKdF5JCIvCHSLyLyR67/QRH5rrrnjjOvF/vvqeo/UdULAbwfwK/PeT3DWE50WH1MdfqMC07svhv1mPH1NIJqI3MUkS6AGwFcDuACAFeJyAVjwy4HsMc99gG4aYpzA+YycajqM/R0F4IkAIZhGO2hoU3CiwEcUtVHAUBEbgewF8AnacxeALdpkdD7XhE5TUTOAnBejXMD5rZBi8hvA3gtgKcBvGLCuH0o/poAwIm79S8emvfem8wZAJ7Y6klMwaLNF7A5hzQhd6qicBHf42+f9wJH8bU779a/OKPG0J0icj8936+q++n52QC+QM8PA7hk7BqxMWfXPDcgu0CLyN0AXhDpul5V36eq1wO4XkTeCOA6AL8Ru457kfvdNe9X1Yty924TizbnRZsvYHPeDBZtvkAx53mvoaqXNTEXADEfxPE/pakxdc4NyC7QqvrK3BjHuwD8FRILtGEYxjbgMIBz6fk5AB6rOWa1xrkB83px7KGnrwHwqXmuZxiG0XLuA7BHRM4XkVUAVwI4MDbmAIDXOm+OSwE8rapHap4bMK8N+gYR+XYUlq7PA7im5nn780Nax6LNedHmC9icN4NFmy/Qojmral9ErgNwJ4AugFtV9WERucb13wzgIIArABwCcAzA1ZPOnXQ/SVYONgzDMLaULcjmbRiGYdTBFmjDMIyWsmkLtIj8pIg8LCJDEbmI2s8TkedcuPjHReTmzZpTjtScXd8bXbjmp0XkVVs1x0mIyG+KyBfpvb1iq+cUY9rw1zYgIp8TkU+493VuN7CNQERuFZHHReQhajtdRO4Skb93P79hK+c4TmLOC/E53gg2U0E/BODHAdwT6fuMql7oHnU3GjeD6JxdeOaVAL4DwGUA/osL42wjf0Dv7cGtnsw4s4S/tohXuPe1rX7Fb0fx+WTeAOBvVHUPgL9xz9vE21GdM9Dyz/FGsWkLtKo+oqqf3qz7NcGEOe8FcLuqnlDVz6LYrb14c2e3bRiFzqrqGgAf/mrMiareA+DJsea9AN7hjt8B4Ec3dVIZEnNeWtpigz5fRP5WRP6PiHzfVk+mBqlQzjZyncuodWvbvs46Fum9ZBTAB0TkAZfGYFF4vvPJhfv5TVs8n7q0/XO8ITS6QIvI3SLyUOQxSREdAfDNqvpSAL8E4F0isrvJeU1ixjlPHbK5UWTmfxOAbwVwIYr3+fe3Yo4ZWvNeTsnLVPW7UJhmrhWR79/qCW1jFuFzvCE0mrB/irBwPucEgBPu+AER+QyAbwOwKRsvs8wZ9cI9N4W68xeRP0aRErZttOa9nAZVfcz9fFxE7kBhqontr7SNL4vIWap6xGVYe3yrJ5RDVb/sj1v8Od4QttzEISJn+g02EfkWFDlUH93aWWU5AOBKEdkhIuejmPNHt3hOFdx/QM+Podj0bBtTh79uNSKyS0RO9ccAfhjtfG9jHADws+74ZwG8bwvnUosF+RxvCJtW8kpEfgzAfwZwJoC/EpGPq+qrAHw/gDeLSB/AAMA1qtqKTYLUnF1o53tQ5HHtA7hWVdtYkfJ3ReRCFCaDzwH4+a2dTpVZwl9bwPMB3CEiQPF/6F2q+tdbO6UqIvJuAC8HcIaIHEaRyOwGAO8RkZ8D8A8AfnLrZlglMeeXt/1zvFFYqLdhGEZL2XITh2EYhhHHFmjDMIyWYgu0YRhGS7EF2jAMo6XYAm0YhtFSbIE2DMNoKbZAG4ZhtJT/D7eCvfUpcS0KAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.pcolormesh(x1, x2, prob)  \n",
    "plt.colorbar()\n",
    "plt.title(\"Theoretical PDF\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Use the same values for $ \\vec{\\mu} $ and $ \\Sigma $ as above to compare the estimated PDF and the theoretical PDF. To generate a sample of a bivariate RV with $ \\vec{\\mu} $ and $ \\Sigma $ proceed as follows:\n",
    "    1. Calculate eigenvectors $ \\vec{e}_i $ and eigenvalues $ \\lambda_{i} $,(i= 1,2).\n",
    "    2. Generate using a function that draws from a univariate normal distribution, a sample of size $ \\textit{m} $ of $ a^{X}_{i} $,(i= 1,2). This forms a matrix $ \\mathbf{A} $ with realizations $ \\vec{a}^{X^{k}} $, $k= 1,\\cdots,m$ as columns, i.g. a $ 2 \\times m $ matrix.\n",
    "    3. The sample of $ \\vec{X} $ is the given as $ \\mathbf{X}$ = $ \\mathbf{E} \\mathbf{A}  + \\vec{\\mu}$\n",
    "    4. Use sample of size $ \\textit{m} $ to plot a 2-dimensional histogram and compare with PDF of 1.\n",
    "    \n",
    " * $ \\mathbf{A} $ = $ \\begin{pmatrix} a_1^{x(1)} & a_1^{x(m)} \\\\ a_2^{x(1)} & a_2^{x(m)} \\end{pmatrix} $\n",
    " * $ \\mathbf{E}$ = $ (\\vec{e}_1, \\cdots, \\vec{e}_q) $ \n",
    " * $ \\Lambda $ = $ (\\lambda_1, \\cdots, \\lambda_q) $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_multivariate_normal_transformed(a_x, a_mu, lbd_mat):\n",
    "    \"\"\"\n",
    "    Gaussian multivariate PDF (transformed on orthonormal basis)\n",
    "    \n",
    "    a_x:        array with factors of elements of random variable (dim: q)\n",
    "    a_mu:       array with factors of mean of random variable (dim: q)\n",
    "    lbd_mat:    matrix with eigenvalues on diagonal (dim: q x q)\n",
    "    \"\"\"\n",
    "    \n",
    "    q = lbd_mat.shape[0]\n",
    "    \n",
    "    Z = np.sqrt((2*np.pi)**q*det(lbd_mat))\n",
    "    \n",
    "    prob = 1/Z * np.exp(-0.5*(a_x - a_mu).T @ inv(lbd_mat) @ (a_x - a_mu))\n",
    "    \n",
    "    #di = np.diag_indices(q)\n",
    "    #np.prod(1/np.sqrt(2*np.pi*lbd_mat[di])*np.exp(-0.5*(a_x - a_mu)**2/lbd_mat[di]))\n",
    "    \n",
    "    return prob\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "L , E = eig(cov)\n",
    "lbd_mat = np.diag(np.array([L[0], L[1]])) #contructing a diagonal matrix with the eigen values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "a_mu = E.T @ mean "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = 20000000 # the number of independent samples we wish to draw"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "a_x1 = np.random.normal(loc=a_mu[0], scale=np.sqrt(L[0]), size=m)\n",
    "a_x2 = np.random.normal(loc=a_mu[1], scale=np.sqrt(L[1]), size=m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([a_x1, a_x2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(2, 20000000)\n"
     ]
    }
   ],
   "source": [
    "X =  E @ A\n",
    "x1 = X[0, :]\n",
    "x2 = X[1, :]\n",
    "print(X.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate probability density using\n",
    "# a: the independent sample on the orthonormal basis\n",
    "# b: the dependent sample in the original data space\n",
    "if False:\n",
    "    prob_a = np.zeros(shape=m)\n",
    "    prob_b = np.zeros(shape=m)\n",
    "\n",
    "    for i in range(m):\n",
    "        prob_a[i] = get_multivariate_normal_transformed(a_x=np.array([a_x1[i], a_x2[i]]), a_mu=a_mu, lbd_mat=lbd_mat)\n",
    "        prob_test[i] = get_multivariate_normal(x=np.array([x1[i], x2[i]]), mean=mean, cov=cov)\n",
    "\n",
    "    assert (np.array(prob_a*1e15, dtype=np.int) == np.array(prob_b*1e15, dtype=np.int)).all  # should be the same (!)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEWCAYAAAByqrw/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dfbitdV3n8fdnrf1w9nkSGRB5UpgkZ9SrkMGHucxSwwQmo5qxsNI0mxMl6cx01WhOZU9zWU1lliNDySg+ZKSRJ6UQmoysUECRQKSOgHLkKAECB845+/E7f9z3zsW6v7+11zp777Me+Ly41sVZ9/1b9/1b9177t+79/f4eFBGYmdnkaA27AmZmtrHcsJuZTRg37GZmE8YNu5nZhHHDbmY2Ydywm5lNGDfsY0DS8yXdNux6ZCS9QNLeI3SuV0n6xGG+tmc9JV0k6eeyspJukfSCwzlvH/V6s6T3bsBxniTpYUntNcr9oKSPrfd8NtrcsG8iSXdKOlj/wq0+fq+P14Wkp6w+j4i/iYinblId3yXpVzbj2OMkIi6IiF8u7Ht6RHwcNq4h3mgR8aWI2B4Ry2uUe19EfMfq8+7Pmk2GqWFX4DHgpRFx9bArMQ4kTUXE0rDrYTbufMc+JJKeIumvJT0o6V5Jf1Rvv6Yu8tn6Dv/7k9DAnZJ+WtJNkh6R9E5Jx0n6c0n7JV0t6fEd5f9Y0lfqc10j6en19l3ADwI/U5/rz+rtJ0j6kKR/lnSHpNd1HGuuvsv/mqTPAc9a432GpNdJur1+n78hqVXve5Wkv5X025LuB94s6XGSLq3P/UVJ/2O1/NcPqd+t38vnJX17x45XS7q1vga3S/qxpD4/W9fjTkk/2LG9+JdLXfYsSWcDPwt8f329PivpZZJu6Cr/U5L+tHCsU+uf+35JVwHHdO1/rqS/k/RAffwXdOz7uKRfrq/Zfkkfk3RMve+U+lpPdVzb2+tyd6y+V3WEswqftZslvbTjnNP19To9ez82oiLCj016AHcCZxX2/SHwJqov1y3At3TsC+ApHc9fAOztOu61wHHAicA9wKeBZwKzwP8DfqGj/I8AO+p9bwVu7Nj3LuBXOp63gBuAnwdmgH8N3A68pN7/FuBvgKOBk4GbO+uWvM8A/qou/yTgH4Efrfe9ClgCfpLqr8c54FLgw3V9T6nLv6ar/H8FpoHvBx4Ejq73/wfgGwAB3wYcAM7ouIZLwG/V1+HbgEeAp3Zfh8L1Pqv+95uB93bsmwXuB/5tx7bPAP+xcD3+vqMO3wrsXz1e/bO8Dzi3/jm8uH5+bL3/48AXgG+sr9XHgbfU+06pr/UUsA14qOO9HQ88veMafqLHZ+1ngD/qeH4e8A/D/l3yY7CH79g335/Wd1+rj/9cb18EngycEBGHImLQpODvRsRXI+LLVA3tJyPiMxExD1xO1cgDEBGXRMT+et+bgW+W9LjCcZ9F1ZD8UkQsRMTtwO8D59f7vw/41Yi4PyLuAt7WR11/rS7/Jaovlpd37Ls7In43qhDMAlVj/ca6vncCvwm8oqP8PcBbI2IxIv4IuI2qQSciPhoRX4jKXwMfA57fVZefi4j5ev9H6/dz2Opr+kfADwHUfw2dAnyku6ykJ1Fd39U6XAP8WUeRHwKuiIgrImIlIq4Crqdq6Ff934j4x4g4CFwGlO6kV4BnSJqLiH0RcUufb+m9wLmSdtbPXwG8p8/X2ohww775vjsijup4/H69/Weo7iw/parXxY8MeNyvdvz7YPJ8O4CktqS3SPqCpIeo7j6hKwTQ4cnACZ1fRlThh+Pq/ScAd3WU/2Ifde0uf0Jh3zFUfyV8sav8iR3PvxwR0bX/BABJ50i6VtL9db3P5dHv82sR8UiPuhyudwM/IElUDeFldYPf7YRCHVY9GXhZ17X/Fqo77lVf6fj3Aeqfc6f6+N8PXADsk/RRSf+mnzcSEXcDfwv8R0lHAecA7+vntTY63LAPSUR8JSL+c0ScAPwY8L+1Ob0TfoDqz+mzgMdR3U1C9aUC1Z/ine4C7uj6MtoREat3jfuoQjCrntRHHbrL393xvPP89/L1v2Q6y3+54/mJdQP6qONJmgU+BPwv4LiIOAq4gq+/T4DHS9rWoy79aEyHGhHXUv218Xyq6126w91XqMOqu4D3dF37bRHxlgHrSERcGREvpvpS+DzVX139ejfVXw8vA/6+/qvQxogb9iGpk24n1U+/RtVgrHZV+ypVbHsj7ADmqWK1W4H/2bW/+1yfAh6S9N/rRGlb0jMkrSZJLwPeKOnxdf1/so86/HRd/mTg9VShi4aouupdBvyqpB2Sngz8N6rwwKonAK+rk3ovA/4tVQM+QxW3/mdgSdI5wHfQ9IuSZiQ9H/hO4I/7qH+nrwKndCV0ocoN/B6wVAqrRcQXqUIrq3X4FuClHUXeC7xU0kvq675FVeL8pOx4JaoS6d9Vf4HMAw/z9c9W9n66P2t/CpxB9bO6dJBz22hww775/kyP7sd+eb39WcAnJT0M7AZeHxF31PveDLy7/nN8XTFgql/ML1Ld9X6OKuna6Z3A0+pz/WnduL6UKnZ7B9Vd9B9Q3e0D/GJ9vDuoYtj9xF8/TJWQvZEqrv3OHmV/kiqpeTvwCeD9wCUd+z8JnFbX61eB/xQR90XEfuB1VF8MX6O6c97ddeyv1PvupgovXBARn++j/p1Wvwjuk/Tpju3vAZ7B2tfjB4DnUCVcf4GOhrPOWZxHFfr6Z6o7+J9m8N/TFvBTVO/zfqpE8U8Uyr6Zrs9aHb//EHAq8CcDnttGgB4drjTbWJICOC0i9gy7LptJ0hxVYveMiPinYddnvST9PPCNEfFDw66LDc4DlMw2xo8D101Io3408Boe3RvJxsjQGnZJW4BrqOKiU8AHI+IXhlUfs8Ml6U6qJO13D7kq61Z3x30rVRL3mrXK22gaWiim7tmwLSIeljRNFU99fd3DwMzMDtPQ7tjrvsgP10+n64cD/mZm6zTUGLuqKUZvAJ4CvD0iPpmU2QXsAmjT/ndb2dldZJATHv5rB+WktNlQ7edr90bEses5xkteuC3uu7/nhJkA3HDT/JURcfZ6zrWRhtqw113rTq9HuF0u6RkRcXNXmYuBiwF26uh47lTWNblPja7HmyhWBii6AV8CA5zP7LHg6vhgP6Oie7rv/mU+deXaY/Dax/9TaST3UIxEP/aIeIBqQqOR+cYzMwtgpY//Rs3QGnZJx9Z36qt9gM+iGvpsZjYSgmAxltd8jJphhmKOpxrx1qb6grksIhoz4pmZDdMo3pGvZZi9Ym6iY2rZwzZA3FytI5g8JV96MounF1epHChO3+d1cCzerG9BsDyGHSE88tTMrIeVMeyF7YbdzKygmnLVDbuZ2UTxHbuZ2QQJYNEx9tFQTJIOMkBpkxKt1Sy2fYpmHUqDmdTqLylaTLI6qWrWEMSGhWIknQ38DlXPij/oXhmrnj/rd6iWdDwAvCoiPl3vu5Nq4fNlqsVczux1rols2M3MNkTA8ga063W37rcDLwb2AtdJ2h0Rn+sodg7VIjKnUS3G8o76/6teGBH39nO+kRh5amY2iqqRp2s/+vBsYE9E3B4RC8AHqFbL6nQecGlUrgWOknR894H64YbdzKxILPfxAI6RdH3HY1fXgU6kWupw1d56W79lAviYpBuSYzeMVyhGasTJ03h6KZaels1j6Wolx9iI2SGTREysFL7zkxh7cTDTcnNYczoYqhCL73uAU1W4/7JmY6xKnvb1e3/vGnHv7CDdv6C9yjwvIu6W9ATgKkmf77UQiu/YzcwKqn7sfd2xr2UvcHLH85OoFhvvq0xErP7/HuByqtBOkRt2M7MeVkJrPvpwHXCapFMlzQDnA7u7yuwGXqnKc4EHI2KfpG2SdgBI2gZ8B3AzPYxXKMbM7AhavWNf93EiliRdCFxJ1d3xkoi4RdIF9f6LgCuoujruoeru+Or65cdRrVcBVZv9/oj4i17nc8NuZlYQiOUNCmxExBVUjXfntos6/h3Aa5PX3Q588yDnmsyGvThAKUtGFrKRgyRPk/OpkMCNJPGo0gpK2Yi3QqI1sveWvT5JstaFk4M60WrWZ6hlpExmw25mtgECsRClrmijyw27mVlBNUBp/PqYuGE3M+thI5KnR9pjqmFPBx1l2wDaze0DxeMLx00/IsUBSlmMvFQ22b641CxWGpCVHbcUW1xvPN5xdxsTEWI5fMduZjZRVnzHbmY2Oark6fg1k+NXYzOzI8TJ03GQxZdLMecsnj5VuFxZPD2J0fc6Xyrtx15YaCOL0081Y+wsJdsgj8cX+rxn8fgodY9PJh3zYh82Tpbdj93MbHJs5MjTI8kNu5lZDyvuFWNmNjmqScDcsJuZTYxALHpKgTE0wIRhxcFMU8kPfnq6/+OWEq2Z4gClJKm6uNg8fZIkBfJE60Lz9UCegC1dxzQBm2daPbmYjZoIxnKA0tBqLOlkSX8l6VZJt0h6/bDqYmaWEyt9PEbNMO/Yl4CfiohP16uD3CDpqoj43BDrZGb2L4LxvGMfWsMeEfuAffW/90u6lWpFbjfsZjYynDw9TJJOAZ4JfLJnwYgktjpAYiOLQ5dkA5RKsfAknh7ThUs71TxGJNuAwQYzLTffm5aaddBSITY9n8Tjpwsx9oWFZFteNrI4f2ERkjQe78nFbIiCvtc0HSlDb9glbQc+BPyXiHgo2b8L2AWwha1HuHZm9lgWwKLnihmMpGmqRv19EfEnWZmIuBi4GGCnjh7gltvMbL3k+dgHoWrJ7XcCt0bEbw2rHmZmJcF4jjwdZo2fB7wCeJGkG+vHuUOsj5lZw3J9197rMWqG2SvmExQWFFq3wgyIg+RZU6VkZpZULSREV2azRGtetphU7ZOWmudKV0oCWoeSsgeTJCmkyWJNzedl55P3UBoktdS8vlEo61kj7UiI0FjesY9fVsDM7AipkqeeUsDMbIJ4zVMzs4lSJU9HL4a+lrFr2KM7fq5kYM5G/OWUrUhUksTeS/HxmGlWbnk2r/DKTDaYaYAPWZJraC3l+YfWluZHoT2XT2TWOpgMOjqQv4d0JapsgBPAfGF7Zr2Tiznubn3yyFMzswnikadmZhPIi1mbmU2QCFgcZJ2AEeGG3cysoArFuGE/8rIkWCEmFknCTcuFTOtUkmQszQ6ZbS+sthTtZt2yJCnA8pbm9uWZ/L2t9PmTbBXGBrWWmtehfShPME4dSBKtM3kFWo8kM0wezMumidbCoLBIErDFSGiSYI88z+qkqjWM4sjStYzfV5GZ2RGy2t1xrUc/JJ0t6TZJeyS9IdkvSW+r998k6Yyu/W1Jn5H0kbXO5YbdzKyoCsWs9VjzKFIbeDtwDvA04OWSntZV7BzgtPqxC3hH1/7XA7f2U2s37GZmPWzQmqfPBvZExO0RsQB8ADivq8x5wKVRuRY4StLxAJJOAv4D8Af9nGz8Y+yJxiCmmpJYa3EgUjZZVmlysbQSednsr7bSF/7SlmbhbBvA8mwySGqAr+0s9t5eyA8wlQxGmp7LcxXTc0k8/qFSPL55DBVi7Fk8PuYLE5GVJh1LFGPvjYKOxT8WVL1i+hrxeIyk6zueX1yvJbHqROCujud7ged0HSMrcyLVEqJvBX4G2NFPZSayYTcz2wgDDFC6NyLO7LE/O0j33V9aRtJ3AvdExA2SXtBPZRyKMTPrYYNCMXuBkzuenwTc3WeZ5wHfJelOqhDOiyS9t9fJ3LCbmRVsYK+Y64DTJJ0qaQY4H9jdVWY38Mq6d8xzgQcjYl9EvDEiToqIU+rX/b+I+KFeJ3Moxsysh40YoBQRS5IuBK6kWvLnkoi4RdIF9f6LgCuAc4E9wAHg1Yd7vvFr2LuSVtnMfdnqOlXhJMGYzhJIOnugSmWXku3LhQRukoDNBi1BnvwsDVBamsvKJscs5IGUXLLWYn6u9vbm9ultednZh5pvYqYwIGsqmeWyNV0YzPTIgea2VuHOSUlSdYCJJLNEqVdremyIEEsbNPI0Iq6garw7t13U8e8AXrvGMT4OfHytc41fw25mdgR5dkczswnihTbMzCaQG/YRlw1cUjYQCfIVeoox9uYAGCXbqvM1VyVSIR6ficJPLIunL29tblvaUjhw8tnN4u4ArSQ+vbS/FPtvxs2X5vKY5WwSe5+ZypMC7eQQxcFM2fbCJG3Z6k6RDnAqfBbk1ZomiRfaMDObQH32Ux8pbtjNzAoiYMkLbZiZTRaHYoZhgD7Gaf/2ARblyGLpAEr6sWshj8G2ku2txVKf96xiadG0f/rybLJtax7vjdnkOhZuVLTQvGaLO/PC00nsfWmuEI+fbX4cS4uQzEw3t09nC3WQTxhWisevV/q5yeLu4Nj7GHCM3cxsAoUbdjOzyTKOydOhZgUkXSLpHkk3D7MeZmaZiI1bGu9IGna6913A2UOug5lZgVheaa35GDVDDcVExDWSTjli58tWQMpWVQKUrYBUWolnajEp2xyIBKCF5jGyhCpAeyFJ+pXGxWR54XYy4ViSJAVobW++h6nZ/GTZSlRLC/lH6dDO5valbXnZpWQiseUteUJ0LpkMLabyO6fpZJK1VmkwU2ng0noUkqRpkt8J1ZHjGPsmkLSLamFXtpAMpTQz2ySeK2aT1OsGXgywU0cPsOiomdk6RXH54pE28g27mdkwjWOvmMls2EtxykEmaEoGmkRpUMtCMz6tqXwlBy00Y+/tQ3nsvn2o+eOZms9jwItLzbppOYlDF0LI01uadXj8toNp2bnp5ntbLhz44fnmKKkHtyWrggAHtzXLLs/kMfaVZIDSynQyExowlyzAMVOIpWex9/X+WueTiOUD5ryAx2iJOnk6bobd3fEPgb8Hnippr6TXDLM+ZmbdItZ+jJph94p5+TDPb2a2FveKMTObINUduRt2M7OJ4u6Oo26Q1eaTFXJUmqUvm/VxMRm0BHBovrGpNZP/GKYONBOt04/kycSF7ckgnCR/qyTJWpIlSQFO2PZQY9u2dl52fqX53u7fkY9HuGv7UY1tD2zZnpZdmWkmSpenCzN1tprXsfS7OpMkT9uDJFQHmDUyS6qmM5DiwUzDNIox9LU8thp2M7MBBGJlDHvFuGE3M+thDG/Y3bCbmRU5eTpZ0gnDstVxAJIBMNmgJQBNNS+5Dubx6fYjzdjw1FweY595OFnVKIm7Lx7K/6xcWmwet/SBnms339uT5+5Ly862mmUPbE2WdgKOm9vf2LZnyzFp2btmj24ed3ZLWjaSn8/KVD5JW1Y2q23+UyjE3gcI0now0wgaw1t2N+xmZj34jt3MbIIEsLLiht3MbHIE5b6xI8wNu5lZD+7HPo42aybI0gClZHZGtfNUXOtg88cz/XBednmuedzFrc2y2SpFAId2NJOJ+3fmSc7FleZxt2ajoYATpr/W2LZjSz5r5JNn721sO372wbTszckxbpl+Ylr2QCsZEFVcKSn7lWgmZfMrA+2kFdiI+708qVpI5tvGcsNuZjZJNJbJ0/EbUmVmdiRFH48+SDpb0m2S9kh6Q7Jfkt5W779J0hn19i2SPiXps5JukfSLa53Ld+xmZiUBsQG9YiS1gbcDLwb2AtdJ2h0Rn+sodg5wWv14DvCO+v/zwIsi4mFJ08AnJP15RFxbOp8b9pL1ThhWOu5CcozSaksPN8u2kwFOADOzSTx9S7MWS0ksHmBxe/O4D+7IVzq6d9u2xrb9W/PBQVtnm5OePbHdnEQM4OSpBxrbTphqxugBjp56pLFtSzJwCuAzOqmx7YCa76GS5VaSQWWRv9/ZJMbeKmTftDLAQKJswFxBZKF3D1pahw0JxTwb2BMRtwNI+gBwHtDZsJ8HXBoRAVwr6ShJx0fEPuDhusx0/ej5gXAoxsysl/5CMcdIur7jsavrKCcCd3U831tv66uMpLakG4F7gKsi4pO9quw7djOzXvr7Y+neiDizx/7str/7yMUyEbEMnC7pKOBySc+IiJtLJ/Mdu5lZyeoApbUea9sLnNzx/CTg7kHLRMQDwMeBs3udzHfsGyCdMEyFr/msP/J8IcaeLNrQSiYGA5iabsbYZ2ea39tZ3B1gaWtz+8G5PI785W2Pa2z7yvadadkHZpt9yJ9UiJufkHTRP7aV92P/V+1mjD2bcKzkRnX/FVw5EMnCHtG8jq3l/FdHK81rNlOIsWexd5VGw2Tx+PxjM1h+yLH3NW3QAKXrgNMknQp8GTgf+IGuMruBC+v4+3OAByNin6RjgcWIeEDSHHAW8Gu9TuaG3cyslw3oFRMRS5IuBK6kmhz0koi4RdIF9f6LgCuAc4E9wAHg1fXLjwfeXfesaQGXRcRHep3PDbuZWQ+lP74HFRFXUDXendsu6vh3AK9NXncT8MxBzuWG3cysZIABSKPEDbuZWVHfydGR4oZ9EBsxYVjyISlOGJYkTzV1KC3anmpmHmdmmtuWZ/Mk2nIywCnbBvDAXDPBeOvccWnZJ8w0V0U6MZkYDOAkmgnRx7fyBO6OVnMUzrT+KS2bWUkSogCfWWkmVQ8uNwcztZbyX3YlSdUsoQows9T8jKiwSleaVB0kq7eUf8bSpKoTqo/mO3Yzswkzht9zbtjNzErGdKGNngOUJO2U9A3J9m/aiJOvNduZmdmwKdZ+jJriHbuk7wPeCtxTzyj2qoi4rt79LuCM9Zy4z9nOxsMmTRiWfV5Ki3Jk8fj2VLMOs1P52Vamm0tHrCQDnABWZpuDpL40+6/SsjdMN0fRPK6dL7Sxo/WPjW1PncoWmIDtSez95MJxnzV3R2Pb/FH5QK+HF2ca225daF7zgwv5BGmtpWQw01J+Li0238N0YWIwLSfx+MLEYJEcQ1H43CSfRw9m6jKCDfdaet2x/yzw7yLidKqO8u+R9L31vo342+RfZjuLiAVgdbYzMzNbh14x9nY9XSQR8SlJLwQ+IukkNuY7LJvJ7DndhepZ0nYBbCFZ4szMbBONYqhlLb3u2Pd3xtfrRv4FVHfVT9+Ac/cz2xkRcXFEnBkRZ04XV5o0M9sEQTWlwFqPEdOrYf9xoCXpaasbImI/1axiP7oB5+5ntjMzs+HaoKXxjqRiKCYiPgsg6WZJ7wF+nWq59l8HzgTes85z9zPb2fgqJJryxFRhUEoy8CnmmysSAaiVDGba33z9VCH5uqXdLLsy3UwkAixPN891YDofhHPb9BMa27ZO5YNltraS9zb3xbTsU5OBWttb+V90p0w1k6r3b9mbln3w6GZS9JEkoXrnwrHp6w8uNMu2FvP7p1Zy3NZi/rmZSgYzpTM+kg9yiuKskdn2/PM40EC8CTJpoZhVz6G6s/47qsb4buB56z1xRCwBq7Od3Uo1Y9kt6z2umdmGmqQ79g6LwEFgjuqO/Y6IjfmazmY7MzMbKSPYcK+lnzv266ga9mcB3wK8XNIHN7VWZmYjoJ/BSaMYqunnjv01EXF9/e+vAOdJesUm1mmyDTKYqTAhVErN+HS6AlMWJwWmk81zhcFM0UoG3CTnAnik3eyi+lmdkJbNHHp8PriHrXsam7JYOsByEl8+eSpfmekrs//c2Hb3445qbHtoPs8p3H+wWbY9n+c12oeaF709n+c1sgFK7aVCbib73GzAMkCRHfexEHcfwV4va1mzYe9o1Du3rTdxamY2FkbxjnwtngTMzKwXN+xmZhNkRGPoa3HDbmbWixt2OyylwUwD5E5pTqKYyhKqAK1kgNJMK0+0RjIYaqWdf5QiOe4BmisSAdxI/0nVxWS2wge2fCkt+8Sp5ipOK4U5to+deqix7clz9zW27duxM339Q0c3k6oLB/I5jqYONOswNZ9fx9Zic/CV5vMPSCtJqqowmClN0BdmjcwT/2PY6g1IY5gL7qe7o5mZjRHfsZuZ9TKGf5S4YTczK3Hy1I6EdJBIyVISRz5UmEQs2dYqxONnku2hfMBO+hGLPAKYxd4/EyfmZZNVie7Zmce9T9/WjL3vaOWDmbLY/WyrOWnZ0bMH0tfv3N6cnOy+nfnkZAuPNK/N1MH8mk8dbNarfagwmGmpuepUtg0oxNgLk4ulqzA9BlZgcsNuZjZh3LCbmU0OMZ69Ytywm5mVOMZuZjaB3LDbhsqSTYXZGQdKqvZJhQFKWQqtvBptM6mqKHzskqTqweV8MNPnFo9vbNu/kNfioaVmHU6YzWd33NFuJj8PLDePO9vKk5E7Z5uvf2BbPkBpcWfzSi48XEqeJqthHchnvmwfbCZVNZ/XN50JMplJEgqrMCUDlNQaYPWwcUioumE3M5ssDsWYmU2aMWzYPaWAmVlJVL1i1nr0Q9LZkm6TtEfSG5L9kvS2ev9Nks6ot58s6a8k3SrpFkmvX+tcvmMfN6WYZJ8r2cRiYaBKesw83pttbRfK5lHvwmCmJPauwuo1hxbmGtvuXDg2LZutdnT89uZkXwDHb823d3tkKR8clJmeza/5/Ezz57O8NV9taXF78zosbi+tzNSs2/RiYbWlbMKwgQYzNW9nY2kM4uaD2IA7dklt4O3Ai4G9wHWSdkfE5zqKnQOcVj+eA7yj/v8S8FMR8WlJO4AbJF3V9dpH8R27mVkPG7Tm6bOBPRFxe0QsAB8Azusqcx5waVSuBY6SdHxE7IuITwNExH7gViAfkl1zw25m1kv08YBjJF3f8djVdZQTgbs6nu+l2TivWUbSKcAzgU/2qrJDMWZmJV9vuNdyb0Sc2WN/FlPsPnLPMpK2Ax8C/ktE9IwZumGfFOkiCM0/yIp9jAeJvSf9mUvruLeTCaW2ZP2hgdZyMxbeWi701V5snvHgQh73vv/gUY1tDx7VjNFDvoBG1je95MBisw5RWNQjkhj70pY8bp50pWdxa6HP+1wyYdjB/Fe9PZ9c34X8mpN9RpKfbz5ZGIzjhGFiw7o77gVO7nh+EnB3v2UkTVM16u+LiD9Z62QOxZiZ9bBBMfbrgNMknSppBjgf2N1VZjfwyrp3zHOBByNin6plz94J3BoRv9XPyXzHbmbWywbcsUfEkqQLgSupBm9fEhG3SLqg3n8RcAVwLrAHOAC8un7584BXAP8g6cZ6289GxBWl87lhNzPrZYMGKNUN8RVd2y7q+HcAr01e9wnK0c7UUEIxkl5Wd7RfkdQr4WBmNjx9hGFGccqBYd2x3zpLEDEAAA4wSURBVAx8L/B/hnT+x4Y+E6qQJ1UHSqgWpCszFcrOZEnZZAAMQCtJ8LUW8yNPHWom8+YfyZOn9z2uedwHtzfLTk3nA34iGVC1tJAnE7WcXJ1W/n5Xpppll2fym7jlueZ1WN6S/6q3DjXfr+b7X5kpG7SUThYG6WCmLKFanWyEJgwbwYZ7LUNp2CPiVgAVRiuamY0KL7SxCeqO/rsAtpBPf2pmtllGMdSylk1r2CVdDTwx2fWmiPhwv8eJiIuBiwF26ugxvMRmNrb6H6A0UjatYY+Iszbr2LYOhTjlpg1mSigZ1ALQShZ4mClMKKXF5mCmdjI4CKA933xv7UN5GHDhkWbMeXFH89dkfkvh2kwnA3ayWDqgpF6lslmyYqUwD1kWe1+ZzeP8KzPN99aeLjQL08nApT7j7gBqJW+iNHirkFsZihGqSr9GPhRjZjYsGzjy9IgaVnfH75G0F/j3wEclXTmMepiZrUUrseZj1AyrV8zlwOXDOLeZWd8cYzczmzzjGIpxw26VIzmYqZA8VTZAqZDsnU4SdK3FQlI2mfVx6lD+0Z860EzmLTySDA5KZlAEWE6SlFEIeGb9o1uLeVmyt1ZocKKdJE9n8krETFLfZBuAkoFWmkquYzv/LESSHC+1mulnbFgzQbphNzObLL5jNzObNG7YzcwmSHhKAZs0AwxmKk7mlBmge1gWdwdQEq+dWi6UTWLv7cJqS1PJSkPTB5rvd3FbPrBmqTluipXCgkRZ7L1VSFVksXeVLnlyzQrjgNJ4fGRxcwrx9DTGXlhBqd2scCmHEgN8nDbTuPZjd8NuZtZLabbKEeaG3cysB9+xm5lNEg9QMjObPE6e2mNDlvDKVrwBojDTX2ohO0Dhdik5bmnWyPZSUnYxr1crWT1o6kAy4+P2PEGYJVVLKx1FcohSkjNLqrYXCqtLJWVL4YTIZlxMEqrVgZOfcbJYjorJ02QwVDZoCdIKF2cbLQ5cyjcPyg27mdkkCZw8NTObNE6emplNGjfs9pg1wERMAw0+KR03i6cXBj4pmzAsibsDaKEZoG4dmm1sax/KRx1NPdKMIy/P5THgLPaeDRiqdjTfW2kwUxZ7zwZpAYPNJd5O3ke/2yCP0WcxfkBJsmEYg5Y8QMnMbNLEaC6ksRY37GZmvYxfu+6G3cysF4dizPoxyMRPhY7d2dYodUvLFvAoLfax1Axct5NFRLL+7gDtuWbsfXlL/mu2kizKsTJV6kOebC/lFJK31l4ovd9k+yDd+5J+7Ok2SGPvpT7vkb23ARblAPLFSQYVDDRp3ahww25m1sv4tetu2M3MenEoxsxswrhXjJnZJPHsjmbrlCRVixM8ZSs2le6ssu2lycmyycUWm8sXabE5aAmgvdDc3ioMZlqZaf76xXRhMrWp5vbShGHZHaYGWF0qTahC/0nVbCDSmKoGKG1Myy7pbOB3gDbwBxHxlq79qvefCxwAXhURn673XQJ8J3BPRDxjrXNNzk/AzGwzrPTxWIOkNvB24BzgacDLJT2tq9g5wGn1Yxfwjo597wLO7rfKbtjNzHpQxJqPPjwb2BMRt0fEAvAB4LyuMucBl0blWuAoSccDRMQ1wP391tkNu5lZSfT5gGMkXd/x2NV1pBOBuzqe7623DVqmL0OJsUv6DeClVEsrfAF4dUQ8MIy62IgrDWZKYu/FhRgGmTwqG7iULAZRGuBEMphJpcFM00mMfaYwYGcq+VUt3ZZlA4QKd5Vp7L0UY88WxciOW7o2aQX6H5CVTQwGmz05WN9zxdwbEWf22J+OqTuMMn0Z1h37VcAzIuKbgH8E3jikepiZ9Rax9mNte4GTO56fBNx9GGX6MpSGPSI+FhGrtzbXUr0BM7PREtUUDWs9+nAdcJqkUyXNAOcDu7vK7AZeqcpzgQcjYt/hVHsUYuw/Avz5sCthZpbagDv2+kb2QuBK4Fbgsoi4RdIFki6oi10B3A7sAX4f+InV10v6Q+DvgadK2ivpNb3Ot2kxdklXA09Mdr0pIj5cl3kTsAS8r8dxdlF1/WELWzehpmZmPWzQAKWIuIKq8e7cdlHHvwN4beG1Lx/kXJvWsEfEWb32S/phqg733x7FafkgIi4GLgbYqaPHcAyYbYr1DmYqJGUVzeRl+vEcZIBTMmMkANPNgUtayJOn6SyIhZkRiysYZfpNiAJkq05lr58wxUT5CBtWr5izgf8OfFtEHBhGHczM1hRszPS/R9iwphT4PWAWuKoaRcu1EXFB75eYmR1Zou8BSCNlKA17RDxlGOc1MxuYG3azISrEzfOi65xcrDjhWLZ8USHGnsTelQ1EgjyeXoqxZ6stlQYCZUoNWRZPz95vKSa9WdPfapM797lhNzObII6xm5lNHveKMTObKH1PGTBS3LCbmZUEbtjNxsY6Z42MwgyIfQ9wgnyQUynRmswEqVLSMBugVJxFcZDBTNlAr+Z7i9Lgrew6jEOjOX6RGDfsZma9uB+7mdmkccNuZjZBIsZyPhw37GadBppcLNPnACdAyUCiweLxpQFKg8TYBxi4lMneW+k9JN0GYxy6EvqO3cxswrhhNzObIMHmTYWwidywm5kVxUBzEI0KN+xmZiWBk6dmE2nds0aWBuwkictBEq2lBicpm67ABIPN+tivQkw6TZQWE60jFP5wjN3MbMK4YTczmySeBMzMbLIE5YVDRpgbdrONtM54fDbhGEBkYXoV4vFJ7D5KMev1DlBar1K90oFi/ZfdUL5jNzObJJ5SwMxssgSE+7GbmU2YUep62Sc37GbDss4JxzYiHp91sc/6zA+ktADIQPmHEWpMHWM3M5sgEe4VY2Y2cXzHbmY2SaK8husIc8NuZlbiaXvNbN3WPeFYrpRoTY9bukEtJUUbNuAOd5QSrWPY3XGQNb82jKRflnSTpBslfUzSCcOoh5lZL0H1xbHWox+SzpZ0m6Q9kt6Q7Jekt9X7b5J0Rr+v7TaUhh34jYj4pog4HfgI8PNDqoeZWVnUC22s9ViDpDbwduAc4GnAyyU9ravYOcBp9WMX8I4BXvsoQwnFRMRDHU+3UX0xmpmNnA1Knj4b2BMRtwNI+gBwHvC5jjLnAZdGtaL5tZKOknQ8cEofr32UocXYJf0q8ErgQeCFPcrtovr2Api/Oj548xGo3kY6Brh32JUYwLjVFx6rdR7kdmj9YeJxvMZPXe8B9vO1K6+ODx7TR9Etkq7veH5xRFzc8fxE4K6O53uB53QdIytzYp+vfZRNa9glXQ08Mdn1poj4cES8CXiTpDcCFwK/kB2nvjgX18e8PiLO3Kw6b4Zxq/O41Rdc5yNh3OoLVZ3Xe4yIOHsj6gJkw3m7v5pLZfp57aNsWsMeEWf1WfT9wEcpNOxmZhNgL3Byx/OTgLv7LDPTx2sfZVi9Yk7rePpdwOeHUQ8zsyPkOuA0SadKmgHOB3Z3ldkNvLLuHfNc4MGI2Nfnax9lWDH2t0h6KlXk74vABX2+7uK1i4yccavzuNUXXOcjYdzqCyNU54hYknQhcCXQBi6JiFskXVDvvwi4AjgX2AMcAF7d67W9zqcYw3kQzMysbFj92M3MbJO4YTczmzAj37BLepmkWyStSDqzY/spkg7W0xLcKOmiYdazU6nO9b431sOCb5P0kmHVsRdJb5b05Y5re+6w65QZdJj1KJB0p6R/qK/rurvjbQZJl0i6R9LNHduOlnSVpH+q///4YdaxW6HOY/E53gwj37ADNwPfC1yT7PtCRJxeP/pNwB4JaZ3rYcDnA08Hzgb+dz1ceBT9dse1vWLYlel2OMOsR8gL6+s6qv3C30X1+ez0BuAvI+I04C/r56PkXTTrDCP+Od4sI9+wR8StEXHbsOsxiB51Pg/4QETMR8QdVNnvZx/Z2k2MfxmiHRELwOowa1uniLgGuL9r83nAu+t/vxv47iNaqTUU6vyYNfIN+xpOlfQZSX8t6fnDrkwfSkOGR9GF9Qxzl4zan921cbqWnQL4mKQb6ukyxsVxdZ9q6v8/Ycj16deof443xUg07JKulnRz8uh1B7YPeFJEPBP4b8D7Je08MjU+7DoPPDR4s6xR/3cA3wCcTnWdf3MYdVzDyFzLAT0vIs6gCiG9VtK3DrtCE2wcPsebYiQW2hhg+oHO18wD8/W/b5D0BeAbgSOSkDqcOtPfsOIjot/6S/p9qqmVR83IXMtBRMTd9f/vkXQ5VUgpyx+Nmq9KOj4i9tUzDt4z7AqtJSK+uvrvEf4cb4qRuGM/HJKOXU08SvrXVHMY3z7cWq1pN3C+pFlJp1LV+VNDrlND/Yu76nuoksGjZuBh1sMmaZukHav/Br6D0by2md3AD9f//mHgw0OsS1/G5HO8KUbijr0XSd8D/C5wLPBRSTdGxEuAbwV+SdIS1VpcF0TESCRPSnWuhxBfRjWP8hLw2ojiQmTD9OuSTqcKbdwJ/Nhwq9N0OMOsR8BxwOWSoPrde39E/MVwq9Qk6Q+BFwDHSNpLNUHfW4DLJL0G+BLwsuHVsKlQ5xeM+ud4s3hKATOzCTO2oRgzM8u5YTczmzBu2M3MJowbdjOzCeOG3cxswrhht7En6S8kPSDpMTMAxawXN+w2CX4DeMWwK2E2Ktyw29iQ9Kx6Qqct9SjOWyQ9IyL+Etg/7PqZjYqRH3lqtioirpO0G/gVYA54b0Q8ZoaJm/XLDbuNm1+imifmEPC6IdfFbCQ5FGPj5mhgO7AD2DLkupiNJDfsNm4uBn4OeB/wa0Oui9lIcijGxoakVwJLEfH+esrmv5P0IuAXgX8DbK9n9ntNRFw5zLqaDZNndzQzmzAOxZiZTRg37GZmE8YNu5nZhHHDbmY2Ydywm5lNGDfsZmYTxg27mdmE+f8DuPXwKYi7JAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot in 2d     \n",
    "X_grid = np.array([np.linspace(-3*sigma_1, 3*sigma_1, 50),  # x1\n",
    "                   np.linspace(-3*sigma_2, 3*sigma_2, 50),  # x2\n",
    "                  ])\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "hs = ax.hist2d(X[0, :], X[1, :], bins=X_grid, cmap='viridis', density=True)\n",
    "ax.set_xlabel('x1')\n",
    "ax.set_ylabel('x2')\n",
    "ax.set_title('Estimated probability denisity')\n",
    "fig.colorbar(hs[3])\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "py3-general [conda env:py3-general]",
   "language": "python",
   "name": "conda-env-py3-general-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}