{ "metadata": { "name": "Gauss_Markov" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "Introduction\n", "-------------------\n", "\n", "In this section, we consider the famous Gauss-Markov problem which will give us an opportunity to use all the material we have so far developed. The Gauss-Markov is the fundamental model for noisy parameter estimation because it estimates unobservable parameters given a noisy indirect measurement. Incarnations of the same model appear in all studies of Gaussian models. This case is an excellent opportunity to use everything we have so far learned about projection and conditional expectation." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Following Luenberger (1997), let's consider the following problem:\n", "\n", "$$\\mathbf{y} = \\mathbf{W} \\boldsymbol{\\beta} + \\boldsymbol{\\epsilon}$$\n", "\n", "where $\\mathbf{W}$ is a $n \\times m$ matrix, and $\\mathbf{y}$ is a $n \\times 1$ vector. Also, $\\boldsymbol{\\epsilon}$ is a $n$-dimensional random vector with zero-mean and covariance\n", "\n", "$$\\mathbb{E}( \\boldsymbol{\\epsilon} \\boldsymbol{\\epsilon}^T) = \\mathbf{Q}$$\n", "\n", "Note that real systems usually provide a *calibration mode* where you can estimate $\\mathbf{Q}$ so it's not fantastical to assume you have some knowledge of the noise statistics. The problem is to find a matrix $\\mathbf{K}$ so that $\\boldsymbol{\\hat{\\beta}} = \\mathbf{K} \\mathbf{y}$ approximates $\\boldsymbol{\\beta}$. Note that we only have knowledge of $\\boldsymbol{\\beta}$ via $\\mathbf{y}$ so we can't measure it directly. Further, note that $\\mathbf{K}$ is a matrix, not a vector, so there are $m \\times n$ entries to compute. \n", "\n", "We can approach this problem the usual way by trying to solve the MMSE problem:\n", "\n", "$$\\min_K \\mathbb{E}(|| \\boldsymbol{\\hat{\\beta}}- \\boldsymbol{\\beta} ||^2)$$\n", "\n", "which we can write out as\n", "\n", "$$\\min_K \\mathbb{E}(|| \\boldsymbol{\\hat{\\beta}}- \\boldsymbol{\\beta} ||^2)\n", "= \\min_K\\mathbb{E}(|| \\mathbf{K}\\mathbf{y}- \\boldsymbol{\\beta} ||^2)\n", "= \\min_K\\mathbb{E}(|| \\mathbf{K}\\mathbf{W}\\mathbf{\\boldsymbol{\\beta}}+\\mathbf{K}\\boldsymbol{\\epsilon}- \\boldsymbol{\\beta} ||^2)$$\n", "\n", "and since $\\boldsymbol{\\epsilon}$ is the only random variable here, this simplifies to\n", "\n", "$$\\min_K || \\mathbf{K}\\mathbf{W}\\mathbf{\\boldsymbol{\\beta}}- \\boldsymbol{\\beta} ||^2 + \\mathbb{E}(||\\mathbf{K}\\boldsymbol{\\epsilon} ||^2)\n", "$$\n", "\n", "The next step is to compute\n", "\n", "$\\DeclareMathOperator{\\Tr}{Trace}$\n", "$$\\mathbb{E}(||\\mathbf{K}\\boldsymbol{\\epsilon} ||^2) = \\mathbb{E}(\\boldsymbol{\\epsilon}^T \\mathbf{K}^T \\mathbf{K}^T \\boldsymbol{\\epsilon})=\\Tr(\\mathbf{K \\mathbb{E}(\\boldsymbol{\\epsilon}\\boldsymbol{\\epsilon}^T) K}^T)=\\Tr(\\mathbf{K Q K}^T)$$\n", "\n", "using the properties of the trace of a matrix. We can assemble everything as\n", "\n", "$$\\min_K || \\mathbf{K W} \\boldsymbol{\\beta} - \\boldsymbol{\\beta}||^2 + \\Tr(\\mathbf{K Q K}^T)$$\n", "\n", "Now, if we were to solve this for $\\mathbf{K}$, it would be a function of $\\boldsymbol{\\beta}$, which is the same thing as saying that the estimator, $\\boldsymbol{\\hat{\\beta}}$, is a function of what we are trying to estimate, $\\boldsymbol{\\beta}$, which makes no sense. However, writing this out tells us that if we had $\\mathbf{K W}= \\mathbf{I}$, then the first term vanishes and the problem simplifies to\n", "\n", "$$\\min_K \\Tr(\\mathbf{K Q K}^T)$$\n", "\n", "with\n", "\n", "$$\\mathbf{KW} = \\mathbf{I}$$\n", "\n", "This requirement is the same as asserting that the estimator is unbiased,\n", "\n", "$$\\mathbb{E}( \\boldsymbol{\\hat{\\beta}}) = \\mathbf{KW} \\boldsymbol{\\beta} = \\boldsymbol{\\beta}$$ \n", "\n", "To line this problem up with our earlier work, let's consider the $i^{th}$ column of $\\mathbf{K}$, $\\mathbf{k}_i$. Now, we can re-write the problem as\n", "\n", "$$\\min_k (\\mathbf{k}_i^T \\mathbf{Q} \\mathbf{k}_i)$$\n", "\n", "with\n", "\n", "$$\\mathbf{k}_i^T \\mathbf{W} = \\mathbf{e}_i$$\n", "\n", "and from our previous work on contrained optimization, we know the solution to this:\n", "\n", "$$\\mathbf{k}_i = \\mathbf{Q}^{-1} \\mathbf{W}(\\mathbf{W}^T \\mathbf{Q^{-1} W})^{-1}\\mathbf{e}_i$$\n", "\n", "Now all we have to do is stack these together for the general solution:\n", "\n", "$$\\mathbf{K} = (\\mathbf{W}^T \\mathbf{Q^{-1} W})^{-1} \\mathbf{W}^T\\mathbf{Q}^{-1}$$\n", "\n", "It's easy when you have all of the concepts lined up! For completeness, the covariance of the error is\n", "\n", "$$\\mathbb{E}(\\hat{\\boldsymbol{\\beta}}-\\boldsymbol{\\beta}) (\\hat{\\boldsymbol{\\beta}}-\\boldsymbol{\\beta})^T\n", "= \\mathbb{E}(\\mathbf{K} \\boldsymbol{\\epsilon} \\boldsymbol{\\epsilon}^T \\mathbf{K}^T)=\\mathbf{K}\\mathbf{Q}\\mathbf{K}^T =(\\mathbf{W}^T \\mathbf{Q}^{-1} \\mathbf{W})^{-1}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following simulation illustrates these results." ] }, { "cell_type": "code", "collapsed": false, "input": [ "from mpl_toolkits.mplot3d import proj3d\n", "from numpy.linalg import inv\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "from numpy import matrix, linalg, ones, array\n", "\n", "Q = np.eye(3)*.1 # error covariance matrix\n", "#Q[0,0]=1\n", "\n", "beta = matrix(ones((2,1))) # this is what we are trying estimate\n", "W = matrix([[1,2],\n", " [2,3],\n", " [1,1]])\n", "\n", "ntrials = 50 \n", "epsilon = np.random.multivariate_normal((0,0,0),Q,ntrials).T \n", "y=W*beta+epsilon\n", "\n", "K=inv(W.T*inv(Q)*W)*matrix(W.T)*inv(Q) \n", "b=K*y #estimated beta from data\n", "\n", "fig = plt.figure()\n", "fig.set_size_inches([6,6])\n", "\n", "# some convenience definitions for plotting\n", "bb = array(b)\n", "bm = bb.mean(1)\n", "yy = array(y)\n", "ax = fig.add_subplot(111, projection='3d')\n", "\n", "ax.plot3D(yy[0,:],yy[1,:],yy[2,:],'mo',label='y',alpha=0.3)\n", "ax.plot3D([beta[0,0],0],[beta[1,0],0],[0,0],'r-',label=r'$\\beta$')\n", "ax.plot3D([bm[0],0],[bm[1],0],[0,0],'g-',lw=1,label=r'$\\hat{\\beta}_m$')\n", "ax.plot3D(bb[0,:],bb[1,:],0*bb[1,:],'.g',alpha=0.5,lw=3,label=r'$\\hat{\\beta}$')\n", "ax.legend(loc=0,fontsize=18)\n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAAFdCAYAAACgiL63AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXuUHGWZ/7/V99tck7lkLkkwiWQSQhIIm+MSBFlYBRaM\nGyJwjCggRolyVBTUn7vCQeXoYQ+Lhz2aJYAHwRADZAlLQDZKBIIhCgRyJ7dJ5kJmksytr9V1+/0x\nvpXqmuru6q5r97yfczySmZ7ut7urvvXU8z7P92EkSZJAoVAoFFvwOL0ACoVCmUxQ0aVQKBQboaJL\noVAoNkJFl0KhUGyEii6FQqHYCBVdCoVCsREquhQKhWIjVHQpFArFRqjoUigUio1Q0aVQKBQboaJL\noVAoNkJFl0KhUGyEii6FQqHYCBVdCoVCsRGfkT9ubGzE8PCwWWuhqGhoaMDQ0JDTy6BQKCbCGPHT\nZRgG1I7XOujnS6FUHzS9QKFQKDZCRZdCoVBshIouhUKh2AgVXQqFQrERKroUCoViI1R0KRQKxUao\n6FIoFIqNUNGlUCgUG6GiS6FQKDZCRZdCoVBsxDbRPT1wGru27sL7f3gfu7buwumB0xXzGps2bYLH\n48G6des0fz9//nzMmTPHlNeiUCjVjS2ie3rgNLr/3I2OdAfas+3oSHeg+8/dpgqvla9x3XXXobW1\nFY8//viE3+3YsQP79+/HbbfdZvh1KhFJkpBKpZDJZCAIAvWKoFCKYIvo9u7uxczwzJyfzQzPRP/e\n/op4Da/Xi1tuuUUWWCWPPfYYfD4fvvzlLxt+nUpDkiRks1kIgoBUKoV4PI6xsTEkk0mwLEtFmELR\nwBbRZQRG8+cSZ94JafVr3H777WAYBo899pj8s2QyiQ0bNuCqq65Ca2urKa9TKQiCAJZlIYoiPB4P\nfD4fvF4vGIYBz/NIp9NUhCkUDWwRXcmrfaIxfm2hdONrzJw5E1dccQV++9vfgud5AMDvf/97JBIJ\nfOUrXzHlNSoBSZLAcRyy2SwAwOPJPYQYhoHH44HX65VFmOM4KsIUyt+xRXQ7FnSgO92d87PudDfa\n5rdV1Gt89atfxalTp7B582YA46mFadOm4ZprrjHtNdyMKIrIZrPgeR4Mw4Bhil/QGIaRBbiQCJM0\nBRVhSrVjaHKEXqa2TAUuBfr39kPiJDB+BjP/Yeb4zyvoNT772c+iubkZjz32GObPn4+33noL3//+\n9ydEe9WGJEkQBAEcxwGAbsHVgoiw8rk5jst5br/fL6crPB5P2a9FobgROjmiRL7//e/jwQcfxOc/\n/3k888wzOHToEGbNmmXJa7nh8yWiKAhCXrHleV7+vRmvJ4oiJEmSX4+KMKWaoKJbIkeOHMHHP/5x\nSJKEyy67DH/6058sey2nP1+STlAKoBaCIMgpB7MhIszzPCRJQiAQkEXY5/MZiropFCeo7vtiC5g1\naxY+9alPAUDV1uZKkgSe58GyLAA4Gl2SdARJ4TAMg2w2i2QyibGxMcTjcaRSKXAcJ0fIFIqbsSWn\nW20EAgE0NDTg+uuvd3oppkOaHUgVgtuiSGVOmAhsNpvNqabw+Xzw+/3yxp3b3gNlckNFt0QOHz6M\nP/zhD/jGN76BYDDo9HJMhWyWxeNx1NbWul6syPpKEeFq3/SkuB8qujp5++23sX//fvzyl79EKBTC\nXXfd5fSSTIOkE0opBXMjVIQplQAVXZ38+te/xpNPPolZs2bh6aefxvTp051ekimIoijnQytZcLWg\nIkxxI7R6wcVY+fkWqr0dGRlBTU1NTj1tIaysXiCQsrVQKGTac0qSJP+PQEWYYjU00p2EFEsnTJaL\nqfq9EwFWR8LqOmEKxQhUdCcZemtvJyP5RJhl2ZzyOSrCFCNQ0Z0kKNMJpBzM7OevNgHPJ8KpVAqC\nICAQCFARppQMFd1JALllrsbNMjshn53H44EgCPB4PBMiYa/XK+eDSccchaKEim6VU046odScLjGs\n8fl8rmyosAqtSFgURWQyGflnRISVfsOUyQ0V3SpFvVlmxW0vudXOZrPweDzIZDKQJEkWmMnWEUZF\nmKIHKrpViB3pBEEQkEgk4PF4UFNTA57n4fF4IIoiBEGAIAhyBYDST3cy5TypCFO0oKJbRSinOlgZ\nZbIsi1QqhXA4jGAwOKHOlWwukY0ns0S40svYqAhTACq6VQNJJ6TTaXAch5qamrKfK19Ol6QTyPP7\nfD755/meh6Q2iAiTSJi4mDEMMyEdke+53IJZlRpUhCcnVHSrAHUrrxWQdILX60VtbW1ZaQLl6B4A\nOSLMcRwymUzOfLXJJjJUhCcHVHQNMDg4iHXr1gEY7+cfHR3FL37xC/j9flteX93Ka1XlgDqdYNZr\nlCLC5PeTiXwinE6nc3wlqAhXFlR0y+TYsWPYuHEj7rnnHlkUrr/+ejz++ONYvXq15a+fb4yOmS28\nkiQhmUyC5/mcdIJVaIkwyQeTyRHpdDonHzyZREaZrgHOinA8Hs9J01ARdjeTZyvZRDiOw4svvoi7\n7747xxTm4MGDiEajlr++KIryCHMrNssYhoEgCBgbGwMA1NXVWS64+dbh8/kQDAYRDAZzcsMsyyKZ\nTCKdTk/aScJKo3ny3yQSTiQSGB0dRSKRAMuy8kWL4jw00i2DDRs2YNWqVRN+5vP5sHLlSste167a\nW1JpEI1GEQgEXBMtERFWbuCRSJhlWYiiOKEywi1rtwOtSFjZ+g0gx0GNRsLOQEW3DM6cOYPGxkY8\n+uijOHjwIHbu3ImjR4/igw8+sGyahN7aWyPpBZJOEEVRzt+6GbUIK2uEOY6b1I0aQH4R5nlefgwR\n4cnWTegk9qcXGMae/1lEb28v2tvbAQCtra3w+/1YtGgRMpkMHn74YUtek5RXWdnswPM8RkdHZSGr\nxCYGkn4IhUKIRqOIRCLwer0QBAHpdBqpVAqZTEau9KgW9JawEQFW3gmQzyYej2N0dBTJZFJOXdF0\nhDVQE/MSWb9+PT7zmc+goaEh5+ff+c53sHXrVnzwwQemvRaZfFvKGB2e55FMJlFXV6frNUh+NJ1O\nIxKJIBgMIplMwuv16jYMJyVrVkZJPM+D4ziEw+Gy/l7dqCEIAoDyGjWsMFQ3QiaTkTfQjKA2dVfe\nSdBI2DxoeqFE4vH4BMEFgNHRUUsuQKXOLSvlQiiKomxTWFtbq3tSRCViZqNGtaJVokYuduT3VISN\nQ0W3RMgBqGbHjh24/PLLTX89q27zeZ5HIpGA3++viMm/ZmOkUaMavYO1oCJsDVR0S6CnpwdHjx6d\n8PNt27ahu7sb3/3udx1YVWko0wmkOoFSmgi7LaVm13qoCJsDFd0S2L59O9LpNAYGBtDS0gIA6Ovr\nw+23346nnnoKM2bMcHiFhdMLoijK1QmF0gmTMVevJp8I8zwvbzKlUqmciRFOCowTr60lwsRwieM4\nBIPBCUM+qQhT0S2J0dFRPPjgg3jggQfkW8yTJ09i48aNWLRokdPLKwjHcUgmk/D7/YjFYqYf/NV+\ny60UYYZhIIoifD4frRFWQD4jSZKQyWTAMIxscE9+T0WYim5JsCyLSCSC+++/3+mlFEQZpZITIJPJ\nWJpOmGwnD23UyI9ySonybopEwkoRVs+XmwyfERVdnfT09KC1tdXpZRRFedCSdIIkSZZWJ0yGE6UY\neho1yNy0ydioAUBThLPZrFw5MllEmIquTrZt24bLLrvM6WXohqQTAoEAwuFwSQcvzekaR2nmDuSK\nsNkTNSo1tVOKCJMhn5X4PtVQ0dXJ8ePH0dnZ6fQyikLEMpFI0OoEF2HlRA23Ue5FoJgIezyenJxw\npYow7UhzMQzD5BhYF0MURSQSCfA8j7q6urLTCel0GpIkIRKJ6Ho8OTnc3JFmJmS6slneFESESWUE\ncY/T26iRSqUQDAZd09xC8tp6jx89EJ1Rtm9XqgjTSLdK4DgOiUQCwWBQHhJJqQyIWJC7EjpRYyJK\n03bgrAhns1n5TsHj8WBoaAhjY2M4//zznVmoDqjoVjjE2JtlWcRiMfj9fnkUerknJimJomhjdQ61\n1G45K8c0lYMdOeZ8IvzOO+9gz549VHQp1kDSCcC40TiJbt10AlKMU6hRg0R5mUym6nf9C0HeLzFu\ncjNUdCuUbDaLZDKJUCiEUCg06U6yyYxahImHBpko4nSNsJP7PFR0KaZD0gnZbFZOJ6ihG5yTCyLC\n5FhwQ6OGU0FAKpXSbWvqFFR0KwhBEJBMJsEwTNlj0PVQqmgTP91qqqWsJNQ51HzdcqQCRKtRY/jU\nMPr39oMRGUgeCW3z29DY3OjUWyqbdDqNtrY2p5dRECq6FYJb0wmkasLj8cjtndVSb1otFOqWy2az\nGDo1hJNvn8THIh8bj4I9DLrf7AaWoSzhdbJZI5VKuaKssBBUdF0OcbPiOC5vOkGNHekFtacDeT3l\n8EyWZSd9qZMbUTdqHDl6BB+LfGz8joUfv3C2e9vRs7sH9Z+qr6gLZyqVQiwWc3oZBaGi63LGxsbg\n8XgsTSeUitoi0uPxyM0ReupNJ/Muu9tgGAZeyQuvzwsvvPBJPkACBFGAwApIpVJyzliZjtBiaHAI\n/Xv7IbACRI+IGYtm2J6ioBtpFMMEg0EEg0FbxalQpKycOEEsIvM9VqvUyekNnmrDjDsayXP2ORiG\nARjA5/EhGA4iGo3qatQYGhzCiTdPYGZ4pjwv78SbJ8pOUZRLKpWioksxRjnDD61IL5BWX3JQl9MC\nW2yDB6D54HIxcrFqm9+G7je7MTM8U/5Zd7ob0y+crrtR4/iu45gemJ5z3M0Mz0Tfvj4quiqo6JbJ\nwYMHce655zq9DFuQJAnJZNL0AZbFNngAyLe0biqBkySpqi4Ijc2NwDKgb18fIADwAtMvnK4plvlE\nmBEZ8DwPSZQABvAwHogeERJv7/eWTqdpTrcaWbt2LdauXYs77rgDX/nKV5xejqUIgoBEIgGv12v5\nAEv1Bo8yoiIj01mWpZtyFtDY3FhWRCrne4M+BKWgfEcEjFe2ZLwZpNNp21JIlZDTrZ7LtU08+uij\niMfjePfdd3Hq1Ck89thjTi9pAkbTC+Tvs9ksxsbGEAyO5/bszit7vV7ZDzgYDMrRJSmfI00iZGYZ\nxTna5rehO909bsHIjOd8PxI/wjmLz5EvoizLWv69pdNpWjJWbXzyk5+U0wo/+MEPcODAAYdXZD4k\n15pKpVBTUyPf/jsJiZBILjnfptxkm8zgFgNzZYqClVj4Qj5MP+9siqJYowb57sz43tye+nH+bKow\n1HncuXPnOrQSaxBFEalUCgBcVaYG5O7U68kH0005eyEpinQ6LU98UFPoe0un0wAwoTLCDRcVM6Gi\nW4WUm14g3WV+vx+CIFSUUGlNZqBNGu7HzIka5O/dDhVdAwwODmLdunUAxvOMo6Oj+MUvfqGra8xN\nKLvLSO1tMpl0elllQ5s0nMXIuB4yliffxVNPo4bbv0squmVy7NgxbNy4Effcc49cPnP99dfj8ccf\nx+rVqx1enX6UE4OJJy/P804vy1SquUnDLTldKyhlosbQ0BBCoVDBSLinpwc333wzBgcHwTAMvvrV\nr+LOO+/Mecy2bdvw2c9+Fh/72McAACtWrMCPfvQjU98XFd0y4DgOL774Iu6+++6cnx88eBDLly93\naFVn0ZteIN1l5UwM1nrNSkGrSYPMJ6NNGu6lUKPGCy+8gPvuuw81NTX4/ve/j8svvxwXX3wxotGo\n/Pd+vx8PPfQQFi1ahEQigQsvvBBXXnklurq6cl7n0ksvxebNmy17H/RoKoMNGzZg1apVE37m8/mw\ncuVKh1alH5JOiMfjiEQiiEQiFSWaZkNGfYdCIUSjUYTDYXi9XrmCI5lMgmXZ8eL/CsgZOo1d0bey\nrPBrX/sa9uzZgxkzZiAUCuEnP/kJ/u3f/i3n8a2trVi0aBEAIBaLoaurC/39/ZrrtxIa6ZbBmTNn\n0NjYiEcffRQHDx7Ezp07cfToUXzwwQemTYi1Cj3dZZPdBL1QkwaZzkyEhW7KuQdRFNHa2op7770X\n9957b8FjuLu7G++99x6WLl2a83OGYfDWW29h4cKFaG9vx4MPPoh58+aZuk7bRZe5z54DVPqxNaLR\n29uL9vZ2AONXzqNHj2LRokXYt28fHn74Ydx3332WvG4p5BNNpVmN1d1l1YLWLS0pbcpmsxWdD7YK\npy7Yat+FfN9DIpHA9ddfj4cffnhCy/AFF1yAnp4eRCIRvPzyy1i+fDk+/PBDU9dpu+haJYZ28cYb\nb+Azn/kMAODaa6/FtddeC2C8+HvTpk2uEF0tWJY1ZFZDGYds7pA6VOWmHJnCTHbW7cgHG72VJ3aM\nWhMjCv1OyZG9R7D/1f1gBAaSV8L0ZdOxYMmCstdULnrMbjiOw4oVK7Bq1SrN/Zeamhr5v6+66irc\ncccdGBoaQmOjeaY9NL1QIvF4HA0NDRN+Pjo66spbcqUJupXdZdW8i14I5aZcMBgs2KRBOq7cwpG9\nR7Dv+X2Y4Z8BiZEwZfoUnBget2MEIFs1ErSmSRzZewS7f7sbF8Uukn/21vq3EAlHMPu82Ta9k3GK\n+S5IkoTbbrsN8+bNw7e+9S3NxwwMDKC5uRkMw2Dnzp2QJMlUwQWo6JYM2d1Ws2PHDlx++eU2r0Yb\nkl4gZjUejwd1dXW6T/jJntM1QqU0aQwNDmH387txEXMR8PcKwcF9g2ia14SP9n0ESZJyBBfQtmrc\n/+r+HMGVJAlLokuw949784qu3gi6VIpFutu3b8dTTz2F888/H4sXLwYA/OxnP8OJEycAAKtXr8az\nzz6LX/3qV/D5fIhEInjmmWcMr0sNFd0S6OnpwdGjRyf8fNu2beju7sZ3v/tdB1aljSAIGBsbk81i\n3BRhTRYK1Zk6nQ/e/+Z+jB0aw64zuwAGiNZHMa1jGoZ6huBt9oJBnnUIuf9kBO3HMVz+6RJ6Iuhy\nSCaTBUV32bJlEEWx4HOsWbMGa9asMbSOYlDRLYHt27cjnU5jYGAALS0tAIC+vj7cfvvteOqppzBj\nxgwcOnQIW7ZsQU9PDy6++GLE43G88soruOeee7Bnzx4cOXIEM2bMwJe+9CVL1ihJEjiOA8/zqKmp\nqbjuuGpGuSkXCAQKNmnoNX8pJ60zNDiEnrd7wJxkMIuZBQAYOTmC/kw/wnPDqPHW5L/TURW7SN48\nj8tz2PXv7dcVQZdDJdg6AlR0S2J0dBQPPvggHnjgAflgP3nyJDZu3CjX//3lL3/BjTfeiPPPPx8/\n+MEPMGXKFLz88svYunUr7rrrLuzevRs//OEPLRFdURSRSCQgiiICgYBtgluNpiR2UKr5i1mbcv17\n+xEeDWNO4xwcP30cM7wzUO+thyflwd6Te/HZeZ8FgLzTJJR0/XMX/vrbv55NMUjA31J/wwVXXKD9\nnkV9EXQ5pNPpnGYIt0JFtwRYlkUkEsH999+f9zErVqzAq6++ipUrV2LKlCkAgP379+O//uu/AAC7\nd++2xJmMmNWQ8T7FbqMKQXO6zpDP/EXpO6CsjCj3QseIDOqm1oGJM5g6dSp6xnqQZbPoHu1GdH4U\n/Xv70Ta/DdOXTS86TWLW/FnAF4FdW3cBHCD5Jcy9eu74zzVQzmPLwYRhJKlUyvRNLyugoquTnp4e\ntLa2Fn1cNBrFG2+8gUsuuQQAMDw8DI7j5IPhxRdfxOrVqzE2Noba2lrD61KPQg8EAnLpEqVy0TJ/\n0fIdIMJbSppB8kgI1YZQ11mH+Ok4aj21SJ9J4/zZ5yPQFEBHpgPdb3Zj+rLpmH/Z/KLPN2v+LFlk\nRVGUo3QtCs1jM0oqlUJHR4fh57Ea2gask23btuETn/iErsdu375dFt033ngDF198MYDxza0dO3Zg\n2bJlpuyKknRCNptFbW2tvGFDb/WtxYnyOPUkDXKBJWJcykSGtvltYJtYnAycRNPMJgQiAQSbgki3\npNHc2YzESALMIQZ/efov2PPaHgwNDmFocAh7XtuDvX/cK/9Mi2KfTWNz43gEHe5DX6APfeE+TF+m\nPY+tVGh6oco4fvw4Ojs7iz5OFEX4/X60tbUBGE8tXH311QDG83P/8A//gEcffRQ33HCDofVojUI3\nm8lae1sJkFSDKIqQJAmBQEB3k0ZjcyMWXLMAB7YfwK7uXRgIDGB623TM6ZoDP/w4te8UmoPNyIgZ\nRE9GseWVLQgiiI83fhxTpk9BrD5mqOKg3HlsxaiEScAAwEgG7kNp7s9aGIaRe/0JZNYU2anV6i5j\nWRbZbDanu6ZUhoaG0NDQoEt0eZ5HMpnMyTmaLdZkh98NJ1UqlUIwGDRtKrIRstksJEmacBwoN+XI\nUM98TRp7XtuDjsz4bfnxD45jKjcVALAvuw8RRMD1c6hHPZpmNmGQHa/ljdXH0Bfum5B+IGN4nJhT\ndvfdd+PWW2/FRRddVPzBDkLTCxUEMathWRa1tbV523ntjE5ZlkU8HofP54PH4wHHcUgmk0ilUmBZ\nlg6NdAiSCw6FQohEIrLXLLlAku+H53lMmzcN3eluAAAjjR873dnxf88MzBz/2d+/wuZgM4Z6/p5a\nMKHiwExoyRjFVOwcha4HdXuxslrCjPpTinkU8qHNZrMIRoNourAJ3Ye60ePrQcqTQtusNgx2DwIC\nIDGqiyb5qp0P9HNIpVITDGzcCBXdCoCMHNfbXWZ12ods4DEMIw+vJD4D5PXzmYQr60+tSkVMJsrJ\nu2s1aYRCIUxtnYqpc6aid3svouEoeImHJEpALTDEDKEJTeNP4MlfcaBej1Utv1rQSJdiCslk0nKz\nGi2IcKtPaJI+CAaDCIVCuk544sql9iNQj1whKQoqwsYpReyUF8mOGR0IBUPo29sH7hwOfz3wVyw4\nfwFERsSJj05gMDOI5vOace4/nltUPK1s+dUilUrR6gWKcURRdMUodOUGHilXKod8fgSkAYCmIoxj\nROyGBodwcv9J+ODD1Nap6FjQgfTpNERWRHBaEIvnLMbUlqnwer3geb7gnYqVLb9asCxbEbalVHRd\nTjnlYGanF/RMmygXrXyjVirC4/HQDTmdlCt2mmL9YXdOHW2h4ZBjQ2Po3d0LRmTgC/owdmYM0Ao8\nLdyAczo40QMVXZfjdJSndwPPLKHPl4og88lSqRRNRSiQJGmC0JTrb6BHrPNNVj710Sn0bO/B9MB0\ngAF8nA9H9h3B0OwhNDaohN5lG3B2Q0WXognDMOA4DqlUCqFQSHf+1uw1kFSE1+tFJpNBMBh0PBXh\n9oi7XH8DvWKtlS8+c/gMPl73cXAcB0mUIIkS5k+fj3cPvItLL7xU9pU4zh43peVXjdu/EyVUdKsQ\no1EniTBJCY5b7CH1piLsmNLg5gi7XH8DPWKdL1+ckTLA3wsHPF4PfD4fprZMRVu0Df2xfohZEaJH\nRPP8ZkRqI0XzweXi5u+FQEWXkoMkSUgkEpAkqWQ/XrvbhgulIpRTGiZbKqKxuRFYhqIOYWr0iHW+\nFMSbB98Ezp34nHVNdXLX2pmBM+jd04uhA0PgwaPp3CZ5U86MSRqV8v1S0aXICIKAeDwOv98vO1xV\nCpO1KiLfha4cfwM9Yp0vBdHY0YjudDdqE7UY7R2F1+PFce445v3r+PjyocEh9GzvyRX097oR/kQY\ndVPqDE/SEEWxYr5PKrplcvDgQZx7rsal3QWUk15QNmCEQiF546pScVMqwklKbU4oJNZDg0M49sEx\nJDIJSF4JLZ0taKgfH9JaN7UOoaYQdm3chWn+afB4PDhv9nkYPTiKoaah/Jt0B/vQdFlTwUkayiaa\nfN9TKpVyxO+hHKjolsHatWuxdu1a3HHHHfjKV77i9HIMIUmSbAlodwOGnVRjKmJocAjH3j0Gv8cP\nj98zQVDNbE4gz7WodZHsQta9vxvoAkaDo5h+4XT07+3HJYsvGf8sfeMXuwY0oG9fn65NukKTNNST\nldXOaZXSjQZQw5uSefTRRxGPx/Huu+/i1KlTeOyxx5xeUtmQdl6e51FbW1u1gquGpE6U3rSkqJ5l\n2RxvWmKd6DaICHZmOtGebUdHpgMn3jyR43ObL7r8aN9HJb8eea5YfQxN85pwJnAGNZEa7BvcJ9fx\nMuLf77DU+iqUV1GhNu0Jh8NyU0YqlZJNe44ePYrTp08XFN2enh586lOfwvz583Heeefhl7/8pebj\n7rzzTsyZMwcLFy7Ee++9V+RTKY/JcZaZyCc/+Uk5rfCDH/wABw4ccHhFE9GTXiB+vER4KiGys4pS\nUxFugIggy7Lyz+rYOrz1u7dwzoJzIHmkvM0Jo6dHsee1PSX5ISgj1Vh9DLH6cWMZb8B7tnGigLC2\nzTM2MaLQJI0nnngC69atQ2NjI3784x/jiiuuwNKlS3O6Jv1+Px566CEsWrQIiUQCF154Ia688kp0\ndXXJj9myZQsOHz6MQ4cO4e2338bXv/517NixQ9f6SoFGuiWizuNaMe/MaogdYzgcRiQS0RRct3kl\n27kekopQRlhKW0RJknRNaLB0jQoRZBgGwyPD6N/fj+mZ6XLke2b/GQyPDOf83fDIMIYODKEj05E3\nQtZCT6TaNr8NxzPHc37dne7GtHnTTJ8YoZykcf/99+O5557DkiVLwLIsvv3tb+Ppp5/OeXxra6s8\nPDYWi6Grqwv9/f05j9m8ebM8MHbp0qUYGRnBwMBAWesrhDsu2xRbUNsxuiVqczNaVRHJZBIAnG3Q\n8EgYHhlG75Fe+D1+nOg5gXl18yDGzlpsLp69GO8ceQeXXnip/LP3D7+Puc1zcfyD42AkBhIjYer0\nqfho30cFBVBPOVljcyNSn0ih/2g/PKJnQvWDVRMjgPHKm3PPPRcPPPBA0cd2d3fjvffew9KlS3N+\n3tfXlzMdpqOjA729vWhpaTF1rfSsM8Dg4CDWrVsHYHz3f3R0FL/4xS8cbybQGlaoZcdIKR3yeRKL\nTaeqIsLNYbz/h/exKLgIPq8PHtaDw8cOY/Y5s+XHxOpjaOlqQV/4bAlYtC2KdHcazcFm+XGD+wbB\nBsbTFPmqHfTW/jY2NWJaxzR4vV75uT7a/ZHlto7JZFLXRloikcD111+Phx9+WNN7V33nYsX3R0W3\nTI4dO4Yoz6fjAAAgAElEQVSNGzfinnvukXOB119/PR5//HGsXr3a4dXlohzPblU772TNCTtVFZEe\nTGPJwiUYODYAL7wYDYxiQfMCnB49nfO4mik1OSN1jr13LEdwgfFpELt6d02odhgeGcbrb76O1nmt\niDXG0Da/reh0YHKht9vWUc9QSo7jsGLFCqxatQrLly+f8Pv29nb09PTI/+7t7UV7e7vpa6XhThlw\nHIcXX3wRd999d47j1sGDB13l50nGsycSCUSj0ZI2zNyW060E7KyKYEQGsfoY2rvaMWPhDHRd2oUx\n/1hOCRbJpypp7GyUR/HIj8t2o6G9Afvf3I/M4QwO7jqIv/7lrzj4t4O4iLkIkSMR3blfgpmVE3oo\nVjImSRJuu+02zJs3D9/61rc0H3PdddfhySefBADs2LED9fX1pqcWABrplsWGDRuwatWqCT/z+XxY\nuXKlQ6uaSDKZlP143TBEcbKRz5GL53nDqQj1xlasPgbMA7Yf2Y5T+0+BAYOGcxom/F3tlFpEu6Lo\n7emV0wRts9rwkfcjDO8exmLPYgBAb38vkpkkhv3DQP3435bihVuu01m5JJNJNDU15f399u3b8dRT\nT+H888/H4sXj7/FnP/sZTpw4AQBYvXo1rr76amzZsgWzZ89GNBrFE088YclaqeiWwZkzZ9DY2IhH\nH30UBw8exM6dO3H06FF88MEHrjBRJtNfAbhinhplHGXxv9FUBNnYamVa5Z8dzxxHa3srFrYulH+m\nvqVvm9+GE8MncO6Cs1U4ZChlZ7AT4P6+VpHBDN8M9J7pRUOjQryLiCZJL5TrdFYuxdILy5Yty5nj\nl49HHnnEzGVpYmt6YePejfjPHf+JX//t18jwmeJ/4MLXUOZ5Wltb4ff7sWjRImQyGTz88MOmv16p\ncByHsbExAMhbDkZxHqOpCLkEK9KH3kDv+GZZHXIEF5h4S5+vdKuutg5Tpk/BIDsI4GwkPcKOoLFT\nEdnqFM22+W2ymBO00h1mUUkdabZGun3xPoxkRpDm0vj93t/j5oU3V9xrvPHGG/jMZz4DALj22mtx\n7bXXAgB8Ph82bdqE++67z9TX0wvJ32YyGcRiMbmsqVxoTtdeyklFNDY3InBxANFoFAzDYO8f9wJZ\njSdXRadapVv9nn45RXGm5wykVglHe4+ibmad3AhRrJlBebyU63RWLqlUioquFiFfCGkujYZwAz4/\n//MV+RrxeBwNDRNzZaOjo46JFLFjFEURdXV1tBysCtCbigDMuaWX63DrZ8oiu/fkXjD1DPoCfSWJ\nJrm7srIuVw0V3Tx8edGX8fu9v8fn538eIV+oIl+D4zjNn+/YsQOXX3656a9XDKUdYznz1Cjup5Bt\nJQB5hNHUOVNx7O1jmBmeKR8HeltttSLTrmu6bBNNo+gpGXMLtke6VqQU7HqNnp4eHD16dMLPt23b\nhu7ubnz3u9+15HXzQewYI5HIhA08O9MDZFKwKIoV5dJVqShTERzHIRKJQBRF1E+tB3chh6MHj8Ij\neuAJeNB5Qadu4bQzMjWbShm/DtDqhZLYvn070uk0BgYG5Pq9vr4+3H777XjqqacwY8YMW9Zhhx2j\nXtEmnW5kQKKWD2o1pTvclOcma1GmIjpmdKB9ervcJScIghwJmzGdodh6nLrY0vRClTI6OooHH3wQ\nDzzwgHyAnTx5Ehs3bpTNNKyGiByAou28VguE0qnM7/fL7v1aPqhEgK086e3ETe9BuZbhU8M5bbzT\n5k1D/dR68Dw/YTpDNU3QSKfTmm29boSKbgmwLItIJIL777/fkdcvxY7R6hNJndog+UUA8uRXtQVf\nNZ/0bkCz9XZ7N5hlzFn7RRMbNNQ4Gelms1nHPU/0QkVXJz09PWhtbS3+QItgWVbOWyl9Qu0mX2oj\n38mmzD8qR7KoT3rlSBZKeeQdiaPoItNbFVGJdyWVksaioquTbdu24bLLLrP9dZ2yY9TK6SonBZfr\nVJbvpOc4DplMBh6PJycXXEknvdOU2npbqCrCrruSUme4VQNUdHVy/PjxHK9Nu4jH42XZMZpdvaAs\nTTOr003rpCe54ExmvJtQax4WZRz17Typ0x0eGcZAzwAYgYHklSDOKd7+CuTelQwNDqFvTx8kXgIv\njY9Lb2xqLJiKKDW9YJYTmZs2N/VARVcnP/rRjxx5XTLBwMmIj1hDkknBVqGMgoPBoJwLJjvxwHia\nhZaladM2vw3vv/Q+fMd8mBmYCQAYTA0iNZLC0OCQbiHTFMN3uxH+xzBqG2vzpiJKRU86pBQq5Xig\noutyjIyVNhoBkPwtaS22e6NCuSGnNAmv9rK0cmlsboS/frxJ5ox4BvAATbOaEKuPFRQy9S3+2MgY\n5oXn5TxmZngm+g70YeplUwFopyI8Ho98t6LnomimE1mlCC5ARbdqMeMgzGazcmrDaWtIkoogTSCT\noSytHOpq69C+QMN4O4+QaUW1r3/wOhIfT8jtwFrPoeUVwbIsBEFAf08/Tu4/CR98YHwMOs/vxJSW\nKRNe2ywnMiLylULlrJRiG4IgyFNm3SC4WqjHc5MUDCllc/sIdTPQyqGWKmRat/jN4WYM9WiYlRc4\nDIgIj54ZxZl3z2C2OBvT+eloS7bh0B8Pofd4L1iWBc/z8vdhlhNZKpUydEdoNzTSrVLK3Ugj+dtK\nihjtLEtzu4AXGiCpVSmgdYvf0tmCvYf2YjqmT3iOQkiShIEDAzgnfA4AwOvxwuvz4tzAueg92oum\n1qacqohYfQztn2hH78FeMAJTthNZJdk6AlR0KQoymYxsHOLkeHGjWF2W5uayqXyWigA0KwWSSKI9\nnJuOaKhvQOOCxpyBlrodxjREnGEYeESPnBpSXhQjtRGcc9E5hjbkKsl3AaCiS0FuLTBJJ2QyGV1O\n++rncVtkXEpZmhMbclYMcNQyrtnz2h7NSoH97H50pydGxnOXzdX1+soLxukzp9G7txeesAcSI2HK\n9Cln88IKLdW6KBIRLqdBg6YXKK6gVMMahmFQV1fnOtE0m0JlaVZO782HkbKpUi5y+SoFamtqMW3B\ntLLMxpUXjOGRYTC7GUxLTcPpxGmc23AuBvcNAvOA08HTeVMT5KKobhvXatDId2dSSWY3ABXdSU0p\nXg7VipZPBBFgO8rS7BrgWGiDTR0ZDw0OYc9re4qmO5QXjIGeAcwIzAATZHBcOI7eQC/gBd4ffB+f\nuOkTuqP2QhM0iJe1WoSLie6tt96Kl156Cc3Nzdi9e/eE32/btg2f/exn8bGPfQwAsGLFCkvr8g2J\nbkNDw6Q8Ue1Ca0KFWbjFy8FNqE94ZVkay7Ly78ktsRnHvl0DHAttsCkpJd2hvGAwwtn/rovUYfqC\n8eeNBWKG2nqLpSL+/d//HSzLIhgMIpPJaDbv3HLLLfjmN7+Jm2/O77N96aWXYvPmzWWvsxQMie7Q\nkEZZiU2QLqXR0VFEo1H4fD6IoohkMglJkhCNRktKymezWbAsi5qaGkPrymQyEAQhJ7GvtGOMxWJF\nIybScltfX1/2OvKlF/R68ZZa/VCNF998UTAwPvLbDF8CvWJoFL0zy0pJdygvGJJXAqS/HwfKw9vE\ni4dWKuKGG27A7373O2zduhVNTU24+OKL8fDDD+Pcc89OO77kkkvQ3d1d8Lnt3DSu+PQCEQdyq1yu\nN4CZXgXK53HTLbzyolSuYc1khUS5DMOM77pHIqaUpdk5wFHPZIhS0h3KC0ZLZwuO7z6OKBNFy8fH\nDf6tuHjkrJVhsHTpUhw+fBhLly7FDTfcgNdeew1NTU0lP89bb72FhQsXor29HQ8++CDmzZtX/A/L\npGJFV3lQZ7NZZLNZzbE1pTyfGaKrXJebbuGtMKyZzJhZllbumBwrqkVKSXfkXDACgLhIRNqbxmj9\nKEa9o5ZO/1WSTCbR2tqK+vp6fO5znyv57y+44AL09PQgEong5ZdfxvLly/Hhhx9asNJxKlZ0gbNJ\ndlEUXdM5RSYnJJPJnBKsUp/D6AWArAMoPEuNYhy3l6WVQqnpDuUFI51Ow+/322Y/SjDaHKFMKV51\n1VW44447MDQ0hMZGay4YFSu6PM9jbGwMDMMgFAoZFlyzIl2y+eKGEiwjhjV2DrasJPR8Jm4rSyuF\nUtId6qaOho81oKW9xfY1Gy0ZGxgYQHNzMxiGwc6dOyFJkmWCC1Sw6EqSNGFMjBHMEBmO45BOp8Ew\njKFx6GashdzuCoKAuro6V0dXenHLRaDU79XusjSjHW560h1aVQ4f/uVDBC4NoGlaaTlVoxQbv37T\nTTfhz3/+M06fPo3Ozk7cd999cvnZ6tWr8eyzz+JXv/oVfD4fIpEInnnmGUvXW7GiGwgEwDAMUqmU\n4ycjcVhKp9MIhULgOM7R6EVtWOOmSKpcquE9APrK0vS6pWnldK3ocNNCq8phRmgGTh44abvoFot0\n169fX/Dv16xZgzVr1pi9rLxUrOiaTbnRpSRJSCaTEAQBtbW1kCRJtho0SjkbJcSwhkRV1WzuUg2Y\nPcTTqDG43ig5X5WDxNl/zFDDG5sxO/dYilAJgiA7cpGI0ox0RzlCqYy2Y7EYRFGUb6HsWgPFGKW6\npWk+h4EOt1KiZK0qB0mSwATsP26KpRfcRsWKLhEF5S69Gc+nF1IREA6HEQwGc9Zjd4RIDGt4nper\nJbLZLI1UK5xiZWkk8lVOajDS4VZKlKxV5XA8cxyz5s6S/23X0MlkMklF107MFDnyXMXyaJlMxvIR\nNnrWAuQa1jiVvyUpFp7nCw4upJSPVllaJpOR/x8Yj4KbPt6E7rfL63ArJUrWqnJo6zorqnbllgEa\n6VY1yo6ufBUBdka6pNstGAw6NrxSFEXE43F4PB6EQqEJpVF0nLo1kHZYIsQkF1zbWAt2MYsjh47A\nK3rhDXrReUGnLqErNUpWVzmQiz9g/tDJQnAcZ/v8PiNUvOhaEelqYbTN2My1AMW73Yx+Lnr+Xin6\ngUBAjnSB/A0CZkxvoExEuSHXObMT7dPb5ZJBURSRTqeLlqWZ6QNhl3ua/HoVdDxVrOhakUPN91xE\n4PR0dFkd6eo1rLEaktMmoq/Oq2s1CCjzkcU8UinlU25Zmt7GCK1cbUNTriOeXe5plUjFii7BSpEz\nInBGy7W03pfSrcwpwxqSQySObHo/E4/HM6FNVu2RWilRsNvXp6aUsrSGpgY0Xpb/9j9frla6WEIw\nenZD2S73tErcLK5o0TX74FcKXbkCZ9UJWapbmRUXI3VNsvoz0XuhKbYrb4ZlolW46SSXJKnkC2+p\nZWnqzz5frrZ3fy9mLjn7czvd09x2jBSjokUXsCa9YNSOUW/lgZ7nAJwxrFF/rmTDTFmTrPU35byO\nllmMWgTcYGZUjagvgCQKzuuWVqApQv39l+ueVu1Q0VXBcRyy2awr7BhJeqPUW3mz0VMlYVakUUgE\ngPHyIKVZDMU81Llgrc1Qlmch8AI8XlUe3qHiAeXmbaVQ0aJLoiQzRFcZYRm1iTTLsIYc6OXkb836\nXNQbZnaiFAG/349kMgm/358TBev1KaCUjtZmaPt57Tj6xlF0BjrBeMbL1nqyPehc1OnIGittEjBQ\n4aKrxMjtPLl1liTJFJtIo5DowufzoaamxrGGB9Lp5mSUrUYrCtbaEJoMUbDdI+89Hg9a2lvgv9yP\n/n39ELMiREZE88JmhGvCkCQJHMfZ6hlMyiYrCXecSQYwetARgxhS2G/GQWwkyiTr8Xg8Oe3FdkLS\nGoWaQJxGa0OI1KXSKNhatHK1pCGGlKURz2CrP/9KM7sBqkB0gfI2rrTaeZ20iVQb1hBrxnIpV/iV\nHWYAXCm4WjAMA7/fX7AsyspR6pWG2b4IJNUXCoVMcUvTC00v2Ey5DRLq0iez0wlmrMeo6JaDusNs\nbGzM9jWYQaEoOJvN5vx+MkbBVvgiKIMeo2VppVBpZjdAhYtuOWjZMRKccAjLV4pl1macXrQ6zNxU\nk2qEao2Cy83p2umLABQvSzPSnUjTCw6hV6Dy2TEqn8csm0g961Hmk802rNH7XPk6zKo1+isUBbt9\nfplZWOGLUEpjTLGytFKc6ozOR3OCSSG6dtkxlkImk5Et6fKVYlkdaRbrMJsMaEXBheaX2V0xYAVu\n8kXQKksrZYhnJUa6FX2WKb+AfAJF2nk5jkNdXV1BwTUrvVDoeYjQZTIZ1NbW5hVcs07sQp8LydlO\nVsFVQ6KwYDCISCSCSCQCn88HQRCQSqVko3hSTleptM1vQ3e6O+dn3eluTJs3zZkFKSAeEeFwGNFo\nVO7AZFkWyWQS6XQaHMfJd6SpVAqxWEzzuW699Va0tLRgwYIFeV/vzjvvxJw5c7Bw4UK899575r8h\nDariTMsnUGRMu8fjQU1NTVFhsTqnS/K3oihasoGnpJBok88lEAggGo0WfGwli4tRiACEQqEcARAE\nQb5wKgXAbsr9bhqbGzF92XT0hfvQF+hDX7gP05cZ80Ww4g6g2EXw85//PF577TUcOXJETk0oueWW\nW/DKK6/kff4tW7bg8OHDOHToEP77v/8bX//6101dfz6qRnTVByDLsojH4/IV085bQq31EKHz+XyI\nxWKOXQCy2Szi8TgikUhBX4lKv4U2GyIA5FY4EonA6/WC53k5CiZ1qnZeqMr9nhqbGzH/svmY/0/z\nMf+y+RXhkaC+CN51113wer3YtGkTmpubcfXVV+PMmTPy4y+55BI0NDTkfb7NmzfjS1/6EgBg6dKl\nGBkZwcDAgOXvo6Jzulo7/UbsGK1KL5Tix2s2JAIp15KRoo2WXaIyF0zywHRskTUwDIOlS5di69at\nuO2227BkyRL86U9/Kiiyavr6+tDZebZ9uaOjA729vWhpabFiyTJVceYRUVHPCysnT2lmlGL0AmD0\ntpWc7HTDzFrUO/Jam0HVPrbIqQ3GTCaDSCSChoYGrFixouS/V5/vdryHqhBd4OzGULl2jIC5Tlnk\nAiBJkqNCR+a6FbJkpJiLOgo2UhJFKYyRkrH29nb09PTI/+7t7UV7e7tZS8tLVYQ8pMaSJNvLPZDN\nSi+IoohsNqt7A8/KtcTjcV0bZlauYTJDSqKCwSCi0SjC4bCcC04mk0ilUshms2XlgquhfM0oRkrG\nrrvuOjz55JMAgB07dqC+vt7y1AJQBZEuGf3t9/sd978FxjeqWJaFz+dztD2RZVlIkoRIJIJQKOTY\nOii56ImCK2VskZJypliYQaGSsZtuugl//vOfcfr0aXR2duK+++6TfZlXr16Nq6++Glu2bMHs2bMR\njUbxxBNP2LLmihZdEh2EQiFbO8m0UDZg2L1ZprUOMoDQDY0g1YSZkb9WY4BybNFkyAUbpVB6Yf36\n9UX//pFHHjF7SUWpaNEl5VekVMco5YqucqOqrq4OHMfJV1Q716LeMCMewUag6YWJWCV+WsM7KyUK\ndirVQbo6K4mKFl2C2bnHUg6gQgY6dqJnhlmplPIcPM/L5iVuHChZaSijYKVHhDIKBsa/98kcBfM8\nX3Hlj5W1WhVmOnIpn08v+QxrzFhPKc/B8zzi8bglxjl6IEZCfr8foihSE3GTIRcwZRRMytEqIQq2\nmkp7vxUtuoB1Y9iLtcYqDcedzJuSxgst4xyrqw+U+eNYLAZRFHNqpqvBPtGNkCiYVOy4YYQ9raTQ\nT8WLLmBvaROZG1ZogKUdka5dHWaFjHuUnwO5EDEMI9/uBgKBvPaJpW4Q0RNaG60oWG0YXq13HETo\nK+09UdEt4bnUHW9Ofdlkw4wY51gVOeZ7f+rPARhPtZD0AqkkEQRBFlUt+8Rybo1pRDVOvmNUr2H4\nZBne6Uao6Op8LjLKRk/Hm5WRt3LDrNikYCvWIQgC4vG4bL9HRFYpmMT6kOy+k78jUTA58bXKpKgo\n6KfYxSefYbgVUTC9GOqn4kXXji+6UN7UKvI5lTm5YUY2DsnkDSK4anEkt3zqJgASdZH6anKi6701\nphgjXxRcqSPsOY6ruMoFoApEF7Au0i3XsMaKCNMJ4QfO3sYqX1+ZRtBzcpKTnTwfOdlJGoLneVmA\nC4kCcPZEqwRRcDPKKNiMEfZORLrJZLLipkYAVSS6gLlfPMlbAqVPVjBzI81Jq0ry9+l0Wt6w83q9\ncqRarqeEliMXEXGSjlAOKySDMlOplCzCtFvLXLTy7m6vPqnEUT1AlYguoK/US+/zkMkAfr/fkIGO\nUSRJctSpjFRIkNcnn42ZffbEi4C8Hkk/KPPBSlElPhJOdWu5JXdp5ToKRcHKuWVOV0Sk02mEw2FH\nXtsIFS+6Zn/hpBRKOZ6lnDWZ1X7LMAxisZjtB7Zyg4xUKPA8L+dqrUCdXlCKMPk8SCqikGdBpeUm\n3Y6eKBhAzn/bATlPK42KF12CUaEjUR3P8wiFQo6a1pANMwCGRg2V+5mQCgUAORtmdtZEKjfjvF6v\nfOfh8Xg0o+B8m3Ekaqvmkep2UiwXbOcIe5rTdRijDmHkNj4QCJgWHZVzC6gc7ZNMJk1ZRykoKxSy\n2awc2TgVMZJUTzAYlC+EWptxyppg9WYcHaNjHcruuGg0WnSEvZnQnK5DGPU7UBrWxGIxpNNpUzbB\nSkVrw8xu0SWCH4vF4PP5wHEcWJZFIBCA3+93pESN5O2UrdbKaEuZhlBuxqlrgoHqGqPjltwykNsZ\npvV5q3PBZkXBlegwBlSB6BpBXXdq5kFcysaeVRtmei9E6pZiUqFAqgaI+Hm9Xvj9fltuHbPZrDz/\nqlDFhjINQd6LejNOWROszE1WknViJaJ3eGe5UTBNLzhMKZFuIcMaO30cgNwOM/WGmVkVGYVQe/Aq\nKxSI2Clzd2TMDLmt9Pv9pooU+W44jkM0Gi15Y0ZrM47n+ZyaYCIGhawT3VgiVcnki4KVdx2lRsFG\n5qM5yaQTXbXIqE9qhjE+hVfvepzuMCO1yB6Pp2iFgjJKJJM6iEAJgiBHwEYqBkiKRRRFRKNRw2Kn\nNoNR5oGVm3HkGFBbJ2rdFrvptt4tlPOZmDHCPpPJYMqUKWa8BVup+Et4KTldMjGY3MY72VrKsizi\n8bg8rFDrwDKruUELQRAwNjYGn8+HSCQiH/R6KhRI1BIKhRCLxeQcNMdxiMfjSCQS8jSPUu4+UqkU\nJEkyRXC1ICc52ZQjFwgiwBzHyWsmFxdl6SDLsnK0z3GcrXdEapx8bbMhx1MwGJSHy+oZ3lmsZOyV\nV17B3LlzMWfOHPz85z+f8Ptt27ahrq4OixcvxuLFi/GTn/zEkvenpmoi3WLojSrNSi/kex4jHWZm\nQXLZkUhEztsaqVAgY2bKTUOQMfE+n8+2iF/dmqzHoAeA3Czihs04t0TcZkf/eoZ3vvbaa/IxrIUg\nCPjGN76BrVu3or29HRdddBGuu+46dHV15Tzu0ksvxebNm01bux6qRnQLiWWpvgVWRRGlbphZ6eFA\nKhSMCq6aUtMQypKwQCDgiJCoDXrIZ6JOQ5DPiFxEtAShWr1rnUKZoyc145lMBo899hi2b9+OV155\nBTfccAOuuuoqLF68WP7Md+7cidmzZ2PmzJkAgBtvvBEvvPDCBNF14o6h4tMLBC2BIvnbdDqN2tpa\nXYJr1omiXg9JbTAMg5qaGls2Z9TmPalUSv4sfD6fnN+0ssOsUBqCpCJIlOwWkSLRayAQQDAYnNCU\nQfK9ZLOR3BaTKJ2ML0qn0+A4zpQ9Aso4Ho8HkUgEzz77LFauXIk77rgDp06dwje/+c2c862vrw+d\nnZ3yvzs6OtDX15fzXAzD4K233sLChQtx9dVXY9++fba8h4qPdJU5XeXBrTba1issZkaXytbVcjfM\nzFgLufgQ03MrPBT0oExDkPlePp8P2WwW2WxWjpDdFCWSC4fH45ErH0gapZBBT6HNuEqrCS6GU5uL\n2WwW//RP/4Svfe1rE36nZz0XXHABenp6EIlE8PLLL2P58uX48MMPrVhqDlUT6SrheV7eJIrFYo6U\n/JAvXc+GWbHnMAKJsAHIpuc8zwOAI5+LsiY4FoshGo2ipqZGbnfOZDIYGxtDMpnMsXR0EvUmH4mA\nyWYcuYhxHKdrMy6VSiGTyRjajKNVFIWbI9rb29HT0yP/u6enBx0dHTmPqampkXPCV111FTiOw9DQ\nkHUL/jsVH+kSSIRKbu0ikUhZ/glmRrpkB9+pDTNJksBxnBxhO+GhoF4P8bdQXgy1ajhJhQBpyiB5\nVLujRGIp6fF4Jlw0i9UE6+2Mc2KQZDVQqE53yZIlOHToELq7u9HW1oYNGzZg/fr1OY8ZGBhAc3Mz\nGIbBzp07IUkSGhsbLV93xYuu8gAlu+ZGRM4M0SUnX6mpDTPXohQs5Vgdpwr9SbQIoKhrWqFqCAA5\nm3FWChSpqiBlZoVeS6smWLkRpzboUXfGVfIgSaei7kKRrs/nwyOPPIJPf/rTEAQBt912G7q6urB2\n7VoAwOrVq/Hss8/iV7/6lVw2+cwzz9iybkaqgoK/TCaDeDwOnudRX19vSFiIw1Z9fX3Zf08qFILB\noCG/z3g8Lu/ql0Imk0E6nUYwGIQgCHIdrlOCWyhaLAWSR+U4Ts6XkgjYbBtHLaOdclEb9BCRIgKs\nbGFWRsGFDGOIGZGTbnhOr+Wqq67C66+/XnEdgxUf6ZKcJbmFM6OLqdzrkHLDjNxa2rkWZQ1wbW2t\nfIvO87xjs6SIeJE8qJHPxK40BM/zSKVSCIVCpoxGKsegp9joeicvolo4FZFXwp2AmooXXY/HI99i\nmOHKVa7oqmuBycaLXahrgMntLsMwSCaTjuRFzRYvNVakIfI5m5lFPoOefEM7842uJ5UnRKidFB8n\n0guVvJFY8aILjJ9cPM+bKnKlOIRZ1WGm9wKgNs0Bzhq7RKNR+cTmOC5HkKwsz7JavNTka8ogF0M9\naQi9zmZmr1vdGZfPoEc5up4I72SellGpwlsVokswq323lNfL12FmRTeZFjzPI5FIyLlHrQoF5S57\noS4xszxzWZaVTa2d8LcoJw3h9JrJuksx6FGmIZzcjLO73ptQiYILVInoWlFmU+wqqjY/d6o4nJTH\n6Vr4/6kAACAASURBVPVQ0CNIJGIs9UTKVxLmNOo0hDrqJxdIJwVXC62hneT7UpelFRpdT20q3UVV\niC5gbmRZTED1dJipO+TKXUe+90QqFIx6KGjlRcltObll12PsrWwgcOoipAd1L386nZZH/SQSCcuq\nIYxC1k3EF4DsEwFo1wSTC3G1TMtQwnGcLWkrK6ga0SWYkecx0zzHbIi4cRwnpzTMauktVD9KfAaI\nICk/Y2VJmJMj60uB5OKBs516JIok79nJpoxCaxZFMefCVsigR+s7NXtahhO51UqdGgFUkejakbcq\nZcPMjMhb/RxaFQpWeSio88DKaEm5MeXxeJBOp221ZTQK8aJQ1w3nS0OQCN6upox8a1YavCtfX52G\n0DO0s9JH11fq1AigSkRXuWFkRaRbqiWj8u/MglQokO4ZIP+UBytQ7pyTXXNiLE1y6iS94Wbh1evd\nq5Un1bro2CFQym4+teBqrVur5ZhEwspqCHKRIa+htRnnVoOeSp0EDFSJ6BLMyusqn6fcDTMzDlLl\nLW+xCgU7UYpsOByGx+OxtRytXMpt1FAKmTJCJDWzJE9qRRoiX1SuF63NuHytyVqbcXrG59D0QmlQ\n0S2AEUtGs9ZC2pKj0WiOwbaTt3+kvEpZz6pVjkZOVjPL0crFzEYNu9IQRgVXjVpY9Rj06NmMcwKa\nXnAYs09mhmFkj1cnN8xICy+Z5+a04EpS4Um9VpajGYEIrhWNGkohA5A3911qGsLqEUZaNcGlGPQo\nN+PIz+xMQxQyu3E7VSG6BDOiS2Vuyym3MpLD43levk13wnRcvaZSJ/WaVY5mBLu7zAqlIZTiVUig\nSnE3MwuShlCWpREhBrQ348j3So7VbDZrW00wFV2XYFR0lRtm4XDYMQ/cRCIBAAiHwzmm2E4KLhky\nWWwjJx/llqMZwekus2JpCK33TCaemOFuVi5arcn5hnYSISbROPles9msfOdjxWZcMpk05ODnJFR0\n/45yw8yMW9By1kLWQCoUstmsfAA7VQhuxW2u3nK0cisDlGkQt3TG6UlDeL1esCyLUCjkCstGYOLQ\nTrVBD/E8IWmvfBUfxTbjSoXmdB3G6BdIRpKTDTO7HcKAs5t24XBY3rxgGAbBYFC+PbU7J2rXpN58\nt+TlNCgoW5H1pkGcQF2CR6Jy4Kw/rdWpl3JQRsGiKCKdTsPj8cg+wTzPy2suNC3DaGdcJpNBc3Oz\n6e/PDqpCdAnlRJdaHWZmRXN610I8FNQVCuoTUysnalW3lNW2jPkoZNdITvh85WiV0oqshvgkkLwz\nSUNYmXoxCrm4KT9rdU2w0qCHpCLMGl1P0wsuoRShK9RhZnbpWaE1ZDIZZDIZ1NTUFKxQ0MqJWtUt\n5YTFoRbK91ysHA2A4byzE2hVVhRLQygjfycg5w4xCVJ29JVaE6wsSSvFoKeSN9Lcee9VJnrFkmxW\nkXIsK4Sl2FpIVJbNZlFXVydXKOgpCSMHbTgclifqMsy4PeHY2Jj8vKVeOJQXgWg06qjgqiG3q6FQ\nCLFYTDb64TgOY2Nj8sTjSmlFBiBfNCORSN6cPbnTicViqK2thd/vl3P/8XhcTqXYlQ5TCm4hnw1y\njJINwUAgIAcVJGAg6yb1wJFIRD7uBEFAKpVCKpWSB7wq32Mx0X3llVcwd+5czJkzBz//+c81H3Pn\nnXdizpw5WLhwId577z1jH0wJuOesMkApJ5meDjNyq2QGWt06ZIeaYcYHV5J1lVOhUKw2Vm9kRE4m\nQRBcs/lUCJKG8Pl88vfJMIwc7dpRjmaEcu4mSF2tuhrCrjSEsh25FGMj9WZcKQY96s04r9eLP/zh\nDwU30gRBwDe+8Q1s3boV7e3tuOiii3Ddddehq6tLfsyWLVtw+PBhHDp0CG+//Ta+/vWvY8eOHQY/\nIX1UhegSiomlesOs0FXaDLMaLUiHmd/vRyQSkU8eszwUCtXG5mtXVZ5MlZYLVQ+PtKsczQhEcI2U\nsuWrhshms5akIcoVXC1KMehRBxQjIyP47W9/i+3bt2P37t1YuXIlrrnmmhxB3blzJ2bPno2ZM2cC\nAG688Ua88MILOY/ZvHkzvvSlLwEAli5dipGREQwMDKClpaXs96X7/Vv+CjZSSCxZlkUikUA0GjWl\npbIcyK1wKBSSJ/QqDWPMhkQNkUgENTU1CIVCcs1tIpGQc9rxeBwMw1SMLSMAeXNNXV5FxCgUCqGm\npkYWNpJ6SSaTcu7QCViWNSy4WpA0RDQaRW1tLQKBgGlpCDMFVw0RVuKHQdIQykietCp7PB40Njbi\nueeew+WXX441a9bgyJEj+N73vpfznH19fejs7JT/3dHRgb6+vqKP6e3tNe19FaIqIl1yEGiJbqmW\njPmep9x1kfSCskrCCQ8FdW2sKIrIZrOyoxQAeWqw24W3lPlrZpajGcHO2mEzG1GsFFytdSvv+LQ2\n40gOOJvNYvny5bj55ps1n0cP6nPcruO+KkQ3H+VaMpopuqIoyjWYxSoU7IK8P47j5M67Uoc4OoWR\nygoj5WhGcLJ2WOtiS2w5i6Uh7BTcYmtXGvQkEgm8//77eT/H9vZ29PT0yP/u6elBR0dHwcf09vai\nvb3dmjeiompEl1wljVoymonSr2Br/1b0J/oR8ATwhflfQCTgXDcNES5lpEjyokqjHaVJDdmUchIz\n23pLKUczcuwQwRUEwRXNGh6PJ+e7zueHQczpGYZxLB2nhJzfgiDgtttuw8MPPyxPvlazZMkSHDp0\nCN3d3Whra8OGDRuwfv36nMdcd911eOSRR3DjjTdix44dqK+vtyWfC1SR6AK5EZyeDbNiz2ME0qFD\nouz+g/0YyYwgxaWw6cNN+MJ5XzD0/OUgSRKy2WxB4VLvkKujQSeqApSRohW35lobNhzHlVwBorXu\nfNMe3EChNATpiAyFQk4vUyaTyeCLX/wivvzlL+OGG27I+zifz4dHHnkEn/70p2WR7urqwtq1awEA\nq1evxtVXX40tW7Zg9uzZiEajeOKJJ+x6G2Aku/tdLYL4FJB6TSOWjOQWpr6+vqy/JxUKpFvH6/Vi\n7TtrsbV7KwRRwBXnXIFbFt6CkM++A9roLa7ypOQ4zvSGjEKvS0QgEonYHikqo0EyqUPPhcfpW/Ny\nUZob+Xw++aLrdFMGy7K4+eabsXLlSnzxi1+smM9Ti6oRXZZlkUwmwbKs4YYHIprliC6JssPhsGxe\n4vF4wAos7njlDsysnwlO5NA1pcu2aFfZHmtWxKUUYEEQ5BPSTLNytwmXekc934XHbevWCxFctWl6\nuRces8hms7j11lvxL//yL7jlllsq5vPMR1WkF8iGGSmyNtpJVW56gWxEkW4pYoQeDAYR8oVwcefF\n6B7tRn2wHivmrjC0Rr1YNalXqyrATLNy5brdkFME9DmFke/d6/W6Zt16yCe4gDO2nASO43D77bfj\n05/+dFUILlBFkW4qlQLDMBgZGUFjY6Oh5yJF2Hqfh9y6sywrpxNIXpCIEcMwEBgBm49sxsp5KxH2\nW2/WQVop7TTDzhcVlVKWZfXUBCsg3zfLsnJZk13laEYpJLjFUB7nxGHMrI1Xnufx1a9+FcuWLcOa\nNWtc/RmWQtWILrnNHR4eRkNDg+Ed5+HhYV2iSw5YQRBQU1Mj77AqW3rz5UP/99j/YiA1gJAvhFXn\nrcKWI1vQH++X/20k50uMVJw0w853O16oLMsuO0mzUU57II0JyguuVeVoRiHHL/G1MHreqC+45b5v\nQRBwxx134IILLsC3vvUtV31mRqka0eV5HoIgYGhoyDTRLfY8Sg8FUr6ix0OBnJCP/O0RjLKjyIpZ\nzGuahzFuDKPsKNJ82lDOt5TmAbsg7Z7KDiP1balTdpJGUQquerc/3/s2oxzNKGYKrtZzqwMNve9b\nEATceeed6Orqwve+972qElygSnK6SpRdYEaeoxhksy0QCCAcDueMRS+WxyT50LpoHd4eeBtZPovW\nSCtEUUQ8E0dDuAGf+/jnylq71qReN5CvLIsU6RMj7EoTXC3/ByVWlaMZxUrBBfI3Zajz/uo0hCiK\nuOuuuzB79uyqFFygCiPdkZERufPLCMPDw6irq9M8EUiFQiQSkf1AieCWcpBk+Ay+/vLX5YqGWXWz\n4GW8uGbmNfDBV9KGlJNdT0YhfgQkF57PmMdtEMEt90LhVFUA2aS0SnCLofW+X3zxRUybNg3/+7//\ni+bmZtx7772u/d6N4p5QyCTM9k1Qo65QMNLSq65ouOm8mxDyhfD8wefRH++Hn/Fj5ZyV8LG+ghtS\n6iL8ShNc5QakMg9sdXuuEcxIhWhVBSitGq2og3bDJqXW+x4YGMB//Md/4Pjx41i+fDmeeeYZXHfd\ndRVrVF6Iyjk7i2D1wUNqL9PptFwHbIaHwqrzVqFrShe+seQb8sZZf7wfo+woehI9ePnEy6ipqUE4\nHJ7gEEZyhMlkUq7BrRTBJZF5NpuVBRfINWivqamRy9zS6TTi8ThSqZScI3QK5bQHs1Ih6vdNvksz\n3dHUuWc3XMBIOm54eBjXXHMNDh48iEsvvRRPP/203OhUbVRNeoEk7cfGxkzZQBodHZVd7InYiaIo\n+zhYORZ9zR/W4NjwMdQEa7D2qrV44C8PoHukG2FfGA9d8RBqA7VyfoyIfigUqgiHMKD8VIgyL6hs\nyLDTmIdEonbmzJV10OqyrFLL8OwsH9SDJEn46U9/imQyiYceeqhiggYjVN07NDu9IIqifMXNVxJm\nNrPrZ6MmWINz6s/BS4dfQvdIN0bZURwbPYYf/flH8kknSZLsk2B0VI9dkDuGciZUELOWWCyGmpoa\n+P1+cByHeDyORCIhj3WxCicEFzjrjka8cklDivquJ9937mbB/fnPf46RkZFJI7hAFeZ0gYk+meVC\nNkrKqVAwQk2wBnMa58ida/937P9wMnkSdYE6/PSyn2qWhGl1hrnNolHZ12+0HbnQhAwrNqTMmPZg\nBnrc0ZRleIXK2ZxEkiQ89NBD+Oijj7B27VpXHJ92UTXpBXLwJRIJ+YpuhNHRUdlkxUiFQjlk+Aye\nO/AcVsxdgZAvhJHMCP7ftv+Hn172U0Q8EV1+skqLRo7jTGnNNYJdGzhWGPOYaSlpJVrdYaIowu/3\nu2pcuSRJeOSRR3DgwAGsW7fO1Z+pFVSd6CprD8slk8nIGyUkmnDSdJxUM3glL1bMWoEpdVNKOlDN\naM01gpO3t2pjHiLAehsTMpkMOI6rqE1K4OxdmrJu3Q3DOiVJwtq1a/Hee+/hN7/5zaQTXKAK0wtG\ncrqk9CqbzconpdOCC4xXM5xOnEaaS+PVvlexqnFVSX+frzSJOI8Va801QrHmAasp15hHUozXqTTB\nJXW4SrNyO8rRiiFJEh5//HHs3LkTTz311KQUXKCKRJccOOUeQMSpjJiOk80er9fr6AknSRIYgUGG\nz6CppgnXd11v6PmUnULKFlUrHKPc1tarJw9MLrYsy1Zkowm5qyCDHgF97mhW5/4lScKTTz6J119/\nHb/73e9c1S1pN1WTXiBTEYh4RCL6x+GIooh4PA6v1ysXY2cyGcfH1ZATiAePl7pfwvVd11tqfJ6v\nJKscjwA3+j/kQxkJkjwwwzDy/Di37PYXg3iBlHJXYUY5WjEkScL69evx0ksvYcOGDa64ADtJ1Ymu\nch6VHsiUiGAwmJO/JQcb2ZQgblF25UKVlRNburegP2GO+5he1CdjKR4BRoZHOgkpZyObT/mMedxI\nOYKrhkT/ZrujbdiwAc8//zyeffZZxxzv3ETVia5yzHUxstksksmkrgoFdTRk5caE+racuJHpdR8j\nG29mibTyVlxZCaGO/pXfgdt3+tXkm/agrghwWxkeYE3eXJl6UrqElXrxee6557B+/Xo8//zzripZ\nc5LKCUOKoPYiKAZJH+j1UMjnmqTcmDBjM0orSgz5QjiZPFlw4oRSaOPZOFJcCieTJ/HcgecMjwVS\nb8TlG1aZzWbLanpwmkIm3urpuepJycpI0AmMmu7kQ+mOVurodsLmzZvx1FNPYdOmTaYI7syZM1Fb\nWytf9Hfu3Gn4OZ2gaiJdYLyekkRaNTU1mo8hEQ3HcaipqZEtBY1UKKjzgeVEBIWiRHXdLjAxml23\na50cDQ8mB9EcbUZ9sD7H08FslNF/NpsFANnE203mNIUod2qCVbfipWCV4BajUAkiaRzasmULfv3r\nX+N//ud/dN116uGcc87BO++8Y3gyjNNUTaQLQE4L5LuOqCsUzGrp1YoIyM4wORgLCbDSi0ArSgz5\nQhOiVWKKQ6JZZTR8buO5+GP3H9He1l72e9IDiYZYlpXfP8/zebuj3IaRhg11Z5idc8MA5wQXyF+C\nODg4iMsvvxwXXXQRDh8+jFdffdU0wSVUQ4xYOfeAOsknusRDgWEYSz0U1P4AXq8X2Ww2xylKub5y\nvQhCvhDSfFpOOSjdykbZUcxpnIO+RB+eO/DchL99/uDzeORvj2DdrnXI8Jmy3ysRLdLWS8QrFovJ\nzmFmumSZiZmOWyTKDYVCqKmpsfy9Oym4apTuaNOmTcOPf/xjnDx5EpFIBHPmzMFNN91kmlAyDIMr\nrrgCS5YswaOPPmrKczpBVUW6+dCqUCAHgpW5R2U+UMsXwev1yhtTpU7qXXXeqgkpBxINF8sBq6Pk\ncnK+xaLEQu/d6ZZkqxs2rHzvbhJcNa+//jp+85vf4MUXX0RDQwNOnz6Nd955x7RIf/v27Zg2bRpO\nnTqFK6+8EnPnzsUll1xiynPbSVXldLPZLHieRzweR319vfyzZDKJaDQKv99vq4dCPpTlbQBMFyGt\nHDDh+YPP4/n9z2MsO4Z/7PhHfOsfvlVyzteIaKnzgcopEXZsRjnZIWd0UoSbBXf79u2477778MIL\nL2DKlCmWv959992HWCyGu+66y/LXMpuqSi8oc7okT5pMJmUbQOWGmZM5RkEQwLIswuGwbNUnCAIS\niQQSiYRca1wuJAesJab98X7MmTIHXo8XrdHWkgWXVC2EQqGyRIsITSQSQU1NDUKhkLyZFY/H5dy2\nFbGAUrScqBdVv3diTE8M2tPpdF6DdjcL7ttvv417770XmzZtskxwU6kU4vE4ACCZTOLVV1/FggUL\nLHktq6m69AIR3VQqBZ7nUVtba0qFgllodWppbUooy7HMbMYI+ULgRA51wToMs8NYt2ud7lpes7vM\n1GV46s0oM/0BlNMe3NAhV6w1V2nMQ1I5blm7kr/97W/44Q9/iE2bNqGpqcmy1xkYGMDnPjc+rJXn\neXzhC1/AP//zP1v2elZSVekFcts2MjICn89ny5SHUijFIlApwCTyM6MWmKQeBlODSHEp3Q0XdneZ\nKd3BjFZCVFJLMjCxGxCAvOHn9DGsZNeuXfjOd76D559/Hm1tbU4vp2KoqkiXeCgAcJXgFisJ06KY\nMU25AkxSD+t2rcNf+v4CQRDQHmtHhs/kjXad8JNVu4OpC/P1doU5Ne3BCMSYx+PxgOd52aAnkUi4\nZlLy7t278e1vfxvPPfccFdwSqapId3h4WB5iSMrCAOuHVhaC5OyIIboZ4m9GFJjhM/jay1/DOfXn\nYO/pvYj4IljWuQyxQAynU6fHxXn+F8AIjKvsDUlXmNqgRWsT0i3THspBKx2ibkUHYMrdT6ns27cP\na9aswYYNGzBz5kxbXrOaqCrRJe5YIyMj8qBGJ4WC+JqW2u1U6msoncGKmXQrO9k4gUNfog/7T+3H\n+S3ngxM59I71ghd5DGWG0BpuxS8++Qs01ja6QnDVFOqMIvnRahFcNWZ5I5TKgQMH8LWvfQ3r16/H\nrFmzLHmNaqfqRJfMi1KOarE7EgDsG0+jfk2lMU1CSGCYG8aZzBmcSp/CQHIALx95GUOZIYxlxsBL\nPOKZOMa4MVw95RPYm+0BK/LgRA5z6udgYdNCLG5bjFULSjNNd4JqaEkGyt/ws8OY59ChQ7j99tvx\n9NNPY86cOaY852SkMpJcOhBFEf/5n/+J2bNn41Of+hRisZhlpjTFsKoWVJIkjLKjGEwNYiA5gMHk\n3/8/NYiB+EkMjvZhIP4RBtOnMZgdQoQJoFmMoJkLopX1oSXlQQ1GERSyWDQq4Pa3s5ieCaIm0oY1\nl+2Fr6senMeDFJvCKDuKaXXTDJum2wVpSSYRbzgchiAIjk5KKBUiuOXkn4sZ8xj1hD569Chuv/12\nPPnkk1RwDVI1ka4kSXjnnXfw7LPPYtu2bZg9ezaWL1+Oyy+/PGczxuoIuNSdckmSEM/Gc4V0uAen\nhk5gYLQPg/GTGMicxmB2GAPiGAKSB618CC0ZL1oSDFrGBLQOZdEyxKJFiqLF34Dm4BQ0R5sRamiC\nNGUKpIYGCPX14GprkaiNYJPwAVbMW4notE54I5HxhokDz+PNnjfhY3yY3TAb/zr3X3Hz+Tfb4t1r\nBsrNSnX+2eicNDswIriF0DLmKdWS9MSJE7j55pvxxBNPYP78+aasSxAELFmyBB0dHXjxxRdNec5K\noWpEV4koiti9ezc2btyIP/7xj5g+fTqWL1+OK6+8EpFIxBRXMC2UZVUZMYPBxEkMDBzFqVPHMDB0\nYjwSTQxgMHMaA9wwBqU4BjxpeCSgNe1Fy/9v786jojrvP46/h002lUUZFLRWY1yIS1HUWIygsS4g\nMy1pUYMSpdHG2lRz0hq15xDbqIkxicdoUpXYNDXHGIYlBiNFrVq1VhOtORqOrVlQBAzBiDiB2e/v\nD38zGcZBB5gVntc5/oFH5z5Owoc79z73+1FD7C0DcjXITaHI/XoQExiBPLgXvcPkyHv2JSQ69k6Q\nRkVBVBTS//8iIgIcOIux902YfzGfBl0DZ2rOoDfqKcgsIDIkst3vg7uZA9doNN73ZqW94eyeno/r\nqsC1Za8p+X7/71dXV5Odnc3OnTsZOXKk09by6quvcvbsWW7fvs2+ffuc9rq+oFOGrjVJkqioqECl\nUlFeXk6fPn3IyMhgxowZLS5BtHUnwL/LdlDx5b+pa6qjrrme6/qbfC2pqQto5utAHSYkYtUgb/Yn\nVt+NGCkUuV935EFR9A7pjTxMTkxEHDFR/Qnr1RcpOvpOiEZHQ2gouOEMzPxNuP3sdq7cukLPbj35\nQc8fcKz6GM36Zn7c78csHLXQq892rXeHhIWFtemHpvXHcE/V1HtyS5vtJwDrrWgBAQHU1tby+OOP\n88Ybb5CYmOi04167do0nnniCNWvW8Oqrr4oz3c5MkiQuX76MSqWirKyMyMhIFAoFM2fOpGfPnnZ3\nArQWwFteyeSLhi/oHRRJ76AookNjiOv9Q+S9BhAjH0hY73hk0dHg5ZvxDQYD3zZ+S1lVGZlDM9lx\nbgeHrxzmauNV/P38GSMfw+szXickMMTTS71La20P7X0td9fUe9MeYutPALm5udTX19PY2MgLL7yA\nUql06rF+/vOfs3r1ahobG9m0aZMI3a5CkiS++uorVCoV+/fvJzw8HIVCQVpaGhERES32g7Z2HdCZ\n3/SeYO+bPv98Pjv/s5MadQ2RwZGkxqcyovcI5g2f51WzcV353rtjP6w3Ba6tq1ev8uyzzwJw5swZ\n+vTpQ35+PklJSR1+7dLSUg4cOMC2bds4evQor7zyigjdrkiSJKqqqigsLOTDDz8kKCiIjIwM0tPT\niY6OvussyHwXXKvVOmUeqye09uCAxqBh8UeLuXTjEsEBwfTv0Z8dM3dQ+r9Sqm5VESgL5MvGL6lp\nqiE0MJTN0zYTERzh1rW3t+2hvcdy9n5Ybw7cGzdukJWVxUsvvcSkSZMwGo3861//YujQoU6ZrbB6\n9Wr+9re/ERAQgEajobGxkczMTN555x0nrN43iNC1IUkSNTU1FBcX88EHHwCQnp5ORkYGMTExwJ3h\nGyEhdz5ue8ONmLaQJMlSa2TvwYGi/xbx+sevc6P5BqGBoZQ8VkJseKylHLNJ38SRyiOEBYbRrG8m\noVcCr//kdbftBPDE/mfb47flYRRb3hy4N2/eJCsriz/96U+kpqa6/HjHjh3rkpcXvD8l3EwmkxEX\nF8eyZcsoLy9n9+7ddOvWjaVLlzJ79mzy8vIYN24cV65coUePHgQGBqLX67l9+7bXtSPYMt/lN7cl\n29uzWXO7hvCgOxUr8jA5hysPA983VUQGRzIochA6SYcePQMiB7DzPzup+7bO5f9+Z7Y9tJdtM0hA\nQAB6vd7SEKHValv995sD19yy4U0aGhqYO3cueXl5bglcM1/7hOgM4kzXQZIksX37dlauXMmUKVO4\nefMmM2bMQKFQ0L9/fwC7NeWeakew5ehd/vzz+ez73z40Rg2pP0hlxbgVBAcEozFoWHV0Ff2798df\n5k/FjQoGRgxEb9LTbGhmaNRQsoZmuezfbxu43uZ+NfXWgettjyU3NjYyZ84cVq5cycyZMz29nE5P\nhK6DTp8+TVZWFvv37ychIYGGhgb27dtHUVERN27cYPr06SgUCgYOHAh4VwC35aaTxqDhvc/eA2BO\nwpwW7cN7K/ZiMBkYHDWYEb1H0GxopvJW5V2tw85uh/Bk20N72O6FNv+e+YEZbzq7U6vVzJ07l+XL\nlzN79mxPL6dLEKHrIEmSuHXrlqUGyFpjYyOlpaUUFhZy/fp1Hn30UZRKJQ8++CDAXQHkzgB21tCd\nrZ9s5ciVI9Q31xPfPZ630t4CoPBSIf5+/pbJZLYD0W13ArR1K5Y3NyY4QqvVotFoCAwMtIwZ9dQ8\nEFvfffcd8+bNY+nSpZYB4YLridB1MrVazYEDBygsLOTKlSukpqby05/+lGHDhiGTyew+kumqfjBn\n3nTKP5/PFze/oKqxiq3Tt7bYsWC+yXa/gehtHczubW0PbWW7Q8R6J4TBYPBoTX1zczPz5s0jNzeX\nX/ziF247riBC16WampooLy9HpVJx+fJlUlJSUCqVjBgxwjJgvb1ngPfj7I/k9yq7zD+fT/mX5RhN\nRh794aNEhkS2euZrZm8rlvVQGqPR2KkC1x7bnRDmAHb1JQiNRkN2djbZ2dnMmzfPZccR7BOhmvx5\npAAADxxJREFU6yZarZZDhw5RWFjIxYsXSU5ORqlUkpiY6PQANp8huusjucag4akDTzEgYoBlJu/g\nqMEOVwHB3Y+kApYfGN50DdQR7RmebjsTwlX3AbRaLTk5OTz22GPMnz/f597bzkCErgfo9XqOHDlC\nQUEBn376KePHj0epVDJu3Dj8/PxafAQHWtyEut83iaf6wPLP51N5q5JrjdfQGrSo9Womxk+07H5w\nlF6vp6mpiaCgIIxGo9dOBWuNM9oqXFVTr9PpWLhwIenp6SxatMjr38vOSoSuhxkMBo4fP45KpeLj\njz9m7NixKJVKJkyYgL+/f5tGUrq7PNKaxqAhqziL6sZqJCT6du9L5pBMnhj1BNCysaK1Sw72Assb\np4K1xhVdcq3diGzLaEa488PsySefZMqUKSxZsqRDgavRaJg8ebLlIRuFQsGGDRva/XpdjQhdL2J+\n5FKlUnHq1ClGjRqFUqkkOTmZgICAewawTqfzeD1NVnEWl25cQq1T0ze8LwfnHbSE6/1utjkSWG3p\nR3M38/odLR5tD3ujGR3ZCWEwGFiyZAkTJ05k2bJlTjnDNY+iNBgMJCcns2nTJpKTkzv8ul2Bdz0W\n08X5+/szadIkJk2ahMlk4syZMxQUFLB27VoSEhJQKBSkpKTQvXt3y9mP+YEHgJCQEI+GT0hACD26\n9SDIL4iCnxUQHBBM0X+LKPu8jIvfXCQqJIrkfslkDs1s8fesn5K71/plMhlBQUGWdlxz+Gi1WstW\nvI60I7SXOwIXWjZEBwcHW/4f0Gg0re6EMBqN/PrXvyYpKclpgQt39nvDnU8nRqORqKgop7xuVyBC\n10v5+fkxYcIEJkyYgMlk4ty5c6hUKl588UUeeOABFAoFjzzyCHl5eeTk5DBkyBDLx3NPbUPaPG0z\na46uYV3KOsuWsprbNXzd9DVqvZqvv/saeZjc8ufNcyDa0zZsfbPR+mGE7777zi1jGc3MH7FdHbj2\nmGvq4fudEOaa+l27dhEREcGZM2d46KGHWLFihVPfB5PJRGJiIl988QVPPfUUw4cPd9prd3YidH2A\nn58fY8eOZezYsZZWjPfee4/f/va3xMbG8vDDDzN48OAWQ9m1Wi1NTU33nAnsbBHBEWybsa3F7wUH\nBGOSTJgkEw9GP0hMWAyFlwqZlzCv1XqdtrIOWeszQOsAbus1UEdYDw7y9OUN2460+Ph48vPz+eST\nT5gwYQIhISHMnTvXKZPCzMc7f/48t27dYvr06Rw9epSUlBSnvHZn5//8888/7+lFCI4zly6+/PLL\nJCQksHHjRk6fPs369espLy9HkiQGDhxIeHi45WO4+QzYfBnC1Wd/1ob3Go7JZMJoMjIoahDRIdHk\njMhBr73T1OHsM0SZTGa51BAUFGS5Fm6+BGAymZDJZJZf7WW+JOINgWtLkiTefvttEhIS+PDDD5HL\n5Rw5coRhw4bRp08fpx4rODiY2tpavvnmGyZOnOjU1+6sxI00H1RfX8/27dt57rnnLB8vHWnFsLcL\nwF3bsMwPV/xsyM8w6e6Ev7sHvzurG8+bA9dkMrFmzRpCQ0NZv369S97f+vp6AgICiIiIoLm5menT\np5OXl8fUqVOdfqzOSISuHWVlZSxfvhyj0cgvf/lLVq5c6ekltUlrrRizZs0iMjLS7lB2d2zD8qam\njbZUM1nz9sBdu3YtJpOJl19+2WXru3DhAjk5OZhMJkwmE/Pnz+d3v/udS47VGYnQtWE0GhkyZAiH\nDh0iLi6OpKQk9uzZw7Bhwzy9tHax14oxe/Zs0tPT6dWrF0CLbViuCmB3tj20lfU8hNY+BXTkpp87\nSJLEunXrUKvVbN682evWJ3xPhK6NU6dOsXbtWsrKygB48cUXAXjuuec8uSynsG3FkMlkpKWltWjF\ncMVISvOkM39/f6+vNrJ3GcZ8XdgZN/1cQZIkNm7cSF1dHdu2bfO69Qktif86Nqqrq+nXr5/l6/j4\neKqrqz24Iuex14oRHBxsacX485//TF1dHSEhIfTo0YNu3bphNBq5ffs2arX6nq0IrTFPOvOFwIU7\nNxmDgoIICwuzNIPodDp0Oh0ymQy9Xu9VzSCSJPHaa69RXV0tAtdHiC1jNrw9FJxFJpMRExPDkiVL\nWLx4Md9++y0lJSWsWLECtVrdohUjJCTE7oMI9zsDtm578MXBNYBl+I7tdjxnzUPoCEmS2Lp1K59/\n/jlvvfWWCFwfIULXRlxcHFVVVZavq6qqiI+P9+CKXE8mkxEdHU1ubi65ubmWVoxVq1a12ophMBhQ\nq9WtzgT2tbYHW+Y+OetLCubLLdbzENz9MIb1+rZv386FCxd4++23va4CSGiduKZrw2AwMGTIEA4f\nPkzfvn0ZN26cT99I66jGxkb2799PYWEhNTU1TJs2rUUrhr1hLH5+fpZJZ+0dLenIgBxXMQeu0Wi8\nZ5+c+c+2ZTC7s9a3a9cuTp48ye7du72u5FK4NxG6dhw4cMCyZSw3N5dVq1Z5ekle4X6tGObB45Ik\nWeYktPfsz9E2CmdrS+Da+7v3GszujACWJIl33nmHw4cPs2fPHp8c8N7VidAV2sVeK8bgwYN54YUX\nOHHiBD169Gj3TGD4fj6vbemlK3UkcO2xngjmjGoeSZLYs2cPpaWlvP/++z7ZGSeI0BWcQKvVsnHj\nRjZs2EBKSgpDhw5t0YrRlpnAZveqB3IFRyvq28vewxhtHcy+d+9eioqKKCgo6HANfVVVFQsWLKCu\nrg6ZTMbixYt5+umnO/SagmNE6Aod9o9//IM5c+ZQXFzMuHHj7tmK0Z4AdjVXB66t9gxmLyoq4t13\n36WoqIiQkJAOr+H69etcv36d0aNHo1arGTNmDCUlJV323oU7idD1cosWLWL//v3ExMRw4cIFTy/H\nroaGBq5evcrIkSNb/P69WjHM5ZOuvP7pCHcHrr3j328w+759+9i1axclJSWWObbOplQq+c1vfiPm\nJ7iBCF0vd/z4ccLDw1mwYIHXhq4j7tWKERgY2OIM2F3V5J4OXHvrMZ8BV1RU8Mwzz5CUlMSnn35K\neXk54eHhLjluZWUlkydP5rPPPnPZMYTvidD1AZWVlcyePdunQ9eadSvGiRMnGDZsGEqlkpSUFIKC\ngto9jKYtzIErSZLHh+/Yo9Pp2LlzJ++++y61tbXExsaSlZXF6tWrnXoctVpNSkoKf/jDH1AqlU59\nbcE+scFPcDtHWjGmTp3qsqHs3h64ACdPnqS0tJSjR4/SvXt3Tp065fQfunq9nszMTLKzs0XgupE4\n0/UBne1MtzUmk4mLFy9SUFDAoUOH6N+/P0qlkmnTphEaGuqUmcDeNF6yNcePH2fdunWUlJS4rHtM\nkiRycnKIjo7mtddec8kxBPtE6PqArhK61iRJoqKiApVKxcGDB5HL5SgUCmbMmGE5A25rAPtC4J48\neZK1a9dSUlJiGb3pCidOnOCRRx5h5MiRlvdhw4YNzJgxw2XHFO4QoesDumLoWrPXipGRkcGsWbPo\n2bPnXTsA7G3B8oXAPX36NGvWrKGkpMQyalPofEToerm5c+dy7Ngxbty4QUxMDH/84x9ZuHChp5fl\nMeZWjMLCQkpLS+9qxQAs4Ws9E9g8G8JbA/fs2bP8/ve/p7i4mNjYWE8vR3AhEbqCz7JuxSgtLSUw\nMPCuVoympiZ0Ol2LKWEdHcrubOfPn+eZZ56hsLCQuLg4Ty9HcDERukKnYK8VY9q0aXz00UdMnjyZ\nZ599tsU0MEdnArvahQsXePrpp1GpVC2G5wudlwhdodORJIkrV64wa9YsZDIZcrncUkvUt29fgBY3\n4VqbCexqFRUVLF26lPfff58BAwa47biCZ4nQFTodSZJIS0sjKiqKv/zlLzQ2NlJSUkJxcfFdrRhg\nfyawqwP40qVL/OpXv2LPnj0MGjTIZccRvI8IXcFhvjSZ6uOPPyYxMfGu4DS3YhQVFdltxWgtgJ3Z\nCHH58mWefPJJdu/ebRkGL3QdInQFh3W2yVT2WjEUCgVDhgwBaHENGNo+E9ieL7/8kkWLFvHXv/7V\nZ983oWNE6Art1pkmU5lbMYqKivjqq6+YMmVKi1YMZ4ykvHr1KgsWLGDXrl089NBDHV6zL0ygE+4m\nQldol848maq5uZm///3vllaMyZMno1QqLU9vtSeAq6uryc7OZseOHYwaNcop6+wsE+i6GhG6Qpt1\npclUWq2Ww4cPo1KpuHjxIsnJyZZWDD8/P4dmAtfW1vL444/zxhtvkJiY6NT1dfWnFX2RCF2hTfR6\nPenp6cycOZPly5d7ejlupdfrW23F8Pf3b3EG3NTURF5eHqmpqbz55pts2bKF8ePHO31NInR9jwhd\nwWFiMtX3bFsxxowZg1Kp5OGHHyYgIIDbt2/z5ptvsnv3burq6lAoFGRmZpKWlubUBl8Rur5HhK7g\nMDGZyj6j0cipU6coKCiwtGJMnTqVLVu2WGYEFxcXW+rru3Xr5rRji9D1PSJ0BcGJzK0Y69evJz09\nncWLF7v0eCJ0fY8IXUHwUWICnW8SoSsIguBG3jPfThAEoQsQoSv4NI1Gw/jx4xk9ejTDhw9n1apV\nnl6SINyTuLwg+LympiZCQ0MxGAwkJyezadMmkpOTPb0sQbBLnOkKPi80NBQAnU6H0Wh0WYOuIDiD\nCF3B55lMJkaPHo1cLic1NZXhw4d7ekmC0CoRuoLP8/Pz4/z581y7do1//vOfHD161NNLEoRWidAV\nOo2ePXuSlpbGJ5984umlCEKrROgKPq2+vp6GhgbgzkjGgwcP8qMf/cjDq3JcWVkZQ4cOZfDgwbz0\n0kueXo7gBgGeXoAgdERtbS05OTmYTCZMJhPz58/3maHqRqORZcuWcejQIeLi4khKSiIjI0M0SnRy\nInQFnzZixAjOnTvn6WW0y5kzZ3jggQcsTcBz5szhgw8+EKHbyYnLC4LgIdXV1fTr18/ydXx8PNXV\n1R5ckeAOInQFwUOc1S4s+BYRuoLgIXFxcVRVVVm+rqqqIj4+3oMrEtxBhK4geMjYsWO5fPkylZWV\n6HQ69u7dS0ZGhqeXJbiYuJEmCB4SEBDA1q1bmT59OkajkdzcXHETrQsQA28EQRDcSFxeEARBcCMR\nuoIgCG4kQlcQBMGNROgKgiC4kQhdQRAENxKhKwiC4EYidAVBENxIhK4gCIIbidAVBEFwo/8Dy73/\nzpdc584AAAAASUVORK5CYII=\n", "text": [ "" ] } ], "prompt_number": 18 }, { "cell_type": "markdown", "metadata": {}, "source": [ "The figure above show the simulated $\\mathbf{y}$ data as magenta circles. The green dots show the corresponding estimates, $\\boldsymbol{\\hat{\\beta}}$ for each sample. The red and green lines show the true value of $\\boldsymbol{\\beta}$ versus the average of the estimated $\\boldsymbol{\\beta}$-values, $\\boldsymbol{\\hat{\\beta_m}}$. The matrix $\\mathbf{K}$ maps the magenta circles in the corresponding green dots. Note there are many possible ways to map the magenta circles to the plane, but the $\\mathbf{K}$ is the ones that minimizes the MSE for $\\boldsymbol{\\beta}$. \n", "\n", "The figure below shows more detail in the horizontal *xy*-plane above." ] }, { "cell_type": "code", "collapsed": false, "input": [ "from matplotlib.patches import Ellipse\n", "\n", "fig, ax = plt.subplots()\n", "fig.set_size_inches((6,6))\n", "ax.set_aspect(1)\n", "ax.plot(bb[0,:],bb[1,:],'g.')\n", "ax.plot([beta[0,0],0],[beta[1,0],0],'r--',label=r'$\\beta$',lw=4.)\n", "ax.plot([bm[0],0],[bm[1],0],'g-',lw=1,label=r'$\\hat{\\beta}_m$')\n", "ax.legend(loc=0,fontsize=18)\n", "ax.grid()\n", "\n", "bm_cov = inv(W.T*Q*W)\n", "U,S,V = linalg.svd(bm_cov) \n", "\n", "err = np.sqrt((matrix(bm))*(bm_cov)*(matrix(bm).T))\n", "theta = np.arccos(U[0,1])/np.pi*180\n", "\n", "ax.add_patch(Ellipse(bm,err*2/np.sqrt(S[0]),err*2/np.sqrt(S[1])\n", " ,angle=theta,color='pink',alpha=0.5))\n", "\n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAENCAYAAAD0eSVZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXl8VNX5/z93liyTBLKQBRJKVFABJQGluIGhSoEoSAtY\nUVQUhLpUaVXQ/mzVr61VSmuxUgWkIlIFRK1QQgQsg0hMUxUqFRoQCCSBJJB9n+38/jgzQ9bJbHc5\nM8/79ZqX3pnr3Mczk+ee+ZzPeR6JMcZAEARBhBQ6tQMgCIIggg8ld4IgiBCEkjtBEEQIQsmdIAgi\nBKHkThAEEYJQcicIgghBPCb30tJSTJw4ESNHjsQVV1yBV199tds5ZrMZ/fv3x+jRozF69Gj85je/\nkS1YgiAIwjsMnl40Go145ZVXkJ2djaamJlx11VWYNGkShg8f3um8G2+8EVu3bpU1UIIgCMJ7PM7c\n09LSkJ2dDQCIjY3F8OHDcebMmW7n0T4ogiAIbeFx5t6RkpISHDhwAOPGjev0vCRJKCgoQFZWFtLT\n07F8+XKMGDGi2zkEQRCE7/g9eWZe0NjYyK666ir20UcfdXutoaGBNTc3M8YYy8vLY8OGDet2DgDG\njhxn7EQpY43NjDkc3lxWMzz77LNqh+A3IsfOGMWvNhS/uniZonukT7eM1WrFzJkzMXfuXMyYMaPb\n63FxcTCZTACAqVOnwmq1oqampvsbRUUBdgdQWgmUnAGaWwFB5JySkhK1Q/AbkWMHKH61ofjFxWNy\nZ4xh/vz5GDFiBBYvXtzjOZWVle6fDUVFRWCMITExsec3NBiAqAjAZgdKzwKnxEryBEEQouBRc9+/\nfz82bNiAUaNGYfTo0QCAF198EadPnwYALFq0CFu2bMHrr78Og8EAk8mEjRs3er6iJAFGA2DQA1Y7\ncLoCMEUCyQlAdBR/XWPMmzdP7RD8RuTYAYpfbSh+cZEYk3/aLEkS2Mnynl9kjM/k7XbAFMWTfFSk\nJpM8QRCEkkiS5PeCqvo7VF0z+cgIoM3C9fjyKqDdonZkbsxms9oh+I3IsQMUv9pQ/OLitRVSdiQJ\niDACzMB1+KYWoH8skBTPnycIgiC8Rn1ZpjcYAyxW/s+EfkBifz7DJwiCCBMCkWW0my0liUs1jAF1\njUBdA5AYDyT2A/R6taMjiJAnMTERtbW1aocRsiQkJPRsGw8S6mvufeFK8kYjUF0HHC8FauoBh0Ox\nEETW7USOHaD41aS2thaMMXrI9JD7xqn95O5Cp+MuGoMBqKoBTpQB9U3kkScIgugB7WrufWG3A1Yb\nX2xNSQRiosk+SRBBJBC9l+gbb8Y3NDX3vtDr+cNmA8oq+AaolET+T4IgiDBHHFmmNwwGIDISaLc6\nPfKVQffIi6ybihw7QPEThL+IO3PvSEePfFMr0NgCxMdxjzzZJwmCCEPE1dw9wRhgcc7eyT5JEH5B\nmru8kObuD5LEpRqHg9snaxt4zZr+sdx1QxAEEeKEdqZz2Sf1eqDiPLdPNvhunxRZNxU5doDiJwh/\nCc2Ze1f0Ou6isdmB8nNAVD131pi0WWKYIAj1qaqqwptvvgkAsFgsqK+vx7Jly2A0ilHrKjQ1976w\n2niij3HaJ6Mi1Y6IIDRHOGvuJ0+exPvvv4/HH38ceud63axZszBp0iQsWrQoKNeQW3MPbVmmN4zO\njlBtFqCkHKg4xxM+QRBhj9VqxbZt27BkyRJ3YgeA4uJixMTEqBiZb4Rncgcu2CcjI3kZgxNlfPG1\nh5o1IuumIscOUPzCIUk9P4J5vsxs2rQJc+fO7facwWDA7NmzZb9+sAgPzd0THZ0152q5syYlCYgz\nkR5PEGFIdXU1EhMTsWbNGhQXF6OoqAgnTpzAN998g8hIcSTc8J25d8XlrIHEO0GdPgu0tQMAcnJy\nVA0tEESOHaD4CWUpKytDeno6ACAtLQ1GoxHZ2dloa2vDihUrVI7ON8JzQbUvGOMavMPBvfEDEmin\nKxF2+LWY19uv3d7ex5/zZUxZ7733HqZMmYKEhIROz//iF7/A7t278c033wTtWrSgqgZuPT4CaGiC\neeOWXvV4rSO65kvxCwZjPT+Ceb6MNDY2dkvsAFBfXy+cc4iSuydcerxez/V4PzdBEQQhBlartcfn\nCwsLhZPYSJbxBZsdsFr55qeUJCBanMUVgvCVcPO5l5aW4k9/+hP+8Ic/dHrebDbjlltuweHDhzFk\nyJCgXY9kGS1h0PNF13YrcOoMcJb88QQRKuzfvx+tra2orKx0P1deXo4HHngAGzZsCGpiVwJaJfQC\nc2EBcq65jh+4ywszLtE0NPEF14R+mixKZjabhfs52RGKn1CK+vp6LF++HL/73e/AGIMkSaioqMD7\n77+P7OxstcPzGUru/tKjPz4RiIshfzxBCEh7eztMJhNeeOEFtUMJCqS5Bwu7U4+PigJSqd0fIT7h\npLmXlpbiiy++wO23367YNUlzFwW9ns/krU49/gzp8QQhCmazGddee63aYQQVSu5eYC4s8O5ESQKM\nzno1jc3AiVLuj7er548X3WdN8RNKcOrUKQwePFjtMIIKJXc5kCS+AcpovOCPryd/PEFolWeeeUbt\nEIIOae5KYLcDFivX4UmPJwQhnDR3NSDNPRTQO/3xVhtQcgY4U8W1eYIgCJmg5O4FXmvunpAkZ5OQ\nSKCxhUs152pl1+NF13wpfoLwD/K5K41Lj3c4+GJrXSP3x/cjfzxBEMEjbDT3hQVP4mjDSZgMUXh3\nwkrER/RXNR43bj0+kterMZEeT2gD0tzlhTT3IHG04ST2Vn6BHeV7sLBgidrhXMCtx9ud/vgqnuwJ\ngiACQLnkbrcrdqmeMBn4jHhsUhZWX7fMp/82KJq7J7rq8SddenzgYya65kvxE4R/KJjcHarOSN+d\nsBKzh9yKnT98TzuSTFc6+uOr65z++EbyxxME4TPKae7tFr4lv62dJzBaPOwbux2w2ICoCCCV9HhC\nWUhzlxdVNffS0lJMnDgRI0eOxBVXXIFXX321x/MeffRRDBs2DFlZWThw4EDPbxZhBL6XBsTH8QSv\n4pZ8YdDreWK32XnDbqpXQxCEl3hM7kajEa+88gq+/fZbFBYWYuXKlThy5Einc/Ly8vDdd9/h2LFj\nWL16NR588EEPV9PxGeigFMBmEyZRya65e8Klx0dGAI1NXKqpqfe6n6vomi/FTxD+4dHnnpaWhrS0\nNABAbGwshg8fjjNnzmD48OHuc7Zu3Yp7770XADBu3DjU1dWhsrISqampnd5r3rx5yMzMBADEx8cj\ne8RI5FxyOdDWDvOBrwAJ7oYYrmSqleODh7/VTjwOB8z/yAMMBuTMmAZER8K8dy9/3dkUwpVQ6JiO\nAzkOZ4qLi3HZZZcpcq2O4282m7Fu3ToAcOdLf/Facy8pKcGNN96Ib7/9FrGxse7np02bhqeffhrX\nXccT0c0334yXX34ZV1111YWL9KYb2e1AxXmgoZnPTDXYyUizWG18/PrFAMmJfHZPEEEkXDX3VatW\nYdWqVXjooYewYMEC2a6jCZ97U1MTZs2ahRUrVnRK7C66XlzydrFUr+cSTWoSYLFwbZnwDrdU0+Kz\nVEMQRM+sWbMGjY2N+Prrr3Hu3DmsXbtW7ZD8ps/kbrVaMXPmTMydOxczZszo9np6ejpKS0vdx2Vl\nZUhPT/c+AkkCEvsD3xsEMAdP8hqbLaiquXvCbZ00AFU1vChZS1un8RP9JzbFTyjJhAkT8MQTTwAA\nnn76aVx//fUqR+Q/HpM7Ywzz58/HiBEjsHjx4h7PmT59OtavXw8AKCwsRHx8fDe93StMUUBmOm90\n0daueIJfWPAkcvJnIXf3XNRZ6hW9dsDodHwDlN3BXTVnyVVDEP7QVWe//PLLVYokcDxq7p9//jkm\nTJiAUaNGuaWWF198EadPnwYALFq0CADwyCOPID8/HzExMXjrrbcwZsyYzhfxRTdijO/OrK7j9km9\n3p//L5/JyZ+FvZVfAABmD7kVm3NWKXLdoMPYhc1iyQncekprGYQfhKvmrhRya+7KbWLy9TKNzdzX\nLUk8yctM7u652FG+B2OTsrS9i9VbHM4dwRFGIC2JNwihjWOED4R7cq+qqsKbb74JALBYLKivr8ey\nZctgNAYnH2liQVUV4mKAi9IBo14RmcZTeQLNau6ecEo15oL9wOkKYaUa0TVr0eMPV06ePIl169Zh\n6dKl+OUvf4nnnnsOpaWl+Otf/6p2aF6j3eQO8FnnkEFA/1jZd7XGR/TH5pxV4s/Yu2LQk6uGIHzA\narVi27ZtWLJkCfQdZOHi4mLExMSoGJlvaFeW6Qhj3AtfcQ6QdIrINHKgek15kmoIHwhXWWbDhg3I\nzc1FYmKi+7lNmzbhpZdeQmFhISIjI4NyHbllGTF2vkgSn71HR/J654IWH3PVlAeAhQVLlF+0dblq\nrDYu1dAGKILoRnV1NRITE7FmzRoUFxejqKgIJ06cwDfffBO0xK4EYv1VRxiB7w3kbpqaesXcNObC\nAncpgEAIpKa8v/QYu9HA5ZrGFv7QsKvGbDa7t8WLiOjx+4r0vDITLvasPL8oOu7TSUtLw4kTJ5Cd\nnY3Dhw9jxYoVeP7552W5rhyIldyBC8XHTFHA2fN8V2uEUYhZ/LsTVmJhwRKsvm6Z+tp+x16uVTW8\nlytJNUSAyJV0lWLfvn2YMmUKAF5aZdq0aQAAg8GAjz76SKjkLobm3htWK7dLtoop02gKm43fKEmq\nIZyEo+a+evVqLFy4sNvz8+fPR1FREQ4dOhS0a4WvFdIbjE6ZJikeaKfaNAFhoFo1BGG19twtrrCw\nUDh5TezkDvDZenICbwTCHDzJB3m24a3PXYslDHzy6PdYq6ZV1Vo/ovvERY8/nCgtLcWJEye6PW82\nm1FSUuKuOSMKofPb2xTNNz2dPQ80twARypcQVt0NEyxcrhobuWqI8GH//v1obW3t1I+ivLwcDzzw\nADZs2IAhQ4bg2LFjyMvLQ2lpKa6//no0NjYiPz8fS5cuxX//+18cP34cQ4YMcfe4UBOxNfeeYAyo\nbeAzT71e0YQUciUMAKpVE8aEm+a+atUq3H333fjd734HxhgkSUJFRQUefvhhZGdnAwDWr1+PyZMn\nY9SoUTh8+DCSkpIwZ84cXH311Xj88cdx6NAh/PKXv8S2bdv6vB753H3FVUI4OsrpibcAkcq4aTTl\nhgkW5KohwoT29naYTCa88MILvZ4zc+ZM7Ny5E7Nnz0ZSUhIA4MiRI1i5ciUA4NChQ5qpJBm6U7Do\nSCBzENDPxDc9BbA46K1urcUSBkGri+OSahwO4JRyZYVF16xFjz9cKC0tdbcU9URMTAz27duH8ePH\nAwBqa2thtVrdu1m3bduGW265BQ0NDbLG6w2hm9wBLssMTOYPi5VbJ4nAMBh4kidXDRFCmM1mXHvt\ntV6du3//fndy37dvn7uhh91uR2FhIW644QZs3LhRtli9JbSTO8Clg/g4vtiq96/CZDB2p6qFLLF3\ndNVUVnNXTbM8rhrR7GddET3+cOHUqVMYPHhwn+c5HA4YjUYMGjQIAJdkcnNzAQB6vR7f//73sWbN\nGsyaNUvWeL0h9BZUPeFwAFXVQG0jd9PoQ//epgg2G5do+seSqyaE0MzfbYhCm5iCiU4HpA4A0lN4\nQrJ4J9OoXc89EP+8IrF3kmpKgZq6oEk1omvWosdPiEt4JXeASwr9YrlME2FUrF9rIAna5Z/fUb4H\nCwuWyBRhgLilGiNQWSOrVEMQRN+ElyzTFYcDOF+nSL/WQHq0CumfJ6lGeDT7dxsikCwjJzodkJII\nDE7jXZ5kKF3gIpByv55aAGqWnlw1lCgIQjHCO7m7iDVxmSY6km96cnROQsHQrQNJ0IH451VdL+jk\nqqkBTp3hMpgPiK5Zix4/IS70W9mF0cBn8DX1vBmIQc9nn0HClaCVwtXSr/V4Kz4Z8666M36djt84\nrVagpJzvIE6KV6TRCkGEK+GtufdGSxtQXgU47NwyKeA2+0A0fllhjMtfBj1vuhJrEnJ8wwHh/m4F\ngzR3NTBFcZkmJvDSBWqhRks/r5AkrsVDAsoq+U1UgTIGBBFuUHLvDYOe++FTk2D+fJ9wCSg5MgnJ\nkYlwnNTojcmg50m+pdXjgqvomrXI8SckJECSJHrI9EhISJD186Pk7glXhcm0AYBOAtqV8cQHg1PN\n5TjXXoOvag5p2xsf0aGMgR8LroR81NTUYM+ePWCMCfvQcvw1NTWyfn6kuXuL3c4dH/VilC4Q0htv\ntXJLamI/WnAlCASWOym5+wJjQEMTUFHtnHUa1Y6oV+os9WLWlqcFV4JwQwuqMuPWTSUJ6B/H68Qb\nDYqVLugNTyUNXNbLg19/q1J0ftJpwbUK5g/+Ltx6R0dE1twBil9kKLn7Q2QEMGQgLyXc3s4lGxUQ\nouaMvxj0QFQE0NrOF1xrG4RZ7yAILUCyTCAwxrfXV5zn/x6hTDs/F0Lq6v7gcJaGiI7ki9tRkWpH\nRBCKQJq72lhtvO1cSytfbFWogbSwuro/MMaLkdkdHXa40g9PIrQhzV1m+tTtXKULUpK440Ohdn7e\n1JxRuxZ9oLjjlyReTjgygnviT5YBjc2al2pE13wpfnGh5B4sXJ74IelcL1Z5sTVk6bLgijO0w5Ug\neoJkGTlwOIDztUB1PZ9tGsivLQuMXeimlZLIF7jJNkmEEKS5a5XmVuDMOZ7sFV5sDSvcC65RQFoS\nLbgSIQNp7jLjt24XE80LkMU5C5DZla/zEjKauyd0Op7QrVbe3q+qRpWx7gnRNV+KX1woucuNQQ8M\nTAYGOZtyy9jtKawRdMGVIOTCoyxz//33Y/v27UhJScGhQ4e6vW42m3Hbbbfh4osvBgDMnDkTzzzz\nTPeLhKss0xWL1WmZbAcijYpZJsMSm53fTONM3MVEPVwJAZFNc9+3bx9iY2Nxzz339Jrc//jHP2Lr\n1q2yBRhyMHah25NeT0lHTmjBlRAc2TT38ePH91lzOBySdlB1O0niG3AyB/FNODJbJj1p1p5q02iF\ngNYMOvZwrTgPnDqreElh0TVfil9cApo2SpKEgoICZGVlIT09HcuXL8eIESN6PHfevHnIzMwEAMTH\nxyM7Oxs5OTkALnwAWj0+ePCgPO8/YQJwrhbmT3YCej1yrr+Bv+5MaDnXXCfrsas2DUqAGd/Nh/mn\nWxS9vmLHRYUAA3KuGguUnIH5yCGgfyxyfvAD/rrGvm90HL7HZrMZ69atAwB3vvSXPq2QJSUlmDZt\nWo+yTGNjI/R6PUwmE3bs2IHHHnsMR48e7X4RkmU809TCtXiHQ9GerWFTm6YjVFKYEAjVrJBxcXEw\nmUwAgKlTp8JqtcreXSQkiTWp0rP13QkrMXvIreGT2IHOO1zLK2mHKxGyBJTcKysr3XeVoqIiMMaQ\nmJgYlMC0hCK6ncHAe7YOTOZ+bYs1KFq8J83am9o0aiObT9+gByIjgSZnD9c6eUoKi675Uvzi4lFz\nnzNnDvbu3Yvz589j8ODBeP7552F1FsVatGgRtmzZgtdffx0GgwEmkwkbN25UJOiQRZK4o8MUxXe2\ntrU7q0ySbCALrgVXh4MvuNY10Q5XImSg8gNaxW2ZrAH0BrJMyg1jXJ5xUElhQjtQbZlQprX9gi4c\nqdxiq1osLHgSRxtOwmSIwrsTViovGXVccE0bwEtIhPiYE9qFasvIjKq6XXQkkJkOxMc669P41tJP\ntNoyXVsHKh5/p5LCFQEvuIqu+VL84kLJXQT0OiB1AJCRygtihXB9GpMhCgAwNikLq69bpl4gCi24\nEoRckCwjGlYbUFnNi2JFKtfSTyk02TrQ4QAsrpLCA/i4E4QCkOYebjAG1DcClTVcRogwqh1R6NNx\nwTU5AUjoF3I3VkJ7kOYuM5rT7SQJiO/H69NEGDzWpxFNc++KZuJ33UQjjLxe/KmzfLG7DzT33fER\nil9cKLmLTGQE8L2BQFJ/rsPbaKel7Oh0XJ6x2YFTzsYgCu0oJghfIFkmVGhp4xuf7DZF69OENS7b\npNEADBwAmKLVjogIMUhzJzh2O9fh6xu5fKDXfmNu1X3twcBm43p8QhyQnCjEuBNiQJq7zAij2+n1\nfAaZnuK2TJq/0Ihm3Qtdfe1dkUtzD2ote4OBe+PrmoAT5Z3a+wnz3ekFil9cKLmHGpIE9IvlVSaj\nIwGrRdOasFq+9r5uKj7j2vwkSUAZVZsk1IdkmVCGMaC2gden0el4A2mNoZavXdZa9oxxX7ykA1IT\n+c2W1kAIPyDNnfBMm4U3A2m3+FWfJiR08S4oclOxO/gvpxgTbwxC+xEIHyHNXWZE1u3MZjMQFQEM\nGcirHba1cxufDwRdwvABuTR3RWrZ63UwH/iaO5lOlvEqn4JNckT+7gPixx8IlNzDBZ0OSEnkSR7M\np8bcmqn3IiISnE26jbxsxKmz/BcUQcgMyTLhiN3ON9/UNwIGIy+S5QFN1nsRkY4lDAbE819SVMKA\n8ABp7oTvMMYbc1ecBxyM68G06KcMDgdvoxhh5NbV6Ci1IyI0CmnuMiOybtdr7JIExMUAF2UAcSa/\nasUrQVfNPaj+dAXocc1Ap+O2SbuDyzRV1fzfNYjI331A/PgDgXq3hTsGPW/KHRvDZ/E2ecoXBMtx\n41rc5e+5BJtzVmkiLr8wGvj41zQADS3OEgZR9AuKCAokyxAXsDlrxTc0B718QU7+LHdSnj3kVr+T\ncrD96cGKK2BsdsBq5Q3SkxP7XAchwgOSZYjgYDAAg1I6lS8IlnUvWI6bdyesxOwhtwZt45FmnEAG\nPZdq6pu5bbJDCQOC8AdK7l4gsm7nc+wdyxfERDu1+MD1YH+TclfNOtj+9GDfLLrik09fkvieBJ0O\nKKviG89ULuMs8ncfED/+QCDNnegZo4HP4BuauFRjQ0COGldS1hqajEuvB6J0fPbe1AqkJfHFb9Li\nCR8gzZ3oG6sVqKjm1skQ7NvqCdVLL9jtgMXGHU2pSfymS4QNpLkT8mI0Ahmp3FVjtQVVi9c6apZe\nAOCcxUcALa1ci69rDJuxJwKDkrsXiKzbBS12SeJODlcp4fZ2RUoJq91DNdAF16DEL0ncnqo3cB2+\ntIJvglIAkb/7gPjxBwIld8I3IozA4DQgJYnLNRZrSM8k5V5w9Qm9c/NTa7uwhcgI5SDNnfAfi5XP\nJFvawk6LVx2Hg9eMj44C0gbw8SdCDqotQ6iHqyFIlashiIFcHUrRsRBZcgIvREZjH1LQgqrMiKzb\nyR67JPGkclE6l2zag9vWT23NPVBkjV+S+JhHGIFztcCpM7wxSxAR+bsPiB9/IFByJ4JDpLMhSHIC\nl2sUWvAjcKEQmcUGlJQD52s13TeXUAaSZYjg42rr19ZOWrzSOJz9WyONQFoydzYRwkKaO6E9HA7u\n5jhfx73atPnGzcKCJ7GtbDcsdgvGJF2J9+Vo92e1AXYbkBjPG4PQDVZISHOXGZF1O9Vi1+mAAQlA\n5iDAoPOprV9HQlFzP9pwEhWtVaix1GH32X3ybI4yGoDISH6DPVnOHU1+IPJ3HxA//kCg6RQhL1GR\nwJBBQHUdf+gNITGLD6QsgWtjFACMTrxCvmqUksTH32YDTp8BEvrzWXwQSzkT2oVkGUI5WtuAs+f5\nYmtk8BuCKEkgdeDrLPWY9/nPIQF464ZXlNkcxRh3MhkMvClITLT81yQChjR3QhwcDu7mqGngNcwN\nYs7ivWkaonrRsZ5wNQVJ7M9lMz0ps1qGNHeZEVm301zsOh0vXfC9gXzm3ocWr1XN3ZuyBEcbTmLv\nv1QsOtYTrqYgtQ3cNtmHFq+574+PiB5/IFByJ9TBFAVkpgMJ/XiCV7kpha940zREM12euuLS4hkD\nTjsbdJMvPuTwKMvcf//92L59O1JSUnDo0KEez3n00UexY8cOmEwmrFu3DqNHj+5+EZJlCE80tzq7\nDtmF1+I7Umepx8KCJVh93TJtSDI94dLiI4xci4+O6vu/IRRDNlnmvvvuQ35+fq+v5+Xl4bvvvsOx\nY8ewevVqPPjgg34FQYQ5MdG8fEF8HE80NrvaEflNbXsdTjaeBuB5dr+w4Enk5M9C7u65qLPUKx3m\nBVyzeLuDly84V0Oz+BDBY3IfP348EhISen1969atuPfeewEA48aNQ11dHSorK4MboQYQWbcTJna9\nnlc3HJwGgLm1eK1q7j3RZG1G7qd3Y/3xLfyJ2tpe41e9CUhXXL746npnjZp2AAJ9f3pB9PgDISCr\nQnl5OQYPHuw+zsjIQFlZGVJTU7udO2/ePGRmZgIA4uPjkZ2djZycHAAXPgCtHh88eFBT8YT0cUw0\nzKe/A2obkDNiFGB3uBNkzjXX8fM1eGxxWLCs6XWMjL8ME1qvgfnDD5Dz8m+B8RNgBrD8v2+gKaMF\nJkMUHoq8F63HW4FYrsffo58Fc2GBNv5/oiJh/nwf4HAgJ3cK4HBo6/sR4sdmsxnr1q0DAHe+9Jc+\nrZAlJSWYNm1aj5r7tGnT8NRTT+H6668HANx8881YtmwZxowZ0/kipLkTvsIY79lacZ7XSwmgObfc\n2Bw2zDIvRITOiPcm/AX6klPAnFlAlfNX7GO/QM7lX3TyxUfro7Cj/J/IShyB93NWa0+Td9WoiYrg\nNWqiqF68GqhmhUxPT0dpaan7uKysDOnp6YG8JUFwJAmIiwEuyuDNodvaebNojeFgDty//3G0OyzY\nMP7P3RM7AKz4Iyb+twnABefMqeZynGuvwe6zn2tDlumKzqnFW52VJqvrqOuTYASU3KdPn47169cD\nAAoLCxEfH9+jJCM6Iut2IscOgEsEA5OB9FS+6NfuX40aOWCM4dGiX6GkuRQf5KxBxKmybondDAA/\nmYPHHnmvky9eszbJLpi/+neHevFn+YK3QIj+/Q8Ej5r7nDlzsHfvXpw/fx6DBw/G888/D6uV1+le\ntGgRcnNzkZeXh6FDhyImJgZvvfWWIkETYYYkAf1iAFMkUFkNNDYDRqOqNVIWFjyJT87sRZ2lHt9M\n/xQmQzRvYh3RRb6Y+APgxWWI1+k6lSh4d8JK7dskXeh03KJqsfIiZCmJfH+CRmUygkPlBwixYIwn\n94pqgDkQF/JmAAAc6klEQVR4MlUhyVzywbU40cQtj51qy5Q5Z+9lpcBP5gAvLgutcrsOB5+9m6L4\nL6oIo9oRhTRUW4YIP2w23re1vpnb+AzKzeJXH92AxUXPotXe1nNtmbIyYNPfgJ8/GVqJ3YWrdytj\nvPMWzeJlg2rLyIzIup3IsQMe4jcY+MwxI4XP4P2sF+8r7534O57/zyv4fOpHvdeWycgAHl8K6HRC\n+fR7osf4Xb1bjQZeuqC0QrNtFUX//gcCJXdCXDo6avrHOmvUyOeo+Ufpbiz+97PIv3kDxtTHYPOR\nKxFv7Cfb9TSPTsc3PrW2AyfLgLpGzSx2EyTLEKECY7zCoUw1aswVBbh970/xjx+8je839gPmzOau\nmId+BjyxlGQJlxYfawLSkviCNxEwpLkThAu7nfdtra3n0k0Q6sX/+/xB3PLpPdg04XVMbE27kNhd\nUILnMMblGQlAahLQL5bGJEBIc5cZkXU7kWMH/Ihfr+eJ5XuDvKoX3xff1hVj2qfzsPa6P/Sc2AGg\nurrXa4Sk5t4bksR/MekNwJlzQHkVX3hVEdG//4FAyZ0ITUxRvNJkYn++8cmPJHOi8RQm77oLfxz7\nLKZl3Aw8vaR7Yv/JncCLL4emK8Zf9Dq+u7W5lfvi65tIi1cBkmWI0Ke1nWvxFiv3xev6lgrKm89i\nfP6PseSKB/HTy+7hT549A9x5O1Bykh8HmNg12YYv2NjtfNzjYrgWL2hbRbUgzZ0g+sLhAGrquR6v\n03ncfHO+rQYT8n+Mey+ZjaVXPtz5RVeCH3dtwDP2QJpsCwVzFiGTdLysc5yJtHgvIc1dZkTW7USO\nHQhi/DodbwidOQiIMHAtvoemFA2WRkzZfRduGzy5e2IHgIGDgA+2ep3YPWnWItSXCcqagSRxy6Re\nB5RXXnA0KYDo3/9AoN9IRHgRFQkMGcRn8edqeYI2GgBJQoutFdP+OQ/fH5CNFwfP5zPOnmaYiYlB\nCaWv+jJdZZslX/5GbBlHrweidEBjC9fj0wZw6yTN4mWBZBkifGm38HrxLW2wGIAZex9AYmQ81g96\nFLo7fwJMmw78v2dVSz5dZZuqturQkXFsdsBmBfrFAamJqhaB0zIkyxCEP0RGAN8bCHtyPO7e9yiM\n0OGttId5Yq+qBNauAX77vGpOj66yTdfjYPVhVaWfq0HPpZqGJuBEOW/MQhPAoELJ3QtE1u1Ejh2Q\nP34G4KdfLMV5XQs2XbQUxrvmdLY7rl0D5P3D7/cPRLN+d8LKTvVruh4Hqw/rtrLd7veZ9/nPgxZ/\nn7iac0sSr09TeT7oDVlE//4HAmnuRNjCGMOTu57EocpD2HXtXxA16VbgXFXnk35yJzD1FlXii4/o\n30l66XocrAVZi/1CAw5VBCiDni+21jUBTa28IFxMtBqRhBSkuRNhy28++w02f7sZ5nlmJNZbgIkT\ngf/978IJs37CXTEarZNSZ6n3qeFHb776STvvwO6z+zA68Qr8c/JmdRdqbTa+4SyxP3c36cNbXCCf\nO0H4yJ//9We8WvQq9t23D2mxafzJiooLCX7BAuD3fwCqarkWrOEG3d7Sm6/e001ClY1WjPHFbqOz\nrLMpSv5rahRaUJUZkXU7kWMH5In/7YNv4/cFv8euu3ddSOwAkJYG7NkDPPccsGoVEN8PuDiD7670\ns5ywlmrL9CbjuOSenhJ30ZcHg6Lr+4RLi2cMOH2W14zvYU+CN4j+/Q8ESu5EWPHhkQ/x1KdPYefd\nO5EZn9n9hLQ04NlnL2xQMuiBgQOAwWkAmGJNQeSg64KsN0TqeU9YVTZaGQzc0VTTAJScAVrblL2+\n4JAsQ4QNu47vwl0f3oX8G97AmO1fAy+84JvU4i4n3MB92cbQ9yP4quvLhtUG2G1AUjx/hEmhNtLc\nCaIPCkoLcNvG2/DRtStww+1PAGfPAg8+CKxc6buW3toGnD3vUyEyOQiLwmMdcWnxkRHAoGT+zxCH\nNHeZEVm3Ezl2IDjx/6fiP/jRph/hnbG/u5DYAeD114GHH/ZdZomO4jVqkuJ5QSyLVZV67sHyuXtC\nS2sGbi3eauOlhGvq+/zsRP/+BwIldyKkOVp9FFP/NhUrs/8fptz56wuJ3YW198TsEZ0OSE4AMtOd\nhcgsfi/6+YsIhcdkIcLIH5XVQFml6g1BtArJMkTIcrr+NMa/NR7P3vgs7n96M/DJJ51PWLCAu2IC\n1W8Z4zr8uVp+rJBtUjN6uFowxm/OkPiidwgWISPNnSC6UNlUifFvjcdDYx/C4msWA+fOATfdBBw6\nxE8IVmLviMXKC5E1t3ItPsw34CiG3c6TfP84ICW0ipCR5i4zIut2IscO+Bd/bWstJm+YjLuuvIsn\ndgBITgY+/RS48kp5EjvAZ+yD0/hin90OtFtg/kJDmrUfaEpz7w29swhZfVM3y6To3/9ACH0vFxFW\nNFmacMu7t2DiRRPx6xt/3fnF5GRg3z4gLk4+K50k8RmkKZprwhYLT/QKzCbDzj3TEddiq80GnDoD\nJCUASWH0/98DJMsQIUO7rR23vncrBhuS8ObtG6DTqzx3YYyXsq04D9gd3LonoyasZts+Td1YHIw3\nRY+O4lq8wJZJkmWIsMfmsGHOB3MQb9VjzRN7oVv0U8XdK92QJF664KIMID6We7Rt8jk7AnXPBFLX\nXQlbptfonLN4i5XLNHWNwu4qDgRK7l4gsm4ncuyAd/E7mAPzt85HS905bPj1f6A/UwGsXQssXKh6\ngjebzbyEQVoy1+MlSbYSBv6UF+hITwnaW81dc7ZMSQIijDB/9W/es7W8StYbqxah5E4IDWMMi/MX\n40T5t/jwt98hsrziwotr1wJ/+5t6wXUlJpr74hP7c9kgiP7shQVPYsY/56PJ1uz3ewSSoAO9sciG\naxbf3Mo7PjU2h80snjR3Qmh+tedX2H7oQ+z5YzX6n6rs/OL8+cDq1dqsQ9LaDlSc45ufIiMCjjEY\nenvI++ZtTstkQj8gOVEIqypp7kRYsrxgObYc3oJPZn6E/omDOr+o5cQOANGRwJBB3JdttXE9PoAJ\nUDBkEU+lf0MCg57P4usauaOmtV3tiGRFo998bSGybi1y7EDv8a/5ag1W/nsldt29C8nplwK7dwOj\nR/MXNZTYPY6/Tsfr01yUzhtStFn87iEqlywip89dicbc3eJ3WSbtDp7gq+tCVqYhnzshHJv+uwnP\n7X0Oe+ftRUa/DP5kYiJP8G+8ATz1lCYSu9dEGIGMVK4HV1TzmbyPtsmu/VW1RG82SdcCLj9nibLx\nGw1cljlXy+2qA5P55xBCkOZOCMX2o9tx/9b7sfvu3bgy9Uq1wwk+NjtwvhaoawD0BkVrxsvlVe9t\nPSB391zsKN+DsUlZ6i3EMsZvpowBaUlAv1hN1achzZ0IC/aW7MW8D+/Gx5U/wJUDRqgdjjwY9EDa\nAOB7g/jMsrVNMTunXF713tYDNOGwcVomYTAAZ84DZ6r8aqeoRSi5e4HIurXIsQMX4v/yzJeYvfHH\n2LRFwjXLNwL33ee3Pq0kfo+/yVkzPkgLrl5dsockHAzNvbckrsQCrtfx63VAVATQ1AqcLOPWScGh\n5E5onm+rvsWt70zFmo8ZfvBVDX/ynXeESfB+E8QFV2+QayYtjAtHkpxrHTremLvS/8bcWqBPzT0/\nPx+LFy+G3W7HggULsHTp0k6vm81m3Hbbbbj44osBADNnzsQzzzzT+SKkuRN+cqL2BG5ccz1e+kcb\n7tpX1/lFuao7ahHGLiy4OuSvUyMKstW0cbX0izDyKp9RkcF5Xx8JJHd6XK2x2+145JFHsHv3bqSn\np2Ps2LGYPn06hg8f3um8G2+8EVu3bvUrAILojTONZzDpnUl4+mgK7tr3TecXwymxAzyR94vl1SZV\nWnDVIrI5bjq29Cs5A6QkAAn9hbqhevzLKCoqwtChQ5GZmQmj0Yg77rgDH3/8cbfzQn1WLrJuLWrs\n1S3VmPTOJNwk3YSHlu8Fvv/9Cy8KlNiDPv4KL7hqvZ57X5u3Ao7faHC29KsBTlfwYmSC4PG2X15e\njsGDB7uPMzIy8K9//avTOZIkoaCgAFlZWUhPT8fy5csxYkR3J8O8efOQmZkJAIiPj0d2djZycnIA\nXPgD0OrxwYMHNRVPqB9v37Udj3/yOGZMmYEphikwHzwI/OpXyHnhBWDUKJjnzAE++0wz8ap2PGEC\nUNsAc14+P75hPCBJ7oSWc811/PwQPn53wkrM+Ot8PNF/kVuSCfr1igoBBuSMuRooKYf52BHAFIWc\niRP560H8fM1mM9atWwcA7nzpLx419w8++AD5+flYs2YNAGDDhg3417/+hT//+c/ucxobG6HX62Ey\nmbBjxw489thjOHr0aOeLkOZOeEmrtRVT/zYVw5OH4y+5f4HU8WdwUxNgMgkxY1cUi5Uv/jW18mbd\nIdRmTnPY7Xy8+8cCqUmyj7VsPvf09HSUlpa6j0tLS5GRkdHpnLi4OJhMJgDA1KlTYbVaUVNT41cw\nRHhjsVswe/2tSI8dhJW5KzsndgCIjaXE3hOuHa7pyXxbvUwlheVGiXIEAaN31qdpaAFOlgMt2rVM\nevxLufrqq3Hs2DGUlJTAYrFg06ZNmD59eqdzKisr3XeWoqIiMMaQmJgoX8QqIKpuDYgTu91hxz3r\nfwTd5/ux7n0rdHauI4sSf28oFr9rwfXiDCA+LmglhZXU3OXYRCVL/JLEPfGQuGWyqkaTlkmPmrvB\nYMBrr72GyZMnw263Y/78+Rg+fDhWreIr0osWLcKWLVvw+uuvw2AwwGQyYePGjYoEToQOjDE8+O6d\nqPpiN/I2WGC0bQGg01YtdlFwLbj2i+Xt/dra+cxew794XHbGb+uKAWio4UdfGPR8Ubumnm96GpSs\nqZZ+VFuGUBXGGJZsXoDP9r6D3WutiLN0eHHNGu6MIfzD4QBqG3hxLNc2ew1a+TrWnskwDcSh2z7V\n/oanrlisfLxTEnm9+CCNs2w+d4KQm99tfRL5+9dj7zpb58S+YAFw//2qxRUSuHa4xsU4F1xbeILX\n2IKry84Ya4jB5f0vUTkaP4kw8uReWcNn8WkDVN+DoN3fahpCZN1Xy7G/VvQa/lryEXb+dzQSO65L\ndfCxazl+b9BE/O4F1xTngqv3dWqU0NzfnbASyZGJaLI1Y/fZz4NatExRn77OWZ+mtY0vtjao29KP\nZu6EKqz/z3q8vP9l7LtvHwbenwRMnQrs3y/UBiWh6LrDtbaBV0LUwA7X+Ij+uHpAlrv8ry96u0uv\nP95YgiGxGehnjA1uGQJfkSQgIoJbJsurnJbJRFV+LZHmTijOR0c+wkN5D+Gf9/wTw5OdpSwaG3lS\n/8UvKLErQUsbX3C1WDWx4Opv/9aOer0Lf3vIBh3GAIuF30QHpfDWij4SSO6k5E4oyu4Tu3HnB3ci\nf24+xgwco3Y44Y0gC66ecDX86G+MQ721sVvjD9kKi/mC1cZn8n4stlKzDpnRhG7qJ1qK/Yv9GzFn\n3TR88OONXid2LcXvD5qO37XgenEGn1W2tXcrKaz12jKuMsX/mb6rx3LFRV8elKUBiU+469NUc6lG\noWYglNwJRfhP4ceYsfUurH+vDeOf+DP/uUpogwgjMDiNSwc+LriqjatW/JDYwT3WjI/Uc9+56t55\nnY7vbG127Wxtk/2SJMsQsnPsy09w4+ZcrNjuwOzDzidnzAA2bwaModWUWHhsduB8DVDbqJkF10Dw\nV8uXFZuNSzXJCfyXkweZhnzuhGY5/fUeTNp4C174Z4fEDgADBmjOb03AucM1GegXJ8wOV0+4Zvaa\nwmAAdHrgXB2fwQ9MluUmKuYnpjCa1k37QM3Yq5qrMOmD2/DYfjvmH+jwgg92R5HHHhA4fmcPV/PR\nb7mjRoEernKg2TUDnbM+TWs7l2kam4N/iaC/I0EAqGurw+QNk3HHjY/g59ETL7xAPnZx0OkuFCPr\nZcGVCABXz1a9DiirBKqC27OVNHci6DRbmvHDDT/E2EFj8crkVyC1tADTpgGXXEKJXVQY4zsuK6sB\nBwMixbNNahpXz9aoSF6ALIKvRZHPndAM7bZ2THtvGtL7pWPt9LXQSc5E3toKREZSYhedEFtw1RwW\nK0/0AwcAcTGQdDryucuJsLoplI3ddvwY7nz/DvSL7Ic109ZcSOwAEB3tV2IXeeyBEIzfteA6xNnD\nta1dk7XMXWhWc++NCCO/aZZXAZXnA3oruu0SQcHxvyNY8H9XoSktBlv/7zsYdPTVCmmcC66obQDO\n1/HnBNzhqkn0Tk98XVNAb0OyDBEw7H//w2PPXIWv+7fgk3eAmEm5wAcfAFFRaodGKIHFyrsRNbZw\nmcZAFtdgIV2UTrIMoRLFxXjuybHYN6AF/3gXiLECyMsD1q5VOzJCKSKMvJzw4FRAgualmnCBkrsX\niKybyhr7yZP446NjsSmzCZ9sAOJdO6oXLAAefDAolxB57IEwil+SgFgTcFE633lptWnCGy+c5h5E\nKLkTfvNm5Q68OsaKXeuBFNceDPKxhzeuYmQXpQMx0bxOjUKFsojOkOZO+MXmbzfj55/8HOaf5GPY\n/U8AO3dSYic6wxhvOVdZzWfyApcxUItANHdK7oTP5B3Lw30f34ddd+/CqNRR3MO+Zg3wyCP0x0t0\nx1U3nlw1PkMLqjIjsm4a7Ng/O/UZ5v19Hj6+42Oe2AHuYX/0UVkSu8hjD1D8ALpLNe3KSTXhrLmT\nGZnwjuJifPX7n2PWZV/ivZnv4ZqMa9SOiBANl6vGJdUIXnFS65AsQ/RNcTEO/+gG/OCW83jj9CjM\nWFfIZ+sE4S9uqaYWgJgt/pSAZBlCPoqLcXL6eEyech6/3wnM2PwNMH060CZ/JxkihHFLNRmKSzXh\nAiV3LxBZNw0o9uJinL1lAm6ecg5PfQ7c/Y3z+cxMICIiCNH1jchjD1D8feLeAJUmywaocNbcKbkT\nvVK9/P8waXIV5h8AHv6380myOxLBRpL47N29AUrc5iBagjR3okca2xtx09s/wMSva/DSqhOQAErs\nhDK4atU0tfA6NYbw9X2Qz50IKq3WVuS+m4vLki7D6zf/CdKsWcDAgZTYCeVgjPcXrTgf1hugaEFV\nZkTWTX2N3Wq34vYtt2Ng7ECszF0JKSoK+PBD1RK7yGMPUPx+EySpRhTNfWHBk8jJn4Xc3XNRZ6kP\nyntScic4x4/D3tiAe/9+LwDg7RlvQ69zlm6NiAjLWROhAXp01djUjiroHG04ib2VX2BH+R4sLFgS\nlPckWYYAiovBJubgwVslFF8zFHl3f4JoI/nYCY0RwlJN7u652FG+B2OTsrDzh+8hPqI/AJJliEBw\nJvanRlbga5zF1r85EG2hWtyEBglhV827E1Zi9pBbOyX2QKHk7gUi66YeYy8uBiZOxEuXVGD7pcCO\nDUDcP/cDq1crFl9fiDz2AMUvCz5INaJo7vER/bE5Z1XQEjtAyT18OXUKmDgRf8k4izfHADvfAZJa\nwe2Ojz2mdnQE0TedNkBJ1AGqC6S5hytWKzYsuhZPx3+Fz94CLqoD+dgJcQnRWjWkuRM+8/HxPDwx\nogyfnL2JEjshPt2kmvaQdNX4Av0le4EmdUcv6Sn2T098ige2PYDtd27HiPU7gDfe0GxiF3nsAYpf\ncdxSzUBAkmD+bG/YSjXa+2vWIAcPHlQ7BL/pGnthWSHmfDAHW27fgqsGXQUYjcCiRZpM7IDYYw9Q\n/KrQwVVzsLQkpFw1vtDnX3R+fj4uv/xyDBs2DC+//HKP5zz66KMYNmwYsrKycODAgaAHqTZ1dXVq\nh+A3dXV13BUzaxa+OV6A2zbehnUz1mHCkAlqh+YVIo89QPGrik6HOmt72Eo1HpO73W7HI488gvz8\nfBw+fBjvvfcejhw50umcvLw8fPfddzh27BhWr16NBx98UNaACR+prgYmTsSxPR9g6psT8eecZcgd\nlqt2VAShHF2kmnBx1XhM7kVFRRg6dCgyMzNhNBpxxx134OOPP+50ztatW3HvvXzL+rhx41BXV4fK\nykr5IlaBkpIStUPwm5Jt21DafBaT7gGe32nB7Y+tBpqb1Q7La0Qee4DiVxt3/CG8AapXmAfef/99\ntmDBAvfxO++8wx555JFO59x6661s//797uObbrqJffnll53OAUAPetCDHvTw4+EvHgslS176RLv6\nMLv+d+RxJwiCUBaPskx6ejpKS0vdx6WlpcjIyPB4TllZGdLT04McJkEQBOELHpP71VdfjWPHjqGk\npAQWiwWbNm3C9OnTO50zffp0rF+/HgBQWFiI+Ph4pKamyhcxQRAE0SceZRmDwYDXXnsNkydPht1u\nx/z58zF8+HCsWrUKALBo0SLk5uYiLy8PQ4cORUxMDN566y1FAicIgiA84Lda74HNmzezESNGMJ1O\nx7766qtez9uxYwe77LLL2NChQ9lLL70kRyh+UV1dzW6++WY2bNgwNmnSJFZbW9vjeUOGDGFXXnkl\ny87OZmPHjlU4yu54M54/+9nP2NChQ9moUaPY119/rXCEnukr/j179rB+/fqx7Oxslp2dzV544QUV\nouyZ++67j6WkpLArrrii13O0OvZ9xa7lcWeMsdOnT7OcnBw2YsQINnLkSLZixYoez9Pq+HsTvz+f\ngSzJ/ciRI6y4uJjl5OT0mtxtNhu75JJL2MmTJ5nFYmFZWVns8OHDcoTjM08++SR7+eWXGWOMvfTS\nS2zp0qU9npeZmcmqq6uVDK1XvBnP7du3s6lTpzLGGCssLGTjxo1TI9Qe8Sb+PXv2sGnTpqkUoWc+\n++wz9vXXX/eaILU89n3FruVxZ4yxs2fPsgMHDjDGGGtsbGSXXnqpUN99b+L35zOQZc/55Zdfjksv\nvdTjOd546NWio3f/3nvvxd///vdez2UacQKJvifB2++DVsa7K+PHj0dCQkKvr2t57PuKHdDuuANA\nWloasrOzAQCxsbEYPnw4zpw50+kcLY+/N/EDvn8GqhUUKS8vx+DBg93HGRkZKC8vVyucTlRWVroX\nhVNTU3v9EkiShJtvvhlXX3011qxZo2SI3fBmPHs6p6ysTLEYPeFN/JIkoaCgAFlZWcjNzcXhw4eV\nDtNvtDz2fSHSuJeUlODAgQMYN25cp+dFGf/e4vfnM/C4oOqJSZMmoaKiotvzL774IqZNm9bnf++t\nh14ueov/t7/9badjSZJ6jXX//v0YOHAgzp07h0mTJuHyyy/H+PHjZYm3L4K1J0EtvIljzJgxKC0t\nhclkwo4dOzBjxgwcPXpUgeiCg1bHvi9EGfempibMmjULK1asQGxsbLfXtT7+nuL35zPwO7nv2rXL\n3/8UgHceejnxFH9qaioqKiqQlpaGs2fPIiUlpcfzBg4cCABITk7Gj370IxQVFamW3EXfk+BN/HFx\nce5/nzp1Kh566CHU1NQgMTFRsTj9Rctj3xcijLvVasXMmTMxd+5czJgxo9vrWh//vuL35zOQXZbp\nTSfyxkOvFtOnT8fbb78NAHj77bd7HOyWlhY0NjYCAJqbm7Fz505ceeWVisbZEdH3JHgTf2Vlpfv7\nVFRUBMaYphKMJ7Q89n2h9XFnjGH+/PkYMWIEFi9e3OM5Wh5/b+L36zPwe4nXAx9++CHLyMhgUVFR\nLDU1lU2ZMoUxxlh5eTnLzc11n5eXl8cuvfRSdskll7AXX3xRjlD8orq6mt10003drJAd4z9+/DjL\nyspiWVlZbOTIkZqIv6fxfOONN9gbb7zhPufhhx9ml1xyCRs1apRHm6oa9BX/a6+9xkaOHMmysrLY\ntddey7744gs1w+3EHXfcwQYOHMiMRiPLyMhga9euFWbs+4pdy+POGGP79u1jkiSxrKwst1UwLy9P\nmPH3Jn5/PgNFeqgSBEEQyqLN9jsEQRBEQFByJwiCCEEouRMEQYQglNwJgiBCEEruBEEQIQgld4Ig\niBDk/wOkm9bO2t6uJAAAAABJRU5ErkJggg==\n", "text": [ "" ] } ], "prompt_number": 19 }, { "cell_type": "markdown", "metadata": {}, "source": [ "The figure above shows the green dots, which are individual estimates of $\\boldsymbol{\\hat{\\beta}}$ from the corresponding simulated $\\mathbf{y}$ data. The red dashed line is the true value for $\\boldsymbol{\\beta}$ and the green line ($\\boldsymbol{\\hat{\\beta_m}}$ ) is the average of all the green dots. Note there is hardly a visual difference between them. The pink ellipse provides some scale as to the covariance of the estimated $\\boldsymbol{\\beta}$ values. I invite you to download this [IPython notebook](https://github.com/unpingco/Python-for-Signal-Processing/blob/master/Gauss_Markov.ipynb) and tweak the values for the error covariance or any other parameter. All of the plots should update and scale correctly.\n" ] }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "Summary" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The Gauss-Markov problem is the cornerstone of all Gaussian modeling and is thus one of the most powerful models used in signal processing. We showed how to estimate the unobservable parameters given noisy measurements using our previous work on projection. We also coded up a short example to illustrate how this works in a simulation. For a much more detailed approach, see Luenberger (1997).\n", "\n", "As usual, the corresponding IPython notebook for this post is available for download [here](https://github.com/unpingco/Python-for-Signal-Processing/blob/master/Gauss_Markov.ipynb). \n", "\n", "Comments and corrections welcome!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "References\n", "---------------\n", "\n", "* Luenberger, David G. *Optimization by vector space methods*. Wiley-Interscience, 1997." ] } ], "metadata": {} } ] }