{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <span style=\"color:red\">Losowanie liczb pseudolosowych z wielowymiarowych rozkładów normalnych w Python na podstawie wykładu 8</span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Importowanie niebędnych pakietów"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from numpy.random import randn\n",
    "from numpy.random import randint"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Przykład obliczania wartości i wektorów własnych (Przykład 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Definiowanie macierzy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[3 2]\n",
      " [4 1]]\n"
     ]
    }
   ],
   "source": [
    "A = np.array([[3, 2], [4, 1]])\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " #### obliczanie wartości i wektorów własnych"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(array([ 5., -1.]), array([[ 0.70710678, -0.4472136 ],\n",
      "       [ 0.70710678,  0.89442719]]))\n"
     ]
    }
   ],
   "source": [
    "results = np.linalg.eig(A)\n",
    "print(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Tylko wartości własne"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 5. -1.]\n"
     ]
    }
   ],
   "source": [
    "print(results[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Tylko wektory własne"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.70710678 -0.4472136 ]\n",
      " [ 0.70710678  0.89442719]]\n"
     ]
    }
   ],
   "source": [
    "print(results[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### .... lub"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 5. -1.]\n",
      "[[ 0.70710678 -0.4472136 ]\n",
      " [ 0.70710678  0.89442719]]\n"
     ]
    }
   ],
   "source": [
    "eigvals, eigvecs = np.linalg.eig(A)\n",
    "print(eigvals)\n",
    "print(eigvecs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Jeśli jesteśmy pewni, że wartości własne są liczbami rzeczywistymi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 5. -1.]\n"
     ]
    }
   ],
   "source": [
    "eigvals = eigvals.real\n",
    "print(eigvals)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Sprawdzanie ortogonalności"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.70710678 0.70710678]\n",
      "[-0.4472136   0.89442719]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.3162277660168379"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v1 = eigvecs[:,0]\n",
    "print(v1)\n",
    "v2 = eigvecs[:,1]\n",
    "print(v2)\n",
    "\n",
    "v1 @ v2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mnożenie macierzy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Tworzenie macierzy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[5 0 6 5]\n",
      " [0 0 7 9]\n",
      " [0 3 6 1]\n",
      " [5 7 8 4]]\n"
     ]
    }
   ],
   "source": [
    "n = 4\n",
    "P = randint(0,10,(n,n))\n",
    "print(P)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Transpozycja macierzy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[5 0 0 5]\n",
      " [0 0 3 7]\n",
      " [6 7 6 8]\n",
      " [5 9 1 4]]\n"
     ]
    }
   ],
   "source": [
    "print(P.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Mnożenie macierzy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 86  87  41  93]\n",
      " [ 87 130  51  92]\n",
      " [ 41  51  46  73]\n",
      " [ 93  92  73 154]]\n"
     ]
    }
   ],
   "source": [
    "S = P @ P.T\n",
    "print(S)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Wyznacznik macierzy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1612900.0000000033\n"
     ]
    }
   ],
   "source": [
    "detS = np.linalg.det(S)\n",
    "print(detS)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Macierz odwrotna"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.07643623 -0.03638787  0.04428297 -0.04541261]\n",
      " [-0.03638787  0.03155806 -0.03031186  0.01749023]\n",
      " [ 0.04428297 -0.03031186  0.11938992 -0.06522785]\n",
      " [-0.04541261  0.01749023 -0.06522785  0.05438899]]\n"
     ]
    }
   ],
   "source": [
    "invS = np.linalg.inv(S)\n",
    "print(invS)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Algorytm generowania liczb z wielowymiarowego rozkładu normalnego"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Algorytm generowania próby $n$ - elementowej z $d$ - wymiarowego rozkładu normalnego:\n",
    "1. Wygeneruj macierz Z rozmiaru $n\\times d$ i wypełnij ją zmiennymi losowymi z rozkładu $N(0,1)$.\n",
    "2. Wyznacz faktoryzację (rozłóż na czynniki) $\\Sigma = QQ^{T}$.\n",
    "3. Użyj transformacji $X=ZQ+J\\mu^{T}$,gdzie $J$ - kolumna złożona z jedynek. \n",
    "\n",
    "Każdy wiersz macierzy $X$ jest zmienną losową z rozkładu $N_{d}(\\mu,\\Sigma)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Faktoryzacja - Metoda spektralnej dekompozycji macierzy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### W metodzie tej mamy $\\Sigma^{\\frac{1}{2}}=P\\Lambda^{\\frac{1}{2}}P^{−1}$,gdzie\n",
    "1. $\\Lambda$ - diagonalna macierz z wartościami własnymi macierzy $\\Sigma$ na przekątnej;\n",
    "2. $P$ - macierz, której kolumnami są wektory własne macierzy $\\Sigma$ odpowiadające wartościom własnymz $\\Lambda$.\n",
    "\n",
    "W metodzie tej $P^{−1}=P^{T}$.Macierz $Q = \\Sigma^{\\frac{1}{2}}$ jest rozkładem macierzy $\\Sigma$ na czynnik takim, że $Q^{T}Q= \\Sigma$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Przykład 2 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Wylosujemy próbę $1000$ elementową z dwuwymiarowego rozkładu normalnego o średniej $\\mu= (0,0)$ i macierzy kowariancji\n",
    "$\\Sigma=\\left[\\begin{array}{cc}1.0 & \\ \\ 0.9 \\\\ 0.9 & \\ \\ 1\\end{array}\\right]$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0. 0.]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x485e088>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dbZBU53Un8P/pngvqwY4b4tmy1dYI7DiwKyMxq4nNLqndwDpCMbY8RisRl6hNVbZKlQ+pDVpnyiNLZZAjryZFyVIq2Q9RrV0blxQHWeAxDnYhuZArCRsUD55BGAuyliWQWqo1CbQcmAZ6Zs5+6LnN7dv3uS99b7/c7v+vymUx03P79iCdfvo85zlHVBVERJRemU7fABERxcNATkSUcgzkREQpx0BORJRyDORERCk30Iknfe9736urV6/uxFMTEaXW8ePH/0lVh9xf70ggX716Naanpzvx1EREqSUiZ72+ztQKEVHKMZATEaUcAzkRUcoxkBMRpRwDORFRynWkaoWIqBOmZorYe/gM3iqVcWM+h/GtazE2Uuj0bcXGQE5EfWFqpogHD5xEubIAACiWynjwwEkA8AzmaQr6TK0QUV/Ye/hMLYjbypUF7D18puGxdtAvlspQXA/6UzPFNt1tNAzkRNQX3iqVQ389StDvBgzkRNQXbsznQn89StDvBgzkRNQXxreuRc7K1n0tZ2UxvnVtw2OjBP1uwEBORH1hbKSAx7avRyGfgwAo5HN4bPt6zw3MKEG/G7BqhYj6xthIIVTlif2YtFStMJATEXkIG/S7AVMrREQpx0BORJRyDORERCkXO5CLyA0i8g8ickJETonII0ncGBERhZPEZudVAFtU9ZKIWAD+TkS+p6rHErg2EREFiB3IVVUBXFr6o7X0P417XSIiCieRHLmIZEVkFsDPAbygqi95POZ+EZkWkenz588n8bRERISEArmqLqjqBgAfAPBREfmIx2OeUtVRVR0dGhpK4mmJiAgJHwhS1ZKI/ADAnQB+nOS1iag3panvd7eKHchFZAhAZSmI5wB8HMAfx74zIuoZpmAdddhDmrXyDSuJFfn7AfyFiGRRTdU8q6p/ncB1iagH+AVrv77f7iCX5pV7q9+wkqhaeRnASOw7IaKe5Besw/b9TvvKPcobVjN4spOIWsovWIft+522iT1urR5UwUBORImZmili0+QRrJk4hE2TRzA1U/QN1mH7fqdtYo9bqwdVMJATUSJMA4s3rxsyBuuwwx7SNrHHrdWDKtiPnIgSYUp/vHj6PB7bvt64URmm7/f41rV1OXKguyf2uLV6UAUDORElwi/9EXdIQ9om9nhp5aAKBnIiSsSN+RyKHsE8qfRHmib2tBtz5ESUiG4aWOy16drLuCInokR0S/oj7TXnzWAgJ6LEdEP6o5WHb6KcLm3nSVQGciJqqXYfrW9VzXmUlX67PxUwR05ELWOqLW9lzrpVNedRTpe2+yQqAzkRtUzUgJbEJmWrNl2jrPTbfRKVgZyIWiZKQEtq9R72tGhUUVb67T6Jyhw5EbVMlNpy0+r9ke+cipxjb8Wma5TTpe0+icpATkSJszc4i6UyBPXT2E0BzbR6vzhXwcW5CoDOlhJGKa9sdymmqLZ/4P3o6KhOT0+3/XmJqPXcFRsAasG84BPQNk0e8Vy9eynkczg6sSWhO04PETmuqqPurzNHTkSJ2nPwVEOKxBnE9x4+47mZ6bVJaZKW9rXtwtQKESVmaqaIUrni+T07LWKqrfZKR1y+Ou95vbS0r20XBnIiSoxfnXRGEHji0r1J6ZWmSVP72nZhICeixPilPBYN23F+P+NcpRdLZWRF6urQO90OoFswR05EiWkm5RH0M2MjhVr+fGGpOKMdJ0TThCtyIkrE1EwRc9fmjd/P5yxcnV8MnSZx9mjJiNSCuC3JKfRpx0BORLF55bKdclYWe+66BUC42mr39dxB3JZU9Uq7G3slLXYgF5GbAHwdwPsALAJ4SlX/JO51iSg9vE5l2ty142ECpN/1nJKoXumF/uVJrMjnAXxOVX8kIu8GcFxEXlDVnyRwbSJKAdPKWADjwR2/VXCYlXZS1Sut7F/eLrEDuaq+DeDtpX/+FxF5BUABAAM5UZdKOpUQdV5n0CrYdL2sCBZVE01/tLtTYSskWrUiIqsBjAB4yeN794vItIhMnz9/PsmnJaIIWtEjPGrr2KD2tqbrPX7vbXhtchuOTmxJbLXc7k6FrZDYZqeIvAvAfgC7VPUX7u+r6lMAngKqvVaSel4iisYURHftm8Xew2ewed0QXjx9PnK3QfvaYX4uaBXczqZT7e5U2AqJBHIRsVAN4s+o6oEkrklEreGXMiiWynj62Lm6P4fd+IvSOjZMKqZd8z+7ZWh0HElUrQiArwJ4RVW/Ev+WiKiVTEHUxLnxl1RuvdtWwd0wNDqO2G1sReTXAfwtgJOolh8CwBdU9bumn2EbW6LOCar5NsnnLFy+No/KwvWYkbOyTU/fSXvtdieY2tiyHzlRl2lHgHMOfojL2Rs87r0zuPtjICdKAVO3vyRmToZ9vqgEwBM7NuCR75yqTfJxWjloYfenbgm8/3a/9jTiYAmiFIg6dT4Mv8n07kHFzcgPWnjwwEnPIA5UR7WFKW9sxWvvF+y1QtRFkjqcYpqZaVehTJ+94FliGGXcGlBdMas29hl3C3NSshcO5nQKV+REXSSJwynOAz9A/eBjoBpUnz52zvNAUJRxa4V8Do9tX493DBOB3IqlsuenAlsvHMzpFAZyoi4S9YSkl7ANp5ycK2ZnqmXFssagnrOyeHLHhtrpyiiB1u8kaZTX7pcu6kcM5ERdZGykgLtvLyAr1Yx1VgR33x5c4+wMbM1WohRLZWyaPIIH9s0CqG5gnvrSnXhyx4ZaYLdX4c77ibKKt3nlvt1vIl7PZb/WpFsMpB2rVoi6SFDlhld5HoDYlScA6nLptnzOwp67wlWcOO/LeczfFGEEwGuT2yLfpymP7yyD7FUsPyRKAb8gtXndEJ45dq4uMOasLJYPZIyT68PyCuLu77n7ioeVVOANqn1v9o0hTUyBnFUrRG0UdODFVKFRLJUbgjhQTVHEXYnnc5bvG4G74gVA6OP6UzNFXL7aOP7Nykrt00TY6wR96ujnTVEGcqKQkji1GDSJxq8Pt2ncWVxX5xcDg7nNmdv2ei3Ossb8oIVLV+ZRWfS476Uvef1OHtg3i137Zus+AQRt4KatW2HSmFohCiHuqcOpmSI+9+wJz2DsHJawed0Q9h8v1j2PX9ojKtMbwspBC1cqi6FX96brRLnXrAjefcOA7xuI/Tt+YN+s8brNpnzSiCc7iWKIc+rQfhMwragXVGvVF/uPF3H37QUUltIEcYK4+6Rmzsoa76E0V8Fj29dj5aAV6tqm60S51wXVwE8B9u/YlDax8+z9EMT9MJATheCXuw6qZY5S112uLODF0+dxdGILCvlcrJW4opr/tl2d988vj40UMPPFO2rlht3irVI5kfr6XsYcOVEIfj28nbXMQOMAhqhHzO3H+/2cCBCUFR20Mrg6v1j7s1eqGqjfeATqe3OvnjgU8q6TTQE52W8yQLqHP7QSAzlRCF6DENxM/UTyg5axoZQXO41gevPI5yxcvjqPSkAkn6ss+n7ftmLZgGdAnJopBm6yBuX3raxgxTL/PPjKQQuDywYa+sIA9avutA9/aCUGcqIQnCtCv5OT7lX01EwRl654l9/t+LWbPDc27ROWnoExI/jFlYpxdd0Mr14pD0+d9Cx3dHv83tvqguvozauMq2bThrGzxa3pwNOmySNciftgICcKyV4Rrpk4ZAxw7k25vYfPeJbfrVg2gEfH1tcCn1eXQnvj89DLb9dW9J6lfDG573lqplg3t9Mkn7MaAqrfqjlMesT982FKNonlh0ShhZ2q4yyHMwV9exiDHdQyCZUFNmPFsiy+/JnrZZQjX3o+VCpo0MrU0jdhh0dE1c/H8b3wZCdRDFEm6RRLZezaN4uHvnXS+Bh7GIN9PVMeOmxuPSPAe3IWSnOVyBuOl68tYPybJwBUV7lhn9OZg784V8H4c9evkRT2KA+H5YdEIew5eMoYxO1OhW6Xry14BtWwwxiiyIpg96duwWuT27Bz43Dkn68sKh488DI2TR5p+h4qC1qrq0+qzSx7lIfDQE4UYGqmaKy6EJhX06bH32DFb3LlVlm8HkQfHVuPnRuHYXh/MSpXFmMPY36rVE60zSzrx8NhICcK4Hd688Z8zrgi96IIny6JyhmER29e1Zqi7gA35nOJzt4M26O83zFHThTALx87vnUtdi0NYugG9qp3/Jsn2h7H7YNFDxh+H83mtVk/HoyBnMjBWZliH4bxazRld+aLm5JwinNC8sEDJ7F8INOSMkUnKyOwsuJZtWL6fTCv3TqJBHIR+RqATwL4uap+JIlrErWbuzLFDt6mLn/bbn0/gHCnPqNQNN+2Non+5EGCug16/T6Y126tpFbk/xvAnwH4ekLXI2q7KM2tFMD+40WM3rwKYyMFTJ+9EOoQTVit6j0e1+shJvC4D/7kBy2oAg/sm8Xew2d4MrMFEtnsVNW/AXAhiWsRdUrUHK5zA+/F0+dbcUtNWTlowcpELFkJIUpHxLGRAo5ObMETOzbgSmURpXKFg5JbqG1VKyJyv4hMi8j0+fPd8y89ka2ZHK4d/JPMkcdhZav15HvvuS3R6wqAzeuGIv+cqYJl177ZWPXlVK9tgVxVn1LVUVUdHRqK/i8EUav5BSrT+jYj0lXBqLKg+MKBl/HId04lel07lRT1tfp9yuHqPDmsIydaYkqPZEVw38bhhoMpQDWXbTdx6hZzlcWW1KqXKwuR3yCCPuU0W19O9RjIiZaYVo8Lqnj62DlkBJ6nJZOsEslZGQiqee5u/I/z4lwl0gra62SmG/umxJfIvysi8g0Afw9grYi8KSL/NYnrErVT0Orx8rWFwKk8YeWsLHZuHEbOavxP8IkdG7D7U7cgm01+wzIJn3v2ROgeKs6TmSasL4+PbWyJlkTpcBiHXYftV7IYZpRbHPmchavz8dvj2lPuw5QTmgZL8Mh9eKY2tt346Y2oI+zVYys5p75/46U3jI9rZRAXAKVyJZE3rCg5bvZNaR0e0SdyaMWReydnPridh37sk6KtGJAcJcfNvimtwRU5da2kelpHNb51LawW5aed+eAoXRObtXPjMF6f3IZXH/sECvlcU0E86D6Z4+48rsipK4WZ1eg1qDdsrtZrwO+eg6dqfcIHrQyWZTO4fO16+iEjiDX02MoI5q7NY/XEoaZ7qUQ1evOq2j83Wx2yqIpCPuf5KUUA9lDpAlyRU1cK6mnd7PCCqZkixr95ou7nHtg3i137ZuuGPcxVFmtBvJDP4ckdG/CVezcEltL5qSxqrb67XWkVZ/662ZWz/Wbnfu0C4L6Nw0yVdAEGcupKQbMaww4vcKdnvnDg5YYWr0Eh1flp4LHt69uSEkmK8/cYpqbbze5a6LVR+cSODXh0rLWbwxQOUyvUlW40fJS3V5VhhvJ6pWeaVa4sYM/BU9hz1y1d25nQi3MV7uxKWCyV60ocB60MlltZXJyr1NI+7na13KjsXlyRU9eZmili7tp8w9edPa3DDOWN0pY2jFK50lXTgILYE3ucxkYKGN+6tnpAR6sr650bh6GQurSPcyVO3Y+BnLqKvYp29wrJ56y6muMwQ3m7pSNhJwiAFcsG8ICry6DX3sIzx855pqk+9+wJzz2HTlUTkRlTK9RVTKvoFcsH6laH7jRBVqQuRz42UmhbZUi7BdWCW1kBFLXNW2eO3+v3a7qWsyGYs1IoqJqI2o8rcuoK9irPtIr2yonbaYKcla0FbGf1Si8GcaA+8OasDHZuHK7bhFyxbKBhQ9deYUf9lOLeQA67yUztxUBOHef8uG9iyombAotfu1UBGtIyLRioE0lWUAvE+ZwV+ufKS8OP7Wk8AOrKKJ383tj8Xr7zTTTMJjO1H1Mr1HFBm5J+g3tNAcSvH7cCuPv2Al48fb42U/KdFvTvjuLdN1iY3X0HAGDNxKFIP/v0sXM49PLbuHRlvmElHkbOyuLu26u9X7yCvfNNNKiaiDqDK3LqOL/VnKmxkp2KaTZ5sv94EeNb1+KJHRtQmqtgscnrJMW5im4mKF6cqzQVxO3f76Nj6/H4vbcFbiCH2WSm9uOKnDrOtMqzOwUC9cfq84NW06tPm10XfnV+MfEmUs2aminWtQloNefvF6jfQDa1PQjzGGo/9iOnjvPqU21XZhTyOaz+5Rz+z6sXuibg9gL2AU8nUz9yrsipJR6eOlnLuWZF8NmP3WQ8zt1w4hDXKzOKpXJf14MnJZ+zsGL5AFfRPYqBnBL38NTJusk39sxLAL7BfGyk4FuCSM3JWVnsuesWBu4exs1OSpxp8o3fRBwbg3g4YaslOYWnP3BFTrF49fY21Sv71TFPzRR9a7/purDzPPM5q24zk3oXAzk1zXRc23SEPCviGfinz17AM8fOcTMzpLD1Ce2qfqHOYyCn0NxBeO7avOepyhXLsnWTdWwbP7iyIfCPP3cClQWGcKdWzNWk3sYcOYXi1TXPdHpy7toCdm4crg1gyIpUZ0f+c7kh8DOI18uKJBbEVw6GP+pP6ZbIilxE7gTwJwCyAP6Xqk4mcV3qHlF6e9+Yz+HRsfUNFSqrIx4970dJNfqysoLdn7olkWtR94sdyEUkC+B/AvhNAG8C+KGIHFTVn8S9NnWWM5USNrzYA4bXTBxqqFfu1baySYryO3LWhr8nZ0EEKM1VWCfeh5JYkX8UwE9V9WcAICJ/BeDTABjIU8zrtKUXZz530MrUDRh296pmEPdnZQU7fu0m7D9eDPV7Z2042ZII5AUAzgLhNwF8LIHrUgeFTaU4Q3O50ti3xNmrmpt4/lYsG8CjY+sxevMq35FynF5PbkkEcq+zCQ3/vYrI/QDuB4Dh4eEEnpZaya8joQDIeKQATEHaXpkziPsrlSvYNHkEby1NPPL6BJMVweP33sYgTnWSqFp5E8BNjj9/AMBb7gep6lOqOqqqo0NDQwk8LbWSqZVqPmfhtcltWIyYJgla3ecsFlAJUKsK8griOSvLIE6ekviv54cAPiwia0RkGYDfBnAwgetSB41vXQvLY2zO5WvzmJopJjpIoJDP4ZU/+i1s+tCqxK6ZNn6HqOzJQTxqTyaxA7mqzgP4fQCHAbwC4FlV5VnrlBsbKeBdNzRm3ioLir2Hz/gOErCDj11HHmR861pMzRTx9z+70Oztplohn/MdgPza5DYcndjCIE5GiXyeVdXvquqvquqHVPXLSVyTOq9kOPATNJ9xcSn4hEm/CIDpsxfw4IGTiDEnIhXcn3ByVhZP7tiAoxNbjG96Yd8Mqb8xMUlGpvRJRsS3qsL+uTDpFwXwzLFzoQ8bpdnee26rm3bvTJU002iMyMZeK2Q0vnWtZy15UHCx0y7jW9eG6qXSD6GqkM/Veq6bvm8ad0cUhCtyMhobKeCx7etrq8gwH/NXLMti7+EzWDNxCHsPn8GAx4ZpP1r9y/4BmUONKQ4GcvI1NlLA0YktoXPe1+YX6xprlSudnk/fWmHfpo6+egFTM0Xj951vmkD1TdM+TOX3c0QAAzlFECbnHWeyfRo9sWND6C6DQYMzxkYKtZW5nb6yD1MxmJMfBnIKzevjf9JWDlq1VWlakjJXQn7quDhXCQzIXq0RnG0OiLwwkFModifEcmUhckncykEr9M9crSzUjqh3+9p+5aAVqb0vgMDVtam0M6jkk/obAzkFcg6VAKpVK1FW5ttufT8ev/e2UD8zt9R4q9vL7ux+31EDbNDq2pS+SvIkLfUelh8SgGqw3nPwVG3O48pBC7s/VW2Tavq4H9bTx87h0Mtv498OvwdHX+2N05srllX/07nRUDboxy/4e5V8snqFgjCQE6Zmihj/5om6jcqLcxWMP3cCQDIf6y/OVXomiAPVToUPHjiJu28vhOof7uS3urbrzN0Dqnk8n/wwkPcBr8n1zsCw9/AZz2oTu69KM6vOflCuLODF0+fx2Pb12Hv4TKjfUZjVtd/BISIvzJH3OK+hye4NN78V91ulcluqVdLqrVK5Vmv/+uQ235OYKwctdjCklmAg76CpmSI2TR7BmolD2DR5JFatsOlaYcrZ/D7q55uozOgn7t/d+Na1xrJJVdROvcb9+yZyYiDvkDAr5SSuFaaczdR7PJsRXLoy3/dpFdOBHysrDWmSsZGCsWyyVK4k8vdN5MZA3iFJHvzwu1aYcraxkQL23nMb8rnrAWvloIV3Lx/ou5OaboV8DoPLvLeSViwb8EyThG10xYM+lBQG8g5J8uCH37W8VttWxnslObv7Drw+uQ2vT27DzBfvwDtl737k/cLemDT9fk2/nyh7CjzoQ0lg1UqHmCpBmjn4EXgtd9ZEqsMcvCpZnBUuXgOWe5099Ljg+J2YKlJMf1deJYRz1+Zx0WNQBw/6UBIYyDskyYMfm9cN4Zlj5+pys/a19h4+09APvLKgdY+387XTZy/U1UQHBXErA/Rac0P71KqzRLOZvyt3CaG9j8GDPtQKDORt5lzx5gctLB/I4J1ypemDH1MzRew/XmzYYBMoHtg3a9x4c3+9XFnAN156I9IK3C+Ib/rQKvzo3DuprHaxc9f230USh3R40IdaSbQDH51HR0d1enq67c/baaZVWZza4k2TR7qyqsQ08SYtBNUWte7ACzAYU+eIyHFVHW34OgN5+5iCbiGfw9GJLb4/azqduWbiUNd3CUyjfM7C1fnFujddKytYWFA4P4hYGcHee25jMKe2MAVyplbaKEylilfABlC3krdz2kBzTZu85KxsKtMgSXC/9pyVhUhjYzCv2aOVRcWeg6d8A3lQiwSiuFh+2EZBNd2mgz17Dp4y1onHPT7vnOjej4N+na/d+bsoeVSYmJR8yjSTPPhFZMIVeRsFVT+YDvaYVsp2nw/7Z+0N1EtX5kMd5LFTOvaKsVgqQ9CeqfY5K4NrC4qFFhw4GrQymAtRTuOsTnGvkMM2wQrid1iLq3JKSqwVuYjcIyKnRGRRRBryNlTPPZXeXv3Z/0FHPRxir+SdA5JnvngH9t5zW+DP2kHMPTRC0Z4Ra+XKYkuCeM7K4n9svxU7Nw4bvp/x/N27Rfmk4zezkxN/qB3irsh/DGA7gD9P4F76gl+L0ij5br8aZL9DLADqDrtsmjzSsGJUACLVJk9p4h6G4WXViuV1G8um/LVXueDmdUPY98M36nLl9qQgkyQPfhGZxArkqvoKAEjEGY7kzSv14qUQYsNsfOtajD93omGDzj6eH/QpIG1BHKgfguy3EjalkpybyHYwd/+OR29eFWnjkhN/qB0SKT8UkR8A+ENVNdYUisj9AO4HgOHh4dvPnj0b+3l7kTPIeAlTqmjb8Mjznhtxzmt0ax16s+zXZnpdKwctXKks+r5ZRvkdh8GqFUpK0+WHIvJ9AO/z+NZDqvrtsDegqk8BeAqo1pGH/bl+Y68CkzjSbWrq5G5h6x7zlmbFUhkPT53E3LX5hu9ZGUGpXAn8tJF0/poTf6jVAgO5qn68HTfSL8KuzpI40h0mPzs2UsAj3znl2dDJTzYjLdmsTMLTx841fC1nZTC/qNAQvWGYv6a0YflhG7lX2e6crFvclVxQftZ+U4kaxIFquVOajg9dm9dQfWSYv6Y0ilt++BkReRPAvwNwSEQOJ3NbvSnJYRJh+JU7Pjx1Eg/sm206Px43FZPNCPI5C4Jq61gTEdQNvGiWXxC3nz2oJJGoW8WtWvkWgG8ldC89rxM1xc5Vvb0C37VvtmXPF9bComLF8gHM7r4DUzNFY6fGf//BVbhndDhUNY+fjABe7z1ZETx+L3ulULrxiH4bhRm71irugz/dwH4DGxspGE8hHfvZxWSGP2u15tspZ2UZxKknMEfeRu2sKXZvql6+Ot90MMy2aFLQjflc7T5Nl19QDXzzCdNWYBHALy0bwIrlAywDpJ7DQN5G7RguMDVTxJ6Dp+rqx+Oswq2MtKw0cfO6ocCUSWbphKnpDnJWFnffXqibbGTyTrmC2d13xLhjou7EQN5mragpblXTq3zOwp67bkmsgZTbi6fPBwbf5QMZlA0NsOz7Gxsp1J24NM0aZVkh9SoG8i4TVGfu/v7mdUN1q9GkgnhWBHvuqvYQuXD5akJXvS6fs3w3ebMi2PjBlTj66gXjY5yra/emLo/FUz9hIO8iQXXmXt93D11OyoIqxp870TARpxkZoOEapXLFmHu3e8nYr92LX+90zsekfsNA3kWCeld7fb+VZyu9JuI04z2DFlQbBzCYNlA3rxsKrFSZuzaPqZmiMTjzWDz1EwbyLhDUKMtOQTRbb97plrQX5yqRphi9ePp84Gu9OFfxPRVL1E9YR95hYeq77U26ZjbrVg5aeOLeDbHGwcWVFYlU+minQ4K08lQsUZowkHdYUArBuUnXzHzO0lyl7qh+uwn8j8d7cQ6dDsJJO0QM5LFNzRSxafII1kwcwqbJI5GH6gYFonJlAbv2zWL1xCF84cDLuPv2Qq13il+PEpui2nMcAI5ObMGTO5Jdnecs/3+FomZ0nHM0w/RYYUkhEQN5LElMSI8SiOYqi/jLl85hfOtavDa5DY/fe1uooOy8L3t1bnoLEABP7tgQevV+g5Vt+o1h5aCFnRuHjTNM99x1i++1WVJIVJXIhKCoRkdHdXraOEwoNUxTaKJMmPFrGGXivL6zrhwBm5rOn1szccj4nIUIs0OBauC37yHK63h9clvgY5yv7z05CyLVdBFLCqkfNT0hiMyS6GY4NlLA9NkLkerBi6UyNk0eaRgUvGbiUOj7NQ2dEEQ/0j999kLk0XErBy1smjwSasAGgzWRP6ZWYkiqm+GjY+vxRIR0BuCdxgl6XjtfPjVT9Nw4bfZ4/zdeeqP2z2E2ZK2s4NKV+VgpKSK6joE8Bq+g1WzedmykgPGtayMNUXCX34UJos7TonffXqhtmGZFmj5c5KxKsXPwKwe9X8fKQQsrlg00NOJiKSFR85haiaGZo+CmXipe/UHCcKZLnPdTLJWNR+DLlQU88p1TuFJZrH0/Tptad/WMfQrVa4Tc4LKBjgzYIP2xrpMAAApnSURBVOplDOQxRcnh+vVSaXZ4gjud4r4f06Zm1DmdAuBX/tUK/N+fX2743mc/dlPD1/yCdZih0EQUHlMrS+LWg4fh10ulmdVomDROnOBor7QL+Rye2LEBL/z338DOjcN16ZidG4fx6Nj60M9rfwpJKiVFRFyRA4g+3b5ZzaxSTQoBaRy/HuU5K4vlA5mGJlZOAuDVxz5Ru84D+2ax9/AZjG9d6xm43fymIbE7IVGyWEeOZOrB4z6PV+AzVZHYjzcFQq98u30t+2enz17A08fOGe/V3qx0p2ByVjb0pPmg3upEFA3ryH20YvPNK4gFzey8wcrUvmcHXq/VtHtEmvsTxJ6Dpzzb3TrfmPwqRDJizqE72+oGYQ04UXswR47kp9ubju4DqGteZXcF3HPwFMafO1EXPNXx/3ZNiH2E3WtEmh1gp2aKxpSJ843J700qaEQnq0uIugsDOZKtBweCB0TYz2eX/JXKFd8hDs7V9NhIwRhIi6Uydu2bNV7H+cZkepMK04iL1SVE3SVWIBeRvSJyWkReFpFviUg+qRtrJ2ebV3fzpmaqWYJSNc2UGrqP1zdj87qh2j+b3ryC6slZXULUfeKuyF8A8BFVvRXAPwJ4MP4tdcbYSAFHJ7bgtclttZVvs90Ng1I1zaQmnNdspi85ADxz7BxWL70hAfB88/JrE5DPWaE3OomofWJtdqrq844/HgPwn+PdTncJSpF4mZop4vLV+YavO1eyUUsNBaj9rL2JWq4sGE9umtiPtN+QHtu+3rMqx6vi5T5DvTgRdV6SOfLfBfA90zdF5H4RmRaR6fPnzyf4tK0TtZrl4amTeGDfbMNm48rB+pXs+Na1xn7gXhRo+IQAVI/V56xspP4sNlNvE6800xM7NiQSxNtx6IqoHwWuyEXk+wDe5/Gth1T120uPeQjAPIBnTNdR1acAPAVU68ibuts2M62cMyJYM3GooVeKqRXt4LKBuhX82EjBd1PSzU53mD4h3GBlkLOysfLuTq0oG2zXoSuifhS4IlfVj6vqRzz+Zwfx3wHwSQD3aSdOF7WQKRe9oFrLmY8/dwIbHnkeu3yGQ3gFzLAta50pGVPgLc1VGlbROzcOB+bR21l94pemIqJ4YuXIReROAJ8H8B9VdS6ZW+oe7qPkGY+cdGVBfY+6A9cDpnvajUlWBIuqtYk49vH4/KDleVDnxnzOcxU9evMqfO7ZE555dGfevR3Y8ZCodeKe7PwzAMsBvCDV+uNjqvp7se+qi0SZwGMyd20eD0+dxP7jxdqqtFSuIANg0fVY+wg8gIZUhJURWFmpqzn3Kwe079u0ednOlAY7HhK1TtyqlV9J6kbSIGq1ie3iXMWzr8kiqiV9K5YPNPQj2TR5pCEVUVlUz8cDMI5N65YGVUHtCYioeey1EoFXMIrrnXIFs7vvaPi6KeXgfnyYTcRu6HnSLW8oRL2IgTwCdzDKD1q4dGW+bmxZ1OoRe46mu3uhVz4eaExFNFPr3ind8IZC1IsYyCNyByOvLod2H/CwnKtooJrT9griXqkIbiISEQN5TF6rTL9e335zNO1SPK8VfVbE83g8NxGJiN0PW+DF094nVwXe8y1tb5XKxpX0oqpnWoJj04iIK/IWMAVjBbD/uPlYur2KjrLC5iYiETGQt4Ap3WEPkvAiqLaZHb15VeQyvTCbiBy7RtS7mFppgWZ6fdur9emzF3CDdf2vJYnWsc224yWidGAgbwHToIqg6TvlygKeOXau7hj+1Xn32c/o2OeEqLcxtRKDX7rCK90RpuOhe82eRE04SxSJelvfBfKkcsXNtGUtNHnEP27AZYkiUW/rq9RKkrniZtIVQSPaTImXuAGXJYpEva2vAnmSueJm0hXu3Hk+Z2HloFXLo9/n0UM8iYDrN1yaiNKvr1IrSeaKm01XBJUKjt68qiVlguxzQtS7UhPIk8htRxndFqRVbVkZcIkoqlSkVpLKbYcZ3Rb2ukxXEFG3kE6M2RwdHdXp6enQj980ecRzJV3I53B0Ykuk53au7E2tYpu5LhFRq4nIcVUddX89FamVJHPbYUa3sb6aiNIkFakV0wZi3LK8Vl2XiKidUhHIW1UHzfpqIuoFqUittKpVK1vAElEvSMVmJxERmTc7U5FaISIis1iBXET+SEReFpFZEXleRG5M6saIiCicuCvyvap6q6puAPDXAL6YwD0REVEEsQK5qv7C8ccVaGynTURELRa7akVEvgzgvwB4B8Bmn8fdD+B+ABgeHo77tEREtCSwakVEvg/gfR7fekhVv+143IMAblDV3YFPKnIewNmI92ryXgD/lNC1OiHt9w+k/zWk/f6B9L+GtN8/0J7XcLOqDrm/mFj5oYjcDOCQqn4kkQuGf95pr3KctEj7/QPpfw1pv38g/a8h7fcPdPY1xK1a+bDjj3cBOB3vdoiIKKq4OfJJEVkLYBHVVMnvxb8lIiKKIlYgV9W7k7qRGJ7q9A3ElPb7B9L/GtJ+/0D6X0Pa7x/o4GvoyBF9IiJKDo/oExGlHAM5EVHKpT6Q90K/FxHZKyKnl17Ht0Qk3+l7ikJE7hGRUyKyKCKpKiETkTtF5IyI/FREJjp9P1GJyNdE5Oci8uNO30szROQmEXlRRF5Z+nfoDzp9T1GIyA0i8g8icmLp/h/pyH2kPUcuIr9ktwoQkf8G4N+oaqqqZ0TkDgBHVHVeRP4YAFT18x2+rdBE5F+jWrn05wD+UFVT0aNYRLIA/hHAbwJ4E8APAXxWVX/S0RuLQET+A4BLAL7e7jMcSRCR9wN4v6r+SETeDeA4gLG0/B2IiABYoaqXRMQC8HcA/kBVj7XzPlK/Iu+Ffi+q+ryqzi/98RiAD3TyfqJS1VdU9Uyn76MJHwXwU1X9mapeA/BXAD7d4XuKRFX/BsCFTt9Hs1T1bVX90dI//wuAVwCkZrKLVl1a+qO19L+2x6DUB3Kg2u9FRN4AcB/S34HxdwF8r9M30ScKAN5w/PlNpCiI9BoRWQ1gBMBLnb2TaEQkKyKzAH4O4AVVbfv9pyKQi8j3ReTHHv/7NACo6kOqehOAZwD8fmfv1lvQa1h6zEMA5lF9HV0lzP2nkHh8LXWf6HqBiLwLwH4Au1yfsrueqi4stfL+AICPikjbU1ypmNmpqh8P+dC/BHAIQGDjrnYLeg0i8jsAPgngP2kXblxE+DtIkzcB3OT48wcAvNWhe+lbS7nl/QCeUdUDnb6fZqlqSUR+AOBOAG3dfE7FitxPL/R7EZE7AXwewF2qOtfp++kjPwTwYRFZIyLLAPw2gIMdvqe+srRZ+FUAr6jqVzp9P1GJyJBdZSYiOQAfRwdiUC9UrewHUNfvRVWLnb2raETkpwCWA/jnpS8dS1PljYh8BsCfAhgCUAIwq6pbO3tX4YjIJwA8CSAL4Guq+uUO31IkIvINAL+BagvV/wdgt6p+taM3FYGI/DqAvwVwEtX/hgHgC6r63c7dVXgiciuAv0D1358MgGdV9Uttv4+0B3Iion6X+tQKEVG/YyAnIko5BnIiopRjICciSjkGciKilGMgJyJKOQZyIqKU+/9l6od4tD270QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n = 1000\n",
    "d = 2\n",
    "\n",
    "Sigma = np.array([[1,0.9],[0.9,1]])\n",
    "mu = np.zeros(2)\n",
    "print(mu)\n",
    "\n",
    "eigvals, eigvecs = np.linalg.eig(Sigma)\n",
    "eigvals = eigvals.real\n",
    "\n",
    "Lambda = np.diag(eigvals)\n",
    "Lambda = np.sqrt(Lambda)\n",
    "\n",
    "P = eigvecs\n",
    "\n",
    "Q = P @ Lambda @ P.T\n",
    "Z = randn(n,d)\n",
    "X = Z @ Q\n",
    "\n",
    "plt.scatter(X[:,0], X[:,1], marker='o')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Funkcja"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x57a0f10>"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAc1ElEQVR4nO3dYYxc1XUH8P/Z8QO/NYnHLhulDF7s0mhdHAevWIGT/ZCapFkSwNmaUIeGKFKkWJUaCVOyqV1QMFIir7Rq4EMiVbSpFAk3MWCyceq0DpHdD0U1zS67i+VgJ0DAMFDhCC9JvAOe3T39MPvWM2/em3kz7828d+f9f1ISPLueuR7is3fOPedcUVUQEZG5uuJeABERhcNATkRkOAZyIiLDMZATERmOgZyIyHAr4njRK6+8UtevXx/HSxMRGWtycvK3qtrjfjyWQL5+/XpMTEzE8dJERMYSkVe9HmdqhYjIcAzkRESGYyAnIjIcAzkRkeEYyImIDBdL1QoRkQnGp/IYO3oGb8wWcFXWxshQH4b7c3EvqwoDORGRh/GpPPY+dRKF4gIAID9bwN6nTgJA4oI5UytERB7Gjp5ZDuKOQnEBY0fPxLQif9yRExF5eGO24Pl4fraAwdFjDadbWpmmYSAnIvJwVdZG3iOYC7D8uJNumXj1bRw/fc43SLc6TcPUChGRh5GhPlgZqXrcfadaobiAAyfOIj9bgOJSkB6fyi9/T6vTNAzkRER+At6E6RXcy4O0X5rG7/FGMZATEXkYO3oGxcXm7zQuD9JXZW3P7/F7vFEM5EREHsLulsuD9MhQH2wrU/F128pgZKgv1Gs4GMiJiDyE2S27g/Rwfw77d2xGLmtDAOSyNvbv2MyqFSKiVhoZ6quoNAFKFSteyZasbWHV5StqlhYO9+da1kgUWSAXkQyACQB5Vb0tquclImo3p+a7UFxARgQLqshlbWzb2INDk/mK4G5bGezbvinWbs8od+T3AHgBwPsjfE4iosgEacpx13wvqC6nSob7cxi4Zm3i5q+IavOnsstPInI1gO8D+BaAv6u3Ix8YGFBe9UZE7TQ+lcfIEzMVlShWl2Dnjesqmnlm5y7iwsWFqt+fy9p4Zs/N7VxyFRGZVNUB9+NR7cgfAfB1AO+rsYBdAHYBQG9vb0QvS0QUzL7Dp6rKCYuLisdOnF3+tVcnZ5CvxS101YqI3AbgLVWdrPV9qvqoqg6o6kBPT9Ul0EREgYxP5TE4egwb9hzB4Oixig7KWmYLxVCvK0uvnURRlB8OAtguIq8A+CGAm0XksQiel4iogpO/rtUO3yoKJHLyIRBBIFfVvap6taquB/B5AMdU9e7QKyMicml0Zkn57l2qx6Y0LKqW+qixjpyIjFFrZom7IqWqVNCjriPTJegCArfir7atJlfeWpF2dqrqf7GGnIhaxa/bUgHsPjhdkXI5cOJs1e4dADIiy92V/3jn9Ri78/qKjsu7t/b67t6j2NW3AnfkRGQMr25LP3577AVV3L21F8dPn8O9B6c9a8EPlFWylJudC3dg2ioM5ERkDCfYOimUZrtg3CWH7kse/C6ViGpaYdQYyIkoVrW6Lf2+5nx9/Z4jkazBOTB1ntdr5x/ltMKoMZATUWxqXYEGoO71aH5DrJrhHKT6zVlJQiu+n0ha9BvFFn2idHOCpV+3ZG4pheH19TXdFlTDN/i4OUHb/cPBtjKRjpwNo9Ut+kREvspTJKttCxcuzqO44L+JrFWvfb5FB44LS5tav2vbkhDI/TCQE1FLudMnQXbSiks75CRIaiOQgzcEEVFLeXVjBpGUIA4kt1rFwR05EQUWZJ63W1S72VzWxoX35iPPjdc7ME1ytYqDgZyIAqlXYeIO8EBpdGxU++o4xsgmvVrFwaoVIgpkcPSYZzD12tFaGcHCgmKxLStrDQHwm9Fb415GBVatEFEofjtir61grYqUdoiivjzpefFyPOwkokASOi+qitUl+MLWcLeQmZAXL8cdOVFKNXJwOT6VjyzX3WoLqhi4Zm3FPJVGmJIXL8dATpRCtQ4uveacdCV1fquHRS2NtG1GEi5YbgZTK0QpVO+mHfeVas3UdGe6zAn+jqQ3/vhhICdKoVo37QDNN/GUWwh4606SmHTAWY6pFaIOVSsHvtq2PBtrFMCGvUeQoKbKtjHtgLMcAzlRB6rXvHPh4rzv701LEM/aFlZdvqKhLtWkYiAn6kD1cuBx13nHzbYy2Ld9U0XgHp/KY3D0WFVgb2YsQbsxkBMlVJgAUi8HnmarLsvgW3+5uSqIe32CmXj1bRyazNes7kkCtugTJZA7sAClXeQdN+Rw/PS5usHdr52eSj70gVV4+dwcFlSREcFKqwsXLlYf7vqN0o2rTJEt+kQG8UuNHDhxdrkxx293OD6Vx/kL77VrqUb69VsXlv95QdUziDtf85K0TzYsPyRKIL9A4Xd7jcPZyc8VTR5XlRwZn0aopJUpMpATJVAjgaI86EdR/00ltpXBXTetg21lqh5PWpkiAzlRAo0M9VUFEL8+yfKgn7SP/CYRlHLfzv/u37EZ3xzejP07Nlc9nqSDToA5cqJEcgJFedXKto09FRUUQPXu8KqszUNOl0yXBOoy/di1a3HgKx+teny4P5e4wO3GQE6UUO4AMj6Vx5Hn36wI5CutruWvjR09g/xsIZJZ3J3krhvXYeCatVVVQG7PnX0H41P5xAdtLwzkRAbwKkcEgPNzRYw8MQPIpSYfBvFKx0+fwzeHNwNAxTRHd0WKc3DMQE5ELVHrELNo4HCqdnLODco/4WzYc6Tm95qGh51ECffA+EnmvUPwqgDyqwpKWllhUAzkRG3kzPPYsOcIBkePYXwqX/P7v/DP/9P0TTdUugTaq1TQqyooiWWFQTG1QtQmQW7lcX//My+93dY1mi7rGs+76jLvEOdVFZTEYVhBcdYKUQiNDLbym3/iHqe6/o9snHj5fFO38qTdIzu3eM6oSWLtdzP8Zq0wtULUJPd1aM4O2y9d4neQNlsoVjzHMy+9zSDeBBFg71PP1xzf26kYyImaVG/mt5upB2mmUAUKPjNmTK1GCYqBnKhJjc78Hhnqg2XghcRJNHjtWt+BVl46/YcoAzlRk5oqYWMcj0SjZwimVqMEFTqQi8g6ETkuIi+IyCkRuSeKhRElXaMlbGNHz6T+irWoNBLE13RbHXHQWUsU5YfzAO5T1edE5H0AJkXkaVX9ZQTPTZRYQUrYymegUPvZVgYP3r4p7mW0XOhArqpvAnhz6Z9/LyIvAMgBYCCnjldrMp7ffBRqjdzShMggV+F1mkgbgkRkPYB+AM96fG0XgF0A0NvbG+XLEiWCu6Z87uI8g3ibxHWHZlJEdtgpIlcAOARgt6r+zv11VX1UVQdUdaCnpyeqlyVKBK+a8vNzxbq/j8IzubU+KpEEchGxUAriB1T1qSiek8gk+w6f4u67jZzSw4zIcu1+vbk1nSyKqhUB8D0AL6jqt8Mvicgs41P5ivke1Hp/0tMN28osV6/U66rtdFHsyAcBfBHAzSIyvfSfz0TwvERG6PT27yT69VsXUtmK7yeKqpX/BtscKMU6vf3bJGn9d8ExtpQ6jUwsDPI8bPFJjk5vxffDQE6p0uhMcPfv3Xf4FPPhCZXm6hXOWqFUaXRioWN8Ko+RJ2YYxBNmTbcFQamOvFNmjjeDO3JKlUYnFgKlIH7f4zOcEZ5A7xYX8fDOLakN4A7uyClVGp1Y+MD4Sdx7cJpBPAZBKijSXKlSjoGcUsVrYqGglCt3X4bsXHzMEN5+uayNh3duQS5rL6dO/KS1UqUcUyuUKuUTC52JhE6gzs8WcO/BaUy8WrrwmBcfx8M5tHQPJPO78zStlSrluCOn1Bnuz2FkqM/zo7sCeOzEWTx24my7l0WofWjZ6Pz3NOGOnFKJ9d/Jk7WtmhMMg8x/TyvRGA5xBgYGdGJiou2vS8SLHpJNAAboGkRkUlUH3I9zR04drzx4C8CdeII5I4CDNmlRCXPk1NHK54QDDOKmYFlhY7gjJ2MFmZni1clJZmBZYXAM5GSkejNTmAs3H8sKg2MgJyPVm5nCS4/N4j67YFlhYxjIyUh+H7vzswXsPjjd5tVQWIpLwTzHqpWGMZCTka7K2kybdBgniNeqJSdvrFohI3l1+ZH5eMDZHAZyMtJwfw77d2yuOUyJkse2MnhkaRiWFx5wNoeBnIw13J/Dto09cS+DGuDMUeHclGgxR07GGp/K4wCHWxnj7q29yweYnJsSLQZyMhYHXyVblwCLCmREcNdN6/DN4c0VX3ePqaXmMZCTMR4YP4kfPPsab+tJOJYPth8DOSVGrZb7B8ZPckZ4wnUJ8O2/4v2ZceAYW0oEd8t9uVzWxhvvFMCNePKxDry1OMaWEinITBQ2/piD/67iwUBObeNOnWzb2INDk3nOROkgGfG6QI9ajYGc2sJrWuEB3lDfcXgQHQ82BFFbeE0r5F/5zsNO23gwkFNbcIZG52NnZnwYyKktVtuW5+PMqHaGXNZebr+n9mOOnFpufCqP37837/m1j127Fs+89HabV0RREoAlhzHjjpxabuzoGSwsemfEGcTNx4mF8WMgp5ZjftxcVkaQtS0IgKxtwcpUJsOYF08Gplao5bLdFs7PFeNeBjXIa2ZKrTEKFB8Gcmqp8ak83mEQN5JX3psTC5OJqRVqmfGpPO57fAaLcS+EGsYOTbNEEshF5BYROSMiL4rIniiek8zmdHKy089Md920Lu4lUANCB3IRyQD4LoBPA7gOwF0icl3Y5yWzeXVykhnu3tpbdQkEJVsUOfIbAbyoqi8DgIj8EMBnAfwygucmg5QfhHEfbh7byrCpx1BRBPIcgNfKfv06gJsieF4ySK154pR8WdvCvu2bGMQNFUUg9zoVqdqQicguALsAoLe3N4KXpTj4lZ/tO3yKQdxgqy5fwSBusCgOO18HUH4ycjWAN9zfpKqPquqAqg709PRE8LLUbs6uO7+UOsnPFrD74DQ2feM/MVtgiaHJ2LRltigC+S8AfEhENojIZQA+D+BwBM9LCeN3gHnhInfipmObvdlCp1ZUdV5EvgrgKIAMgH9V1VOhV0axc6dReI2X2TIiWFCFoDL3yTZ780XS2amqPwXw0yiei5LB60YfdwAgc5RXpLDNvvOwRZ88+d3ow2BuHvfMFLbZdx4GcvLkd/ilKJWq8XAz+XJZm3PCU4KzVsiT3+FXLmtj3/ZNbV4NNYp573ThjpyWjU/lse/wKd/dtm1lsG1jD3YfnG7zyiiIjAgWVZn3TiEGcgJQCuIjT8yg6HOTTy5rY9vGHvzbs2fbvDIKgu316cbUCgEoHW76BXEAmLs4j4P/+xpqfAvF6I4beICZZtyRE4D6nX284SfZjp8+F/cSKEbckROA0nVsZC622KcbAzlhfCqPP7w7H/cyKAS22KcbUysGi6pDr15+nJJNAJYaphwDuaG8Wuj3PnUSABoO5vxYboY13Rb+8O58xQ9dAfCFrb086Ew5BnJDebXQF4oLGDt6JtBfamc3z0FYyScAHt65hXNSyBcDuaH8dtFBdte8zccs5TtuzkkhLzzsNJTf4VaQQy9ejGwOXoRMQXBHbqiRob6qXXX5fA2vj+AAmE4xBO/QpEYwkBvK+QvulS/1OggdeXIGCwuKxTgXTXUJSjX9s3NFjB09A6Dxw2tKHwZyg/nlS71SJ8UFlhcmndUlgFzqog1TiUTpwhx5B2I5oXlyWRtXrFxR9QPXqUQiqoWBvIOMT+Wx5aGf8QYfg1gZwSM7t+CZPTdj1meeDX8wUz0M5B3CGUPLm3vMsabbwtjnrl9Om4SpRKJ0Y468Q7DN3gxrui1MfeNTnl8bGerDyJMzFekVKyNsv6e6uCPvEPz4nXxWRvDg7XWuyXP/LObPZgqAgbxDrLY5hjbJVl2WqUijePH6VFVcVB52Ul1MrXQIkbhXQF5yDcxDCTN2gdKNgbxD8AafZCkfdBXUVVnbs+uWh51UDwN5AtWacOf1tYlX3455xeSmQOBJlI56YxeI/DCQJ0ytOeMAqlvvn5hhtUpCNZoSqTV2gagWBvKEqTVn3PnncgziydVMSoRjaqkZDOQJ47eLy88WwPNMc7D+m9qJ5YcJ47eLE7DEMKlEANu69FfJ3bFJ1GrckSfMyFAf7j047dkXwvb75LG6BGN3MmhTvLgjT5jh/hyb+QyRtS0GcUoE7sgTZnwqH/cSKIBc1sYze26OexlEALgjT5x9h0/FvQRysTKVx8ys7aakYSBPGObBE0hLB5iC0k58/47NTKdQojCQE9VRXFR0X7YCD+/cAgC49+A0BkeP1UyDjU/lMTh6DBv2HKn7vURhMUeeMGu6Lc5NSSCnw9ar49a9O6/VncudPLUCd+QJc+tH/jjuJaROF0o/QGvJiNTsuC1XrzuXKGqhArmIjInIaRF5XkR+JCLZqBaWVsdPn4t7CamzutvCg7dvwiujt+KRnVtgW5mKr9tWBgvqXRTq1YnLcbTUbmF35E8D+LCqfgTArwDsDb8ks4XNjfIve/udnyti71MnMT6Vx3B/Dvt3bEYua1ccbuYauE+Td29Su4XKkavqz8p+eQLA58Itx2xR5EZX2xYrV2LgpD6coVVe/76CjpjlOFpqtygPO78M4KDfF0VkF4BdANDb2xvhyyaHX270vsdncO/B6YqxpF5zxQHgwsX5OJZOqP1pqJERsxxHS+0m6pP7W/4GkZ8D+KDHl+5X1R8vfc/9AAYA7NB6TwhgYGBAJyYmmlhusm3Yc6Rue71tZXDHDTkcmsxXBX2K1uC1a/Hc2XeqdsYrrS7PyiB2a1LSicikqg64H6+7I1fVT9Z54i8BuA3AJ4IE8U7md1VXuUJxAT949jXfwzOKxoc+sAoHvvJR308+TH1QJwmVWhGRWwD8PYCPq+pcNEsyl1du1AuDeOvNXVwEUPuiBqY+qFOEzZF/B8DlAJ6W0jXuJ1T1b0KvylBOIHjoJ6fY1BOzetU/vImHOknYqpU/jWohneTd4mLcS0iNjIjnJxyW+lGasLMzYl6VK9QatpXBXTet82zgYb6b0oSzViLGhp72yJXltQeuWct8N6UaA3nEglSuUDgCVJQJMt9NacfUSsRGhvqqPupTtJj/JqrEQB6Se7YKAOzfsTnmVXUu5r+JqjG1EoLXbJWRJ2ew6jK+ra2QtS3s276JaRQiF0acELwqVIoLyqFXEVuzNGaWAZzIGwN5CKxQiZagNP1RBJidK7IChSggBvIQWKESHQ6sImoeDztDGBnqg8S9iA7AA0yicBjIA/C79We4P1d3bC3V5tzAw/QJUfMYyOtwKlPyswUoLt364wRzvyvAKBjmwInCYyCvo96N6Ns29jC9EgJvlicKj4eddfhVpuRnC1i/5wgEYHoFQLfVhbkmpj6y8ocoPO7I66jXDs4gXjJXXESXAF0Nfjxhuz1ReAzkdXB2SnCLCqxc0YVc1oagtEuvFddZrUIUDaZW6ii/EZ014/XNFRfxS1c9uHNvZn62sHwRRI7NPkSRYSAPwBmTOjh6jMG8CRwzS9RaDOQevG5eH+7P8WAugDXdVtxLIEodBvIl5R//yytRnLpxgC359WS6BA/evinuZRClTqoDuV/wdleiOHXjI0N92H1wus2rjI8zNtbvfKD8PeOEQqL4pDaQu2eJ1ysjzM8WcO/BaYgAHpe2dwwB8JvRW6seL3+vgFLFCVvriZIhtYG8mdvudfm/OpdXXXd55Q4vOCZKntQGch5cVhOURg54YeUJUXKltiGoXkdhp85PcapKMlL9J1QAhybzywPBiMgMqQ3kXh2bTmjLZW08vHMLXhm9taOmG+ayNqa+8Sm8MnorXtr/Gc8/W/lAMCIyQ2pTK0HzviNDfRh5YgbFxXDJcSdoxlW+6NUO75deYtqJyCypDeRAsLzvcH8OD/3kFM7PhbtQOT9bwCM7t8RSvui0wwPA4Oix5R9c2W7L88/FQVZEZkl1IPfj7uwMG8SBUtrmoZ+cCr84ALbVhXeLi4EKaJy7MN3llvnZAqwugZURFBcuPRMHWRGZh4HcZXwqX5FKiSoVokAkPxAA4I4brgYAHDhxtmYwtzKyHJS9yi2Li4qsbWHV5StYVkhkMGMCud/8k6jtO3wqUD7c6hJcsXIFzs8Vlyf6Of/bao+dOIu7t/bi4Z1blt+T1baFi/MLy5c7uDst/fLe7xSKmH7wUy1fMxG1jhGB3Cst4Mw/iTqYzxbq75oFwM4b12HgmrUV61pQhZURQBH6cLSex06cxcA1a/GMa2SsH785McyHE5nPiPLDevdmtpsCOH76nHe6YkFxxcoVbZkC2Mif36vckvlwos5gRCBvZ5lc0AD8xmzB9/Vn54roviz4hx0BMHjt2tJuvgGN/PmH+3PYv2Pz8u09uazNWSlEHcKI1Eo70wIP3r4JI0/OVFRy+K0J8D4MvSprNxRkFcBzZ9+B1SV1X9drDUGxzZ6oMxmxI29nWmC4P4exz12/vHNd023B8rhR+MJ789i2scd3XbWCrG1Vv+2F4kJDt9BbXcK0CBEBMGRH3u7pe+6d6/hUvqopaLZQxKHJPO64IYfjp895rsuv+afQQMAGSnPBRS6VLwpKh6lOjpy7bKJ0E41huPbAwIBOTEy0/XXD8Luv02m48bJ+z5GGXiNrW3hvftF37re7esf9dSLqbCIyqaoD7scjSa2IyNdEREXkyiieL4lafeBqdQn2bd9U80AyadU7RJQMoVMrIrIOwF8AOBt+Ocm12rY8a8xX2/5VLlmf3+PlipUrlgO23+6aQ66IyEsUO/KHAXwdHX53jsf47pqPA8C+7Zs8D0q9zAZo3/c7QGVTD1G6hQrkIrIdQF5VZwJ87y4RmRCRiXPnzoV52Vj4BdpaAXi4P4exO6+vSJVkfXbwQYIxm3qIyEvd1IqI/BzABz2+dD+AfwAQaFCHqj4K4FGgdNjZwBoTodladq8KGK8DyyDBmHdnEpGXuoFcVT/p9biIbAawAcCMlPILVwN4TkRuVNX/i3SVCTAy1Nd0AC4XNhizqYeI3Jo+7FTVkwA+4PxaRF4BMKCqv41gXYkT5W6YwZiIomREQ1BSMAATURJFFshVdX1Uz0VERMEZMWuFiIj8MZATERmOgZyIyHAM5EREhotl+qGInAPwattfOJgrAXRkCWUT+F6U8H0o4ftQEuf7cI2q9rgfjCWQJ5mITHiNiUwjvhclfB9K+D6UJPF9YGqFiMhwDORERIZjIK/2aNwLSBC+FyV8H0r4PpQk7n1gjpyIyHDckRMRGY6BnIjIcAzkNaThUulaRGRMRE6LyPMi8iMRyca9pnYSkVtE5IyIvCgie+JeTxxEZJ2IHBeRF0TklIjcE/ea4iQiGRGZEpF/j3st5RjIfaTlUuk6ngbwYVX9CIBfAdgb83raRkQyAL4L4NMArgNwl4hcF++qYjEP4D5V/TMAWwH8bUrfB8c9AF6IexFuDOT+UnGpdC2q+jNVnV/65QmUboFKixsBvKiqL6vqRQA/BPDZmNfUdqr6pqo+t/TPv0cpiKVyKL+IXA3gVgD/Evda3BjIPTRyqXSKfBnAf8S9iDbKAXit7NevI6UBzCEi6wH0A3g23pXE5hGUNneLcS/ELbU3BEV1qbTpar0Pqvrjpe+5H6WP2AfaubaYicdjqf10JiJXADgEYLeq/i7u9bSbiNwG4C1VnRSRP497PW6pDeS8VLrE731wiMiXANwG4BOarqaD1wGsK/v11QDeiGktsRIRC6UgfkBVn4p7PTEZBLBdRD4DYCWA94vIY6p6d8zrAsCGoLo6/VLpWkTkFgDfBvBxVT0X93raSURWoHTA+wkAeQC/APDXqnoq1oW1mZR2M98H8Laq7o57PUmwtCP/mqreFvdaHMyRUy3fAfA+AE+LyLSI/FPcC2qXpUPerwI4itIB3+NpC+JLBgF8EcDNS/8fmF7alVKCcEdORGQ47siJiAzHQE5EZDgGciIiwzGQExEZjoGciMhwDORERIZjICciMtz/A1Bk0Db6qqJjAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def rmvn_eigen(n, mu, Sigma):\n",
    "    \n",
    "    d = mu.shape[0]\n",
    "    \n",
    "    eigvals, eigvecs = np.linalg.eig(Sigma)\n",
    "    eigvals = eigvals.real\n",
    "\n",
    "    Lambda = np.diag(eigvals)\n",
    "    Lambda = np.sqrt(Lambda)\n",
    "\n",
    "    P = eigvecs\n",
    "\n",
    "    Q = P @ Lambda @ P.T\n",
    "    Z = randn(n, d)\n",
    "    X = Z @ Q\n",
    "    \n",
    "    return(X)\n",
    "\n",
    "Sigma = np.array([[1,0.9],[0.9,1]])\n",
    "mu = np.zeros(2)\n",
    "\n",
    "X = rmvn_eigen(100000, mu, Sigma)\n",
    "plt.scatter(X[:,0], X[:,1], marker='o')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Faktoryzacja - Metoda SVD"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### W metodzie SVD dla każdej macierzy rzeczywistej $\\Sigma$ można znaleźć jej rozkład na trzy macierze $\\Sigma^{\\frac{1}{2}}=UD^{\\frac{1}{2}}V^{T}$, gdzie\n",
    "1. $D$ - wektor wartości osobliwych macierzy $\\Sigma$;\n",
    "2. $U$ i $V$ - macierze ortonormalne tzn. macierze te są ortogonalne ($U^{−1}=U^{T}$) i długość każdego wektora w macierzy ma długość 1.\n",
    "\n",
    "Ta metoda nie wykorzystuje informacji, że macierz $\\Sigma$ jest symetryczna. Dla macierzy symetrycznej $U=V$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Przykład 4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Wylosujemy próbę $1000$ elementową z dwuwymiarowego rozkładu normalnego o średniej $\\mu= (0,0)$ i macierzy kowariancji\n",
    "$\\Sigma=\\left[\\begin{array}{cc}1.0 & \\ \\ 0.9 \\\\ 0.9 & \\ \\ 1\\end{array}\\right]$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-0.70710678 -0.70710678]\n",
      " [-0.70710678  0.70710678]]\n",
      "[1.9 0.1]\n",
      "[[-0.70710678 -0.70710678]\n",
      " [-0.70710678  0.70710678]]\n"
     ]
    }
   ],
   "source": [
    "Sigma = np.array([[1,0.9],[0.9,1]])\n",
    "\n",
    "U, D, V = np.linalg.svd(Sigma)\n",
    "\n",
    "print(U)\n",
    "print(D)\n",
    "print(V)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Funkcja"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0xcd598b0>"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAaXElEQVR4nO3dbYxcV3kH8P8z42t8l5SMIa6ox97ape0agoktVslW+4E6QDaQxFgGYpIYISHVqgQSScxSm0TEaUNtdZWYDyBVVlupUlyw88LgEJAJsqlUkCPWGS/WEhsFSJxMQCyyJ0XxFI93n36YvevZ8b1z5+W+nXv/vy/Ozu7OnF3H/z37nOecI6oKIiIyVy7uARARUX8Y5EREhmOQExEZjkFORGQ4BjkRkeGWxPGi1113na5ZsyaOlyYiMtbJkyd/r6orWh+PJcjXrFmDycnJOF6aiMhYIvKK2+MsrRARGY5BTkRkOAY5EZHhGORERIZjkBMRGS6WrhUiojiUyhVMHD2L16s1rCzYGB8bwpaNxcA/J2oMciLKhFK5gt1Pn0atPgsAqFRr2P30aQDwDOZePicOLK0QUSZMHD27EMiOWn0W9x8+hbW7nsXovmMolSsdfc7E0bOhj7cbnJETUSa8Xq25Pj43fyWD22zb63O8Ho8Lg5yIMmFlwUbFJ4Br9VnsOTK9UBPPiWDW5fKdlQU7rGH2hKUVIjJGqVzB6L5jnqWQdsbHhjr6uGqtjkq1BgVcQ9y28h0/V1QY5ERkBGfh0QlZpxTSaZhv2VhEwbb6GkOxYGPv1vU9LXT280PID0srRJR4pXIFOw9PXTVDrtVn8fAz04uC1a1dEGgsXFZrdQiAXm4qdmbivYZ4mN0vEsfly8PDw8rTD4nIT6lcwcPPTOPCxXrbj9s+MojjZ2Zca+A5AfIiqM9dnXUigL0kh1p9DisLNi5euuz7WnkR3HXTajyyZX3HX8fovmOuYysWbPx4180dP4+InFTV4dbHOSMnolD1sqHmwdJpHHz+HDqdZx48cc5zlj2nwJzHE6kCCsH+bRuwZWPxqpmzm1lVPH7iHAB0HOZhd7+wRk5Eoemlrv1g6TQeP9F5iAO9lUoctfos7jt8CqVyBVs2FrF363rkRXw/75vPv9rxa3h1uQTV/RJYkItIXkTKIvLdoJ6TiMzmtaFm5+Ep10W/UrmyMNuNkiow/sTUQpjfddNq389x62jxMj42BNvKL3osyO6XwGrkInI/gGEAb1PV29t9LGvkRNmwdtezvrNlKye4ZtkSXLjY+0JkUJYPWBhYusS33xxo1Mp/ufejHT93EGe2hFojF5FVAG4D8FUA9wfxnERkvk424dTndGGBMc4QB4ALF+u+i52Okb9Y3tVzb9lYDO18lqBKK18D8CUAcwE9HxGlwPjYEPyrzWZ64dwbgfaC96PvIBeR2wH8TlVP+nzcDhGZFJHJmZmZfl+WiBLM2fxy36FTsc+yw5Kkw7OCKK2MAtgsIh8FsAzA20TkcVXd3vxBqnoAwAGgUSMP4HWJKGFK5Qr2HJlGtdZZecJ0STk8q+8ZuaruVtVVqroGwKcAHGsNcSJKP6fVMCshDgA5kUSUV9hHTkSBcGs1TLtZ1a7OewlLoDs7VfVHAH4U5HMSUfDCuL4sKWWGoHTaCunUyuO8MYhb9IkyptcDnBqf9zPU6o3mtJwAd980uLBN/VrbSk1ZpdhB22SzuH+IMciJMqbd9WXt7q68/9CpRf3Fcwo8fuJcLDsxw7bmHTZ++8b/dbx7M+6LJhjkRBnT7QFOzhGyWdok8pNfnu+4bTIJF01wsZMoY7o5wMkpw3RzrkgaeH21ywcsbB8ZRLFgQ9DfRRNB4oycKGPGx4auOqq1eVbZvBDqdWdlFgmA8lduiXsYrhjkRBnihHStPov8fEgXCzY2rVuBiaNnce+hU4s+niF+Rdx18HYY5EQZUSpXMP7E1MJNOU5IV6q1VC5YBikJdfB2GOREKeTWJ77nyLTrdWfU3vIBCw/dcX3sdfB2GOREKeGEd6VaW7SZpVKtpfrwqrANLF2S6BAHGOREqdC6yac1tBnivYt7s08n2H5IlAJZPOckKkle5HRwRk5kEK8zUkyYNZoo6YucDgY5kSHanZHSyZVq1J2CbWHP5mQvcjoY5EQJ0u5UQq8zUriQGYy8COZUAzsNMkoMcqKE8DuV0GvGzRDvn23lE7HVvlcMcqKE8Jpx7zw8hScmuWEnaCbPwFsxyIkSwmvBclYVP/7l+YhHk24C4NE7bzA6vJsxyIli1FwTp+jcMzKYmhAHGOREsXmwdBoHT5xjjTtio+96+8KtRmnBICeKQalcYYhHLC+Cu25anboQBxjkRKFzaymcOHqWIR6hvEiqauKtRGM4b3h4eFgnJycjf12iqLW2FAKNVjdup4+eoFEbN3lGLiInVXW49XGetUIUIq+WQoqeonFZdKlciXsogWOQE4WI3SjJ8/Az03EPIXCskRMFwGtrPc9ASZ4LF+txDyFwDHKiPrW2ETZvrd+0bgWvUaPQMciJ+uDVRlirz2LPkWm8UUvf7M90BduKewiBY5AT+fA7kdCr76vKEI9dDsBc09tWTrBn8/VxDSc0DHKiNvxOJORiZrLZS/MoDCx1/SGcJgxyypx2M+xWXu2DE0fPcjHTAG9emsX0P94c9zBCx/ZDyhRnhl2p1qC4MsN26y0ulSueIe3MxDetWxHmcIk6whk5ZYrfDBtoBPjDz0y3bVNbWbBRKlfwX8+zIyXJ0riw6abvGbmIrBaR4yLyoohMi8gXghgYURi8atrO486M3a/XuFKt4d5DpzDHA1MSK60Lm26CKK1cBrBTVd8NYATA50TkPQE8L1HgVhbsto+7zdjJDDkBlg9YEADFgo2JT6b3kKxWfZdWVPU3AH4z/99/EJEXARQB/Lzf5yYK2vjYkOshVpvWrcDovmNcuDSIlRNcs2wJqhfrqe5I6USgNXIRWQNgI4DnXd63A8AOABgcHAzyZYk61tz/7XStbFq3Ak+drHAmbpgszbj9BHaMrYhcA+C/AXxVVZ9u97E8xpaShDNxM72877a4hxA5r2NsA5mRi4gF4CkAB/1CnCgO7XrHuanHPFnpRulUEF0rAuDfAbyoqo/1PySiYPn1jnstgFIyZakbpVNBzMhHAXwawGkROTX/2JdV9XsBPDdR39r1jgPA+Tf/GMewqAfFjC9qegmia+V/0LhFiSgRmsso19qW5+FVlWoN409OoT7LZvCkWz5gofyVW+IeRmJxZyelSushV34nEDLEk8+28njoDpZS2mGQU2qUyhXsPDyF2RguFKdwiAB7t65nKcUHD82iVHBm4gzx9LBygv13bmCId4AzcjJSazvhxUuXuaEnJQTI/E7NbjHIyThulz1Qeuzfxll4t1haIePwYKt0c9pCqXMMcjIOd2KmG/9+u8fSCiWW17Z6Xq+Wbtxp2z3OyCmR2m2rHx8bgm3l4x4ihcC28hgfG4p7GMZhkFMi+W2rz3EvceoUCzZ7xnvE0golkled1LlijdKBZ6cEg0FOidBaDy8MWL73ZpIZrBwAkUXHIdhWnrPvADHIKXKtod16Q0+lWoOVE1h54Vkohts+MohHtqxvex489Y9BTpFy28xz8MQ5tMZ1fU5RsC2IgDNzgx0/MwOgccUegzs8DHKKlNsiptecu1qrw8pzVdNk7AmPBrtWKFLd/sNmacVs7AmPBoOcIlUY4F2LWcKe8GgwyClSPGU2O966NM+6eEQY5BQpvxt7KD2sPOMlKvxOU2ScW+spG97gD+3IsGuFQtPaO3yBt9VnChc6o8Mgp1Dw8ods4+FX0WJphQL3YOk07j10ipc/ZISVF2wfGUSxYEPAw6/iwBk5BcIpo3DmnS0F28KezdcztGPGIKe+lMoV7DkyzW6UlHt53208LyXBGOTUs9Y6OKUbz0tJLgY59aRUrmDn4SnMcodP6uWF590kHRc7qWvOTJwhng38e04+Bjl1ze0EQ0qv5TwfJ/EY5NQ1Hk2aLZyQJx9r5NQxp2uB/66zhVvtky+QGbmI3CoiZ0XkJRHZFcRzUrI8WDqN+w6dYp94BnGrffL1PSMXkTyAbwD4MIDXAPxURI6o6s/7fW6KH/vEs+OtS/OYUyxa/+BWezMEMSO/EcBLqvorVb0E4FsAPhbA81LMnO4Uhng2XLw0i71b13OrvYGCqJEXAbza9PZrAG5q/SAR2QFgBwAMDg4G8LIUJLdde+xOyZaVBZubfgwVxIzcbbfAVethqnpAVYdVdXjFihUBvCwFxZl5V6o1KBonFY4/OcV6eAqNvuvtGLDc/9lvWsd/l6YKIshfA7C66e1VAF4P4HkpIm4zb156nD4DVg4H/+5vsPytb3F9//EzMxGPiIISRGnlpwD+SkTWAqgA+BSAuwN4XupBLwcbsS88G2r1OQDef9/8/8BcfQe5ql4Wkc8DOAogD+A/VHW675FR19wuc9j99GkAWAhzt6BfWbBZRskAp43Q6++bbYbmCqSPXFW/p6p/rarvUtWvBvGc1D23EkmtPouJo2cBeNfCz/MKttSxcouXrprbCMfHhmBbec/3k3m4szNFvH41rlRrWLPrWdf31WeV9fAUmvjkDZ4lNudPni2eHgzyFGGJhJr9eNfNnu9jm2G68NCsFHH7lZmyySmnUTZwRp4irb8ys2CSXexAyRbOyFOktSOlYPMc6bTzuruHHSjZwhl5Sri1Hlp5QU6AOU7NU0vR6DjhQVfZxhl5SnjtzmSIp5tzsBUPuso2zshTgjXR7HFm3uxAIc7IU4I10fRamm9Uwpsvsy/YFmfetIBBnhLjY0Ow8l5LX2SqASuHf/nEDbCt/KK7M/94eS6+QVHiMMhTYvKV89yhmTK2lcc/b32f79ELRAzyFCiVKzh44lzcw6AANZdOeFoh+eFiZwo8/Mw0N/+khADYv23Doto3TyskPwxyw7Ru+lnzDhsXLvJOzTTIAXisJcSBxvpH8x4BgL3itBiD3CClcgXjT0yhPt8cXqnWeEhWShRsC3s2X+/ahcLTCskPg9wge45ML4Q4pYNt5TtqI2SvOLXDxU6DVGssoaRJXoS94BQIzsgTrLUeTuYSYNGCdKczcaJOcEaeUG7XspG59m/bwPNQKDSckSeU2yYQMpOz37bdjT1E/eCMPIFK5Qpn4Cmi4I09FC7OyEPWWuf2axtzSiqULtyFSWFikIfI7bIHJ6S9wpwllXTiYjWFiaWVEPVy2BFLKmYr2NZVp1ByFyaFjTPyELU77Mhtq/1PfnU+4hFSkJYPWCh/5Zauy2lE/WKQh8jrsKNrbeuqkgtn4uZ76I7rAXAXJkWPQR4CZ0ZWqdZcN4KIgHXwFBEA94wMMrwpNqyRB6x5Iw/QCHGnYupsBKnytELjNW/u2b9tAx7Zsj7uIVGGcUYeMLcFTkXjH7yzIcSZrZOZmv8uiZKAQR4wrwXOSrWG0X3H8Hq1hsKAFfGoKCjsQKEkYmklYF79wgIsnJvCiyDMxNMKKan6CnIRmRCRMyLyMxH5togUghqYqcbHhmBb+UWPtS54knlsK49H77yBIU6J1O+M/DkA71XV9wH4BYDd/Q/JbFs2FrF36/pFi2EMcbPxtEJKur5q5Kr6g6Y3TwD4RH/DSYfWPuLRfce4uGkAnhlOpgqyRv5ZAN/3eqeI7BCRSRGZnJmZCfBlk61UruDNP16Oexjkw7byuGdkkGeGk5F8Z+Qi8kMA73R51wOq+p35j3kAwGUAB72eR1UPADgAAMPDw6mrNrhtywZw1e3nlAw5Af7sWpvb6CkVfINcVT/U7v0i8hkAtwP4oKqmLqA74XbK4fiTU6jPZvLbYYS7bxrkJh5Kjb5q5CJyK4B/APABVb0YzJDM47YJiCGeXNtHGOKULv1uCPo6gLcAeE5EAOCEqv5936MyDC8NMIOVE0x8ki2ElD79dq38ZVADMZnXKYcUP6cTpcg6OKUYt+gHYHxsiIuaCZQX4SYeygQGeQCcoLj30KmYR0LN5lQZ4pQJDPIOdXLriwiQzb6dZOI9mZQVDPIO+F2iXCpXMP7kFEM8QXhKIWWJxNH6PTw8rJOTk5G/bq+8ttjnRTCnipwIZpniicGFTUorETmpqsOtj3NG3gGv9kInvBniycALHyireB55B661eRFEnGzL/39TllIoyxjkPkrlCt68xEOv4rJ9ZBAv/tNHUGjzw5QHXFHWsbTiY+LoWW63j0E+J3i0aRfmGzX3W5UEYDmFMo8zch/cfh+9YsFeFOKAdyshWwyJOCO/itMvXqnWGn3hcQ8oY7wWLN12z7IuTtTAIG/S2i/OZpRotQtmZ3butymLKIsY5E3cjqOlaAjgu2DZeoUeETWwRj6vVK7wBMOYWDnB/m0bGNJEPWKQ40pJhcIhLX82Wz5g8Yxwoj4xyMGSSpiWD1jYv20DigXbdeF4YOkShjhRn1gjB1sMw9B6Fvh9Hkf88ntP1L/Mz8hL5Qpy4vZLP/Wj9Sxw9oEThSfTQe7UxnnoVfeKPgHcGtDjY0Owrfyix9gHThSMTJZWmjf9UG+cHu7W3nvAPaDZB04UnswFuVvwUPeaL9YAOgto9oEThSNzQc4OlWDU6rOYOHp2IZwZ0ETxyVyNnF0SweH3kigZMhfk7JLwZ+UEb1ni/78Gv5dEyZC5IHfrnmDz4RXFgo1tN672bclkxwlRcmQuyLdsLOLj7y8iPx9UOWGQO5wjZI+fmXFdR8iLQMAbeYiSJnOLnaVyBU+drCz0js+xhRxA44eZM8P2qn3PqeLX+26LcFRE1InMzcj3HJlm10oLAXDPyODCDJu7MInMkqkgL5UrqHrc/ZgVThnJKS0VCzb2b9uAR7asX/gY7sIkMkvqSivOrk23zSkTR8/GPLr4KbyvU3NwFyaRWVIV5K27NivV2qIdiOx7bujk+8BNPkTmCKS0IiJfFBEVkeuCeL5eue3adHYgAqzxOvh9IEqXvoNcRFYD+DCAc/0Px1upXMHovmNYu+tZjO47hlK5ctXHeM00ncfdar8msPKC7SODC3XtfrDWTZQ+QczI9wP4EuB6AUwgnJJJpVqD4krJpDXM/botnB5y0/rG67OK42dm8OidN/S0malYsNn/TZRifdXIRWQzgIqqTonPbFFEdgDYAQCDg4NdvU67kklzKI2PDfkeqXr8zEx4P3FC9Hq11nYRcnTfMddjef0WNonIfL5BLiI/BPBOl3c9AODLAG7p5IVU9QCAAwAwPDzcVZb6lUwcnXRbJHnBc/mAhYGlS1wDufm3CrcZdSc/xIgonXyDXFU/5Pa4iKwHsBaAMxtfBeAFEblRVX8b5CBXFuy24dbMr9vC67nciAC9Xh5UsC2IABcudta3LgAeuuN6AOgpkNkySJRdPZdWVPU0gD913haRlwEMq+rvAxjXIkHONt2ey00OwGN3bmh8zhNTqHexl3/7yCCeOlnpagdp885KoLdAZssgUTYZ0Ufe72yzdZPQx99fxPEzM21n5vm8YPKV8/jm8692fadnL59z/MwM1u56duFrY12biDolGsPFw8PDwzo5ORnJa3ndKbl363rcd+hUIhY+8znBbNOM3xkfZ9dE1ExETqrqcOvjqT9rpV3HS1gbY/zavaXpTwEWhTiweBMTEZGf1Ad5u46XTetWhPKa7XK8WLDx63234WvbNmCZlff8jSDJ3TVElCypD3KvWbeiUcsOw5wCA1buqkBvXqD1uwSa2+iJqFOpD/J22/K7XZDsRq0+h/3bNnjuqmw342b/NxF1w4iulX40d7x02j/eTNCYvefkym1CBdsCgLZnm68s2G3bAb362fMiXOgkoq6kPsiBK/3Va3c969ulYuUE1yxbgurFets2x1K54tlfbuXEd0bt1RvPECeibmUiyB3tZsFzql1vvnn4mWnXnZvXLFvi+xzciUlEQclUkAc9C656bL/3erwVd2ISURBSv9jZbMvGIvZuXR/Ysa68pJiIkiBTM3Ig2FkwTxwkoiTIXJAHiXVuIkoCBnmfWOcmorhlqkZORJRGDHIiIsMxyImIDMcgJyIyHIOciMhwRnattF7dxpY/Isoy44K89eq2SrWG3U+fBgCGORFlknGllXZXtxERZZFxQd7u6jYioiwyLsh5UBUR0WLGBbnb1W08qIqIssy4xU4eVEVEtJhxQQ7woCoiombGlVaIiGgxBjkRkeEY5EREhmOQExEZjkFORGQ4UdXoX1RkBsArkb9w964D8Pu4BxEDft3ZktWvGzDva/9zVV3R+mAsQW4KEZlU1eG4xxE1ft3ZktWvG0jP187SChGR4RjkRESGY5C3dyDuAcSEX3e2ZPXrBlLytbNGTkRkOM7IiYgMxyAnIjIcg7xDIvJFEVERuS7usURBRCZE5IyI/ExEvi0ihbjHFCYRuVVEzorISyKyK+7xREFEVovIcRF5UUSmReQLcY8pSiKSF5GyiHw37rH0i0HeARFZDeDDAM7FPZYIPQfgvar6PgC/ALA75vGERkTyAL4B4CMA3gPgLhF5T7yjisRlADtV9d0ARgB8LiNft+MLAF6MexBBYJB3Zj+ALwHIzMqwqv5AVS/Pv3kCwKo4xxOyGwG8pKq/UtVLAL4F4GMxjyl0qvobVX1h/r//gEaoZeKgfxFZBeA2AP8W91iCwCD3ISKbAVRUdSruscToswC+H/cgQlQE8GrT268hI4HmEJE1ADYCeD7ekUTma2hMzubiHkgQjLwhKGgi8kMA73R51wMAvgzglmhHFI12X7eqfmf+Yx5A41fwg1GOLWLi8lhmfvsSkWsAPAXgXlX937jHEzYRuR3A71T1pIj8bdzjCQKDHICqfsjtcRFZD2AtgCkRARrlhRdE5EZV/W2EQwyF19ftEJHPALgdwAc13RsOXgOwuuntVQBej2kskRIRC40QP6iqT8c9noiMAtgsIh8FsAzA20TkcVXdHvO4esYNQV0QkZcBDKuqSael9UREbgXwGIAPqOpM3OMJk4gsQWNB94MAKgB+CuBuVZ2OdWAhk8bs5D8BnFfVe+MeTxzmZ+RfVNXb4x5LP1gjJy9fB/AnAJ4TkVMi8q9xDygs84u6nwdwFI0Fv8NpD/F5owA+DeDm+b/jU/OzVDIMZ+RERIbjjJyIyHAMciIiwzHIiYgMxyAnIjIcg5yIyHAMciIiwzHIiYgM9/9cTfVwRL0ogAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def rmvn_svd(n, mu, Sigma):\n",
    "    \n",
    "    k = mu.shape[0]\n",
    "    \n",
    "    U, d, V = np.linalg.svd(Sigma)\n",
    "\n",
    "    D = np.diag(d)\n",
    "    D = np.sqrt(D)\n",
    "\n",
    "    Q = U @ D @ V.T\n",
    "    Z = randn(n, k)\n",
    "    X = Z @ Q\n",
    "    \n",
    "    return(X)\n",
    "\n",
    "Sigma = np.array([[1,0.9],[0.9,1]])\n",
    "mu = np.zeros(2)\n",
    "\n",
    "X = rmvn_svd(1000000, mu, Sigma)\n",
    "plt.scatter(X[:,0], X[:,1], marker='o')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Faktoryzacja - Metoda Choleskiego"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Metoda Choleskiego mówi, że dla każdej dodatnio określonej i symetrycznej macierzy rzeczywistej $\\Sigma$ można znaleźć jej rozkład na dwie macierze: \n",
    "$\\Sigma =Q^{T}Q$, gdzie\n",
    "* $Q$ - macierz górotrójkątna."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Przykład 6"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Wylosujemy próbę $1000$ elementową z dwuwymiarowego rozkładu normalnego o średniej $\\mu= (0,0)$ i macierzy kowariancji\n",
    "$\\Sigma=\\left[\\begin{array}{cc}1.0 & \\ \\ 0.9 \\\\ 0.9 & \\ \\ 1\\end{array}\\right]$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1.         0.        ]\n",
      " [0.9        0.43588989]]\n"
     ]
    }
   ],
   "source": [
    "Sigma = np.array([[1,0.9],[0.9,1]])\n",
    "chol = np.linalg.cholesky(Sigma)\n",
    "print(chol)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Funkcja"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0xcd58af0>"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO2df5Ac5Xnnv98djWCkUIwIa0CDhDgHRNDJ0sIekktVPgsbBNjAIuAENhdfclU6J/AHmFNFxBQIHw667CXgGJdt2UfZLggWtmAsgmyBI65sKxFml10hBIhgAkIjKqwNKxs0htHquT9mepmd7e7pmX779/Op2tLOdGv63d2nn37e5ydFBIqiKEr66Yl6AYqiKEo4qMJXFEXJCKrwFUVRMoIqfEVRlIygCl9RFCUjzIh6AW6ceOKJsmDBgqiXoaSU4eHhX4tIb9jXVblWgsRNrmOt8BcsWIChoaGol6GkFJKvRXFdlWslSNzkWl06iqIoGUEVvqIoSkZQha8oipIRVOEriqJkBFX4iqIoGSHWWTqKUh6pYHD7Phwcr2JusYB1qxZioK/U9piixIU4yalvhU9yHoDvATgZwFEAm0TkKy3nEMBXAFwC4DCA/yYiz/i9thJfTAh5eaSCWx7eg2ptAgBQGa/ilof3TB53OqZKX/GDSQXtJsNRyKkJC/8IgJtF5BmSxwEYJvmEiDzfdM7FAM5ofC0D8PXGv0oKMSXkg9v3TX6GRbU2gcHt+ya/tztm6kZSYyZ7mFbQbjLczef5fRj59uGLyBuWgIvI7wC8AKB1BZcD+J7U2QWgSPIUv9dW4kk7Re2Vg+NVx/fdjhnEMmb+GMByANeTPLvlnGZjZi3qxoySUEzJroVJObUeRpXxKgQfPIzKIxXPn2E0aEtyAYA+AE+1HCoBeL3p9QFMfyhYn7GW5BDJobGxMZPLU0LClJDPLRYc33c7Zgo1ZtJNeaSCFRt34PT1j2HFxh0oj1QcZbQyXp08pxM6lVO7NVmYeBgZU/gk/wDAFgA3ishvWw/b/BfbUVsisklE+kWkv7c39DYnigFMKeN1qxaikM9Nea+Qz2HdqoWux4LArzGjhky8cLKWjy/kHf9PNxZ1J3LazoI3YUgZUfgk86gr+wdE5GGbUw4AmNf0+lQAB01cW4kfppTxQF8Jd61ejFKxAAIoFQu4a/ViDPSVXI+ZxoQxo4ZMvHCylklMk93WczqxqDuR03YWvAlDykSWDgH8XwAviMjfOZy2FcANJL+PerD2kIi84ffaSjyxhNlEpoOl3Ds9Zgo1ZtKJk1U8friGu9csxeD2fagYck16ldN2Fvy6VQunBJSBzg0pE1k6KwD8VwB7SI423vsrAPMBQES+AWAb6lkML6OeyfCnBq6rxJgwlHHQqDGTXuYWC7YKvYf1DdvO9edjxcYdtueYjBN5WZN1PROGlG+FLyK/gP22tvkcAXC932spSsioMZNS7KxlAJgQmUzDXLdqIdb9cDdqE1M9dO++dwTlkYpxg8aLBe/XkNJKW0VxQI2Z9GIpzZsf2o0JmarQq7UJ3Lh5FMVCHhNHp+eWjFdrgRRPmXSFOqEKX1GUTDLQV8JNm0cdj49Xa47HqrUJ3PzQ7snXppR00K5QVfiKoqQat+rU4wt5V8XuxoQI1v1gN0BMun0q41XcuHkUdzy6F7dfuih2cSxV+IqipJZ2/Zh+994RX59fs3H5AMDbh53dPlE2U1OFryhKYmnXTdXJR2/lttv56E1RrU3gCw+N4saG26hYyOPTS07BluFKZM3UVOEriaY8UsGGrXsnt+VzZuVjuZVWzOOlm2qrsrc42KhmDZrm58l4tYb7d+2fdo7ppn9uqMJXAsVUm2S7zyiPVLDuB7unbKvfPlzDuh/Wg2mq9NONU2XqjZtHkSMdlT0AHJvvQbV2NOglesZw0z9HVOErgWGi1eyt5T14YNf+SWus+TMGt++z9aHWJiQ0i0kxQzeGgZuSdFP2PUCslD0QXDFXK6rwFSPY3bDd9AJv/hynDArrM9xu+LAsJsU/3RoGTpWpbhD1wQZRQkxtthRk079WdKat4hunLn+d9iJp/Ry3dDnrweJEWBaT4p9u2/52qiTzPa41dKHRrOznzMoH1vTPDlX4im+cbtgc7W8wJ2Vs9zlOWLsIu5s4n2NoFpPin27b/g70ldCJDj8Kh57sEXKoyxqAblGXjuKL8kjF0ZKfEEEhn3PtDVIeqeCOR/fi7cOdCf7Ks3onrSLN0kk27ZqGudFJVmWQKZjdclQQarxJFb7SNZYLxokciSvPLeHJF8cc86TtmlN5YctwBf2nnZCKrpxZx0/b31IXfvy4EWa8SRW+0jXtXDATItgyXHH0Ud7x6N6ulD0wtZeJKvxk003TsLqx8Wzssm26Icx4kyp8pWu8WCbV2gQ2bN077WYG0LEbp5XmVraq9JNNJzu18kgFX9g8Gnm2jQl62Hnw2Q+q8JWu8ZoWN16tTfrYK+NV3LR5FDMMpQuEWaWohI+VplsZr04WU+XIVCj7Wfke/PXqj4Qqu6Zm2t5H8k2Szzkc/zjJQyRHG1+3mbiuEjzlkQpWbNyB09c/hhUbd0wZ4Gw3u9YLAsDkTjwoH6jKdbQ0p+kCHxRTuRVVJYVSsYDn/9fFoRsqpiz87wC4F8D3XM75uYh82tD1lBCwK4i5afPolNL1YiGPY/M9vt0zfgjQB/odqFxHRidpukmjMl6dnJoVZvdMIwpfRH5GcoGJz1Lig90NZ9lWlpU1Xq2hkM+h6KOvuB+CrFJUuY6WtFdL37h5FHdtex6//f1EaN0zwyy8+ijJ3SR/THKR00kk15IcIjk0NjYW4vKUVrzecNXaRCjKvlQs4Lrl81EqFsDG6zCrFB1QuTaAneswC9XS//6797uqMu6WsIK2zwA4TUTeIXkJgDKAM+xOFJFNADYBQH9/f/KddQmgdUu58qxePPniWOyqEgGg/7QTcOfA4qiXYaFybQCnXjrnzD8+tDbGcSOo3U0oFr6I/FZE3ml8vw1AnuSJYVxbcceuD879u/bHspjFUgTNgeMoUbk2g1Nrjp2/eiuTyh4ILi4VisIneTJZb6xC8rzGdX8TxrUVZ6yJQEkKjAW53e0UlWszpN1X3w3vvnckEMPGiEuH5IMAPg7gRJIHANwOIA8AIvINAFcB+HOSRwBUAVwjkoLcqgRjWfZJTHELS0GoXIdDN22O08bsmTm8+/4Hhtd41Xkmrh9MZelc2+b4vaintykxIckpb83bXbvCnJKh1DaV63BYeVav7ei/LPF7m8KUIIoKtT1yRknyNnrBH9YV/q3lPbhp8+i0wpy4+foVdx4ePhD1EiLHaaddGa9OK3j0g7ZWyAB2hR1J3kb/86/emjb6sBVtuRBfWqeaHU5BAzS/uM3gNZmbrwo/pTS7OppHqlnCc+W5JWx++vWuu1VGiQB48KnX22ZwJHkXk1ZaH9RRFOvFjUI+hyvPLWHLcMXRzWrKgFGXTgpp7UHSqhirtQk89uwb8Rv/0wFegs1ZKNxJEuWRiuuuLIv0sH4/PvniGK48t4RiIe94rgkDRi38FOIlIBtl75swCHMwtOKNwe37VNm3YA3hqoxXsfnp112NMBMGjCr8FNDqo0+qb94UprJ0lO6xixupi80dN/eqKQNGFX7CsStLzyr5HmLw6iWq6CPGqVXC8RE12EsDpnpGqQ8/4SQ5n940taMSmyrcLOPUKoFEV/MTsk6pWDBmxKjCTzi6TZ6K/j6ix+lvMH64hivP1d2XG/keTnltOhalCj/haCbKVPT3ET1Of4O5xQIefCrbFbVuzJmVx+DVSwJt/60+/ISzbtXCKf7SLEPUy/SVaFl5Vq9t+mWW40vtKORzuP3SRR0Nc+8GVfgJxxKO1n72WexNIgC2DNdL0J98cSyUkXHKVMojFWwZrmj6ZQeEmVWmCj8FDPSVMPTaW3jwqddRGa/igQwqe4tqbWKKdRn0yDhlalW30jkHx6uTyQZBy6gq/ITgNuj41vKeKRZ91q0ru8pi7asTDK0pmErnWIOHwjBMVOEnALu85ps2j+LGzaMoFQs4eCibllWpgyIzzd4JBk0LNkcYholm6SQAu5uq2WWRxZEbpUaswivdZO+QvI/kmySfczhOkn9P8mWSz5I8p+OLJBx9kJol6N+nKvwEoDfVVIj6CDivgWkfuczfAXCRy/GLUR9afgaAtQC+3s1FkoymwZol6N+nEYWvllCw6E01FYH3trrFQr7rXGYR+RmAt1xOuRzA96TOLgBFkqd0fKEEs27VQq2eNUQYDf9MWfjfgVpCgaFdH7tn9jEzgvSJlgC83vT6QOO9aZBcS3KI5NDY2FhQ6wmdgb4S7lq9GCU1SnxBmuuX44YRha+WUDCURypYsXEHbto8CrY/XbEhYHeY3Z/FNqIiIptEpF9E+nt701Ec1iyfAHDPmqURryjBSDhpw2H58DNvCXVK8xATgaZadkvA7rADAOY1vT4VwMEgLxgXWuWzMl7FjQ3Fr3QOiVBmMIel8DNtCXWDprt1xnXL50/zJYfgE90K4E8aMarlAA6JyBtBXjAuqHya5agAtzy8J3ClH5bCz6wl1C2amdMZ/aedMOlLNtV4iuSDAP4FwEKSB0j+d5KfJ/n5xinbALwC4GUA3wLwFz5/jMSgVbXmsfLwgySswqutAG4g+X0Ay5AhS6hbdHJVZwxu34ed68836gcVkWvbHBcA1xu7YEIIw/WQVRKRh6+WkHk0M6czdEcUHn/18LNRLyG1BJ2CbcTCV0vIHNqIqju0ViEcyiMVHK4djXoZiWdWvmfa7zGMPHztpRMDmpU8oRk5nRLGjaLU5VQzccxQrR3Fig+fgFd/Uw21jbcq/IhpbYymyr5zjpmhHUKCpjxSwRdU2RtDAPzzr97C3WuWhtrFVe+UiNH0Nv+MV2tY94PdGkwMkMHt+6COHLMIEHhWTiuq8CNGg41mqB0VbNi6N+plpBaV02AI+/eqCj9iNNhoDq8N1ZTOOb6Qj3oJqSTs36v68EOmdXLVyrN6sWW4om4dJXTcpqi1nve7945EsML0w5CbZKnCDxG7yVVbhiu48txSJoeOm2bOrA+sJa/KLKvYyaLTiL3B7fswcVTTCYJg/HC4u1J16YSIXYDWGrqt+COfI26/dBEA+8ZeYfQpSRJOsmgXRFT/fXCE7dJVhR8iTjeO2k7dMSvfM9k3Z/CqJZOWaSfKLKs4yaLd+xpnCoYo6kfUpRMi2h/HLALa5jF3osyyipMszi0WUB6pYMPWvZNB8NkzdaKVaUoRuRlV4QdMsy9ZMx3MYlntrTeNmzJT6qxbtXCKDx+oW5wrz+rFuh/sRq3JZ//u+5pQYJJXN34qsmurSydAWn3JmjZoHjur3W7OqrZfmErzaEKiPvv32HwP7t+1f4qyV8xSjNjoU4UfIFpFGzx2VnurMjPRGz+NDPSVsHP9+bh7zVK8d+Qo3g45YySLjFdrWLFxR2QJBOrSCRD1GZulB5hS3u9mtQ/0lYwoeJIXAfgKgByAb4vIxpbjHwfwIwD/1njrYRH5ku8Lh4gaJuHilgIbNKrwA0SDtGYpzMyhOGtmaLn1JHMAvgbgAtSntj1NcquIPN9y6s9F5NOBLcQwrTUKKqPh4xR/ChpV+AGybtVC3LR5VNMuDfHu+xPY+6Xzw7zkeQBeFpFXAKAxse1yAK0KPzHYFVwp0RCFB0B9+AEy0FdSZW+YkH2fJQCvN70+0HivlY+S3E3yxyQX2X0QybUkh0gOjY2NBbFWT6j7Jj5EkTVmasThRST3kXyZ5Hqb4x8neYjkaOPrNhPXjTvlkQpyYTfLSDkhV8za/fFan+HPADhNRJYA+CqAst0HicgmEekXkf7e3l7Dy/SOxpXiQVRZY74VfpOf82IAZwO4luTZNqf+XESWNr4SFdTqBmvrPCFq45sk5IrZAwDmNb0+FcDB5hNE5Lci8k7j+20A8iRPDGuBnaK1CNETZdaYCQt/0s8pIu8DsPycmUa3zsERopX6NIAzSJ5OciaAawBsbT6B5MlkfRtH8jzU76nfhLXATrGrUVDCo1QsYOf68yNLETYRtLXzcy6zOe+jJHejbiH9TxGxnVZBci2AtQAwf/58A8szi9cujBoMC46wrFQROULyBgDbUU/LvE9E9pL8fOP4NwBcBeDPSR4BUAVwjUh8t3WWrOps2vCJQ/GfCYXfiZ/zHZKXoO7nPMPuw0RkE4BNANDf3x+rG8cuw+GmzaMYeu0t9J92wpQHAQnE97ZPNmHeNA03zbaW977R9P29AO4NbUEGGOgr4eaHdqu7MWTiUPxnwqWTOj+nE3ZuGgFw/679WPfD3VPa8eq9FBxDr70V9RISz7XL5rU/STFGqViIXNkDZhR+6vycTrj5jmsTquHD4sGnXm9/kuLKnQOLMTOnGWRhEAdXjoVvl04a/ZxOaFViPFBXhBnUSAmeqNogO2Gk0jaNfk47tHI2Hmhtgxlmzcxp6+MAsTJy4oRW2nbAQF8Jn10+3zZKrYTH8v8wJ+olJJ5by3tU2QdMHIvcVOF3yJ0Di3H3mqUoaQFLZDyz/5DOp/VAeaSCFRt34PT1j01ryfvAUzpHOWiKs+I38EgVfhdYfcTvWbNUi1giQOfTtqfdIHcNgwRPHH/H2i2zC1rHFvZQx8CFTRy3y3HCbZB7XAKIaedQDCfcqcLvkNbiKx1bGA3aE8Ydt0Hu6g4LhzjKqCr8NrRa84d+X4vlVi1LxCmvOa64DXJXd1jwxFVG1Yfvgt0QclX24dGcDWVlYup8Wm84DXJfeVav1pIETLGQj62MqoXvgna8jJa71yyN5U2TBKzf2+D2faiMV5EjUa1N4P5dmp0TNO8dOdr+pIhQC98FDQxGx+yZOVX2PhnoK01a+lqdHB5xziJThe9CHIMuWSGfU9E0ge5SoyGuxqLeVS7osIjoiGNKWxKJq+JJO3E1FlXhuzDQV8JdqxejWIhfxVza6SFtK0TDxsO8ZpL8+8bxZ0meE8U6nYir4kkzcc3QAVTht2Wgr4TZx2hsO2wmRGwrRMPE47zmi1Ef5nMG6pPavh7qIl0oj1Tw7ntHol5G5ohrhg6gCt8Tui2OlgiDYF7mNV8O4HtSZxeAIslTwl5oK1ZKsRYGhktcBp04oQrfA7otDge3tscRPXTt5jW33s1ezgkdDdaGT5xdORaq8D2wbtVC5Hq0KXLQTIg4tp6O6KHrZV6zl3NAci3JIZJDY2NjRhbnhu5Kw4UAjpnRg5s2j0Yed3LDiMJPemCrHQN9JR0HFxKC6Ro0Qsup7bxmj+dARDaJSL+I9Pf29hpfaCu6Kw2XGT2sV+Ij2rhTO3wr/KQHtrxSrcW3ei5tCOq+UCLyVgpt5zU3Xv9Jw6hZDuCQiLwR9kItrB742j4hPAigdnTqpi6uxVcm0k8mA1sAQNIKbD3fdM5kYAvALpJFkqdEeWPY0dwobW7TLMpby3uiXlqmyJGxGA3ncV7zNgCXAHgZwGEAfxrkmpxk1DrW3MlVCQenGuY4utVMKHy7oNUyD+eUAMRG4bfeLNa27AdD+7HzV29FvLpsEac2AB7mNQuA68NYi5OMAnW344ate1XZx4g4utVMKHxjgS2gHtxC3e2D+fPn+1uZR8ojFdz80O5piqZam1BlHwE6PtIet6EmgM5miIpiIY/3jhyd8reJa8aOiaCtscAWEH5wy7Ka4mRVZpm43ihxwMlFUBmv4o5H94a8GgWoy+uGyxbhrtWL4xJ3csWEhT8Z2AJQQT2w9ZmWc7YCuKHh31+GiANbzWi+cnwggCvPLcXyRokDTkNNAODtw2rdh8FJx83EjFzONoaSBLn1beGLyBEAVmDrBQAPWYEtK7iFug/0FdQDW98C8Bd+r2uKOAZWsooAePLF4HPUk4o284ueCxadjHWrFmJusYCD41UMbt8Xy/RLJ4w0iYlTYKtT3KwmJXz0AeyMZUHeuHk04pVkl/t37cfmX74+mYbZGjiPO5mvtLWzmrTEKjrimNkQJwb6ShrUjpik5NzbkXmFb7VAtm4iwjmvVgkWDdh6Q1078SMpO9PMK3ygrvR3rj8fpWJBlX1ExDmzIW60GilK9CRlZ6qN3huURyrqy4+AQj6nir4LBvrq2UwL1j8W9VIyRT5HQKa6dZK0M1ULHx/k4ivBc93y+YnIV04KOo0tPObMymPwqiUYvHpJYmVYLXxoLn6Y3DmwOOolpIbySAXvvq8TrcLguuXzp8huUhR8K6rwkZyAS9JRn7MZrAZq6oIMj7TUh6hLB8kJuCSdpPg544zlflRlHy5pMQpV4QNY8Ieq8INmzqx8YrfBcULdj9GQFqMw8wr/1vIe7YgZMIV8DrdfuijqZaSCtFiaSSMtu9PMK/wHn3q9/UmKL7QhmjnSYmkmiWIhPbvTzCt8bYscPFuGK4lqMAUAJE8g+QTJf238O8fhvFdJ7iE5SnIo6HVplW24EMCGy9KzO828wqc2zgmcJPUaaWI9gH8SkTMA/FPjtRMrRWSpiPQHvajmKlsV3eARJDcF045MK/zySEVvmpBIoO/5cgDfbXz/XQADEa5lClYrkM8uD2ciXJZJWypxphX+4PZ9OKoenVA4PnkVoSdZQ3oa/37I4TwB8DjJ4cZ4TltIriU5RHJobMxMTvcDT+038jmKPUlqmeCVTBdeJdDqTCwxdZ2dSfI5m/e/2MFnrBCRgyQ/BOAJki+KyM9aTxKRTQA2AUB/f78RM0PDT+YpFvI4VK1Nm2aVFjKt8HX4SXi8fbiGFRt3xO0mesnJ707y30meIiJvkDwFwJt254nIwca/b5J8BMB5AKYpfFPcWt6DB3bt166uATF6+4VRLyFQfLl04prJ4JWVZwU/JF35AGs6UEIydrYC+Fzj+88B+FHrCSRnkzzO+h7AhQDsdgxG+Oy3/gX3q7IPjLT56+3w68OPZSaDV9LSHyNJJChjZyOAC0j+K4ALGq9Bci5Ja5znSQB+QXI3gF8CeExEfhLEYsojFS0QDJA0+uvt8OvSuRzAxxvffxfA/wPwlz4/MxS0/32wuE0OS0LsRER+A+ATNu8fBHBJ4/tXACwJYz0JeUgmkjmz8rj90kVxcjUGhl8L32gmAxBMNkMr2v8+WAjgs42+93ZotWjnJOEhmVTGD9cw9Fo2dk9tFT7Jn5J8zubr8g6us0JEzgFwMYDrSX7M6UQR2SQi/SLS39sbjI99w9a92oAqQAR1d5ldVWhWts6m0YdkcAiAB3btT0psyRdtXToi8kmnY3HNZHCjPFLBeLUWxaUzxcHx6uQWeXD7Phwcr6Y21S0oyiMV3PHoXrx9WOU1aAR1OU27bPr14VuZDBvhkskAoEdEfteUyfAln9ftGvWFhoNlkVqzV5XOKI9UsO6Hu1Gb0JycsMiC28yvDz9WmQxeyMIfNWxaa6rUbeOfwe37VNmHTBbcZr4s/LhlMnhBi63MUsjncOW5JTz54pi6bQyihkm4ZMVIyVQvnVvLe/DGIb2RTDFnVh53rV6MOwcWY+f683H3mqUAgJs2j2LFxh2ZCIIFRRaszaixdqalYgF3rV6cCSMlM60Vbi3vwf27tNmUSUZu+6AM3Up1tbKfrKpaIF3tZcNi3aqF6sMPGEFd2e9cf37USwmNzFj4OtnKLMWW7pd2s1artQnc/NButfS7YKCvhMGrlmD2TB12EiRZc51lRuHrZCtz5Hs4bQqQ040zIZKk/jmxYqCvhL1fugj3NFxlinmy5jrLjMLPxbQ/b9IoFQsYvHrJNDeN242ToP45sURdYsGQlUBtM5nx4V+7bJ768H1wz5qlropn3aqFU3z4rWRt62wS3R2Zp5TRbLLMKPw7BxYDqPvy1b3TGfme9lamdfzmh3bb/n6ztnX2S3mkMlmh3KO7U6Pkc8yksgcy5NIB6kr/V3ddglc3fgrX6TxQz9SOejtvoK+Ev/0vS7R/jk+sjKfKeBUCjT+ZpjYhmXUxZkrhN9N/2glRLyGVDPSVcNfqxSgVCyCyleNsCruMJ8UsWXUxZsal00pWn/BhkIb+OSSvBrABwB8DOE9EbCe1kbwIwFcA5AB8W0Q2+r12VpVRmGTVxZhahd/sA7Ur99f2Ct7JaC74cwBWA/im0wkkcwC+hnofqQMAnia5VUSe93Nhbf8RLJYPP4uk0qXT6gO1m6WqaZreyPUQX75icdTLCB0ReUFE2m0DzwPwsoi8IiLvA/g+6lPgfLFu1ULke1Q+TdH8m5wzK4/Bq6anFWeFVFr4TlWfdzy6d/K4BsKc6SEgAm2E1p4SgOYS7gMAltmd2Jj0thYA5s/3kDCg+t4IhXxOY0hNpELht7pvnLbDbx+uaX+SNuRzzJIFdCbJ52ze/6KITJvtYIOdWrYVLhHZBGATAPT397sKoLZG7h4COL6Qx6FqTQ0WGxKv8O2adrmhN5Izs/I9+OvVH8nSDfKSiPT7+P8HAMxren0qgIP+lqRBWz8cX8hj9PYL25+YURKv8DWFzT9zZuVx+6WLsqToTfE0gDNIng6gAuAaAJ/x+6EatO2eQzq+1JXEB23VGvLHdcvnY+S2C1XZt0DyCpIHAHwUwGMktzfen5zmJiJHANwAYDuAFwA8JCJ7/V7bbvi74o2splt6xZfCJ3k1yb0kj5J03BqTvIjkPpIvk1zv55qtdPIHzvcQ+ZxGw5rZMlzRXi02iMgjInKqiBwjIieJyKrG+wdF5JKm87aJyJki8mER+bKJa1vFa5pJ1hla0d0evxa+lav8M6cTmnKVLwZwNoBrSZ7t87qTeLWGioU8Bq9egsGrYjNtMRZoJ8v4YSUhaCaZd4qFvGbjeMDvTNsXAIDulshkrnLjXCtX2VdxioX1B27O0ll5Vq/rjNUbN4+auHRqUH9xfGhNQlC8MfuYGarsPRBG0NZzrjLQRb4yOi/lL2lQbArqOogPmoTQHRrL80ZbhU/ypwBOtjlkPFcZ6Cxf2Suteforz+rV3vhNqOsgeiwZVUOkOzRY6422Cl9EPunzGoHkKnvFLk9/y7AGKZsp6c0SKerG6QxiqsWowVrvhJGWOZmrTHIm6rnKW0O4LgDgjkf32rZZUOrozRI96sbxTqlYwN1rlmr77S7x5cMneQWArwLoRd5h4sEAAAyBSURBVD1XeVREVpGci3qr2EtE5AhJK1c5B+A+E7nKXiiPVPD2YS3EcOPYfOJLMRKP+p+9YyVgqILvDr9ZOo8AeMTm/YMApuQqA9jm51rdoOmGdV7d+KkpPuLmLfHbh2u45eE9AHRYdlRoZa13VEb9kWrzTi2negZOeaSCgb4Sdq4/H6ViYVrEXHPxo0Ura72R0bkMRkm1wtfIfT0Dp3kWgNNDUB+O0WFV1iruHH5f4xx+SbXCV8upTrMF7/QQ1IdjtAz0lTRbqg0qo/5JtcJvHaid5QIjy4K3ewhqpk48UAPFGZVRMyS+PXI7miP6t5b34IFd+52rvlKMZR3ZtaLQIRHxwPobbNi6F+Pa5neSksqoMVKl8N0Gl5dHKtgyXMmksu8BplhHmtYWX6y/TXmkgjse3Zv5tOJSsYCd68+PehmpITUunXaDy7Nc3HI06gUoHTPQV8LIbRfiuuXe+kmlEXXjmCc1Ct9pcLkVrMx6FoqmXXZGB7MeXiW5h+QoySHT63jkmWy2AdEK2mBIjUunXbphFopbcqRjI7SsP/C6wJr18E0P564UkV+bXsCt5T14N4OpiHNm5dWNExCpsfDbpRumPQOikM/h2mXzbFuTAprS1iki8oKIRLYtKo9U8EBGO7q+8/sjOoUtIFKj8NulGw70lXDluencHlrTfu4cWIzPLp8/TemrLzRQBMDjJIcbsxxsIbmW5BDJobGxsbYfOrh9XyYTDACgdlTUBRkQqXHptEs3tFIy08acWXmM3Hbh5Os7Bxaj/7QTNO3SG2eSfM7mfa+zHgBghYgcJPkhAE+QfFFEpo387HTOQ9rdj+1QF2QwpEbhA87phtb2OG0WUz5H3H7pomnva9qlZ14SEceArBcajQIhIm+SfAT1kZ6OM5694haPyQLqggyG1Lh03Ejt9jiVP1RyIDmb5HHW9wAuRD3Y65s0KftOm56pCzI4MqHw3baHc2blQ1yJWdTXGRwkryB5AMBHUZ/1sL3x/lySVqvvkwD8guRuAL8E8JiI/MTE9dPUV8dLplGO1IEmIZAql44TbimZ7/z+SMirMUvWfb1B4WXWg4i8AmBJENdft2phpsYeHhXBv238VNTLSD2pt/DLIxUcfn+6UieAWfke1I4mY+vs1PiNgKawpZDWxn+FBE8mKxbybVOijy8kd6edJHxJUVyqEZ2w2i209iMpFvK4e81SVGvJaDrglmMv0CratCNAYmS1lXwPseGyRZMPLycy3Mg2VPyaDVY1opeshJUistRvVoQT5ZEKVmzcgdPXP4YVG3dMNlKz2xLPPmYGBvpKsbYqLPm3fJp3Dix2jNFqClv6aO4N5ZU4xqP+4Ni619iauOak18cz3iQuLPzOtH0BABjx49m6OSzlbjVOc/J/WgoyaqvCmi1bKhaw8qxePPnimGvufMkhFqEpbOmjk2Z/+R5i8OolGHrtLdwfg1oTt5nJTvE0leFwCCtoa1UjCoBvNopQbGlUK64FgPnzp3YKdGp/7NQ4zSmX2RIur1ZFId8zuaWePTPnmnVwz5qluHHzaNvPzJG4dtk83DngfbSdXSBPU9jSidddW451ZQ8gFsre7p6zmhgO9JVUhiOmrUuH5E9JPmfzdXkH11khIucAuBjA9SQ/5nSiiGwSkX4R6e/t7Z183639sdPNMSHi2m6h6HELfMLsY/Dqxk/hnjVL4RbjJYCbPCh7a21bhisdBVxbA3mawpZevFq81y6bBwCTVnSUFPK5ts37VIajpa2FLyKf9HsRE9WIbu2PnbaJ1qQcpzYDXmtbLGFtt83uNN+n2fLxilbRZgOvaZlbhiv4x91vRJ6+WSzkseGyRRjcvq+ty0ZlODoCz/UyVY3o1v7YrXGaFSy6e81SAHUL3ArqHvI4Rs4K7gYRHNWAq2KHnSVsF5St1iYiH4dYLOQxevuFky4bnZkcX3z58EleAeCrAHpRr0YcFZFVJOcC+LaIXIJ6NeIjjcDuDAD/0E01oluwp13jNKeg7vGFvKebZbxaw4qNOzyfb0GXdTevX1HsaLWET1//WISrcWbDZR/0c9KZyfHGb5ZOaNWI7YI9bttEJ3fQsfkeFPI5T9vhyngVPR1m9VjC7rQ1V8tH6YS4DvFpve/UZRNfElO+5yfY4+Q2GT9cw12rvWfJdFKUS9QfEoPb9+HKc0uTRSdWxawGq5ROieMQnzT1/MkCieql063l0M4d5BRo6oRiIY/Zx8xAZbw6JQ+5Ml7FluGKKnfFN3buknffOxKZDz+fo+5QE0ZiLHw/tAsk+bWcCvkcNly2CDvXn49SsTAtW6d5mLqimOTTS05BPme2grBULOC65fMnrXe7T589M4fBq5aoEZMwEmXhd0u7QJLd8bfefc9T/5IcOcV6bzdMXVG6xS75YMtwBWv+0zw89uwb03pGtTIr3wMBXWNWdgPEnQoeleSRCYUPtHcHdZMRUcjnprlqtHRcCQqn5IMnXxzDyG0XYsXGHY6uyUI+h79uxKvcXJh2tSkahE0PmXDpdINTFW67QQ2ah5wOSA6SfJHksyQfIVl0OO8ikvtIvkxyfZBrard7dNtFWrJq5co74bU2RUkmmbHwO6E8UrEdjJLPsa3fUvOQU8MTAG4RkSMk/zeAWwD8ZfMJJHMAvgbgAgAHADxNcquIPB/EgtrtHt0qzgf6SiiPVHDHo3tdXT+6E003qvBtGNy+z3YwyuyZM6Ypbif/pir4ZCMijze93AXgKpvTzgPwcqPWBCS/D+ByAIEo/Ha1KG7HW/3/duhONP2owm/QrLid0u1bt7tOFbzA9GIUJdH8GYDNNu+XALze9PoAgGV2H+DWBdYr3SQfWMf7vvR42wJDTR1OP6rwMV1xO9G63XVr6KY3TiI4k6RdX6cvisiPAIDkFwEcAfCAzXlOQ8imv1lvCb4JAPr7+7ueq9lp8gFQl+92GTyW20dJN6rw4W3YhN12V1MwE89LbhPYSH4OwKcBfELEtrfqAQDzml6fCuCg2SX6p10NiLpysoNm6cBdQbtl5DgFuDTwlXxIXoR6kPYyETnscNrTAM4geTrJmQCuAbA1rDV6xU2+i4W8unIyhFr4cM9uaC1CaUan96SaewEcA+CJRqfXXSLy+eZOsI0MnhsAbAeQA3CfiOyNbsn2OMm31dZYyQ5q4aP73Hmd3pNeROSPRGSeiCxtfH2+8f7BRttv67xtInKmiHxYRL4c3YqdcZLv5rbGSjZQCx/+cuc1BVOJO1oboliowm+giltJMyrfCuDTpRPH8nNFURTFHr8+/CcA/EcR+QiAl1AvP59CU/n5xQDOBnAtybN9XldRFEXpEF8KX0QeFxGr6cwu1POQW5ksPxeR9wFY5eeKoihKiJjM0vkzAD+2ed+u/FydiYqiKCHTNmhL8qcATrY5ZLz8vPFZvnuOKIqiKNNpq/BF5JNux02Xnzf3HCE5RvK1dmv0yIkAfm3os+JOVn5Wvz/naaYW0gnDw8O/Vrl2RH8e/zjKNe11tDca5ed/B+A/i8iYwzkzUA/ofgJABfVy9M+EXZFIcsitb0qayMrPmpWf0420/Q705wkWvz78ewEch3r5+SjJbwAAybkktwFAI6hrlZ+/AOChOJafK4qipB1fhVci8kcO7x8EMKX8HMA2P9dSFEVR/JGlXjqbol5AiGTlZ83Kz+lG2n4H+vMEiC8fvqIoipIcsmThK4qiZBpV+IqiKBkhUwqf5AaSlUZG0SjJS9r/r+SQpSZ1JF8luafxdxyKej1Rkga5TpvsxlU+M+XDJ7kBwDsi8n+iXotpGk3qXgJwAerFbk8DuFZEno90YQFB8lUA/SKSpiKdrki6XKdRduMqn5my8FOONqlTkorKbkhkUeHf0Ojffx/JOVEvxiBZa1InAB4nOdzov5R1kizXaZTdWMpn6hQ+yZ+SfM7m63IAXwfwYQBLAbwB4G8jXaxZOmpSlwJWiMg5qM9ZuJ7kx6JeUJCkXK7TKLuxlM/UjThs1+zNguS3APxjwMsJk46a1CWdRjU3RORNko+g7hb4WbSrCo6Uy3XqZDeu8pk6C98Nkqc0vbwCwHNRrSUAngZwBsnTSc4EcA2ArRGvKRBIziZ5nPU9gAuRrr9lR6RArlMlu3GWz9RZ+G34G5JLUd8uvgrgf0S7HHOIyBGSVpO6HID7Utyk7iQAj5AE6jL8DyLyk2iXFCmJlusUym5s5TNTaZmKoihZJlMuHUVRlCyjCl9RFCUjqMJXFEXJCKrwFUVRMoIqfEVRlIygCl9RFCUjqMJXFEXJCP8f3OyPSByDoToAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def rmvn_chol(n, mu, Sigma):\n",
    "    \n",
    "    d = mu.shape[0]\n",
    "    \n",
    "    Q = np.linalg.cholesky(Sigma)\n",
    "    Z = randn(n, d)\n",
    "    X = Z @ Q\n",
    "    \n",
    "    return(X)\n",
    "\n",
    "Sigma = np.array([[1,0.9],[0.9,1]])\n",
    "mu = np.zeros(2)\n",
    "\n",
    "X = rmvn_chol(100000, mu, Sigma)\n",
    "plt.subplot(1, 2, 1)\n",
    "plt.scatter(X[:,0], X[:,1], marker='o')\n",
    "\n",
    "X = rmvn_chol(1000000, mu, Sigma)\n",
    "plt.subplot(1, 2, 2)\n",
    "plt.scatter(X[:,0], X[:,1], marker='o')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Algorytm generowania liczb z wielowymiarowego rozkładu jednostajnego"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Algorytm do wyznaczania $n$ zmiennych losowych z $d$ - wymiarowego rozkładu jednostajnego. Dla każdej zmiennej $u_{i}$, $i=1,2,\\ldots,n$ należy powtórzyć rozumowanie:\n",
    "1. Wygnereuj próbę $x_{i1},x_{i2},\\ldots,x_{id}$ z rozkładu $N(0,1)$.\n",
    "2. Wyznacz normę euklidesową $||x_{i}||=\\sqrt{x_{i1}^{2}+x_{i2}^{2}+\\ldots+x_{id}^{2}}$\n",
    "3. Wylicz $u_{ij}=x_{ij}/||x_{i}||$ dla $j=1,2,\\ldots,d$.\n",
    "4. $u_{i}=(u_{i1},u_{i2},\\ldots,u_{id})$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0xbd0a48>"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD4CAYAAADsKpHdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO2dfZRU9Znnvw9FtTTEpCECkVKEEA5GbeyOfQSXczJqQBBWKTUGETPu7owcd+OeRTO90wQ2QBZCRzZCMpOJC1nPMQdG8QUrGIgIOm42DDA2dkNLlIBKkAKho7Qa6EDTPPtH3dsW1ffld+u+1H15PufU6ap7f7fur6qfus99nt/zQswMQRAEQehX6QkIgiAI4UAUgiAIggBAFIIgCIKgIQpBEARBACAKQRAEQdDoX+kJlMMll1zCo0aNqvQ0hJiye/fuPzHz0EqcW2Rb8BM72Y6kQhg1ahRaWloqPQ0hphDRHyt1bpFtwU/sZFtcRoIgCAIAUQiCIAiChigEQRAEAYAoBEEQBEFDFIIgCIIAwKMoIyJ6AsC/B3CCma8x2E8AfgJgOoDTAP4DM7+h7Zum7UsB+AUzN3sxp6SxMNeOdTsPQy9VSAAYQKamGo1TxyFbn0GuNY8lL+7DydPdAICa6jQW3341svWZSk071Ihch4Mpj72GAydOOTpmQIrw9rLpPs0ovpAX1U6J6OsA/gzglyY/nOkA/isKP5wJAH7CzBOIKAXgDwCmADgC4HUAs5n591bna2ho4CSH5uVa85i/YS+6us8rja9Op3DXdRmsf/19dPf0/X8PTPfDD+8cL4pBg4h2M3ND0HINiGyPX/QSPjnT49v79wPw2Ky6xMq6Lttm+z2xEJj5t0Q0ymLITBR+VAxgJxHVENGlAEYBOMjM72qTfVoba/vDSQq51jxWbNmPfGcXUkToKUOBd3X34Kld75see7r7PBqf3QMAvT8U/bxHO7swosjKSBIi1/6yMNeOtTsPB3rO8wDmrW/DvPVtAIBJY4Zg3QM3BDqHMBNUYloGwPtFr49o24y2TzB6AyKaC2AuAIwcOdKfWYaMXGsejc/uQff5woW8HGWgY3ds93nGii37e11L8ze0o6u7cKeW7+zCw+vb0PLHj7A0W1v2HGKIa7kGkifbflsBTtj+zkcY1bQJAPD5i1LYu2RahWdUWYJSCGSwjS22993IvBrAaqBgVns3tfCyeOO+XmXgFhXr4mhnFwBgxZb9vcpAhwGs23kYDVcMSZylYIFruQaSIdtz1uzA9nc+qvQ0LPnkTE+vcliVULdSUArhCIDLi15fBuAogCqT7Yki15rH4o370NlVWOwdPDCNRbdd3fvaLXZrCDojaqoBfKYYSmGg14oQAIhc2/KV+ZtwLoIqTncr3TdxZKKs4qAUwkYAD2m+1AkAPmbmY0TUAWAsEY0GkAdwD4B7A5pTKCh1CwHAydPdaHxujyfvXxxl1HDFECx4oR2nzvY119P9CI1TxwEoKIa8iVIwUxYJReTagFxrvtdHH3XW7jyMtTsPJ8ad5FXY6VMAbgRwCREdAbAIQBoAmPlxAJtRiMQ4iEJ43n/U9p0joocAbEEhPO8JZt7nxZyiwoot+w3dQt09jH4EOPUYZSwWgLP1GaXw08ap4/Dw+jZDH4duReiYWTdxsCJErp1z5YLN+IuFFeoFZgvB5YSnqqK7k+LuSvIk7DRooh6aVxzBY/ftp1PUx80TRP5AaV4DUHA9Lb+z9oJIpFLrRp/zim9eG9kfjl1onp9EWbZ1/7tXeOWu8XoRe+ywQdj6yI2evV+QBBJ2KthTHD6qJ43Zod/tVyL8c2m2Fg1XDLE8t5V1I2sNyWHCsq04/ulZV+9Rne6H5T7lwpS6etwucB84cSq21oJYCAFQGsapQhTuskc3bbJUbHpkU4oIsydcHpnFObEQ1PAicqiSF1Uv5h81a0EshBBgFMZpRVT88FaLz8BnuQ89zL0JSFFRCoI1bvz1YbmzLl6HKHchPG7WghS3CwCVyJxMTTUONc/AoeYZaP3+LZEQrsap45DuZxRyb8xTu963HySEnjlrdpSlDPoTcKh5RihlO1ufwaHmGRg7bFBZx89b34aFuXaPZxU8YiH4QGnZh5qB6d6IHiOq06nekM8oof+wS6OMzD5rDzMmNb+ayDIYcaDcUhNRiuXX3T/lfNa1Ow9j17sfRsqFVIqsIXiM0XpBuh8BhAuihYyqkcaFMfM3W2ZFl0YrhQ1ZQ+hLuQvHUXellPO5+xNwcPkMn2bkDjvZFpeRxxitF3SfZwyq6o9MTTUIBSWwclYdDjXPwPammyP9gzFi9oTLLfd3dfdgxZb9Ac1GcMuVCzY7vijeN3FkaN1DTti1YArum+isvtQ5LmRoRxGxEDxAJa+AALzXHM67Bj9YmGu3rLCqE8YFdLEQPsNpbkGcM3rLsRbC5i4TC8FnFuba8fD6NuRtksxKM3zjztJsLd5ZPh0Zm8998nR3bBbk4oYTZZDuV3APxVUZAOVZC2t3HsaEZVt9mpH3iEJwQa413yeb14ioLhp7QePUcahOp2zHrd15GHPW7AhgRoIKTi5iwy+uwoEfRt89pMLSbC0ONc/A8IurlI85/ulZTHnsNf8m5SGiEMok15rHd5/ZY6kM9PWCMC+g+k22PoPld9baWgpAoTa9WAqVZ86aHY5cI7sWTPFxNuFk14IpONQ8A/0Vo64PnDiFXGve30l5gCiEMsi15tH43B5L/3imphrvxXTR2CnZ+gy2N92spBTW7jwciR9OXJny2GuOsncPJWhdzIiDy9WVQhRco6IQymDJi/ss+woQkFgXkRWqiWzzN7SLUqgAo5s2KSecjR02KPHKQMeJUli783ColYIoBAfkWvOoW/KyZZIZAZgzcWTirQIjsvUZrLj7WgxMW4udhKUGzyibulTFTBozJNLJV37gVCmEFVEIiuilnu26mK2cVReqMLOwka3P4Pf/81ZMGjPEcpw04gmO8YteUh5738SR0pTehIPLZ2BASk0rhDXySBSCIktetO9vXFOdFstAkXUP3GAZwscAJjW/Kq4jn5mzZodyr4BVcrNjy9vLpiuNO/7pWUeKOCg8UQhENI2I9hPRQSJqMtjfSERt2uNNIuohoiHavkNE1K7tC09GThG51rylmwgolKdYfPvVAc0oHizN1mLVrDrTsNR8Z1fF1xPiLtuqC8hRL0ERJKprK5+c6cFoj5sKucW1QiCiFICfAbgVwFUAZhPRVcVjmHkFM9cxcx2A+QD+LzMXS+JN2v6KZIfaseRF6+6HKSKsuDvcvQvCil1Yald3j+337xdxl23VxDNRBs451DwDKs4jBkKVo+CFhXA9gIPM/C4znwXwNICZFuNnA3jKg/MGgp11kE4RfvwtUQZu0MNSzX5AJ093o/4HL1fCUoitbKsqg/skQKJsVEvV+NUHuhy8UAgZAMWF7o9o2/pARAMBTAPwfNFmBvAyEe0morlmJyGiuUTUQkQtHR0dHkzbHj35zHxOCH1XsyhhVd7j5OnuSriPYinbVy7YrDQubHV4osiqWXVK48LiOvJCIRjd2Jmtvt4GYHuJST2Jmb+Ggln+HSL6utGBzLyamRuYuWHo0KHuZqyAXsbaKvls5bfElPYSu9yNCoSjxk62F+ba8ReLHBqdscMGiTLwgGx9Rqn+EQOhKN3ihUI4AqC43vFlAI6ajL0HJSY1Mx/V/p4A8AIKZnrFsWt7KRFF3pOtz6CmOm05Jt/ZFaSVEDvZVomBJ0DyDDxkabZWSSm47e/sBV4ohNcBjCWi0URUhcIPY2PpICL6AoC/AvCrom2DiOhi/TmAWwC86cGcXGMVB1+dTklEkU8svv1q22J4jc/uCUopxEq2Vd0SSSrTHhRLs7VKiWuVdh25VgjMfA7AQwC2AHgLwDPMvI+IHiSiB4uG3gHgZWYuXkEZDuB3RLQHwL8B2MTMFQ3OzbXmMan5VVO/QIoo0cXq/EaPOrKi+zxj/oa9vs8lTrI9Z80OpUxku4RBoXxUuqgx1Nd4/EAa5BRh1P6ymLC3fowTk5pfRd4mW9mvRc84NshRjSqS+kT+o/K/8CvUVxrkOOB7G/aaKoOkl7EOGpU+ClIZVQ2Vdo4DUiTKICDGDhtkO2be+rYAZtIXUQgaC3PtON193nAfAVLGOmB01xHZ+F2DcB1FmVxrHucUnACqJRcE92x95EalmkeVcB2JQtB4atf7pvuS1v4yLGTrM1j5Les47q7u82IlWPDIM/Z3mk7bQgruUVHAf+nhwEtli0LQsMo3kN4GlSNbn8GgKmvXkZTKNmbOmh2wqccIAiTfoEKoJK0FXSo78Qoh15pH/Q9eNt1PBHEVVZhld1hfsKRUtjEqce0SYlo5svUZVCm4joK0gBOtEHKteXz32T2WtYrmTBBzutLYZXv2IxK3UQkqrgYJMa08j37zWtsxQS4wJ1ohLHlxH3osbGqp5RIe9GxPo/upHmY0PhdYslokUHE1SKObyqNa2iKotYREKwS7VpiiDMLF0mwtVs6qg1Fb5u4erliZ7LChohhlITk8qFxnglpLSLRCsEIii8JJtj5julBq18QoKai4GORmJ1wMv7jKdkwQVkKiFYJVITWJLIomSXcbqTRbkQS08LFrwRTbMUFYCYlWCItvvxppA/+DNAUJN1aKvNItNyuNXbMVcRWFF5UwVL+thEQrhGx9BivuvhaZmmoQCuUppJF4+DFT5EChZ8LijclcS1DJbBXZDi8qN6F+Wwn9fX33CJCtz4g1EDH0/5eZr7yzqxu51nzi/q92jW8kzDT83DdxpO1F30/ZTrSFIESXbH0GGYuF/6RlL6usHUiYafhZmq21Vdx+5iUkTiHo/Q5GN23CpOZXE+1vjjpWC/92pbPjht3agUoUixAOVBS3X4XvPFEIRDSNiPYT0UEiajLYfyMRfUxEbdrj+6rHeone7yDf2QVG4aKR9EXIKJOtz2DwQPMF5qALg1UKFflViWIRooNKX+xycK0QiCgF4GcoNBK/CsBsIrrKYOj/Y+Y67fEDh8d6glGf5Ao0bhc8ZNFtVxtmLwPAOg/6JUThZscuIU/CTKNHpaLBvLAQrgdwkJnfZeazAJ4GMDOAYx1jVgRNiqNFl2x9xrQ1JANofLZ8f2tUbnYkIS9+qESDjV/kfUdWLxRCBkBxM4Ej2rZSbiCiPUT0GyLSO9SrHusJZtnHkpUcbawWl7vPF8pAl0nob3bsPtvnL7IuHS6EF7t1n0/OGHd3dIMXCsHIYi+9aXsDwBXMfC2AfwCQc3BsYSDRXCJqIaKWjo6OsiZq1JaxOp2SrOSIY/f/UykDbUIgNztuZNvus+1dMs3R+wnhQWXdRyW6zAleKIQjAC4ven0ZgKPFA5j5E2b+s/Z8M4A0EV2icmzRe6xm5gZmbhg6dKijCeqRRQ+vb8OAdD/UVKd7E9GkT3L08fH/F8jNjhvZtsIkd0+IEHZrCXbRZU7xQiG8DmAsEY0moioA9wDYWDyAiL5EVOiOS0TXa+f9UOVYt5RGFp083Y0z585j5aw66ZMcI6zcRi4I5GanXL4yf5Pl/null0fkCTqz3LVCYOZzAB4CsAXAWwCeYeZ9RPQgET2oDfsmgDeJaA+AnwK4hwsYHut2TsVIZFEysHIbucjQDfXNzjmbyEMpU5EMvAyv9qR0hXZntLlk2+NFz/8RwD+qHuslElmUDHRLr/HZNnSfv3DfoQ+7ykr3Z+ZzRKTfsKQAPKHf7Gj7H0fhZuc/E9E5AF3QbnYAGB7r4iNegN1isniL4oNdOYt1Ow97pvxjX8voC9VpdHb1DcuTyKL4odel0t2EumWoJyDqY5wQ1psdu8XkOVLVNDYszdZaKgQvU9RiXboi15rHqbPn+mxP9yOJLIox4iYUd1HcsLP4vCplEWuFsGLLfnQbpHh/bkB/WUyOMXF3E9r5jKVuUfyws/i8KmURa4VgdgHolMzOWBP3BES78shStyh+BGXxxVohxP3CIBhjlIBIAG660rsYf0EIGruscxcZ+b3EWiFIZnIyydZncNd1mQv8rgzg+d35yFe2jfr8hfKxyzp3kZHfS2wVQq4137u4mCqEiUtmcoL4l7c7+kRfxGFh+Xsb9lrul57JghtiqRCKs5MBoIe51zIQZZAM4rqwfLo0yaIEiS6KN34HDMRSIUjYoWC2TvSFavOGOoIQduwCBty6FGOpEOJ6dyio0zh1HNIG1d1OnT0XWT+8Xbjpqll1Ac1ECCt2zZLsiKVCMLsLlOii5JCtz+BzA/om4nf3cGQtxX/eZR1uKu5QwW2zpNgpBMlOFnTM8k2iaime96eNrhAx/KxTFTuFINnJgo6ZRTiwKn5dxAbF8DMJxqy0cQ26cYnGTiFIdrKg0zh1HFKG6wg9npYMDgK7+S67Q6KLkoLdje289eX3EY+dQpDsZEEnW5/BeRM/i50/PmzYlasQ61fwgtgpBMlOFooxc7ufZ8n6FaKLXwmInigEIppGRPuJ6CARNRnsn0NEe7XHvxLRtUX7DhFROxG1EVGL27lk6zNYfmctMjXV0jdZ6M1SN0Il2ihMsi0IOn4lILpukENEKQA/AzAFhT6yrxPRRmb+fdGw9wD8FTOfJKJbAawGMKFo/03M/Ce3c9HRG6UIwuwJl5u6W+yijcIo20aMHTbIz7cXIsicNTuw7oEbHB/nhYVwPYCDzPwuM58F8DSAmcUDmPlfmfmk9nInCg3HBcF3lmZrUZ02FnOFdaVIyPbWR24M+pRCyCm30J0XCiED4P2i10e0bWb8DYDfFL1mAC8T0W4immt2EBHNJaIWImrp6OhwNWEhWSy/c3y560qhkG1Z6xCM6O9DQoIXCsFoWoZreUR0Ewo/mr8v2jyJmb8G4FYA3yGirxsdy8yrmbmBmRuGDpW69oI6LtaVQiHbf/fsHrt5Cgnk4PIZnr+n6zUEFO6aLi96fRmAo6WDiGg8gF8AuJWZP9S3M/NR7e8JInoBBTP9tx7MSxB6KXNdKRSyfU5SlIWA8EIhvA5gLBGNBpAHcA+Ae4sHENFIABsAfJuZ/1C0fRCAfsz8qfb8FgA/KHcieg+Eo51dGFFTLeWuhQsoQz5CI9tmyIKy4CWuFQIznyOihwBsAZAC8AQz7yOiB7X9jwP4PoAvAvgnKoQBnmPmBgDDAbygbesP4J+Z+aVy5qH3QNDLXuc7uzB/QyG7U5SCUI58hEW2rZAFZcGMciKNvLAQwMybAWwu2fZ40fO/BfC3Bse9C+Da0u3lYNUDQRSCUK58hEG2BaEcyok0ik2msvRAEKwQ+RAEe2KjEMx6IEiHLAGIbo2rKY+9VukpCCHmov7eXsJjoxDMKhRYVC4QEkRUa1wdOHGq0lMQQsyP7hrv6fvFRiGYlbeWstcCIDWuhHjitfx6sqgcBvr3A7rP990edpeAEBxxq3HlV8VLIbnEwkKYs2aHoTIAEHqXgCCUi18VL4XkEguFYBVeFac7QkEQBD+JjctIEFSQbHZBMCcWFoIgqKBnK+c7u8D4LFtZqokKceXKBZvtBxURC4UwacwQR9uFZGKVrSwIceQvPc4KI8ZCIax74IY+F/9JY4aU1TFIiC+SrSzEkUFVKftBisRmDUEu/oIdA6tSOHW2x3C7IESVZXfUYt76Nk/eKxYWgiCocNpAGVhtF4Qo4GVQhCgEITGYeVOl/YwgFBCFICSGlElhK7PtgpA0RCEIiWHilwc72l5pFubaKz0FIWF4ohCIaBoR7Seig0TUZLCfiOin2v69RPQ11WMFwSsOfWgcTWS2HaisbK/bddjpIYLgCtcKgYhSAH4G4FYAVwGYTURXlQy7FcBY7TEXwM8dHCsInuA07LTSss2yuCEEjBdhp9cDOKi1DAQRPQ1gJoDfF42ZCeCXzMwAdhJRDRFdCmCUwrFKSEkCwY4RNdXIG1z8LSrihkK2jZBlD8EPvHAZZQC8X/T6iLZNZYzKsQAAIppLRC1E1NLR0XHBPilJIKhQRpOcisu2GV8ZOkhpnCA4wQuFYHSvUmrsmo1RObawkXk1Mzcwc8PQoUMv2CclCQQVsvUZ3HVdpjeqKEWEu66z7JFQcdk2QzqpCX7ghUI4AuDyoteXATiqOEblWFukJIGgQq41j+d359GjOed7mPH87ryVJVlx2RaEIPFCIbwOYCwRjSaiKgD3ANhYMmYjgL/WIjImAviYmY8pHmtLVBuoC8FShiVZcdkWhCBxrRCY+RyAhwBsAfAWgGeYeR8RPUhED2rDNgN4F8BBAGsA/BerY53OoXHqOKRTF1ro6RRJtzThAowWlK22V1q2V82qczJcEFzjSXE7Zt6Mwg+jeNvjRc8ZwHdUjy1vEjavhcSTIup1F5VuN6OSsp2tz3hWtEwQVIhFpvKKLfvRff7CH3r3eZZFZeECjJSB1XZBSBqxUAhOXQFCMqmpTjvaLghJIxYKQYqWCSqYiYOIiRBlvKx5FQuFIK4AQYXO092OtgtCFFi707uaV7FQCOIKEFSoGWgiJybbBSFpxEIhiCtAUOFMt3FnNLPtghB1hl9c5Wh8LBSCmcl/UlwBQhGnu8872i4IUWfXgimOxsdCIVhlJEuTEUEQBDVioRCsMpKf2vW+6T4hOciNgSDYEwuFYNX3QCKNBMD6xmBwRBeVRckJXhMLhQBILoJgjdWNwaLbrg5wJt4hLTYFr4mNQpg94XJH2wVBJ6qd9cT4FbxuAhYbhbA0W4tJY4ZcsG3SmCFYmq2t0IwEQRD85RGPix/GRiHkWvN44/DHF2x74/DH0kZTiDSlNzmCUIzXAdOxUQjSRlOII+seuKHSUxASRGwUgrTRFARB+IzPX5RyfIwrhUBEQ4hoKxEd0P4ONhhzORH9CxG9RUT7iOi/Fe1bTER5ImrTHtPLnYtZctoXpJ5R4rFyG/YzLXsSHtm24soF7ntLCfFk75Jpjo9xayE0AXiFmccCeEV7Xco5AN9l5q8CmAjgO0R0VdH+lcxcpz3Klu7GqeOQNvh1nzp7TtYREs6SF807V947YaTZrtDIthV/6ZFQI8E73CqEmQCe1J4/CSBbOoCZjzHzG9rzT1HoL+t5nF+2PoPPDejbEbS7h7F4o+M2zUKMsKppZRGFFhrZFgQj/EhMdKsQhjPzMaDw4wAwzGowEY0CUA9gV9Hmh4hoLxE9YWSWFx07l4haiKilo6PDcIxpvfuubrESBKeESrYFoRQv+yDo2CoEItpGRG8aPGY6ORERfQ7A8wDmMfMn2uafAxgDoA7AMQA/NjuemVczcwMzNwwdOtRwjFWRO4k2Eko5/vQCXHPNNX0eAGqcvI/fsi2hp4JTzNbG7OjrYymBmSeb7SOi40R0KTMfI6JLAZwwGZdG4Qezjpk3FL338aIxawD82snkS2mcOg7zTBI1JNoomVhZhlf+p0fR+v1b+mwnok4APWGR7XUP3IBRTZvKPVxIIBZrY5a4dRltBHC/9vx+AL8qHUBEBOD/AHiLmR8r2Xdp0cs7ALzpZjLZ+gwu6m/8kQZWOQ/BEqKPlWVoU8MoVLJtxfhFL/n11kJEKbdCg1uF0AxgChEdADBFew0iGkFEelTFJADfBnCzQQjeo0TUTkR7AdwE4GGX88GZc8a5e6fOSlesJGJlGdrUMAqdbJvxyRmR7aQxZ80OX97X1mVkBTN/COAbBtuPApiuPf8dAEOPFjN/2835BcGOL1Sn0dnVN9jArt922GSbAEiAqaCz/Z2PfHnf2GQqC4IRcem3vXJWneV+6Y0geEHsFMIgi7UCCT1NHqahyBHrt21Xolt6Iwg6900sb0EZiKFCWHaH+WKKhJ4mD7NQZKsQ5SgivRGSg10QgZuS/7FTCFZ3UhJ6mjwap45DdfpCq7E6nbLswy0IYcbPIILYKQQAyCTkrlCwJ1ufwfI7a5GpqQahIBvL76yNZJc0SVAT/CaWCkHuCoVisvUZbG+6Ge81z8D2ppsjqQwA+94IsrAsuFk/AGKqEOJ0VygIqvhR20YIF3ZK323LYFd5CGEmW58RBZBgcq15rNiyH0c7uzCiphqNU8clQh5yrflEfM6k4rfSj6WFUEyuNY9Jza9idNMmTGp+VUJPE0CuNY/5G9qR7+wCA8h3dmH+hvZY/O/HDhtkuX/+hr0BzUQIG4MHum8GFmuFYHRheHh9m/haY06c+2tvfeRGy/1d3V63XRfCgl25CpvaXErEWiEYXRgYwLqdh2NxtygYE/f+2uX0yhWij125Ci9chbFWCGYXAIYkqcWZGhPTOS5hx3a9cqX6qVAusVYIVheAuNwtCheyMNdu2DIznaLEhB1L9dP4YefmHn5xlSfnibVCaJw6zrgUJeJztyh8Rq41bxqFMaiqv0TfCJHFLrpo14Ipnpwn1gohW5/BnIkj+ygFSVKLJ0te3Ge672ODEthRxu6O8MoFmy33C9EhyPVOVwqBiIYQ0VYiOqD9NWwkTkSHtGYhbUTU4vR4NyzN1mLlrLrechYpot6IE1lYjhdGriIdpxZh2GXb7o7wLz1S7S4u/Pfn9lju98pdBLi3EJoAvMLMYwG8or024yZmrmPmhjKPL5tsfaa3nEWPVhYyTrHpgj1lWISRkG0rJLw6Hpy1Ue5euYsA9wphJoAntedPAsgGfLwycY5NF6xjtAem+5WzfhB62barWyOlLKKPXcSY1wUP3SqE4cx8DAC0v8NMxjGAl4loNxHNLeN4ENFcImohopaOjg7HEzWLKspLtFHkybXmLWO0f3jn+HLeNvSyrVK3RizgaGMXMWZX8NAptrWMiGgbgC8Z7Frg4DyTmPkoEQ0DsJWI3mbm3zo4Hsy8GsBqAGhoaHDsIB1RU2148SdI/ZeoY2flmf1vJ0+ejA8++MBoV42D01dUtvsRcN7iiHnr20S2I0ollLmtQmDmyWb7iOg4EV3KzMeI6FIAJ0ze46j29wQRvQDgegC/BaB0vBc0Th2Hh9e39WlUriepyY8muljllKQsmidv27bNcDsRdQLoiYJs3zthpLiGYsrD69ss99vVtSoHty6jjQDu157fD+BXpQOIaBARXaw/B3ALgDdVj/eKbH2mjzLQkSS1aGMVQTR7wuXlvm0kZHtpthYDUuZKD5DF5SiSa82bXq907OpalfWBaO8AABCkSURBVINbhdAMYAoRHQAwRXsNIhpBRHog9HAAvyOiPQD+DcAmZn7J6ni/kE5q8SPXmsepM+cM900aM8RNffjIyPbby6Zb7hcLIno0PmttHfiFq34IzPwhgG8YbD8KYLr2/F0A1zo53i8ap47D/A3tF0QbSZJadNGr2ZZGjw0emMai26525QaMmmzbsTDX7rp5ihAcdkVr/XAXATHPVC6ltJNaTXUaA9L98PD6NumVEEEWb9zXRxkAwMAElqmwu0CIlRAdRjVtstzfn/xxFwEJUwjAZ/11V86qw5lz53HydHfsmqgkgVxrHp0m5SiSuCakcoEQ2Y4HB5fP8O29E6cQdMwS1azq4QjhwSrUNKlrQnY/5nk2UStC5bFLREv7fMVOrEIwu4s8ebpb7qQigFVCYVLXhB6bVWc7RiKOwo1dItqKu+3/x25IrEKwuouUchbhxuqiNnhgOnHrBzrZ+oxtNzVZSwgvdmsHgDdd0axIrEKwuotMog86KizMtZte1Aje9JWNMnbd1ACxEsLIlMdesx1jV7vKCxKrELL1GdRUx7vVYtzIteaxzuIOl+H/HVQUsEtUWys9xUPHgROnbMcEETacWIUAAItvvxrV6QtNbMlLCC+LN+6zzN40SzxMGnaJaoAsMIcJFetglcL6kBckWiGU5iVkaqqx/M5aucsMIQtz7aZhpkDBXSSK3BkTlm2t9BQEqFkHQV2TXGUqx4FsfeaCLzvXmsek5ldxtLMLI2qq0Th1nCiICmPVK1lnzsSR8n8qYtWsOlsr4PinZwOajWCGXZgpEMzagU6iLYRS9FII+c4uSVYLEfM37LXcf9/EkVKWoYRsfUapvIGKu0LwD7sw0/79KFDZFoVQhFmy2uKNkqxWKXKteXRZFHYZPDAtysAElezlAydOSdRRhVAJM/1fdxuWyvINUQhFmIWbdnZ1y4+mQnzPxjpIepipHSqLkZKbEDyjFZRBfwo+ak4UQhFW4aZrdx627NsreE+uNY/TFtbBoKqUrBvYkK3PKPXdFdkOjjlrdtj2OgD8rVlkhiiEIuyiVLa/85FYCgFilzG+7A5xFamg0ndXZDsY7Pp/66gocT8QhVBEtj6DwQONk9V0ntr1fkCzEawyxqvT/cQ6cIDKBUYS1vxHJf9jQIqUlLgfuFIIRDSEiLYS0QHt72CDMeOIqK3o8QkRzdP2LSaifNE++4wan1l029WwyvPsYZYfjc/oob9WZvXyO8f7Ooe4yfa6B26wlGsdSVjzD5UQU0AtsdAv3FoITQBeYeaxAF7RXl8AM+9n5jpmrgNwHYDTAF4oGrJS38/Mm0uPD5psfQZzbOJ+561vE5+rTxSH/hpBKISZBmAdxE6232tW80lLwpr3TFi21TbEFAguI9kMtwphJoAntedPAsjajP8GgHeY+Y8uz+srS7O1tib29nc+EqXgMQtz7Zi3vs2wCxpQyCRfOasuqDDTWMq2ygXn+Kdn8ZX59lEwghpz1uxQSgIckKKKu0HdKoThzHwMALS/w2zG3wPgqZJtDxHRXiJ6wsgs1yGiuUTUQkQtHR0d7matwLoHbrDNENz+zkfiPvKIOWt2WIY/EoDtTTcH+YOJpWxn6zMYfnGV7bhzLJaCV6gsIgOVdRXp2CoEItpGRG8aPGY6ORERVQG4HcCzRZt/DmAMgDoAxwD82Ox4Zl7NzA3M3DB06FAnpy6bpdla24Jp0mHNPSqRF35UoJ08eTKuueaaPg8ANU7eJ2qyvWvBFKVxxz89Kzc8LrlygZqnsNKuIh3bWkbMPNlsHxEdJ6JLmfkYEV0K4ITFW90K4A1mPl703r3PiWgNgF+rTTs4GqeOs1xoO3m6G1Mee823ptdJwC681K8KtNu2bTPcTkSdAHriLNsqtY6AwnpZpd0YUUUlExkIbE1MCbcuo40A7tee3w/gVxZjZ6PEpNZ+aDp3AHjT5Xw8RyWx58CJU1ITxgV2DYkqVIE21rKtWusIUL+wCZ+h+p0Nv7gqVKVX3CqEZgBTiOgAgCnaaxDRCCLqtZWIaKC2f0PJ8Y8SUTsR7QVwE4CHXc7HF9Y9cIOSUhDzujys3EEVvHuKvWxvfeRGpfUEQJSCE5x8V6ruu6AgZpUk6nDR0NDALS0tgZ+3bsnLljX5MzXV2N50c4Azigd6qGlpdNGkMUMqkqBDRLuZuSHwE6Mysj1nzQ7lhc9DiqGrSeUr8zfhnOIltRLfpZ1sS6ayAxbfbl1ILd/ZhdFNmzCp+VWxFhxg1Kho1ay6imVrJg2ViDqdUU2bRLZNGNWkrgyC7HHgBLEQHDLlsdeUOhwRCk1bwuQfrDQLc+14atf76GFGigizJ1weyu8naRaCzvhFLyklTwHSg6IUJ5ZBpSxfQCwEz9n6yI1Ki3EMYJ3UhullYa4da3ceRo92A9LDjLU7D0tBtRCxd8k05bFrdx6WPAUU3J1jHCiD+yaODLXlKwqhDLY+ciNWzarrdXGYwQC++8weUQowLwooxQLDhRO/9vFPzyZaKeiZ9T2KyuDzF6VCb1WJQiiTbH0G25tuxnvNMyyT13qYpQ0n0GsZqG4XKocoBXvGL3rJUWOhz1+UcmSBVQpRCB7QOHWcpaXQ1d2DeevbEr3YnCLjb8hsu1BZnCoF1UqecWBU0ybltRagsGYQBWUAiELwBL1Cqt2lLd/ZhXnr2xKRxLYw144x8zdjVNMmjJm/GV8eOtBw3OwJlwc8M0EVJ0rhkzM9iYhAcpqPUckF5HKQKCMPybXm8d1n9ii5QVIE/PhbdaFJWfeCXGseS17ch5OnjXM1xg4bhHc7TkuUkQVhlG2nF8HhF1eFLuHKLROWbVWqWFpMGCOx7GTbtpaRoI5+cTdKsiqlh4HG5/ZccFyU0aOIrHi34zTeWV75io6CMw41z8Dopk1KfYCBggtpVNMmjB02KPI1vlTDzEtZNSuaN3viMvKY4iQrO7p7OPJrC7nWPL76P36jtMAmC8jR5b3mGRiQcrbec+DEKYxq2hTJviFz1uzAqKZNZSmDQ80zIqkMAHEZ+UquNe+oJeGgqhSW3VGRQm6OWZhrx7qdh5XvGoHCAnIULARxGZlT7h0zEA1XktPfbClhL+0hiWkVxElFSQA4dbYHjc+FP29Bdw85vZWQBeToo+fglIPuSgqjxTDlsdcwqmlT2cpg1ay60CsDFcRCCACnd1UpIpxnxoiaajROHRc6i2HM/M2O3D9EwJwJ4VtgM0MsBDXcWAsA0J+Ag8srexEtZ7G4lCgpAjvZFoUQELnWPBZv3GdZLdWKMEXlOIk6CWOkhR2iEJzhRWnsIBO3nNRssiKKi+aiEEJIrjWPh9e3OXa56GR8shxyrXms2LIfRzu7LK0TFQthYLoffnjn+NBZNyqIQnCOk+JuKngdv++VEgCil1tQjCiEkKISpqlKTXUai2+/2tXF16gnQXU6ZditzGruUVoYN0MUQnm4XZC1Y0CKbBvRO+ntUA5Rcg8Z4atCIKK7ASwG8FUA1zOzoSQT0TQAPwGQAvALZta7Tw0BsB7AKACHAHyLmU/anTfKP5piSt1IRIBb/VxsPeh3/PnOLqSI0MOMmuo0iIDO090XWAGTml9F3qCVpVnTn4W5dqzbdbh3vlG2CEohot0AfgSR7bLw8mYnLEQhQkoFvxXCVwGcB/C/Afyd0Y+GiFIA/oBCm8EjAF4HMJuZf09EjwL4iJmbiagJwGBm/nu788bhR2OEWecwp1SnU7jrugye3523fS/dCjBzYREKMehJQlMI34bItivioBiiuAZmha9hp8z8FjPvtxl2PYCDzPwuM58F8DSAmdq+mQCe1J4/CSDrZj5Rx0lSmxVd3T14atf7Soqlq7sHK7bsN+1rbNXvOM6IbLtnabYWh5pnlB2mWkn0MNI4KQMVgshDyAAoLnp/RNsGAMOZ+RgAaH+Hmb0JEc0lohYiauno6PBtspVGL6t9qHkG7ps4suxqoE7CQo92dqFx6jhUp1MXbK9Op9A4dVxZ508IItsKZOszONQ8o1emw0qKPlMEcXB9loNtLSMi2gbgSwa7FjDzrxTOYXRFc+ynYubVAFYDBbPa6fFRZGm2tvcOpXg9QAV9zUCFETXVvT8AlSijuDB58mR88MEHRrtqFN9CZNshxTLtNo/BK6IcNeQ1tgqBmSe7PMcRAMUpqpcBOKo9P05ElzLzMSK6FMAJl+eKLdn6TO/F2a6qqNM1BN0KKD5HEti2bZvhdiLqVHwLkW0XlMbwe5EkpoJKtFJSCaLa6esAxhLRaAB5APcAuFfbtxHA/QCatb8qFkfiKb1wm+UPNFwxRDnKSCgLkW0PMYricWtFRLXqaKVwG2V0B4B/ADAUQCeANmaeSkQjUAjBm66Nmw5gFQqheU8w8zJt+xcBPANgJIDDAO5mZtsg4jhFYgjhQ4syWgaRbSFmSGKaIDhEEtOEuCLVTgVBEAQlRCEIgiAIAEQhCIIgCBqiEARBEAQAEV1UJqIOAH802X0JgD8FOB0zwjIPQOZihNU8rmDmoUFORsdCtsPyvQEyFyPCMg/AhWxHUiFYQUQtlYoQCeM8AJlLmOehSpjmK3MJ7zwAd3MRl5EgCIIAQBSCIAiCoBFHhbC60hPQCMs8AJmLEWGZhyphmq/MpS9hmQfgYi6xW0MQBEEQyiOOFoIgCIJQBqIQBEEQBAAxUAhEdDcR7SOi80RkXsWPaBoR7Seig1qPW6/nMYSIthLRAe3vYJNxh4ionYjaiMjTKmZ2n5EK/FTbv5eIvubl+R3M40Yi+lj7DtqI6Ps+zeMJIjpBRG+a7A/k+yiHsMi1do6KynZY5FpxLtGWbWaO9APAVwGMA/AagAaTMSkA7wD4MoAqAHsAXOXxPB4F0KQ9bwLwI5NxhwBc4sP3YPsZAUwH8BsUOn1NBLCrQvO4EcCvA5CNrwP4GoA3Tfb7/n24mHso5Fo7T8VkOyxy7WAukZbtyFsI7L4ZuldUuqm6ymecCeCXXGAngBqtm1fQ8wgEZv4tAKseBEF8H2URIrkGKivbYZFr1bkEgl+yHXmFoIhVM3SvUG2qzgBeJqLdRDTXw/OrfMYgvgfVc9xARHuI6DdEdLXHc1AliO/DT4KafyVlOyxy7eQ8kZXtIFpouoaItgH4ksGuBcys0prQk2boVvNw8DaTmPkoEQ0DsJWI3ta0vVtUPqMn34MH83gDhZoqf6ZCx7EcgLEez0OFIL4P85OHRK7t5uLgbfyQ7bDItep5Ii3bkVAIzDzZ5VtYNUP3ZB5EpNRUnZmPan9PENELKJihXigElc/oyffgdh7M/EnR881E9E9EdAkzB10cLIjvw5SwyLXdXCos22GRa6XzRF22k+Iy6m2GTkRVKDRD3+jxOfSm6oBJU3UiGkREF+vPAdwCwDBKoAxUPuNGAH+tRSBMBPCx7grwENt5ENGXiIi059ejIIcfejwPFYL4PvwkCLkGKivbYZFrpblEXrb9Xg33+wHgDhS04RkAxwFs0baPALC5aNx0AH9AIUpggQ/z+CKAVwAc0P4OKZ0HCtEJe7THPq/nYfQZATwI4EHtOQH4mba/HSbRKwHM4yHt8+8BsBPAv/NpHk8BOAagW5ORv6nE9xFluQ6DbIdFrpMg21K6QhAEQQCQHJeRIAiCYIMoBEEQBAGAKARBEARBQxSCIAiCAEAUgiAIgqAhCkEQBEEAIApBEARB0Pj/Du8eZ8lH7LYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def runif_shape(n, d):\n",
    "    M = randn(n, d)\n",
    "    M2 = M*M\n",
    "    L = np.sum(M2, axis = 1)\n",
    "    L = np.sqrt(L)\n",
    "    D = np.diag(1 / L)\n",
    "    U = D @ M\n",
    "    return(U)\n",
    "\n",
    "plt.subplot(1, 2, 1)\n",
    "X = runif_shape(200, 2)\n",
    "plt.scatter(X[:,0], X[:,1], marker='o')\n",
    "\n",
    "plt.subplot(1, 2, 2)\n",
    "X = runif_shape(2000, 2)\n",
    "plt.scatter(X[:,0], X[:,1], marker='o')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Porównanie generatorów dla wielowymiarowego rozkładu normalnego"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Czas dla metody opartej o wartości własne: 0:00:04.348477\n",
      "Czas dla metody svd                      : 0:00:03.518608\n",
      "Czas dla metody Choleskiego              : 0:00:02.761044\n",
      "Czas potrzebny do estymacji macierzy kowariancji: 0:00:02.213640\n",
      "Czas dla metody z numpy                  : 0:00:05.413354\n"
     ]
    }
   ],
   "source": [
    "from datetime import datetime\n",
    "from numpy.random import multivariate_normal\n",
    "\n",
    "n = 100\n",
    "d = 30\n",
    "N = 2000\n",
    "\n",
    "mu = np.zeros(d)\n",
    "\n",
    "np.random.seed(100)\n",
    "now1 = datetime.now()\n",
    "for i in range(N):\n",
    "    C = randn(n, d)\n",
    "    Sigma = np.cov(C.T)\n",
    "    rmvn_eigen(n, mu, Sigma)\n",
    "end1 = datetime.now()\n",
    "\n",
    "np.random.seed(100)\n",
    "now2 = datetime.now()\n",
    "for i in range(N):\n",
    "    C = randn(n, d)\n",
    "    Sigma = np.cov(C.T)\n",
    "    rmvn_svd(n, mu, Sigma)\n",
    "end2 = datetime.now()\n",
    "\n",
    "np.random.seed(100)\n",
    "now3 = datetime.now()\n",
    "for i in range(N):\n",
    "    C = randn(n, d)\n",
    "    Sigma = np.cov(C.T)\n",
    "    rmvn_chol(n, mu, Sigma)\n",
    "end3 = datetime.now()\n",
    "\n",
    "time1 = end1 - now1\n",
    "time2 = end2 - now2\n",
    "time3 = end3 - now3\n",
    "\n",
    "print(f'Czas dla metody opartej o wartości własne: {time1}') \n",
    "print(f'Czas dla metody svd                      : {time2}')\n",
    "print(f'Czas dla metody Choleskiego              : {time3}')\n",
    "\n",
    "np.random.seed(100)\n",
    "now4 = datetime.now()\n",
    "for i in range(N):\n",
    "    C = randn(n, d)\n",
    "    Sigma = np.cov(C.T)\n",
    "end4 = datetime.now()\n",
    "\n",
    "time4 = end4 - now4\n",
    "\n",
    "print(f'Czas potrzebny do estymacji macierzy kowariancji: {time4}')\n",
    "\n",
    "np.random.seed(100)\n",
    "now5 = datetime.now()\n",
    "for i in range(N):\n",
    "    C = randn(n, d)\n",
    "    Sigma = np.cov(C.T)\n",
    "    multivariate_normal(mu, Sigma, n)\n",
    "end5 = datetime.now()\n",
    "\n",
    "time5 = end5 - now5\n",
    "\n",
    "print(f'Czas dla metody z numpy                  : {time5}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sprawdzenie generatora z numpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x6615538>"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3df2xd5Zkn8O9zb07IdShzQ/GoxcSEbauwEwKJ8Jas/McMaUuYEoILbVkGZqvtaqORtlLpUHeSEkHo0iWraIDRttJsdna0uyLL8CvchtJOACWj0WQnqA43xnhIZtspvy6M6i6YaYmXOPazf9jHuT4+5/4659z3fc/5fqSosX1z7lv8nue+53nf93lFVUFERO4qmG4AERHFw0BOROQ4BnIiIscxkBMROY6BnIjIcctMvOlFF12ka9asMfHWlAPHjx//par2mnhv9m1KU1TfNhLI16xZg5GRERNvTTkgIq+bem/2bUpTVN9maoWIyHEM5EREjmMgJyJyHAM5EZHjGMiJiBxnZNUKUacq1Rr2HjqFtyencHG5hOEtazG0sc90s4gaSrvfMpCTMyrVGnYeGMPU9AwAoDY5hZ0HxgCAwZys1Y1+y9QKOWPvoVMLN4NvanoGew+dMtQioua60W8ZyMkZb09OtfV9Iht0o98ykJMzLi6X2vo+kQ260W8ZyMkZw1vWouQVF32v5BUxvGWtoRYRNdeNfsvJTnKGPzHU7VUrIlIEMAKgpqpbU30zypxu9FsGcnLK0MY+EytUvg7gVQAXdPuNKRvS7rdMrRA1ICKXALgBwJ+ZbgtRFAZyosYeBvAtALNRLxCR7SIyIiIjExMT3WsZ0TwGcqIIIrIVwC9U9Xij16nqPlUdUNWB3l4j51lQzjFHThRtEMA2Efk8gBUALhCRR1T1DsPtIseFbdkHOp8QZSAniqCqOwHsBAAR+R0A32QQpyit1lMJ27I//MQoIMD0jC58r51t/AzkRERt8oN2bXIKRRHMqEIA6PzPwwJx/b8Jmp7VJd/zt/EzkBMlRFX/CsBfGW4GWSA4op7RuSAcDMX1gXhXZQz7j72x5DXNtLqNP3YgF5EVAP4awHnz13tSVe+Ne10iIps0GlFHqU1OYd09f4kPzsw0f3GIVrfxJzEi/xDAZlX9tYh4AP5GRH6sqscSuDYRkXHBUXg7Og3i7Wzjjx3IVVUB/Hr+S2/+T7tPEERE1gorRZu2B25e391VK/O1KI4D+CSA76vqi0lcl+zDE3ooL+r7uomRaTv3VSKBXFVnAGwQkTKAp0XkClV9pf41IrIdwHYA6O/vT+Jtqct4Qg9lTdTAJE4qpVUlr4Cp6fANw31tlrhNdGenqk5ibmb/+pCfcfeb43hCD2WJH6xr8yNuf2DiB/c0g/jK5UVMTc9CQn7WSYnb2IFcRHrnR+IQkRKAzwI4Gfe6ZB+e0ENZEjUwuevx0bZWpnTCnwANS9nccnX7lRKTGJF/HMAREXkZwE8APK+qP0zgumQZntBDWRI1APHXhZty5GT7hddiB3JVfVlVN6rqlap6hap+J+41yU48oYeyxNYBSCdPuKx+SC0b2tiHB25ej75yCYK5CZl2lkgR2SRsYGKDTj5guEWf2mLohB6ixNUfwZZ2TrxVXkE6esJlICeiXKpUa9h9cByTU9OmmwIAEAB7v3RVRwMlBnIiyp1KtYbhJ0ZDqw6aUPKKsdKUDORElDt7D52yJoj3JbBDmoGciHJlV2XMipy4V5COUylBDORElFnBLfhrPlrC0Z+9a7pZAIDzVyxLbOEAAzkRZVJYbSAbRuK+ydPJTbIykFMsrIZItjJRerYdwfXice4lBnLqGKshks1sGn0HBXdEx72XuLOTOsZqiGSzooTVFjSvIEsPjYh7LzGQU8eiakLUJqdw2Y5nMbjnMCrVWpdbRTTHdPGrMF5R8OCXNywZZcetLMpATh1rVBMiWN+ZqNtW9Ximm7DErf9idWiqJG5lUQZy6lgrRYeYaiETKtVaoqtCkvK/Xnwj9Gk1bmVRTnZSTM0fX3nwBHWDv+rD5klOfzNpcDKzvoAXV61Q6upvFkErYdzeus/kPheCdxT/adUP1nEqizKQZ0Aaa7nDrglg0RKpVoI4D56gtHTjgOS0JfW0ykDuuDTWckddc4VXaOumSaIYEFEU2zf8tCKpp1UGcsc1Wn/aaQCNumarN03ckpxErXApneIVAIhgeubcc2yST6sM5I7rZP1ps1RMnMe9VT0e7r1xHYM4pa4g5yYPbbesWMQtV/fhyMmJVMpZMJA77uJyKXRkEvXI1igVA8yNxqPujXLJw4dnZxeNzP0JT6ZRqJsq1ZozQRyYe6I9cnICR3dsTuX6DOSOG96ydsmET/0jW3D0/cGHZ0PTJrsPji8J0vUEwNarPo6BSy9kkSwyzsW9CWkuw2Ugd1yj9adho+8ozc4tVABPHa9h4NILUxtVELXKpfy4L81luLEDuYisBvA/AXwMwCyAfar6J3GvS60LW39aqdZw1+OjidabmJqewX3PjHMETka5WvLh2st7U7t2EiPyswDuUtWXROQjAI6LyPOq+ncJXJs64I/E0yga9N7paVSqNQZz6jqXN/8AwJGTE6ldO3atFVV9R1Vfmv/7rwC8CoB3uUFpr691MT9JbvMHJ64GccChHLmIrAGwEcCLIT/bDmA7APT39yf5thSQdm0T1k6hbqifqC+IWFmWth1W58h9InI+gKcA3Kmq/xT8uaruA7APAAYGBtz+jVguakliktePwqPfKK5KtYa7nx7DB2fOPVW6HsTTLlWRSBlbEfEwF8T3q+qBJK5JnWulvGwcaz7aeI16bXKK9cipI5VqDcNPji4K4lmQ9k7n2IFcRATAfwPwqqo+GL9JFNfQxj48cPN69KX0KPe/f/ZuaHDO2tFvIrJaRI6IyKsiMi4iXzfdpqzbe+jUom3sWVAueak/lSaRWhkE8PsAxkTkxPz3vq2qP0rg2tQhf0nihvuea7pGvF2KcxOe9WmUqHSOwzl1rsjqkkq1hvueGcd7Fh4GEYdXEOzeti7194kdyFX1bzC38Y8sU6nWEg/iPj9tUr/ZKKo+uav1yFX1HQDvzP/9VyLir8hiIE+Qn07J2kgcAPZ+6aquzBHxqLcMSzulEUyjKJZ+omelHnmzFVkiMiIiIxMT6a0VzqosplN83Zro5xZ9x9WvEin3eFAF3p+aTn3lShS/gFaWVq1wRVa6HE69NVQude/wZwZyhwVrqdTnF9s5iq0TxYh1vX3lUqZqsXBFVvpMDTrS1o3cuI+pFYc128GZVhD3ioLbrlkd69RvF3BFVnekWYPElMFPXNjVJ1EGcoeZeiRduXwZ7h9av7DEUTD3GLnCK+Abj53A4J7DWVk77q/I2iwiJ+b/fN50o7JkV2UM+4+9YboZiSmK4I5N/dj/7/5lV9+XqRWHmXokfX9+JYy/xDGNc0NtwBVZ6alUa9h9cDy1VVWm3HbNatw/tL7r78sRucPS3sEZJbicMGsbgShd/gd/1oI4ADz64ptG3pcjcocFD5XwV62keYN4RVmSB+/k3FDKr7Src5pkqiYMA7lFOik4FXaoxLp7/jKVWhVRByu3e24o5VuWP+CLYiYTx9SKJZIsOJVWEK/ec13oB0tYiidrK1goOVn+gL/tmtVG3pcjcks0yjP7wdNkidjJBjUwGp0bShTst2s+OrdhLGs7p+7Y1G9kohNgILdGszxzOytDRICkU3UXl0sNP0jCUjxEYf02i5t/BDAWxAGmVqwR9bjpfz9qxH5nyLrtpIO4YO4G/MZjJ1hrnNqS5YnNeqbTRQzklmiWZ240QRQMqklPt2jgf31cYkjNZHli0xe2kqvbmFqxRLM8c7PNP1PTM9h9cBxAelvzw+ThRqXOZbWOim/l8iK++4V0T/9pBQO5RYJ55kq11tbBEJNT09h54OW0mhfK9CMl2atSreHdDz403YxUjX/netNNAMBAbq1KtYbhJ0YxPdve+HpqejalFi3FJYYUZVdlDI9kqIZKmLSOUuwEc+SW2nvoVNtBvJvKJS/1A2XJTXkI4l7BfF68HgO5pWzPPX/w4VnTTSAL5SGIA8D5K5ZZNYhhILeU7bnn6VldmFwlAubSgVkqSdtIow1yJjBHbpHgsW0FAN3LeLcvi9XrqHN7D53K3G7NKLYNtBjILVCp1nDfM+OLjmrz/257MCfy2Z4OTIqNk/xMrRjmb2F+L+JRzeYgvqqne4fLkv1sG6WmYeXyopWT/ImMyEXkzwFsBfALVb0iiWvmhatbmL2i4N4bu3e4LNnHTwWmfdC3Tco9y60L4kByI/L/DsCOlfGOcfFxtK9cwt4vXmVlh6buqC+7DOQjiAP23q+JjMhV9a9FZE0S18obW7Ywf+o3V+L0mdmm5UX7yiUc3bG5a+0iO7n6JBmXremjruXIRWS7iIyIyMjExES33tZalWoNg3sOLzyWmvbLX59Z+HvUKScCWDfJQ2bYOjJNUwH29v+urVpR1X0A9gHAwMBAXp7EFonKKSpgPMf43unphQnXsHMHBcDtm/qZTiEA9jxJdkvJK+CBm6+0tv9z+WGXBAvsB0OljZ9sRRHMqvLEH1pieMta3PnYCdPN6IrX9txguglNMZCnJHiazgcfnnUupzirip870Imp+0Zef9d0E7rC1GHK7UokRy4ijwL4WwBrReQtEfm3SVzXVWEHKbu4C9LWiR0yb/+L+diKH5ZmtFFSq1ZuS+I6WZGFGX0bd6+ReZVqDd8+8HLixwnayqZStY0wtZIC12f0+5gTpxCVag1/+PgJWFxdOVE2HOHWKgbyFETN6Ptb2qO249tAAK4Tp1C7D47nJoiv6vFw743rnBnMsNZKCqIOUr73xnWo3nMdHr51g6GWNVdm/RQKUanWnJzn6US55KF6z3XOBHGAgTwVQxv78MDN69FXLkEwl6qoL7RjcwfJS+6TWucfO5gHBQF2b3OvhhBTKykJHqRcr1Ktdbk1rXs/J6Muap3txw4m6YIVntUDrSgM5CkKriX3J052Hhgz3LJoCmBwz2FOdtIC1yfv2+HqQIaBPCXBnZy1ySnsPDCGFV7B+qWJflsBu9NAlC7/wJN8jMXnuLp3gjnylIStJZ+anrF6xUq9qekZ7D10ynQzyJBKtYbhJ0ed6a9JcHnvBEfkKcnC42gW/j9Q+yrVGu56fNSZXY1JcH3vBAN5SrJQHc7Vx0zqnD8Sz0sQL3l2Ht3WLqZWUnLt5b2mmxCLy4+Z1Ln7nhnH9Ew+grgIMhHEAQby1Bw56dbhGat6vMh173kmIteLyCkR+amI7DDdnrTlKScOzc5kPlMrLQpbStioE7iWX548PY3qPdeZboZVRKQI4PsAPgfgLQA/EZGDqvp3ZltGSchS6pAj8haElaXdeWCs4cYe2zpJX7nUsJKbbe21xKcB/FRV/0FVzwD4CwA3GW4TJSBrxxZyRN6CqKWE/vK8sE0/p8+c7Xo7o9Tnu4efGF2yS8+lKm9d1gfgzbqv3wJwTfBFIrIdwHYA6O/v707LErKrMoZHX3wTM6rOHKKQhKwdW8hA3oKoNEltcgrDT44uTA7VJqfwh4+fQLEg1kwY9XgF/MdAvnv3wfGFAkiuVXnrsrDItuQX6+p5tLsqY3jk2LkDIvKyUqVc8nD/0HrTzUgUA3kLopYSimBJwJ5VYNaSIA4Ap6dnF54c/PovDNotewvA6rqvLwHwtqG2JC4vp/wEuboNvxHmyFsQVZbWlQFMKzl9CvUTAJ8SkctEZDmAfwXgoOE2JcaV/pu0LM4HMZC3IKosrUu45b59qnoWwNcAHALwKoDHVXXcbKsojqzuj2BqpUVhKYn6XLMLXN9paoKq/gjAj0y3I2l5fDorimR2fwRH5DHs3rYOXsGdmf48rUqgaP5y2ryZVc1kEAc4Io/F7xSuFBhyoY2UvrDltHmQxdy4jyPyBFxQcufzcHDP4Vw+VtM5ru06TkLW90okEsjzVo/C5z+iulSfgitYKMsj0zCrejzs/eJVmU2rAAmkVvJcj8LVR1R/BUuWOzZFG96yFnc+dsJ0M1KVlfK0rUoiJ7BQjwIARMSvR5H5QO7yI6rLbafO+IXf8rB66Zar87XxLYlAnvl6FL5gBcRyj2d1WqUoggtKy0LbmLfH67wLbsfPOtfKSMeVRI685XoUqjqgqgO9ve4duhBWAfHX/8+ewlhBJa+IP/7yVbj3xnWhu1KzPPFDi+UtiAP5e+JMYkSe6XoUvrB8eLCKoC3Czh9sp5Y6uSmsZj4A7M9ZEAfy98SZRCBfqEcBoIa5ehS/l8B1reLKJ/xre25Y8j0Wyso+/4nRH2z4q5POW1ZY+niccXl84owdyFX1rIj49SiKAP48a/UoKtUaCiLcUEPWiqqZ7+KqqjiyvA2/kUR2smS1HgVwbqQTFsS9gliVXimXPNNNIENceWJMk1eUzK8Xj8KdnU00Wiu+fJk9//m8gmD3tnWmm0GG5C0nHJSHTT+NuLO33JBGI50Pztjx2Bo2uUn5Mrxl7aIceV6ct6yAU/f/rulmGMdAPi9sxn9oY1/k6UC2EABHd2w23QwyzP8Q9/uwPQm/dH14dtZ0E6xgT27AoLA14n49Ettnv/P+SE3nDG3sw9Edm/HQrRtCN3dQdnFEjugZ//ueGUfPcnv/E+VxmRWFq3+iLIjkZkTOD6w59kapLorKg793etraLfj+MitgrjQtN/vkV3ANeZ6Wyd6+yc1yH0ljagVupidm52/WqJQQ5YerVTjjKIrgjk39uH/IrbNz08IROdyc8b+4XIpMCbFEbb7YPBmfpKII/vjL+V1i2AhH5JibJHrg5vXomx+Z2362pZ8bj0oJcXNIfuTp6SvLZ27GxUA+b2hjH4a3rEXJK1qXY1zV46GvXIJgbs24vwU5KiXkYqqI2pe3Q5TZr6MxkNexNdfoT7g+dOsGHN2xeWFU4n/w1ONKlvywtb+mgf26MebI6zRKSRQLghmDdVX8iUzg3OaP4CYQrlrJl7yk0LhzuTkG8jqNdnGaDOK+sIlMlqjNl+B6cdvSgEl7+NYN7N8tYGqlTliqwjZ5GYXRUsEdyFkP4nds6mcQbxFH5HX8TmPzCeOc8MmvPOXEyyWPa8TbwBF5gM0jAE745FtensZKXpElmduUyxF5VKVD25S8Ai5ceZ717aTusL0SZxJKXiGXJ/zElbtAHnW2IXBuNC6AFUWHbrn6Ej5e0oJrL+/F/mNvWNE302P3Zjxb5S61ErWt/a7HR3HZjmex8TvPGWrZUj8cfcd0E8gSlWoNTx2vZTyIn1uZRe3JXSCPyjPOqEIxt/nGlptlcmo6V1uwKVqeJjrzMheQpNwFctdWfXB0QkC+gptr96gNchfIXVgrXi9PNzBFy0tw48qszsQK5CLyJREZF5FZERlIqlFpqq90KLC/0mFebmBqzLUBSDvCCsJRe+KuWnkFwM0A/ksCbema+m3twVUsNvEKwtEJATi3omr4iROYztB5w+WSx8PDExArkKvqqwAglo9qmzlvWcHKQH7+imUcndAiy4pFTM/a11c75XjosEbX1pGLyHYA2wGgv9/sOXv+hqDa5JQ1a8bDTFp6XmgeiMheADcCOAPgZwD+japOmmxTFleusI8no2mOXEReEJFXQv7c1M4bqeo+VR1Q1YHe3t7OWxxTfeEhwN4gDjA/btjzAK5Q1SsB/D2AnYbbk8mJb/bxZDQdkavqZ7vRkG5xZVTD2XuzVLV+Z9gxAF801Rbfb5Q8TE5lZwTLPp6cXC0/rFRrTtSq4Oy9db4K4MdRPxSR7SIyIiIjExMTqTSgUq05H8Sjjiyk+GLlyEXkCwD+M4BeAM+KyAlV3ZJIyxLmwvmGgrnj3Ni5u0NEXgDwsZAf3a2qP5h/zd0AzgLYH3UdVd0HYB8ADAwMJJqt21UZw6MvvpmJ2uOTp6dRvec6083IpLirVp4G8HRCbUmV7SkVAXA7C+l3VbO0oYh8BcBWAJ9R7X4k3VUZwyPH3uj226aG+fD05KL6oe0plaIIbrtmNSsdWkRErgfwRwB+W1VPm2jDoy++aeJtU+EVuSciTZnPkbuQUplRxVPHayyQZZfvAfgIgOdF5ISI/Gm3G5CFdAowlxvf+8Wr+LSZosyPyHcfHLc6peILO1iZzFHVT5puQ9Hxw5U559M9mR6RuzbTn8V1wtS5Tf9slekmdIxzPt2V6RG5LSVgS14Rt1zdhyMnJ/D25BQKESMtTgZRvdf+r5sf7H08lrDrMh3IbRnhBtfL7qqMLTmyK2xzhCtni1I8Yb9nAFZP0EdhESwzMh3IbTisNlgmN+zILgFwy9V9i4J0K2eLkvvCfs/DT45iZsbN3DiLYJmR6Ry5DTWcZ1Sx88DYwoqUsPXsCuDIycU7AqPOFrUlXUTJCPs9T88oXK1UyyJYZmQ6kNcfImFSfQCOSvcEv9/q68htWft9cp7HjEwHcmAumB/dsRmv7bnBaDv8Gzaqowe/3+rryG0u/z4LgTQKi2CZk/lAXm9Vj2fsvf0bNizdE3YDtPo6cpsN6b9ODH7iQjz45Q0sgmWJTE921qtUazC5t8IPwH5Hb7YapdXXkdv83+edj50w3JL2fGmgf9GRiWRWLgL57f/1b3H0Z++m/j5RO/FW9XiLOnyrNwBvlPyw+aSqMNyFbBfnA3mztdbdCuL+pp+njtcWrUIoeUXce+O61N+f7NTKXoC9h045FcSB7E3Sus7pQN5srXWlWks1iPujqPqdbAOXXsh0CAFofS+Ai0HR5UnaLHI6kDdaaz20sS/VNddR25CZDiFfs/7ps2HjWiMFYNG6dk6628fpVSvN1lqnMdIpCPDanhtwdMdmBmxqqNW9AMNb1sILruWzyIO3cnWK7ZwekUeNZC4ul1Kr7T3rWjKTjGnWP/0U3G+UPKt3co68/i7rp1jO6RF51Frray/vxc4DY6lMIJneJUruaNY/a5NTUACTU9OYsXiE8MixN7CrYvfhLHnndCCv34Jf/9h35OREKodJ+Dfh4J7DuGzHsxjcc5in+lCkbvfPNGXp2LksEgNnymJgYEBHRkZSu/5lO55NfDReLnmYnpnFB2cW34BhK1fILBE5rqoDJt67lb6dRv/sBtNlLii6bzs9Ig+qVGsY3HM40Zuk5BVxx6Z+fHh2aRAHzm3i8JeWcYROUdLon90SLMdMdslMIN9VGcM3HjuR6DKuokhbj8IsM0tR/DXlNi8zbOS2a1abbgI1ECuQi8heETkpIi+LyNMiUk6qYe2oVGtLTtyp1+kE5awqhjb2tbWM0cXNHZS+sDXltvCKgnLJW8jjD37iwoUReFEEd2zqx/1D6802khqKu/zweQA7VfWsiPwnADsB/FH8ZrWn0RZnAXB0x2YM7jnc9mjI373WzoYN7nijMLZ+wHNuJxtijchV9TlVPTv/5TEAl8RvUvsa3STl+dK17e5Ek7p/M7xlLVrJEHLHG0XhBzylKckc+VcB/DjqhyKyXURGRGRkYmIi6mUdaXST+Ityhjb2oVxqrR65ALh9U//CKGVoY1/DCSrueKNmbK07zkn6bGgayEXkBRF5JeTPTXWvuRvAWQD7o66jqvtUdUBVB3p7e5Np/bxGo+D3p86dIbj1qo83vVZfuYTbN/XjyMmJRWvFo/LsfeUSfs4t+9RE8NhBm1aBcJLefU1z5Kr62UY/F5GvANgK4DNqYlE65m6S3QfHMTm19ODX+tF68IDjeiWviAdunpvQCatYF1WilqkUapX/QV/fv2xhaw6fWhN31cr1mJvc3Kaqp5NpUmd2b1sXuR3a34nZaMLST4tEVaw7cnIidJceR+HUDlOrV8olDw/PF78Kwxy+2+KuWvkegPMAPC9zj4rHVPUPYreqA2FHo117ee+SUXSYvnJp4d83qljHErUUl4mRb7nk4cS91y18HXwi4JOl+2IFclX9ZFINaeUklWaCgXZwz+GmQTzYiRtVrCOKy0Tt8d3bzp1QxbNgs8mKMratnqQS9W+jOmWz0U9Blk70nD5zdsnrOGKhpAxvWdvVHHnwvFiAh59kkRWBvNWTVIIafQAAQCHiMGSfXzm0NjmF4SdHAQWmA+VEyyUPu7etY8enRPj96K7HRxv2zU6UvCLPi80pKwJ5qyepBEV9ANz99BhmFW3dKNMz4a9ded4yBnFKjP8EmXQQB4AVXgHnLSvg/alppkxyxopA3mleOirQh1Up7BSXZVFSdlXGGtYEaoVfNjnMe6enUfKKeOjWDQzgOWNF9cOok1Sa5aW7MQHJSU5KQqPCbkURCJpvEvKDdKNlhNzck09WBPKok1SajSrSnoDkJCclpVFht1lVPHTrBnxkReMH5Pp5o6M7NkfW/+FTZP5YkVoBOptJH9rYhzsfO5FYG1b1eOhZvozLsihxjYJrz/JiyytZ6q/DpbLksyaQd2pVj4f3Ti/dmr9yebGtXLk/w8/ATWlotH68nX5aH6TDljLyKTKfrEitxHHvjevgFRc/ZHpFwXe/sL7lAyX8k4AYxClIRL4pIioiF8W5ThLVD4NButOUJGWP8yPyZjvVvvHYiaarBPyTgIjqichqAJ8D8Ebca/n9q9NUYNQBENzcQ0AGAjkQ3ZmHNvZh5PV3my75Yk6RIjwE4FsAfpDExfyibO1s0RcAP+fp9dSE86mVZu4fWo+H6pZrBWf6mVOkMCKyDUBNVUdbeG1Lh6ZUqjV88GF4CYiVy8PTLhxkUCusH5HX11Ip93hQRds71+pH7EkU56JsEJEXAHws5Ed3A/g2gOtCfraEqu4DsA8ABgYGQh/+guUkfKt6vIVt9MGfC+bKRwzuOcx+Sg1ZHciDnb9+dUo7hbXqMadIvqhDU0RkPYDLAIzOl2e+BMBLIvJpVf3HTt4rqg55z/LFJSD81Ev9Ds5O+zrlh9WBvFkR/vpdbBxlU1JUdQzAb/pfi8hrAAZU9ZedXrOVekL+IGNwz+ElefRWishRflkdyFvZoeaPVhqVwGU6hUxrtnmnvo9GTcxzxyZFsXqys5WJnqJIZAlc4Fx6pjZ/g/DUcGqXqq6JMxoHGtcTCvbRKJz4pChWB/JmmyhKXjGyHKg/emlU65yoWxpt3mnlHE+urqJGrE6tBDf7hIQZSDAAAARVSURBVK1aiVqX649eOq11TpS0qIn2Rn1RAKYDqSmrAznQ2iqTRvUmWFiIbBfVR/vKJRzdsdlAi8g1VqdWWtGs3kSntc6JuoV9lOKyfkTeikajdp4aTrZjH6W4MhHIm+EmILId+yjFESu1IiL/QUReFpETIvKciFycVMOIiKg1cXPke1X1SlXdAOCHAO5JoE1ERNSGWIFcVf+p7suViD7gm4iIUhI7Ry4i3wXwrwG8D+DaBq/bDmA7APT398d9WyIimtd0RC4iL4jIKyF/bgIAVb1bVVcD2A/ga1HXUdV9qjqgqgO9vb3J/T8gIso50Ygt7m1fSORSAM+q6hUtvHYCwOuJvHF7LgIQq2ZGimxuG2B3+4Jtu1RVjYwWWujbtvx3tKEdNrQBcKsdoX07VmpFRD6lqv9n/sttAE628u8M3mQjqjpg4r2bsbltgN3ts6ltzfq2LW21oR02tCEr7YibI98jImsBzGJuFPIHMa9HRERtihXIVfWWpBpCRESdcb7WSpv2mW5AAza3DbC7fTa3LciWttrQDhvaAGSgHYlNdhIRkRl5G5ETEWUOAzkRkeNyF8hFZK+InJwv9vW0iJQtaNP1InJKRH4qIjtMt8cnIqtF5IiIvCoi4yLyddNtChKRoohUReSHptvSLhH5poioiFxk6P2N3Qs29Hmb+nfcfpy7QA7geQBXqOqVAP4ewE6TjRGRIoDvA/hdAL8F4DYR+S2TbapzFsBdqvrPAWwC8O8tapvv6wBeNd2IdonIagCfA/CGwWYYuRcs6vM29e9Y/Th3gVxVn1PVs/NfHgNwicn2APg0gJ+q6j+o6hkAfwHgJsNtAgCo6juq+tL833+FuY5mTdFsEbkEwA0A/sx0WzrwEIBvwWChOYP3ghV93pb+nUQ/zl0gD/gqgB8bbkMfgDfrvn4LFgVLn4isAbARwItmW7LIw5gLhrOmG9IOEdkGoKaqo6bbUqeb94J1fd5w/47djzN5QpCIvADgYyE/ultVfzD/mrsx92i1v5ttCyEh37NqTaiInA/gKQB3BkoXGyMiWwH8QlWPi8jvmG5PUKM+CODbAK4z3Q6D94JVfd5k/06qH2cykKvqZxv9XES+AmArgM+o+YX0bwFYXff1JQDeNtSWJUTEw1wn36+qB0y3p84ggG0i8nkAKwBcICKPqOodhtsFILoPish6AJcBGBURYO73/ZKIfFpV/7Fb7ahrj4l7wZo+b0H/TqQf525DkIhcD+BBAL+tqhMWtGcZ5iaaPgOgBuAnAH5PVceNNgyAzEWa/wHgXVW903R7osyPZL6pqltNt6VdIvIagAFV7Xr1PVP3gi193rb+Hacf5zFH/j0AHwHw/PxZo39qsjHzk01fA3AIc5Mtj9sQxOcNAvh9AJvn/1udmB85UDYYuRcs6vOZ6d+5G5ETEWVNHkfkRESZwkBOROQ4BnIiIscxkBMROY6BnIjIcQzkRESOYyAnInLc/wfAQL+KcLHojwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n = 1000\n",
    "\n",
    "Sigma = np.array([[1,0.9],[0.9,1]])\n",
    "mu = np.zeros(2)\n",
    "\n",
    "X = multivariate_normal(mu, Sigma, n)\n",
    "plt.subplot(1, 2, 1)\n",
    "plt.scatter(X[:,0], X[:,1], marker='o')\n",
    "\n",
    "X = multivariate_normal(mu, Sigma, n*100)\n",
    "plt.subplot(1, 2, 2)\n",
    "plt.scatter(X[:,0], X[:,1], marker='o')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}