Conditional_expectation_MSE_Ex.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:9b67a4eac3f1e33cbbfac5318e6bc4298d31615f0abce086d6000ea25b106389"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Introduction\n",
      "-------------------\n",
      "\n",
      "Brzezniak (2000) is a great book because it approaches conditional expectation through a sequence of exercises, which is what we are trying to do here. The main difference is that Brzezniak takes a more abstract measure-theoretic approach to the same problems. Note that you *do* need to grasp the measure-theoretic to move into more advanced areas in stochastic processes, but for what we have covered so far, working the same problems in his text using our methods is illuminating. It always helps to have more than one way to solve *any* problem.  I urge you to get a copy of his book or at least look at some pages on Google Books. I have numbered the examples corresponding to the book. "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Examples\n",
      "-------------\n",
      "\n",
      "This is Example 2.1 from Brzezniak:\n",
      "\n",
      "> Three coins, 10p, 20p and 50p are tossed. The values of those coins that land heads up are added to work out the total amount. What is the expected total amount given that two coins have landed heads up?\n",
      "\n",
      "In this case we have we want to compute $\\mathbb{E}(\\xi|\\eta)$ where\n",
      "\n",
      "$$ \\xi = 10 X_{10} + 20 X_{20} +50 X_{50} $$\n",
      "\n",
      "where $X_i \\in \\{ 0,1\\} $. This represents the sum total value of the heads-up coins. The $\\eta$ represents the fact that only two of the three coins are heads-up. Note \n",
      "\n",
      "$$\\eta = X_{10} X_{20} (1-X_{50})+ (1-X_{10}) X_{20} X_{50}+ X_{10} (1-X_{20}) X_{50} $$\n",
      "\n",
      "is a function that is only non-zero when two of the three coins is heads. Each triple term catches each of these three possibilities. For example, the first term is when the 10p and 20p are heads up  and the 50p is heads down.\n",
      "\n",
      "To compute the conditional expectation, we want to find a function $h$ of $\\eta$ that minimizes the MSE\n",
      "\n",
      "$$ \\sum_{X\\in\\{0,1\\}^3} \\frac{1}{8} (\\xi - h( \\eta ))^2 $$\n",
      "\n",
      "where the sum is taken over all possible triples of outcomes for $ \\{X_{10} , X_{20} ,X_{50}\\}$ and the $\\frac{1}{8} = \\frac{1}{2^3} $ since each coin has a $\\frac{1}{2}$ chance of coming up heads.\n",
      "\n",
      "Now, the question boils down to what function $h(\\eta)$ should we try? Note that $\\eta \\in \\{0,1\\}$ so $h$ takes on only two values. Thus, we only have to try $h(\\eta)=\\alpha \\eta$  and find $\\alpha$. Writing this out gives,\n",
      "\n",
      "$$ \\sum_{X\\in\\{0,1\\}^3} \\frac{1}{8} (\\xi - \\alpha( \\eta ))^2 $$\n",
      "\n",
      "\n",
      "which boils down to solving for $\\alpha$,\n",
      "\n",
      "$$\\langle \\xi , \\eta \\rangle = \\alpha \\langle \\eta,\\eta \\rangle$$\n",
      "\n",
      "where \n",
      "\n",
      "$$ \\langle \\xi , \\eta \\rangle =\\sum_{X\\in\\{0,1\\}^3} \\frac{1}{8} (\\xi  \\eta ) $$\n",
      "\n",
      "This is tedious and a perfect job for `sympy`.\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import sympy as S\n",
      "\n",
      "eta = S.Symbol('eta')\n",
      "xi = S.Symbol('xi')\n",
      "X10 = S.Symbol('X10')\n",
      "X20 = S.Symbol('X20')\n",
      "X50 = S.Symbol('X50')\n",
      "\n",
      "eta = X10 * X20 *(1-X50 )+ X10 * (1-X20) *(X50 )+ (1-X10) * X20 *(X50 )\n",
      "xi = 10*X10 +20* X20+ 50*X50 \n",
      "\n",
      "num=S.summation(xi*eta,(X10,0,1),(X20,0,1),(X50,0,1))\n",
      "den=S.summation(eta*eta,(X10,0,1),(X20,0,1),(X50,0,1))\n",
      "alpha=num/den\n",
      "print alpha"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "160/3\n"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This means that\n",
      "\n",
      "$$ \\mathbb{E}(\\xi|\\eta) = \\frac{160}{3} \\eta $$\n",
      "\n",
      "which we can check with a quick simulation"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "from numpy import array\n",
      "\n",
      "x=np.random.randint(0,2,(3,5000))\n",
      "print (160/3.,np.dot(x[:,x.sum(axis=0)==2].T,array([10,20,50])).mean())\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(53.333333333333336, 53.243528790279981)\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Example\n",
      "--------\n",
      "\n",
      "This is example 2.2:\n",
      "\n",
      "> Three coins, 10p, 20p and 50p are tossed as before. What is the conditional expectation of the total amount shown by the three coins given the total amount shown by the 10p and 20p coins only?\n",
      "\n",
      "For this problem,\n",
      "\n",
      "$$\\eta = 30 X_{10} X_{20} + 20 (1-X_{10}) X_{20}  + 10 X_{10} (1-X_{20})  $$\n",
      "\n",
      "which takes on three values (10,20,30) and only considers the 10p and 20p coins. Here, we'll look for affine functions, $h(\\eta) = a \\eta + b $."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from sympy.abc import a,b\n",
      "\n",
      "eta = X10 * X20 * 30 + X10 * (1-X20) *(10 )+ (1-X10) * X20 *(20 )\n",
      "h = a*eta + b\n",
      "\n",
      "J=S.summation((xi - h)**2 * S.Rational(1,8),(X10,0,1),(X20,0,1),(X50,0,1))\n",
      "\n",
      "sol=S.solve( [S.diff(J,a), S.diff(J,b)],(a,b) )\n",
      "print sol"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "{b: 25, a: 1}\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This means that\n",
      "\n",
      "$$ \\mathbb{E}(\\xi|\\eta) = 25+ \\eta $$\n",
      "\n",
      "since $\\eta$ takes on only four values, $\\{0,10,20,30\\}$, we can write this out as\n",
      "\n",
      "$$ \\mathbb{E}(\\xi|\\eta=0) = 25  $$\n",
      "$$ \\mathbb{E}(\\xi|\\eta=10) = 35  $$\n",
      "$$ \\mathbb{E}(\\xi|\\eta=20) = 45  $$\n",
      "$$ \\mathbb{E}(\\xi|\\eta=30) = 55  $$\n",
      "\n",
      "The following is  a quick simulation to demonstrate this."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x=np.random.randint(0,2,(3,5000))  # random samples for 3 coins tossed\n",
      "eta=np.dot(x[:2,:].T,array([10,20])) # sum of 10p and 20p\n",
      "\n",
      "print np.dot(x[:,eta==0].T,array([10,20,50])).mean() # E(xi|eta=0)\n",
      "print np.dot(x[:,eta==10].T,array([10,20,50])).mean()# E(xi|eta=10)\n",
      "print np.dot(x[:,eta==20].T,array([10,20,50])).mean()# E(xi|eta=20)\n",
      "print np.dot(x[:,eta==30].T,array([10,20,50])).mean()# E(xi|eta=30)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "25.1587301587\n",
        "34.3410852713\n",
        "43.7323358271\n",
        "56.6238973536\n"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Example\n",
      "-----------\n",
      "\n",
      "This is Example 2.3\n",
      "\n",
      "![alt text](files/ex23.jpg)"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Note that \"Lebesgue measure\" on $[0,1]$ just means uniformly distributed on that interval. Also, note the the `Piecewise` object in `sympy` is not complete at this point in its development, so we'll have to work around that in the following."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%pylab inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Populating the interactive namespace from numpy and matplotlib\n"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x=S.Symbol('x')\n",
      "c=S.Symbol('c')\n",
      "\n",
      "xi = 2*x**2\n",
      "\n",
      "eta=S.Piecewise((1,S.And(S.Gt(x,0),S.Lt(x,S.Rational(1,3)))),  #  0 < x < 1/3\n",
      "                (2,S.And(S.Gt(x,S.Rational(1,3)),S.Lt(x,S.Rational(2,3)))), # 1/3 < x < 2/3,\n",
      "                (0,S.And(S.Gt(x,S.Rational(2,3)),S.Lt(x,1))), \n",
      "                )\n",
      "\n",
      "\n",
      "h = a + b*eta + c*eta**2 \n",
      "J=S.integrate((xi - h)**2 ,(x,0,1))\n",
      "\n",
      "sol=S.solve( [S.diff(J,a), \n",
      "              S.diff(J,b),\n",
      "              S.diff(J,c),\n",
      "              ],\n",
      "              (a,b,c) )\n",
      "print sol\n",
      "\n",
      "print S.piecewise_fold(h.subs(sol))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "{c: 8/9, b: -20/9, a: 38/27}\n",
        "Piecewise((2/27, And(x < 1/3, x > 0)), (14/27, And(x < 2/3, x > 1/3)), (38/27, And(x < 1, x > 2/3)))\n"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Thus, collecting this result gives:\n",
      "\n",
      "$$ \\mathbb{E}(\\xi|\\eta) = \\frac{38}{27} - \\frac{20}{9}\\eta + \\frac{8}{9} \\eta^2$$\n",
      "\n",
      "which can be re-written as a piecewise function as\n",
      "\n",
      "$$\\mathbb{E}(\\xi|\\eta) =\\begin{cases} \\frac{2}{27} & \\text{for}\\: 0 < x < \\frac{1}{3} \\\\\\frac{14}{27} & \\text{for}\\: \\frac{1}{3} < x < \\frac{2}{3} \\\\\\frac{38}{27} & \\text{for}\\: \\frac{2}{3}<x < 1 \\end{cases}\n",
      "$$\n",
      "\n",
      "The following is a quick simulation to demonstrate this."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = np.random.rand(1000)\n",
      "f,ax= subplots()\n",
      "ax.hist(2*x**2,bins=array([0,1/3.,2/3.,1])**2*2,normed=True,alpha=.5)\n",
      "ax.vlines([2/27.,14/27.,38/27.],0,ax.get_ylim()[1],linestyles='--')\n",
      "ax.set_xlabel(r'$2 x^2$',fontsize=18);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEaCAYAAADqqhd6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFXZJREFUeJzt3X+s3Xd93/HnC9esRSUKwVumOoYziEeTlRA6ETpo68Ng\nwkSl0RhaZkqnsAJhkxmbNBGGpvpaCCQmVYsgVepmSUQ6EVcNiJkpP4Q0ThYhCI2WnyNmdoM7OzBv\nQMooYao9v/fHPdyenPie873nfq/vOV8/H9KV7vd8P/d73/nok9f9+nu+5/1NVSFJ6pYXbHYBkqT2\nGe6S1EGGuyR1kOEuSR1kuEtSBxnuktRBU8M9yW1JTiZ5fJX925Lcm+SRJE8kua71KiVJa9LkzP12\nYPeE/XuBh6vqSqAP/E6Sn2qhNknSjKaGe1U9ADwzYch3gAuG318AfK+qTrdQmyRpRm2cYd8C/Ock\n3wZeDPzDFo4pSVqHNsL9o8AjVdVP8krgS0leU1U/HB2UxD4HkjSDqspaf6aNu2XeAPzRsIA/Ab4F\nvOpsA6vKr5a+9u3bt+k1dOXLuXQ+5/lrVm2E+2HgLQBJLmY52J9q4biSpBlNvSyT5E5gF7AtyXFg\nH7AVoKoOAJ8Abk/yKMt/LD5cVd/fuJIlSdNMDfeq2jNl/3eBt7dWkRrp9/ubXUJnOJftcj7nQ9Zz\nTWdNvyipc/W7JKkrklCb9IaqJGnOGO6S1EGGuyR1kOEuSR1kuEtSBxnuktRBhrskdZDhLkkdZLhL\nUgcZ7pLUQYa7JHWQ4S5JHWS4S1IHGe6S1EGGuyR10NRwT3JbkpNJHp8wpp/k4SRPJBm0WqEkac2m\nPqwjya8Afw7cUVWvPsv+C4GvAG+tqhNJtg2fzjQ+zod1SNIabdjDOqrqAeCZCUPeBXyuqk4Mxz8v\n2CVJ51Yb19x3Ahcl+XKSh5L8ZgvHlCStw9QHZDewFfhF4M3Ai4CvJvlaVR0ZH7i0tLTyfb/f90G6\nkjRmMBgwGAzWfZxGD8hO0gO+uMo19xuAn6mqpeH2vwfuraq7xsat6Zr70tLSc/4YaHXOldrmmpof\ns15zbyPcfx64CXgr8FeAB4Frq+obY+PWFO7D/6DG489nzpXa5pqaH7OG+9TLMknuBHYB25IcB/ax\nfCmGqjpQVYeT3As8BpwBbhkPdknSudXozL2VX+SZ+4ZxrtQ219T82LBbISVJi8dwl6QOmttw37dv\n32aXsDCcK7XNNbX45vaauyTJa+6SpBGGuyR1kOEuSR1kuEtSB81tuNvXojnnSm1zTS2+ub1bxk/I\nNedcqW2uqfnh3TKSpBWGuyR1kOEuSR1kuEtSB81tuNvbojnnSm1zTS2+ub1bRpK0gXfLJLktyckk\nj08Z97okp5O8Y61FSJLaNfUxe8DtwKeBO1YbkGQL8EngXmDVvzA33vgHa61vQ1155d+g3//lzS5D\nklo3Ndyr6oHhA7In+SBwF/C6SYOeeurvNC5soz3zzLe48MKTm12GJG2IJmfuEyXZDlwD/F2Ww33V\nC+sXXXTpen9da06d+jHwfza7DEnaEOsOd+BG4CNVVUnChMsyg8HSyve9Xp9er7/qQQeDJfr9pVX3\n6y8tLS3ZC0Stck1tnsFgwGAwWPdxGt0tM7ws88WqevVZ9j3FXwb6NuBZ4H1VdWhsXO3b1/xumf37\nw1rGr9XJk4/z+tf/d6677h9s2O84V+wDora5pubHrHfLrPvMvapeMVLE7Sz/ETg04UckSRtsargn\nuRPYBWxLchzYB2wFqKoDG1ueJGkWTe6W2dP0YFX1nvWVI0lqw9y2H5AkzW5uw33XLntbNGUfELXN\nNbX4zmlvmY28+2WtunS3jKTu8klMkqQVhrskdZDhLkkdZLhLUgfNbbiP9qHRZPYAUdtcU4tvbsP9\n/vv3b3YJC2P/fudK7XJNLb65DXdJ0uwMd0nqIMNdkjrIcJekDprbcLe3THP2AVHbXFOLz94y9paR\nNMfsLSNJWjE13JPcluRkksdX2f8bSR5N8liSryS5ov0yJUlr0eTM/XZg94T9TwG/WlVXAB8Dfr+N\nwiRJs5sa7lX1APDMhP1fraofDDcfBC5pqTZJ0ozavub+W8DdbRzI3jLN2QdEbXNNLb5Gd8sk6QFf\nrKpXTxjzJuB3gTdW1fPO9JPU6O2NvV6fXq+/6u/cvz9s5N01XbpbZvhu+maXoQ5xTW2ewWDAYDBY\n2d6/f/9Md8v8VBvFDN9EvQXYfbZg/4l+f6mNXydJndXv9+n3+yvbszZxW/dlmSQvAz4PvLuqjq73\neJKk9Zt65p7kTmAXsC3JcWAfsBWgqg4Avw28BLg5CcCpqrpqwyqWJE01Ndyras+U/e8F3ttaRZKk\ndZvbT6jaW6Y5+4Coba6pxWdvmQ7cLSOpu+wtI0laYbhLUgcZ7pLUQYa7JHXQ3Ia7vWWasw+I2uaa\nWnxzG+733z/bR27PR7N+PFlajWtq8c1tuEuSZme4S1IHGe6S1EGGuyR10NyGu71lmrMPiNrmmlp8\n9paxt4ykOWZvGUnSiqnhnuS2JCeTPD5hzKeSHEnyaJLXtluiJGmtmpy53w7sXm1nkquBS6tqJ/B+\n4OaWapMkzWhquFfVA8CqD70Gfh34zHDsg8CFSS5upzxJ0izauOa+HTg+sn0CuGS9B7W3THP2AVHb\nXFOLr603VMffyV33bTH2lmnOPiBqm2tq8U19QHYDTwM7RrYvGb72PKNn471en16v38Kvl6TuGAwG\nDAaDdR+njXA/BOwFDib5JeDPqurk2Qb2+0st/DpJ6q5+v0+/31/ZnvVfUVPDPcmdwC5gW5LjwD5g\nK0BVHaiqu5NcneQo8CPgPTNVIklqzdRwr6o9DcbsbaccSVIb5vYTqvaWac4+IGqba2rx2VvG3jKS\n5pi9ZSRJKwx3Seogw12SOshwl6QOmttwt7dMc/YBUdtcU4tvbsPd3jLN2QdEbXNNLb65DXdJ0uwM\nd0nqIMNdkjrIcJekDprbcLe3THP2AVHbXFOLz94y9paRNMfsLSNJWmG4S1IHTQ33JLuTHE5yJMkN\nZ9m/Lcm9SR5J8kSS6zakUklSYxPDPckW4CZgN3A5sCfJZWPD9gIPV9WVQB/4nSRtPJtVkjSjaWfu\nVwFHq+pYVZ0CDgLXjI35DnDB8PsLgO9V1en1FmZvmebsA6K2uaYW37Rw3w4cH9k+MXxt1C3A30ry\nbeBR4ENtFGZvmebsA6K2uaYW37TLJ03uXfwo8EhV9ZO8EvhSktdU1Q/HB46ejfd6fXq9/hpKlaTu\nGwwGDAaDdR9nWrg/DewY2d7B8tn7qDcAHweoqj9J8i3gVcBD4wfr95dmLlSSzgf9fp9+v7+yPeu/\noqZdlnkI2Jmkl+SFwLXAobExh4G3ACS5mOVgf2qmaiRJrZh45l5Vp5PsBe4DtgC3VtWTSa4f7j8A\nfAK4PcmjLP+x+HBVfX+D65YkTTD1lsWquge4Z+y1AyPffxd4e9uF2VumOfuAqG2uqcVnbxl7y0ia\nY/aWkSStMNwlqYMMd0nqIMNdkjpobsPd3jLN2QdEbXNNLb65DXd7yzRnHxC1zTW1+OY23CVJszPc\nJamDDHdJ6iDDXZI6aG7D3d4yzdkHRG1zTS0+e8vYW0bSHLO3jCRpheEuSR1kuEtSB00N9yS7kxxO\nciTJDauM6Sd5OMkTSQatVylJWpOJ4Z5kC3ATsBu4HNiT5LKxMRcCvwu8vap+AXhnG4XZW6Y5+4Co\nba6pxTftzP0q4GhVHauqU8BB4JqxMe8CPldVJ2DlsXvrZm+Z5uwDora5phbftHDfDhwf2T4xfG3U\nTuCiJF9O8lCS32yzQEnS2k17QHaTG9O3Ar8IvBl4EfDVJF+rqiPjA0cvtfR6fXq9fuNCJel8MBgM\nGAwG6z7OtHB/Gtgxsr2D5bP3UceB71bVj4EfJ/kvwGuA54V7v780e6WSdB7o9/v0+/2V7VkvkU27\nLPMQsDNJL8kLgWuBQ2Nj/iPwy0m2JHkR8HrgGzNVI0lqxcQz96o6nWQvcB+wBbi1qp5Mcv1w/4Gq\nOpzkXuAx4AxwS1WtO9ztLdOcfUDUNtfU4rO3jL1lJM0xe8tIklYY7pLUQYa7JHWQ4S5JHTS34W5v\nmebsA6K2uaYW39yGu71lmrMPiNrmmlp8cxvukqTZTWs/0Gmf/eznGQwe3+wyWnHddUubXcI5d+GF\ncOONS5tdhjSXzutwf/bZrfR6S5tdRgv2d+S/Y22OHVva7BKkueVlGUnqoLkNd3vLNOdcqW32lll8\ncxvutgduzrlS27wVcvHNbbhLkmZnuEtSBxnuktRBhrskddDUcE+yO8nhJEeS3DBh3OuSnE7yjjYK\ns7dMc86V2uYbqotvYrgn2QLcBOwGLgf2JLlslXGfBO4F1vzEkLOxt0xzzpXaZm+ZxTftzP0q4GhV\nHauqU8BB4JqzjPsgcBfwv1uuT5I0g2nhvh04PrJ9YvjaiiTbWQ78m4cvzc+DUiXpPDWtt0yToL4R\n+EhVVZIw4bLM6LXhXq9Pr9dvcHhJOn8MBgMGg8G6jzMt3J8Gdoxs72D57H3U3wYOLuc624C3JTlV\nVYfGD+YnKSVpsn6/T7/fX9me9f2PaeH+ELAzSQ/4NnAtsGd0QFW94iffJ7kd+OLZgn2t7JfSnHOl\nttlbZvFNDPeqOp1kL3AfsAW4taqeTHL9cP+BjSrMs/zmnCu1zVshF9/Ufu5VdQ9wz9hrZw31qnpP\nS3VJktbBT6hKUgcZ7pLUQYa7JHXQ3Ia7/VKac67UNt9QXXxzG+72S2nOuVLb7C2z+OY23CVJs5t6\nK6Q0r77+9a9x3XVLm11GZzm3i81w18L6i7/4aXq9pc0uo6P2O7dzY7ZLZF6WkaQOmttwt19Kc86V\n2uaaWnxzG+72S2nOuVLbXFOLb27DXZI0O8NdkjrIcJekDjLcJamD5jbc7ZfSnHOltrmmFl+jcE+y\nO8nhJEeS3HCW/b+R5NEkjyX5SpIr1luY/VKac67UNtfU4psa7km2ADcBu4HLgT1JLhsb9hTwq1V1\nBfAx4PfbLlSS1FyTM/ergKNVdayqTgEHgWtGB1TVV6vqB8PNB4FL2i1TkrQWTcJ9O3B8ZPvE8LXV\n/BZw93qKkiStT5PGYdX0YEneBPwT4I1n2z/6Jk2v16fX6zc9tCSdF44dG3Ds2GDdx2kS7k8DO0a2\nd7B89v4cwzdRbwF2V9UzZzvQWj7SbG+L5pwrtc01tXnGT3xnfXO7yWWZh4CdSXpJXghcCxwaHZDk\nZcDngXdX1dGZKhljb4vmnCu1zTW1+KaeuVfV6SR7gfuALcCtVfVkkuuH+w8Avw28BLg5CcCpqrpq\n48qWJE3S6GEdVXUPcM/YawdGvn8v8N52S5MkzWpuP6EqSZqd4S5JHTS34W5vi+acK7XNNbX45jbc\n7W3RnHOltrmmFt/chrskaXaGuyR1kOEuSR1kuEtSB81tuNvbojnnSm1zTS2+uQ13e1s051ypba6p\nxTe34S5Jmp3hLkkdZLhLUgcZ7pLUQXMb7va2aM65UttcU4tvargn2Z3kcJIjSW5YZcynhvsfTfLa\nNgqzt0VzzpXa5ppafBPDPckW4CZgN3A5sCfJZWNjrgYuraqdwPuBmzeoVo1o4wG6WuZctsv5nA/T\nztyvAo5W1bGqOgUcBK4ZG/PrwGcAqupB4MIkF7deqZ7D/4Ha41y2y/mcD9PCfTtwfGT7xPC1aWMu\nWX9pkqRZTXuGajU8Tpr83J/+6ScaHm628Wtx5sz/4wVz+3ayJK1PqlbP7yS/BCxV1e7h9r8GzlTV\nJ0fG/B4wqKqDw+3DwK6qOjl2rKZ/KCRJI6pq/AR6qmln7g8BO5P0gG8D1wJ7xsYcAvYCB4d/DP5s\nPNhnLU6SNJuJ4V5Vp5PsBe4DtgC3VtWTSa4f7j9QVXcnuTrJUeBHwHs2vGpJ0kQTL8tIkhZT628p\nbtaHnrpo2lwm6Sf5QZKHh1//ZjPqXARJbktyMsnjE8a4LhuaNp+uzbVJsiPJl5P8tyRPJPnnq4xr\nvkarqrUvli/dHAV6wFbgEeCysTFXA3cPv3898LU2a+jKV8O57AOHNrvWRfgCfgV4LfD4Kvtdl+3O\np2tzbfP514Erh9//LPDN9WZn22fufuipPU3mEp5/G6rOoqoeAJ6ZMMR1uQYN5hNcm41V1f+sqkeG\n3/858CTwc2PD1rRG2w53P/TUniZzWcAbhv9EuzvJ5eesuu5xXbbLtTmj4d2JrwUeHNu1pjU67VbI\ntWr1Q0/nuSZz8l+BHVX1bJK3AV8A/ubGltVprsv2uDZnkORngbuADw3P4J83ZGx71TXa9pn708CO\nke0dLP91mTTmkuFreq6pc1lVP6yqZ4ff3wNsTXLRuSuxU1yXLXJtrl2SrcDngP9QVV84y5A1rdG2\nw33lQ09JXsjyh54OjY05BPxjWPkE7Fk/9KTpc5nk4iQZfn8Vy7e2fv/cl9oJrssWuTbXZjhXtwLf\nqKobVxm2pjXa6mWZ8kNPrWkyl8A7gX+a5DTwLPCPNq3gOZfkTmAXsC3JcWAfy3chuS5nMG0+cW2u\n1RuBdwOPJXl4+NpHgZfBbGvUDzFJUgfZF1GSOshwl6QOMtwlqYMMd0nqIMNdkjrIcJekDmq7/YDU\nOUneAbyc5U58T1bV/k0uSZrKcJcmSPJK4MKq+ndJfhr4ZpIjVfXZza5NmsTLMtJkvwDsB6iq/wt8\nneVPE0pzzXCXJrsbeNvI9iUs99qW5prtB9RJww57H2Q5jF/O8oMPPlVVd67jmFcCf8jyE3N+3Eqh\n0gbxmru6ah/wB1X1TYAkvwYcSrKtqj691oMl+RmWL8+81WDXIvDMXZ2T5MXA/wI+U1UfGHn968Ar\nq+qlMxzz48DvVdXxJJdW1dH2Kpba5zV3ddEZ4DvAi8defwp4SZK/upaDJfkA8J+AU0m2A29ppUpp\nA3lZRp1TVT8CXnGWXZcC3we+95MXkrwf2Ab8PHAHy9fn/xrLd8ncAPSAm3juidA7N6JuqU2Gu84L\nSa5g+aHD/7Kqzgxfex/wcFX9cZLXAV8CrgP+B/Bx4I6qug//P9ECctGq85K8APg0cFdVfWpk10ur\n6o+H378cOFNVXxi+ebqrqh4417VKbfENVXVekn8LvAR4f62y4JPcBGyvqr9/TouTNohvqKrTkvwL\n4IdV9b6qqiQvGz5wfNybgcG5rU7aOIa7OivJu1i+1PKxkZc/AJxJsiXJ30vygiQ/B7wKuH/kZz98\njsuVWuU1d3VSkrcCHwI+n+QjP3mZ5U+Xnk7yz1i+C+Yy4GrgWeDE8Gd/Dfjmua9aao/X3NU5SV4K\nHANexHKgj/pCVb0jyWuAfwUcAR4FLgDeNPy5Y1V1xzkrWNoAhrskdZDX3CWpgwx3Seogw12SOshw\nl6QOMtwlqYMMd0nqIMNdkjrIcJekDjLcJamDDHdJ6iDDXZI66P8DNocSzQ3ikt4AAAAASUVORK5C\nYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x7f75050>"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This plot shows the intervals that correspond to the respective domains of $\\eta$ with the vertical dotted lines showing the $\\mathbb{E}(\\xi|\\eta) $ for that piece."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Example\n",
      "-----------\n",
      "\n",
      "This is Example 2.4\n",
      "\n",
      "![alt text](files/ex24.jpg)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x,a=S.symbols('x,a')\n",
      "\n",
      "xi = 2*x**2\n",
      "\n",
      "half = S.Rational(1,2)\n",
      "\n",
      "eta_0=S.Piecewise((2, S.And(S.Ge(x,0), S.Lt(x,half))), \n",
      "                  (0, S.And(S.Ge(x,half), S.Le(x,1))))\n",
      "\n",
      "eta_1=S.Piecewise((0, S.Lt(x,half)),\n",
      "                  (x, S.And(S.Ge(x,half),S.Le(x,1))))\n",
      "\n",
      "\n",
      "\n",
      "v=S.var('b:3') # coefficients for quadratic function of eta\n",
      "h = a*eta_0 + (eta_1**np.arange(len(v))*v).sum()\n",
      "J=S.integrate((xi - h)**2 ,(x,0,1))\n",
      "sol=S.solve([J.diff(i) for i in v+(a,)],v+(a,))\n",
      "hsol = h.subs(sol)\n",
      "f=S.lambdify(x,hsol,'numpy')\n",
      "print S.piecewise_fold(h.subs(sol))\n",
      "t = np.linspace(0,1,51,endpoint=False)\n",
      "\n",
      "fig,ax = subplots()\n",
      "ax.plot(t, 2*t**2,label=r'$\\xi=2 x^2$')\n",
      "ax.plot(t,[f(i) for i in t],'-x',label=r'$\\mathbb{E}(\\xi|\\eta)$')\n",
      "ax.plot(t,map(S.lambdify(x,eta_0+eta_1),t),label=r'$\\eta(x)$')\n",
      "ax.set_ylim(ymax = 2.3)\n",
      "ax.grid()\n",
      "ax.legend(loc=0);\n",
      "#ax.plot(t,map(S.lambdify(x,eta),t))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Piecewise((0, x < 1/2), (2*x**2, And(x <= 1, x >= 1/2)))\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 13,
       "text": [
        "<matplotlib.legend.Legend at 0x840a310>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD7CAYAAAB+B7/XAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcTfX/wPHXJ2Rku2Qp65BskS2SIVMkZsoaEmpo8Uuk\nTcn0pU2RlKQkKlokhDQjuykmgzCoLDFj38JMwkxm+fz+ODNzx6x37tx7z7l33s/H4z7Gvffcc9/e\nnXl3vM/nfD5Ka40QQgjfcY3ZAQghhHAtKexCCOFjpLALIYSPkcIuhBA+Rgq7EEL4GCnsQgjhY4p7\n6ouUUjKuUgghnKC1VgXZ3qNn7FpreWjN+PHjTY/BKg/JheRCcpH3wxnSijHBoUOHzA7BMiQXdpIL\nO8lF4UhhF0IIHyOF3QQhISFmh2AZkgs7yYWd5KJwlLM9nAJ/kVLaU98lhBC+QimFtvLFU2GIiIgw\nOwTLkFzYSS7sJBeFI4VdCCF8jLRihBDCwqQVI4QQQgq7GaR/aCe5sJNc2EkuCkcKuxBC+BjpsQsh\nhIVJj10IIYQUdjNI/9BOcmEnubCTXBSOFHYhhPAx0mMXQggLkx67EEIIKexmkP6hneTCTnJhJ7ko\nHCnsQgjhY6THLoQQFiY9diGEEFLYzSD9QzvJhZ3kwk5yUThS2IUQwsfk2WNXStUEvgSqABr4VGs9\nLYftpgHdgMtAiNZ6Rw7bSI9dCCEKyJkee/F83k8CntVaRyulygDblFKrtdZ7Mn1pEFBPa32zUup2\nYAbQtqDBCyGEcI08WzFa61Na6+i0P18E9gDVsmzWHZibts1mwKaUquqGWH2G9A/tJBd2kgs7yUXh\nONxjV0r5Ay2AzVneqg4czfT8GFCjsIEJIYRwTn6tGADS2jCLgFFpZ+7ZNsnyPOdm+u23Fyg4XxVo\ndgAWEmjWF7/2GnTtata35ygwMNDsECxDclE4+RZ2pVQJ4Hvga6310hw2OQ7UzPS8Rtpr2YRUqoT/\njTcCYCtThub16xPYqhUAEdu2Achzee7+5zNmELFoEfj5ZRSQ9H/6y3N5bvbziIgI5syZA4C/vz/O\nyG9UjMLon5/TWj+byzZBwAitdZBSqi0wVWud7eKpjIqxi4iIkDOSNKbkYsIEuHQJ3nrLs9+bDzku\n7CQXdu4YFRMADAJ2KaXShzCOBWoBaK1naq2XK6WClFIHgEvAkALGLYRnlSkDZ86YHYUQbiNzxYii\nZ/Zs2LQJPvvM7EiEyJfMFSOEI8qUgYs5jQEQwnPC94cTnxh/1WvxifGE7w8v9L6lsJtAxujamZKL\nMmWMHrvFyHFhVxRyEVArgNC1oRnFPT4xntC1oQTUCij0vh0a7iiETyldWs7YhelsfjYmdJrAk0uf\no2ftIfxyfj4TOk3A5mcr9L6lxy6Knq1bYfhw46cQJjp1/iJ1JrQjsdxuYkfF4m/zz7aN9NiFcIT0\n2IUFJF5JpsmbvSnpl8TBkTFMjpycrefuLCnsJigK/UNHmZILi7Zi5Liw8/VcpKZqmoQO4ZLfPv4c\n/Qt1K9ZhQqcJV/XcC0MKuyh6LHrxVBQdd7/+GseIInpkJNVslQF7zz3ySGSh9y89dlH0/PcflC0L\nV66YHYkogh75YDbzjr7Njqd+pUmd/CfClR67EI649lrQWgq7cLusY9Vf/Sacr46HMv7ulxwq6s6S\nwm4CX+8fFoQpuVDKku0YOS7sfCUXmceqz129ldd2P0xg3Q6MuLufW79XxrGLoin9AmqFCmZHInxY\net/84flPEbZnNa1uaMvikNkuGaueF+mxi6KpYUNYutT4KYQbRR88yW2ftCGlzLFcx6rnRXrsQjhK\nxrILDzhy5h/aftSFG8tUJ3ZUrEvHqudFCrsJfKV/6Aqm5cKCY9nluLDzhVzEX0yk6dvBlC6t2Dn6\nJ/xt/i4dq54XKeyiaLLgxVPhO64kpXDLuIFcW0Kx58X1VLzOuJbjyrHqeZEeuyia+vWDPn2gf3+z\nIxE+IHx/OAG1ArD52Yy7Ssc8yYmkP5k15Dn63tqzUPuWHrsQjpIzduFCmYc1dnrjdWKTN9GjfQPu\nqR9oSjxS2E3gC/1DVzEtFxa8eCrHhZ235SK9xdL6vR5s+PdzerRpyQfBk90+rDE3Mo5dFE0WvHgq\nvNu4L1dw8J+96LJnmNh1vGlFHeSM3RSy+rqdabmwYCtGjgs7b8vFm/NX8OHBkXS5+S6PDmvMjRR2\nUTRZsBUjvNOM8Ej+Fz2Qjv53Mn/QJx4d1pgbKewm8Lb+oTvJOHY7OS7svCUXizbs4qlfejOo3iiW\nDvkso/3iqWGNuZHCLoomC7ZihLVlnalxXfRB+v14L/fVCOGrx8Zl66nb/GwE1w/2dJiAjGMXRdWP\nP8LMmRAWZnYkwkvEJ8YTujaUCZ0mEHP0Mm1mtqNuRX+2PLfUrRdKnRnHLqNiRNEkZ+yigNLbKyN+\nfJ6FURu5sdwNbi/qzpJWjAm8pX/oCTKO3U6OCzur5uLypWIsjtrGlXL7+Xnkt5Ys6iCFXRRVFrx4\nKqzt/IUEGr3ZjZKlkjg4MoYpm941dUhjXqTHLoqmI0egfXvjpxD5uJhwhTqh95NYKpaDYzZRpez1\nV/XcrdZjlzN2UTRZsBUjrOlKUgqNQweTXOwf9o2OpErZ6wHzhzTmRQq7CazaPzSDjGO3k+PCzsxc\nZB7WmJySSpOxT3BBn+KTwS9RzVb5qm3NHNKYFynsomi69lrQGq5cMTsSYTHpMzWevxxH61ee40Ty\nbnp1aMC99e8yOzSHSY9dFF02G8TGyoLWIpv4xHgaT7yH8/9eot8d7Zh237umjYCRcexCFET6WHYp\n7CKL/u/O4Oy/50kqF8PrnZdbdlhjbqQVYwLppdqZmguLXUCV48LOzFz0nPg+6+Jm0bNVB0vM1OgM\nKeyi6LLgBVRhrgenfEzY2ffp3qI9n/aeaomZGp0hPXZRdHXsCK+/bvwURd6QaZ/z1dFXeaPzeJ7s\n2Oeq9kt8YjyRRyJNGQEj49iFKAg5Yy+yss7U+NQn85h7LJRxga/w8r2PWmqmRmdIYTeB9FLtpMdu\nJ8eFnbtzkXnx6dGff8+MmGfp2iCQpzv1c+v3eoqMihFFl8zwWGSl3zXaacZDbD+9hS43d2LewBle\nN/olN9JjF0XXyJFw883w9NNmRyJMMO7rMN74/REodZ7YUbH42/zNDilHbumxK6U+V0qdVkrtzuX9\nQKXUP0qpHWmPVwoSgBCmkTP2Iuv1b3/ijd0h3OVvjcWnXc2RHvsXQNd8tvlZa90i7fFmQQJQSsmj\nEA9vJz12O+mx27kzFxO+W8n46MEE+ndkcchsrx3SmJd8C7vWegMQl89mhaowWmt5OPEQhSSjYnxe\n1tEv7yxawyvbBxJcrwdLLLT4tKu5YlSMBtoppXYqpZYrpRq7YJ+iiAgMDDTvyy3WijE1Fxbjqlxk\nHv0yZfE6Xtr6IHfWvpOvB0/x+iGNeXHFqJjtQE2t9WWlVDdgKVA/pw1DQkLw9/cHwGaz0bx5cxd8\nvQD7P13TfyHkuQPPDx8mMO2M3RLxyHO3PJ/QaQIdXrqP3//eRfuATvww9HOio6ItE1/W5xEREcyZ\nMwcgo14WlEOjYpRS/sCPWuumDmwbC7TSWp/P8nqOo2LSrvg6Gq/IxBdyFxERYd6Z6o8/wqefGj8t\nwNRcWIwrc/HekvU8v6UP+MVZevRLbky581QpVVWlXcVTSrXB+J/F+Xw+JoT5LHbxVLje5O/X8nxU\nXwKq3+mTo19yk+8Zu1LqW6AjUAk4DYwHSgBorWcqpZ4CngSSgcvAc1rrqBz2I2fsDvrrr7/4/fff\n2bVrF/fffz8tW7bMcTvJXSFt2QJPPQVbt5odiXCDiQtX8/JvA2hfqwM/PvoFNj+bx9YpdSVnzthN\nv0FJilN277//PgEBATRq1Ihhw4Yxb968HLeT3BXSn3/CAw8YP4VPmfDdSv63fTD/d+sLvNXnCctM\n6OUMmQTMoubNm8eUKVPo378/8+fPz3f7Z599ljZt2nD06FHq1KnjgQjN487xyvmyWCvG1FxYTEFy\nkXVI4+vf/sQr2wcyvNmLfDzwRZ8e/ZIbKewO2LNnD/369ePbb78lNja2QJ89cOAA586d4/nnn+ej\njz7iySefdHgfS5YsITQ01JmQhSNkHLtPyDykcdzXYRk3H73Z+zGzQzONFPZ8REVFcc8999CoUSP+\n/fdfSpYsWaDP//HHH7zzzjsAVKpUiXr16rFt27Z8P7ds2TKefvppjh8/7lTc3sLUUSAyjt2yCpKL\nzBN6vfF7CHfVCbzq5qOiSHrs+ejXrx/du3dn0KBBV70eExPDrFmzcv1c27Zt6dGjB0lJSezbt48m\nTZqgtaZmzZqEhYXlOYZ/yZIlvPXWW9hsNgIDA3M9a7d67ixPayhRAhISjJ/Caz332UKm/jUcXeqs\nVw5pzItcPHWDF154gT///JOhQ4dSu3ZtWrdu7fS+wsLCmD17NkuXLgWMs/JixYqxYcMGmjZtyooV\nKwgNDaVhw4YO7c/quXOE6WO3bTY4dMj4aTLTc2EhBcnF/338FZ8eGk2X+oF88sBEJkdO9qpRL/mR\ni6du0LhxYxo0aEBqairFihVzej/x8fHMmTOHr7/+GoAjR47QuHFjgoODWb16NcHBwfTv359atWq5\nKnThCItdQBUF8/DUWXwa+xLdGt3N/EGf+OSEXs6w/Bm7qyYwdOavuWDBAmrUqEG7du2yvedoK8b4\nbs3LL7/MmDFjsNlsHD58mNq1awNw+vRp+vfv79SICF84Yzddgwbwww/g4L+ShHnC94cTUCsg40y8\n7+TpLD4ziYEtezKt1xtePaQxL9KKcbFBgwbx2WefFfiCaVbTpk0jICCA6tWrc+TIERISEqhatSqJ\niYns2LGD2NhYXn/9dcLCwrjvvvsc3q+Vc+c1WrUyphVo1crsSEQ+Mt9cNOj92aw4N53uzTrwed8P\nfabtkhNpxbjY4MGD6dChA2+//TaHDx92ah8bN27k2WefpXXr1lSrVo077riDevXqsWrVKsLDw9Fa\nk5iYyJIlS6hSpYqL/wbWZ/rYbQu1YkzPhYXklAubn4037nqTJpPuZdW5GfRqGejzRd1ZsuZpHu69\n917at2/P4sWLeeSRRwgJCSEkJKRA+2jfvj0pKSnZXn9almOzBhnL7jVSUzXd3prI3wnxJJWLYXLQ\nWinquZAz9nyULl0am81Gv379MvriwnVMHwViobHspufCQrLmIjkllaYvD+f3y6vp27Z9kZrQyxnS\nY/dikjsXGDIEOnSAoUPNjkTk4nJiEo1DQ4hLPUyvDo2YGjTZayf0cob02IXXMb2vbKEzdtNzYQHp\n876k5yL+YiJ1xvYg7pr9zHzkuYyiDr63nJ0rSY9dFG0Wungq7PO+dCvRjVPnL9LwjSCSS51g7+hf\nqVEh++CCojChlzOkFePFJHcu8OabxpQCEyaYHYlIE58Yz8gfX2Bx1DZK+CWyf8xGqpS93uywTONM\nK0bO2EXRVqYMnD1rdhQik0PHE1gQtZEr5fbx14iDRbqoO0t67MJUpveVLTTc0fRcWMC66IPcNuMO\nyp1THBwZw/tRU2TkixOksIuizUIXT4u6Bb/spPO89tStWJsvHppM3Yp1ZN4XJ0lhF6Yyfey2hS6e\nmp4LD8m64hHAu0t/ov9PHele6xG2PPcD93UxptaQkS/OkcIuijYLtWKKiswrHgGM+fI7Rv/Wm3HN\nvmDpiIlFcik7V5PCLkxlel/ZQq0Y03PhIeln4aFrQxn44XtM2j+EaXf8yGsP9srYpqjkwl1kVIwo\n2izUiilKbH429u60sS7peT7ttIrH7+psdkg+Rc7YXaCgC1xndfLkSS5fvuyiaLyL6X1lC7ViTM+F\nhySnpNJ87FNEXJzB4uBIoi8vzdZzLyq5cBcp7IUUExNDVFRUofZRuXLljAWvhYdZqBVTFFxMuELd\nF/vxp17Cb49vo9dt7WTkixtIYS+kmTNnMmDAgDy3mTRpUp7vFy9enODgYL788ktXhuYVTO+lWuiM\n3fRcuFjW0S/H/r5AzbFdiCvxO/ue3UaLOnWAnEe++FouPM2ShT2n4VDxifGE7w/36D7SDRw4kA4d\nOjB48OCMx4gRI9i1axc1atTItv2ePXvo168f3377LbGxsSQmJub7Ha1bt2bNmjUFjk0UUsmSkJIC\nSUlmR+JzMo9+2RVzipvfbo8qfZaDoZHUqXLjVdvKyBfXsmRhzzocKn16zoBaAR7dB8D58+dp164d\ngwYN4quvvsp4TJ8+nR9//JG77777qu2joqK45557aNSoERcvXizQsnqVK1fmwIEDBYrP25neS1XK\nMu0Y03PhYuln4g/Pf4qWM9pwfflSHBi7waEpAnwtFx6ntfbIw/iq7HJ7PS4hTg8PG65j42L18LDh\nOi4hLsft8uKKfSxYsEDv27dPz507N9t7PXr00KmpqVe91rdvX/3VV19d9dqrr77q0HfNnTtXz58/\n3+HYcsudKKBq1bQ+etTsKHzSzOW/avViJc2r6Ni4WLPD8Uppv+cFqreWHe5o87MxOmA0dT4w+nAf\n//ax0/v6+LePiR0V69Rk/DExMdSvX5+vv/6a1atXA9CzZ0/69OnD5cuXUerqSddq1arFvHnz8PPz\no3bt2rRu3Trjve3btxMVFcXx48dp27YtSUlJLFq0iHnz5gFQoUIF9u/f7/Tf0xtFRESYf3ZmkTN2\nS+TChV78YjGT9z5B6xvbsiBkOpMjJzu8KIav5cLTLFvY4xPjmRw5OWMJLGdWSUlvv4wOGO30PooX\nL87mzZtZtWpVtvdyWsu0cePGpKSkkJqaSrFixa5678yZMzRp0oSVK1cyYcIEUlNTeeaZZzLeL1Wq\nFFeuXClQfMIFZCy7y/V+5wOWnp1I53p3s/DhT6+6KcnXVzyyhIKe4jv7oACtmPQWSnrrJOtzR7hi\nHzExMXrZsmX6ww8/zPH9Ll26XPX8u+++05GRkdm2y9yKGTt2rF6wYIHWWutNmzbp3r17Z7y3aNEi\n/cknnzgcX245FQXUoYPWERFmR+GVwvaFXfU79d+VZN10zP/p4qNr6ndWf57t9y0uIU6H7QvzdJhe\nDSdaMZa8eBp5JPKq/6s7MxGQK/axYcMG2rRpQ6lSpXJ8/4YbbuBipjO9ZcuW0apVqzz3uXbtWjp1\n6gTAihUruO+++zLeO3nyJPXq1XM4PuEiFmnFeKPMgxTOX0ig9ks92ccyfhv+C6M7D5F5X0xiycIe\nXD+40AeEK/YRGRnJkCFDWLhw4VVDHdevXw9Ax44d2bJlS8b2gwcPpkOHDrz99tscPnw42/4SEhIo\nW7YsFStWBGD16tUEBQVlvB8dHU1AQMFG7Xg7S4xXtshYdkvkooDST5iG//AcNV8N4GKpffz1wnaa\n+fsXar/emAsrsWyP3QpmzpyZ5/u9e/fm3XffzRjyeO+999K+fXsWL17MI488QkhICCEhIRnblypV\nKuMC7Pnz50lOTqZq1aoAJCYmUq5cOfz8/NzzlxG5kzP2Qvl11ykW/LaGlPJHOTgyhloVq5odUpFn\nyTN2b2Gz2ahUqRJnMy2tVrp0aWw2G/369aN27dq5fjYqKopeveyz2c2fP59hw4a5NV4rssTIB4tc\nPLVELgpo8vdrCf6+A/Ur3UzsqFimbHrXJVMDeGMurETO2Atp1KhRzJ49m8cffzzjtfvvvz/fzwUF\nBWW0YY4ePUqFChVo0KCB2+IUebBIK8bbhEz7jLnHxxBQM4Cwx+bIyBcLkTP2QlJKXVXUc5Lf2UfN\nmjXp0aOHC6PyHpbopVqkFWOJXOQi8xQdySmp3P7KS3xzeAIDb30wo6iD61Y8snIuvIEUdg/o2LGj\n2SGIvFikFWNl6aNfDp4+Qe3RfdlzaQP9br+L6X3ekJEvFqSMYZIe+CKldE7fpZTCUzH4Gsmdi8ya\nBZs3w+zZZkdiaT//sYfO33TkxqT23Nu+CpPvzb6MnXC9tN9zlf+WdvmesSulPldKnVZK7c5jm2lK\nqb+UUjuVUi0KEoAQprNIK8bK5q7eSqc5XbijVAhHyywhtOMYKeoW5kgr5guga25vKqWCgHpa65uB\nJ4AZLopNFAGW6KVa5OKpJXKRg1GzvmPImmBGNppI05aXMqb5cOfCGFbNhbfIt7BrrTcAcXls0h2Y\nm7btZsCmlJKBrMJ7yBk7kH0Ng+SUVNqNf5HpB0byWZfvuVL1VyZ0moC/zV9WPbI4V1w8rQ4czfT8\nGJB99QkhcmCJ8coWuXhqdi4yTw9wJu4StUb35LcrX7Fh6Eaq1LxQ6Ck6CsLsXHg7V41jz9rYlyt6\nwntYpBVjtvRi/diiUYTt+I2SpVI48Fw0ta6vCtTPcXsZ/WJNrijsx4GamZ7XSHstm5CQEPzT5pCw\n2Ww0b97cBV8vwN6TTD/T8Zbn6a+ZGk+ZMkScOweZ5gA3I57o6OiMaZzNysfuC9ewePdP6HN/83nv\neWlF3fPxTJ06lebNm5t+fJrxPCIigjlz5gBk1MsCc2QKSMAf2J3Le0HA8rQ/twWictkurykphRN8\nIXfr1683OwStz57VumJFs6MwNRcpKam6/7sfaV6spG+fel+hVh1zBUscF+4WG6v1xx9rHRys9fjx\nuW6GE9P25juOXSn1LdARqAScBsYDJdIq9cy0baZjjJy5BAzRWm/PYT86p++SsdjOk9y5SGIi2GzG\nzyIgfH84AbUCMvrlFy79R6tXHydWrSf4lo7M7T8dm58tY6EamR7ARZKSYONGWL7cePz9N3TrBkFB\n0KULVKiQ48ecGccuNygVQmxsLHXq1Mlzm5MnT1K+fHmuu+46l3+/N+fOUrSGEiUgIcH46eMyF+xD\nxxNoP70HutRZpnQfx4PNe15VxOMT44k8Eim9dGedPAk//WQU8jVroH59o5AHBcFtt8E1+Y9fccsN\nSiJnMTExREVF5btd5cqVeeeddzwQkXeyxHhlpYwLqCYPefRULtIvkvae8xitZrbEZlMcDv2N/2sb\nYpnpASxxXDgjJQU2bYL//Q9atoTGjWHlSrj/fti3D7ZsgVdfhTZtHCrqzpLC7qSZM2cyYMCAfLcr\nXrw4wcHBfPnllx6ISjitiI1lf/qThaw/vI7U0qfYOOo7KpWuaHZI3uvcOfjmGxg4EKpWhSeeMNou\nH3xgtFu++w4eecR4z0OksDth586d1Kjh+FD91q1bs2bNGjdG5L0sM17ZAmPZPZGL+IuJNBj9GN8d\nmUJw/W4euYvUGZY5LnKSmgrbt8Obb0K7dlC3LixcCB07wo4dsHs3TJwIHTpAcXNmRpfC7oSwsLCM\nVZMcVblyZQ4cOOCmiEShFYGx7JF/HKb6/9pzMeUsD7W/k68HfCR3kTrqn39g0SIYOhSqV4cHHzTO\n1F97Dc6cgaVLjTP1mjXz35cHSGF3wtatW2ncuHGBPtOsWTO2bdvmpoi8l2V6qRZoxbgyF1mnB5i4\ncDUdvmzNLZWa8cmwx3g/6B2P3UXqDNOPC63hzz9h8mS46y6oUQM++wxatIANG2D/fnj/fbjnHihZ\n0txYc2D9FZRUgS4G586J0SM7duxg06ZNHD9+nLZt25KUlMTChQu5fPkyKktcy5Yto1ixYmzYsIGm\nTZuyYsUKQkNDadiwIQAVKlRg//79LvmrCDfwsTP29OkB3rjrTQa89wmrL3xAu1rtCHtsSo5DF+Uu\nUoz/sa9fbx+OqDUEB8PzzxvFvXRpsyN0mPULu4nD+U6fPk2TJk1YuXIlEyZMIDU1lVGjRmVbwu7I\nkSM0btyYevXqMW7cOMaMGUP58uWpVatWxjalSpXiypUrnv4rWJ5leqkWOGN3ZS5sfjaeuGU0Nd9s\nCf/Z6NPyXmb1+cBrxqN77Lg4eBDCw41CHhlpDEEMCjKeN2rkuhNLD7N+YTdR165dCQ0NZdCgQQBs\n2bKFNm3acClLAUgv4KdPn6Zs2bLYbDbuu+++q7b5559/qFhRRh5YlgUunrrS12u3MeSnvtx8bSB7\nys1hcvASrynqbvXff0YrJb2YX7hg3CT0+OPG6JXy5c2O0CWkx56PtWvX0qlTJwBWrlxJcHAwVatW\n5WKmIrB371527tzJ8uXLufPOOwHjAmtmJ0+epF69ep4L3EuY3ktNZ4FWjCtykZqqGTBlBg+v7srj\nDcZxV8B1lh35kheXHhfHjhmrZPXsCVWqwLhxcP31MH8+HD8On38Offr4TFEHKex5SkhIoGzZshln\n2qtWrSI4OJiOHTuyZcuWjO1WrVpFWFgYWmsSExNZsmQJVapUuWpf0dHRBAQEeDR+UQAWaMUUVNYL\npKfOX8R/dD++P/UOC3v8hKq2tWjOn56cbJyVv/wyNGsGzZsbvfO+fY3Wy6+/wiuvGBdC3XiTkKkK\nOrmMsw+8fBKwc+fO6TZt2mittY6Li9OhoaEOfzYhIUE/++yzLo/JW3LnFd54Q+sC/De1griEuIyJ\nun7Y9Icu8Wx9XX5sQx1z+oQO2xeWbQKvuIQ4HbYvzKRo3ez0aa3nztW6f3+tK1TQukUL47/nr79q\nnZxsdnSFghOTgEmP3UFRUVH06tULMKYcrlSpEmfPnqVSpUr5fnb+/PkMGzbM3SGKwihdGs6eNTuK\nArH52XjjrjdpN7UPe/7ZQaMqLfj1me+x+dmoUyX7CBefGvmSmgrbthl98vBwY/hhp07Ghc/33oNq\n1cyO0FQ++u8Q1wsKCmLMmDEZz0eNGsWSJUvy/dzRo0epUKFCtpE0wmCZHrsFWjEFzcWJc/9y2/+e\n5uDfR8AvjuX/95nPXCDNMRdxcfbb82+80fh58SJMmmTcJPT99/Doo0W+qIMUdqcppXj88cfz3a5m\nzZr06NHDAxGJQrHAxdPcZO2lA8xe8wu13mlIsWsUA9vf7ZUXSPOlNezaZdyef+edULs2fPWVMYFW\nVNTVNxBde63Z0VqKTNvrxSR3LrRsGcyebfy0mMzT7Ja7tjy93p3Msn/H8/ANkyjjvy9jvnSfmD/9\n4kVYu9Z+k1CJEsZNQkFBEBgIpUqZHaHHyXzsRYzkzoXWrTMmdVq3zuxIchSfGM/IH19gZdRRzpfc\nzvc9V1JB2VgqAAAWWklEQVTcdvKqBTPSt/Oq+dO1Nvrj6YU8Kgpuv90o5t26QYMGXnuTkKvIfOzC\n61imx26BVkxeufhq1S7mbf2Jv8utYseIjfRo05Lg+sGWmT+9QBISYMUKePppuPlm46Lnnj3w1FNw\n4gSsWUNEixbQsGGRL+rOklExQoAlLp7m5HJiEl3eeo1fE2dxW7XWLAiJZHLkZCZU9LJ2y6FD9rPy\nX34xxpcHBRkXPG+9VQq4i0krxotJ7lzo8GHjAt3hw6Z8fdZ1SAGWbd1BvwUPUjqlJve0rsknvd73\nnl66k+t7iuy8tscunCeF3UXOnjX+6W/SWPasF0iHfPQJX555juAy43iiRxPu9O9g/V56but7dusG\nrVv77l2ebuaVhb0oioiIsM6shiazTC4SE8FmM36aJGxVGAsuLGVd1BlOXruR2Z1+YEinDqbFk6+U\nFNi82X5WfuiQMT95cDDce2+hloKzzHFhAc4UdumxCwHGYglJScY8IyYtZxa2aR/fXFpOatmT7Bi6\nh+Y1G5oSR57OnjUWZ16+3PhZvbpxVv7BB3DHHablTlxNztiFSFe+PBw54vFZ/s5fSKDT2y+zM3kB\nbWq0Yv4jHxoXSK3QQ09Nheho+637f/5p3BAUFGQ8CrD2r3COtGKEKIzq1WHLFuOnm2S9SPrNuu0M\nCRtA2eIV6XpbYz7qPsX8C6QXLsDq1fYWS/ny9kLeoYMll4LzZTKO3UtYZuy2BVgqFx4Yy56+ZN2p\nC2fp8sZbDFrZhQY31OKjwc8wwNbLnHVI09f3fPdduPtu4yx81ixjutsNG2DvXmNirc6dPVbULXVc\neCFpiAmRzgNj2W1+Nu6uOIRabzeh1H830adlMLMfMJasizgXkW1bt416uXzZmKM8fSWh9PU9n3vO\nKO7XXeee7xUeIa0YIdLdeacxrUDaKliudiUphX5TprEs/i26VhjBT4mvEjsqFn+bv1u+L5uDB+3t\nlY0boVUr+zwsjRvLTUIWJa0YIQqjdGmXnLHnNBvjkqht2F5uyvpTS1jaZxV1mpxx/4yM//1njCd/\n7jljzpWAANixw5ja9tgxiIiA0aPhllukqPsYKewmkP6hnaVy4aIFrdP76PGJ8SSnpHL/pIn0Dg+g\nU+WHOThhKSvPzM5xyTqX5OLYMfj0U/v6nv/7n3GX57ffGvOwfP45PPCA5df3tNRx4YWkxy5EOhdd\nPE2/8PnIdyNYv/MvEkseYXGPSHq1bUX4/vCrRrpkvkhamtIF/7LkZNi0yd5iOX7cuDmob19jGmIH\nVvgSvkd67EKkGzHCmFZgxIhC7SbxSjIPTJlK+IUJ4BfPvuEHqF/5JhcFibFa0IoVRiFftQr8/Y3b\n9oODjSlvixVz3XcJ00mPXYjCKGArJqde+tx1G7GFNmLjmWV0bxhM7KhYPtjyXuH66KmpsHUrvPaa\nsXpQ/frGgiBdusDvv8P27TBhArRrJ0VdAFLYTSH9QztL5aKArZjMvfQLl/7jjvGjCVl3D/dVe4yH\nOjdhbv/p2froebkqF1nX9wwJMWKbONE4Y1+0CIYO9dn1PS11XHgh6bELka5MGTh/3uHN0/vjPb94\nlF8P7KREyVQiBm7nYokYAmoNy7GPnuu4dK2N4YhRUcbY8p07jWGXQUHGmbq/vwv+gqKokB67EOk+\n/dRoecya5dDmh0/HE/zeWP5kIfq6sxwcGUPdinUc/76LF43hiOkXPkuWtC8JV0TX9xTZSY9diMLI\n4c7TnPro5y/H0XP6GOpOuYUkncDA23oQOyqWKZvezbvdojXs2wfvv29Mb3vjjfDRR9CokbGA84ED\nMG2aUdilqItCkMJuAukf2lkqFzlcPM3cRwcI/20n1d9szk8xP/BOwGd07ngdH97/bu699PT1PUeO\nhHr17Ot7Dh9ujCtfvRqefRYaNCDi5589+be1NEsdF15IeuxCpMvh4ml6f/ylVS+zb1d5fv5vKneX\neokf/hfKz0dX82it7GPSt/+6mLv3JtrX92ze3OiVL1kCTZvKXZ7C7aTHLkS6zZvh6aeNn5l8uOwX\nXvj5ca6U288390TwULuOV3/uyhWIjJT1PYVbyApKQjgpfH84HUpUpFymM/aNe/fxwBf/x9+pf9H0\nxiYsHrqCKZveJSixGba4BGN9z/Bwoz+evr7nF1/AbbfJ+p7CVPkefUqprkqpvUqpv5RSL+XwfqBS\n6h+l1I60xyvuCdV3SP/Qziq5CKgVwOToj0m9+C+JV5K5f9JEOnzTnBtKNOThO4KIePIb6u47xZSN\n13Ghyc2k3tLYuOuzRw/jguiWLfDqq8YNRE4WdavkwgokF4WT5xm7UqoYMB3oDBwHtiqllmmt92TZ\n9GetdXc3xSiE29n8bLzQ5VUuPN2Yqi/fCtde4Pv2P1Bn72puWXCIa59oANWq4RccTIVP5rKiWjJB\njeWQF9aUZ49dKXUHMF5r3TXt+RgArfXETNsEAs9rre/P84ukxy4s7NgxGDs6ntkLKvB6R3jlYkv8\n9h2Q9T2F6dwxjr06cDTT82Npr2WmgXZKqZ1KqeVKqcYFCUAIMyUkGNOsNGsGxer8xc/1ijG64VBm\nda9O/JH9sHQpPPGEFHXhVfK7eOrIKfZ2oKbW+rJSqhuwFKif04YhISH4p90abbPZaN68OYGBgYC9\np1YUnmfuH1ohHjOfp7/m6e9fvz6CX36BOXMCadkSJr0fxmf7J/DmmAB+HvIZdVaFEfLxMOY8O8dY\nts4D8UVHR/PMM8945O9v9edTp04t0vVhzpw5ABn1sqDya8W0BV7N1Ip5GUjVWk/K4zOxQCut9fks\nr0srJk1ERETGf9Cizoxc7NoFzzxjjEqcOtW4Zyh8fzh7z+5l37l9fHr/pwDEJ8bnPb+Li8lxYSe5\nsHOmFZNfYS8O7AM6ASeALcCAzBdPlVJVgTNaa62UagMs0Fr757AvKezCVCdOGAsKhYXBuHEwbBgU\nz/Rv1hdXv0gFvwq83OFl84IUIguX99i11snACGAl8CfwndZ6j1JqmFJqWNpmDwC7lVLRwFTgwYKH\nLoT7XLpkjERs2hSuv94YnfjUU1cXdYDY+FjqVCjAJF5CWFS+A2611j9prRtoretprd9Oe22m1npm\n2p8/0lo30Vo311q301pHuTtob5e5v1zUuTMXKSnGEp/16xvFfNs2eOcdsNly3j4mLoa6Feq6LZ78\nyHFhJ7koHLnzVPgcrY15t8aMgbJlYfFiY8W4/MTGxVLHJmfswvvJXDHCp0RFGQX91Cl46y3o1cux\nObfiE+Op8V4N/n35X5RM0iUsROZjF0XWnj1GEe/bFwYPNpYC7d3b8YkUY+NiqVuhrhR14ROksJtA\n+od2hc3F0aPw6KPQsaOxlvP+/cbzrBdG82OFC6dyXNhJLgpHCrvwSidPGjPsNmsGVasaBX30aOcX\nHpL+uvAl0mMXXuXvv2HSJGO0yyOPGP30qlULv9+nwp+iQaUGPH3704XfmRAuJD124bPOn4exY6Fh\nQ2N+l927jaVDXVHUIa0VI2fswkdIYTeB9A/t8svFuXPG3aL16xtn69u3G+s/V886FV0hSY/dWiQX\nhSOFXVjSmTNGm6V+fWPo4ubNMGsW1K7t+u9K1akcij8kZ+zCZ0iPXVjKyZMweTLMmQMDBsBLL0Gt\nWu79zhP/nqDFzBacfuG0e79ICCdIj114rdhYGDECbrkFUlONHvpHH7m/qIOMiBG+Rwq7CaR/aDd7\ndgQDBxrrP5cpY9xoNHWq63voeYmNjzV1jph0clzYSS4KR+aKER6nNfzyC0ycaKwB/dJL8PHHUL68\nOfHExMXIGbvwKdJjFx6TkmKsNDd5MsTFGTcUDR4MJUuaG9eQH4YQUDOAx1o+Zm4gQuTAmR67nLEL\nt7twwbih6IMPoFo1o6D37AnFipkdmSE2LpbBtw42OwwhXEZ67CYoKv3DQ4fg+eehTh1j1sX58yEy\nEvr0sRd1K+TCKq0YK+TCKiQXhSOFXbhUev+8b19o1QquuQZ27DCKuiNzonvalZQrnL50mprla5od\nihAuIz124RKXLsE338D06XDlirH0XEiIsdCFlR04f4AuX3UhZlSM2aEIkSPpsQuP++svY0TLl19C\nhw7w3nvQqZPj86CbLSYuxvSpBIRwNWnFmMDb+4dJSfD993DvvRAQAH5+xhwuS5dC584FK+pm58JK\nNyeZnQsrkVwUjpyxC4cdOACzZxu3+zdoAE88AT/8YBR2b2WVm5OEcCXpsYs8/fefUbw//RR27jTm\nQH/sMWP6XF/Qb2E/ejXsxYCmA8wORYgcSY9duITWsG2bcWY+fz7ceqtxdt6rl/k3E7maFabrFcLV\npMduAqv2D9NnVmzaFPr3Nxax+O03WLcOHnzQPUXd7FykL2JtBWbnwkokF4UjZ+xF3MWLRqvlm29g\n0ybo3RtmzID27b1nZIuzLvx3gcTkRCpfV9nsUIRwKemxF0FXrsDKlTBvHixfbhTxgQOhRw8oXdrs\n6Dxn56mdDFw8kN+H/252KELkSnrsIlfJybBhA3z3HSxaBI0awUMPwbRpULmInrDKiBjhq6THbgJP\n9Q+Tk2HNGhg2zD75lr+/cWF0wwZ48knzi7qZvVQrjWEH6StnJrkoHDlj9zH//Qfr1xs3EC1dakzA\n1bevMQlXXTk5vUpMXAw3VbzJ7DCEcDnpsfuA8+eNXvkPP8Dq1cbycr17wwMPuGfxZ19x37z7eKLV\nE3Rv0N3sUITIlfTYi5C//oLwcKOYb9sGd98N3bsb64RWqWJ2dN4hNt5arRghXEV67CZwpn94+bJR\nyEeOhHr1IDAQfv8dnnsOTp0y2i5Dh3pfUTerl6q1NnrsFro5SfrKdpKLwpEzdovS2ijca9bAihXw\n66/G/ObdusHixcZNRL4+ztydTl86TelrS1Pm2jJmhyKEy0mP3UJiYmDtWuOxfr0xl3mnTtC1q/Gz\nXDmzI/Qdm45u4pmVz7D5sc1mhyJEnqTH7kW0hoMHjWGHv/wCERGQkGAU8C5dYNIkufDpTlZZDk8I\nd5Aeu4ckJ0N0NHz4IQQGRlCtmtEnX7UKWreGsDBjrpZvvjF65UWlqJvVS7XihVPpK9tJLgpHztjd\n5MQJY+z45s3GY9s2qFHDuH3/jjvgiy+Mm4WkT26O2LhY2tZoa3YYQriF9Nhd4ORJYwWhHTuMn1u3\nGm2V22+Htm2Nn61bQ4UKZkcq0t019y5CO4TSuW5ns0MRIk/SY3ezK1dg/35jtMru3fZCnpwMLVpA\ny5bQrx+88w7cdJOcjVuZ1aYTEMKV8i3sSqmuwFSgGDBbaz0ph22mAd2Ay0CI1nqHqwP1pIQEYxm4\nffvgzz/hjz+MYn7woNE+adLEeAwbZhT0mjULvs5nYGCgu8L3KmbkIikliZMXT1KrfC2Pfm9+5Liw\nk1wUTp6FXSlVDJgOdAaOA1uVUsu01nsybRME1NNa36yUuh2YAVi+eXnhAhw6ZDxiYowz8fTHmTPG\nvCo33wyNGxt3dI4da6zz6Yr1PaOjo+WgTWNGLo78c4Qby9xIiWIlPPq9+ZHjwk5yUTj5nbG3AQ5o\nrQ8BKKXmAz2APZm26Q7MBdBab1ZK2ZRSVbXWp3PaYXxiPFM3TeWZO57B5me76vXII5EABNQKKNB7\nWfd36RLsPRTP6n2RnDsH5eIDOH/CxuHDaYX8RDyXbp1Kvb+f4abqNvz9jWlsOwXHE182kspV4E7/\n7N8z0Ym4c/q7njp7ivD94S75u1rpM87sb8vBLYTvD/do3JmXw0t/L7h+MGaLj483OwTLkFwUktY6\n1wfwADAr0/NBwIdZtvkRaJfp+RqgVQ770nEJcXp42HB9KO6QHh42XMclxGmtdcbrcQlxV/05IUHr\nfYfj9EPfDNfrfo3TS1fE6a4fDtfTZ8fpyZO1fvrFON3w+eG6c59D+oahw3XdxnG6XDmtS5aP0+UG\nDNd33BWnew2I083GDtdvTI7TCxdqve7XOD1k0XAde96xGDK/52jc+X3mtgG35fkZT8Tgjs84s7+X\nQl/yeNxTN03VQ5cOzbYPs40fP97sECxDcmFnlOnc63ROjzxHxSil+gBdtdaPpz0fBNyutR6ZaZsf\ngYla68i052uAF7XW27PsS/uNqc11lxuiUkuQQhIJpfdy7aWbSCx1kBL/NCQ1uQQpKZBCEqkV9qLi\nbkJdf5CS/zakRLESlCgO15RI4vJ1e6nITfxb4iDVr21IKb8SXFM8iWOJe6lX4SYOXzxIw0oNM/6p\nnZSSxN6ze7mp4k0cPG9/L7fX3f2ZpO+TaDm8ZZ6fsWLcjnymoPv7/dPfafFkC4/GvfnYZh5t8SgJ\nyQlM6DThqrN7M4WEhDBnzhyzw7AEyYWdM6Ni8ivsbYFXtdZd056/DKTqTBdQlVKfABFa6/lpz/cC\nHXWWVoxSyjfHOgohhJsVtLDn12P/DbhZKeUPnAD6AwOybLMMGAHMT/sfQXzWou5MYEIIIZyTZ2HX\nWicrpUYAKzGGO36mtd6jlBqW9v5MrfVypVSQUuoAcAkY4vaohRBC5Mpjd54KIYTwDJdPAqaU6qqU\n2quU+ksp9VIu20xLe3+nUqqFq2OwivxyoZQamJaDXUqpSKXUrWbE6W6OHBNp27VWSiUrpXp7Mj5P\ncvD3I1AptUMp9btSKsLDIXqMA78flZRSK5RS0Wm5CDEhTI9QSn2ulDqtlNqdxzaO182CDqPJ64HR\nrjkA+AMlgGigUZZtgoDlaX++HYhyZQxWeTiYizuA8ml/7uqLuXAkD5m2WweEAX3MjtvEY8IG/AHU\nSHteyey4TczFq8Db6XkAzgHFzY7dTfnoALQAdufyfoHqpqvP2DNuaNJaJwHpNzRldtUNTYBNKVXV\nxXFYQb650Fpv0lr/k/Z0M1DDwzF6giPHBMBIYBHwtyeD8zBHcvEQ8L3W+hiA1vqsh2P0FEdycRJI\nX16mHHBOa53swRg9Rmu9AYjLY5MC1U1XF/bqwNFMz4+lvZbfNr5Y0BzJRWaPAsvdGpE58s2DUqo6\nxi/1jLSXfPXCjyPHxM1ARaXUeqXUb0qpwR6LzrMcycUs4Bal1AlgJzDKQ7FZUYHqpqtnd3T0FzLr\n0Edf/EV2+O+klLoLGAoEuC8c0ziSh6nAGK21Vkopsh8fvsKRXJQAWgKdgOuATUqpKK31X26NzPMc\nycVYIFprHaiUuglYrZRqprX+182xWZXDddPVhf04UDPT85oY/2fJa5saaa/5GkdyQdoF01kYd/jm\n9U8xb+VIHlph3AcBRi+1m1IqSWu9zDMheowjuTgKnNVaJwAJSqlfgGaArxV2R3LRDpgAoLU+qJSK\nBRpg3F9T1BSobrq6FZNxQ5NS6lqMG5qy/nIuAx6GjDtbc7yhyQfkmwulVC1gMTBIa33AhBg9Id88\naK3raq3raK3rYPTZn/TBog6O/X78ALRXShVTSl2HcaHsTw/H6QmO5GIvxsyypPWTGwAxHo3SOgpU\nN116xq7lhqYMjuQCGAdUAGakna0maa3bmBWzOziYhyLBwd+PvUqpFcAuIBVjEj6fK+wOHhdvAV8o\npXZinIS+qLU+b1rQbqSU+hboCFRSSh0FxmO05Zyqm3KDkhBC+BiX36AkhBDCXFLYhRDCx0hhF0II\nHyOFXQghfIwUdiGE8DFS2IUQwsdIYRdCCB8jhV0IIXzM/wMwVLXJbGPMdAAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x830c350>"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The figure shows the $\\mathbb{E}(\\xi|\\eta)$ against $\\xi$ and $\\eta$. Note that $\\xi= \\mathbb{E}(\\xi|\\eta)= 2 x^2$ when $x\\in[0,\\frac{1}{2}]$ . Assembling the solution gives,\n",
      "\n",
      "$$\\mathbb{E}(\\xi|\\eta) =\\begin{cases} \\frac{1}{6} & \\text{for}\\: 0 \\le x < \\frac{1}{2} \\\\ 2 x^2 & \\text{for}\\: \\frac{1}{2} < x \\le 1 \\end{cases}$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This example warrants more a more detailed explanation since $\\eta$ is more complicated. The first question is why did we choose $h(\\eta)$ as a quadratic function? Since $\\xi$ is a squared function of $x$ and since $x$ is part of $\\eta$, we chose a quadratic function so that $h(\\eta)$ would contain a $x^2$ in the domain where $\\eta=x$. The motivation is that we are asking for a function $h(x)$ that most closely approximates $2x^2$. Well, obviously, the exact function is $h(x)=2 x^2$! Thus, we want $h(x)=2 x^2$ over the domain where $\\eta=x$, which is $x\\in[\\frac{1}{2},1]$ and that is exactly what we have.\n",
      "\n",
      "We could have used our inner product by considering two separate functions,\n",
      "\n",
      "$\\eta_1 (x) = 2$ \n",
      "\n",
      "where $x\\in [0,\\frac{1}{2}]$ and\n",
      "\n",
      "$$\\eta_2 (x) = x$$ \n",
      "\n",
      "where $x\\in [\\frac{1}{2},1]$. Thus, at the point of projection, we have\n",
      "\n",
      "$$ \\mathbb{E}((2 x^2 - 2 c) \\cdot 2) = 0$$\n",
      "\n",
      "which leads to\n",
      "\n",
      "$$\\int_0^{\\frac{1}{2}} 2 x^2 \\cdot 2 dx = \\int_0^{\\frac{1}{2}} c 2 \\cdot 2   dx $$\n",
      "\n",
      "and a solution for $c$,\n",
      "\n",
      "$$ c = \\frac{1}{12} $$\n",
      "\n",
      "Assembling the solution for $x\\in[0,\\frac{1}{2}]$ gives\n",
      "\n",
      "$$ \\mathbb{E}(\\xi|\\eta) =  \\frac{2}{12}$$\n",
      "\n",
      "We can do the same thing for the other piece, $\\eta_2$,\n",
      "\n",
      "$$ \\mathbb{E}((2 x^2 - c x^2) \\cdot x) = 0$$\n",
      "\n",
      "which, by inspection, gives $c=2$. Thus, for $x\\in[\\frac{1}{2},1]$ , we have\n",
      "\n",
      "$$ \\mathbb{E}(\\xi|\\eta)= 2 x^2$$  \n",
      "\n",
      "which is what we had before."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Example\n",
      "-----------\n",
      "\n",
      "This is Exercise 2.6\n",
      "\n",
      "![alt text](files/ex26.jpg)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x,a=S.symbols('x,a')\n",
      "\n",
      "xi = 2*x**2\n",
      "eta = 1 - abs(2*x-1)\n",
      "\n",
      "half = S.Rational(1,2)\n",
      "\n",
      "eta=S.Piecewise((1+(2*x-1), S.And(S.Ge(x,0),S.Lt(x,half))),\n",
      "                (1-(2*x-1), S.And(S.Ge(x,half),S.Lt(x,1))))\n",
      "\n",
      "v=S.var('b:3') # assume h is quadratic in eta\n",
      "\n",
      "h = (eta**np.arange(len(v))*v).sum()\n",
      "\n",
      "J=S.integrate((xi - h)**2 ,(x,0,1))\n",
      "sol=S.solve([J.diff(i) for i in v],v)\n",
      "hsol = h.subs(sol)\n",
      "\n",
      "print S.piecewise_fold(h.subs(sol))\n",
      "\n",
      "t = np.linspace(0,1,51,endpoint=False)\n",
      "\n",
      "fig,ax = subplots()\n",
      "fig.set_size_inches(5,5)\n",
      "\n",
      "ax.plot(t, 2*t**2,label=r'$\\xi=2 x^2$')\n",
      "ax.plot(t,[hsol.subs(x,i) for i in t],'-x',label=r'$\\mathbb{E}(\\xi|\\eta)$')\n",
      "ax.plot(t,map(S.lambdify(x,eta),t),label=r'$\\eta(x)$')\n",
      "ax.legend(loc=0)\n",
      "ax.grid()\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Piecewise((2*x**2 - 2*x + 1, And(x < 1/2, x >= 0)), (2*x + (-2*x + 2)**2/2 - 1, And(x < 1, x >= 1/2)))\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAT8AAAE4CAYAAAAto/QTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXd4FOX2x7/vT8SgAgFBhECIUgSUC1IDAUGRIhEQuMLl\ngtLUACoIgghRyFUCUkUEKQGkCaEoLSBFTCChE0go0mJCkU6yS0tCSc7vj0liym62zeyUPZ/nmQdm\n592Zc2Zmv3nPeZsgIjAMw3ga/6e2AQzDMGrA4scwjEfC4scwjEfC4scwjEfC4scwjEfC4scwjEdS\nqPgJISoJISKFECeEEMeFEIOtlJshhDgrhIgXQryijKkMwzDyUcTG8YcAhhJRnBDiaQCxQojtRHQy\nu4AQoj2AqkRUTQjRGMBsAP7KmcwwDOM6hdb8iOgqEcVl/f8ugJMAKuQr1hHA4qwy+wF4CyHKKWAr\nwzCMbNid8xNC+AF4BcD+fId8AFzMtf83gIquGsYwDKMkdolfVsi7BsCQrBpggSL59nnMHMMwmsZW\nzg9CiMcB/AJgGRGts1DkEoBKufYrZn2W/zwsiAzDKAIR5a+A2cRWa68AsADAn0Q03UqxDQDeyyrv\nD8BMRNesGGjYbezYsarbwL6xf57on7PYqvkFAOgF4KgQ4kjWZ6MB+GaJ2Vwi2iyEaC+ESABwD0Bf\np63RMefOnVPbBMUwsm8A++epFCp+RBQDO/KCRPSxbBYxDMO4AR7hIRN9+vRR2wTFMLJvAPvnqQhX\nYmaHLiQEuetaDMN4DkIIkNwNHu5ACMGbC5s7iIqKcst11IL980xsdnVxB1wjdA53iR/DGBHVw96s\nKqtbbDAafO8YRsdhL8MwjD3I/XeexY+xidFzRuyf9jl1CmjRQl4BZPFjGEbzhIQA7dsDcqa5Oeen\nQc6ePYvjx4/j6NGj6NChA+rVq2exHN87xhM4ehRo2xZISACeeqrgcc75GYiIiAj4+Phg2LBhmDJl\nitrmMIyqjB0LjBxpWfhcgcXPDSxfvhxTp05F9+7dER4ebrP80KFD0ahRI1y8eBHPP/+8GywsHCPk\njAqD/dMuhw4BBw8CAwbIf25N9PPTOidPnsTYsWPRuXNn+Pv7OyRICQkJSE5OxmeffYabN2+iWrVq\naNy4sV3nWLt2LYKDg10xnWF0zVdfAcHBgJeX/Ofmmp8N9u3bh9atW6NmzZq4c+cOnnjiCYe+f+LE\nCUyaNAkAUKZMGVStWhWxsbE2v7dhwwYMHjwYly4VmBrR7bRs2VJtExSF/dMmMTFSK2///sqcnxs8\nbNCtWzd07NgRvXr1yvN5YmIiwsLCrH7P398fnTp1wsOHD3H69Gm8/PLLICJUqlQJERERqFu3rtXv\nrl27FuPHj4e3tzdatmxptfan9XvHMM5CBLz2GtC7N9DXxiR5zjZ4sPjZYPjw4fjzzz/Rr18/VK5c\nGQ0bNnT6XBEREZg/fz7WrZMmxN6wYQMee+wxREdHo3bt2tiyZQuCg4NRo0YNu87nrnsXFRWl29qD\nPbB/2mPHDmDgQODPP4EiNpJzzoof5/xsUKtWLWRkZCAzMxOPPfaY0+cxm81YtGgRli1bBgC4cOEC\natWqhapVq2LMmDH44osvULJkSfj6+splOsPoEiLgyy+B//3PtvC5guZrfnJ1anTGzVWrVqFixYpo\n2rRpgWP2hr3StQmjRo3CF198AW9vb5w/fx6VK1cGAFy7dg3du3d3qkVO67VmhnGGTZuAL74A4uOB\n/7OjVYLDXgXo1asXFixY4HAjR35mzJiBgIAA+Pj44MKFC0hLS0O5cuWQnp6OI0eOICkpCV9//TUi\nIiLw1ltv2X1eLd87hnGGzEygfn2plbdLF/u+w52cFeDdd99F8+bNMWHCBJw/f96pc8TExGDo0KFo\n2LAhKlSogCZNmqBq1arYtm0bNm3aBCJCeno61q5di2effVZmD+RBz/3E7IH90w6rV0uhbufOyl+L\nc36F0LZtWzRr1gy//vorevfujT59+jg8JXizZs2QkZFR4PPBgwfLZCXDGINHj6Qa36xZ8o7htQaH\nvXawceNGXLx4ETVr1sRrr72mtjk56OHeMYy9zJ8PLF8utfQ6In6c8/NA+N4xRiE9HahWTQp7/f0d\n+y7n/BjF0FPOyBnYP/WZPRuoV89x4XMFzvkxDKMqd+4A334L/P67e6/LYa+O4XvHGIGvvwbOnAGy\n+v87DOf8PBC+d4zeuXkTqFED2L8fqFLFuXNwzo9RDD3kjFyB/VOPb78FunVzXvhcgXN+DMOowt9/\nAwsXAsePq3N9Dnt1DN87Rs+8/z5QpoxU+3MFntWFYRjdcPIksH691NChFpzzY2yi5ZyRHLB/7mf0\naODzz4FSpdSzgWt+DMO4lb17pYWJli9X1w7O+clAUlKSS6usXblyBSVLlsSTTz7p0PeMcO8Yz4II\naNEC6NMH6NdPnnNyVxeVSExMxL59+1w6R9myZXMWOWIYI7N5M5CcDLz3ntqWsPi5zNy5c9GjR49C\ny0ycOLHQ40WKFEFgYCCWLFkip2myocWckZywf+4hIwMYNQoYP17Z6entRZPit+nMJpjTzXk+M6eb\nsenMJreeI5uePXuiefPmePfdd3O2jz/+GEePHkXFihULlD958iS6deuGFStWICkpCenp6Tav0bBh\nQ/zu7sGNDONGfv4ZKF4c6NhRbUuyICK3bNKlCmLpc1OaiQZFDCJTmsnivj3IcQ4iouTkZJo5cybN\nmTOnwLFx48bR8ePH83y2d+9e8vHxoTFjxtC8efPo0qVLFBISYte1hg0bRmfPnrXbNmv3lGG0Rloa\nUeXKRNHR8p8763fguCY58yWnLuSA+BH9I1ZJpiSnREuuc6xatYpOnz5NixcvLnCsU6dOlJmZmeez\nd955h5YuXZrnM3vFb/HixRQeHm63bSx+jF6YNo2oQwdlzu2s+Gkg8raMt5c3RgSMwPPfS62oPx76\n0elz/XjoRyQNSYK3l7fD301MTET16tWxbNkybN++HQDw9ttvo2vXrkhNTYXIN+Wsr68vli9fDi8v\nrwLr/B4+fBj79u3DpUuX4O/vj4cPH2LNmjVYntXmX6pUKZxRs9enFfS47qsjsH/KYjYDEyYAkZGq\nmWARzYqfOd2MybsnI2lIEibvnozQVqEOi5c53YzgHcEYETDC6XMUKVIE+/fvx7Zt2wocs7Q2R2Hr\n/F6/fh0vv/wytm7ditDQUGRmZuLTTz/NOV6sWDE8ePDAIfsYRutMmAB06gS89JLaluTDmeqiMxt0\nmPNLTEykDRs20A8//GDxeJs2bfLsr1y5knbv3l2gXO6wd/To0bRq1SoikvKDXbp0yTm2Zs0ai7lF\na1i7pwyjFc6fJypdmujSJeWuASfDXk229u6+sDtPLc3byxuhrUKx+8Jut54jOjoajRo1QrFixSwe\nf+6553D37t2c/Q0bNqB+/fqFnnPHjh1o1aoVAGDLli151um9cuUKqlatard9DKN1vvwS+OgjoEIF\ntS0piCbFL7B6YIHw1NvLG4HVA916jt27d6Nv375YvXp1nm4ukVnJixYtWuDAgQM55W2t85uWlobi\nxYujdOnSAIDt27ejffv2Ocfj4uIQEBBgt33uQiv9xJSC/VOGI0eAbduAESNUubxNNJvz0wJz584t\n9HiXLl0wZcoUvP766wBsr/NbrFixnEaTlJQUPHr0COXKlQMApKeno0SJEvDy8lLGGYZxI0TSxAVj\nxkh9+7SIJmt+esHb2xtlypTBzZs3cz576qmn4O3tjW7duqFy5cpWv7tv3z50zrUsfXh4OIKCghS1\n11mM3BIKsH9KsG0bcOEC8MEHbr+03XDNz0WGDBmC+fPn44NcT7lDhw42v9e+ffuckPfixYsoVaoU\nXnzxRcXsZBh3kZEh1fq+/RZ4/HG1rbEO1/xcRAiRR/gsYesvb6VKldCpUycZrZIXzonpG3f7t3Sp\nFOq+/bZbL+swXPNzAy1atFDbBIZxC6mpwFdfAatWAcLhSabcC8/np2P43jFaIzQUiIsDVq923zV5\n3V4PhO8doyWuXgVeftm1NXidgSczZRSDc2L6xl3+jRkjzdCsxhq8zsA5P4ZhXObYMWk1ttOn1bbE\nfjjs1TF87xgtQAS0bQt06AB88on7r89hL8MwqrBlC3D+PDBggNqWOAaLH2MTzonpGyX9e/QIGD4c\nmDxZ2x2aLcHixzCM08yfD5QrJ4W8eoNzfjqG7x2jJrdvA9WrA7/9Brzyinp2cM5PBZKSkmyWuXLl\nClJTU91gDcO4lwkTgHbt1BU+V2DxcxJ7Fys3woLknBPTN0r4l5gIhIVJa/DqFRY/J7FnsXJA+wuS\nM4wzjBgBDB2qzRma7YXFzwni4+MtLlZuDb0vSM7z3ekbuf2LigJiY4Fhw2Q9rduxKX5CiIVCiGtC\niGNWjrcUQtwSQhzJ2r6U30xtERERkTN7s72ULVsWCQkJClnEMO4hIwP49FNg0iTAytI2usGemt9P\nANrZKLOTiF7J2sbJYJemOXjwIGrVquXQd+rUqYPY2FiFLFIWzonpGzn9W7gQKFECeOcd2U6pGjbH\n9hJRtBDCz0Yx5WbukmtSMCe6hBw5cgR79+7Ns8j46tWrLS5WvmHDBjz22GOIjo5G7dq1sWXLFgQH\nB6NGjRoAtLsgOcPYy61b0lx9mzZpf64+e5Aj50cAmgoh4oUQm4UQjlWJbJ6d5Nmc4Nq1a3j55Zdx\n/PhxdOjQAW+//TZiYmKQmZmZp9yFCxdQq1YtBAYGYvv27QgMDET37t3h6+ubU0bPC5JzTkzfyOVf\naCjQvj1gY3VW3SDHrC6HAVQiolQhxJsA1gGoLsN5Vaddu3YIDg5Gr169AAAHDhxAo0aNcO/evTzl\nskXu2rVrKF68OLy9vfOsxwsAt27dylmykmH0RkKCFPIeP662JfLhsvgR0Z1c//9NCPGjEKI0EaXk\nL9unTx/4+fkBkFY+q1u3rquXV5wdO3bgs88+AwBs3boVgYGB2LlzJ+7evYunn34aAHDq1Cncv38f\nhw8fxquvvgpAahTJvyB5zZo1FbMzO6+T/Vdezv3cOSMlzq/2Pvtne79Pnyh06QI895w2/Fm0aBEA\n5OiJM9g1vC0r57eRiGpbOFYOwHUiIiFEIwCriKiARXoc3paWloaOHTvmrLUbEBCAX3/9FZs2bYKf\nn19Oi++MGTNw584dlC9fHqdOnUKTJk3g4+ODRo0a5Zzr/fffx8yZM2Vdl9dd9y4qKsrQoSH7Vzjb\ntwNBQcCffwJaXFba2eFtIKJCNwArAFwG8ADARQD9AAQBCMo6/hGA4wDiAOwB4G/lPGQJa59rjeTk\nZGrUqBEREZlMJgoODrb7u2lpaTR06FDZbdLLvWP0y4MHRDVrEq1bp7Yl1sn6HdjUsvybPa29hQ5j\nIKJZAGY5rLo6I/ci47kXKy9TpozN72p5QXKGKYxZs4BKlYCOHdW2RH54hIedtG/fHl988UXO/pAh\nQ7B27Vqb3zPCguTcD07fOOvf9etSC+/06cbo2pIfXsPDSexZrByQFiSvVKmSGyxiGHkZNQp47z1A\nwXY6VeH5/HQM3ztGKQ4cADp1Ak6dAkqWVNuawuH5/BiGkYXMTGDwYGm6Kq0Lnyuw+DE24ZyYvnHU\nv6VLpUFRvXsrY49W4JwfwzA53L4t5frWrQP+z+BVI8756Ri+d4zcDB8OJCcDP/2ktiX242zOTxM1\nv/wzpDAM435OnAAWL5b+9QRUr9g60zNbi1tkZKQq13UHnBPTN/b4RwR8/DEwdizw7LPK26QFVBc/\nhmHUJzwcMJuBAQPUtsR9qJ7zYxhGXe7ckToyr1oFNG2qtjWO42zOj8WPYTwcPTZy5IY7OauMkfNG\nRvYN8Gz/TpwAliwBJk50nz1agcWPYTwUT2zkyA2HvQzjoYSHSzW+Q4eAxx5T2xrn4ZwfwzB2c/s2\nUKuWfhs5csM5P5Uxct7IyL4BnunfmDFAu3b6Fz5X0MQID4Zh3MeRI1LI6ykjOazBYS/DeBAZGVJt\nLygI6NdPbWvkgcNehmFsEhYGFC0K9OmjtiXqw+InE0bOGxnZN8Bz/Lt2Tcr1zZ5t/Omq7IFvAcN4\nCCNGSDW+l19W2xJtwDk/hvEAIiMl4TtxAnj6abWtkRfO+TEMY5EHD4BBg4AZM4wnfK7A4icTRs4b\nGdk3wPj+DRwYhWrVpNXYmH/gfn4MY2DOnAHWrAGOHVPbEu3BOT+GMShEQKtWQIcOwNChalujHJzz\nYxgmD0uWALduAZ98orYl2oTFTyaMnDcysm+AMf27cQMYOVLq1BwTE6W2OZqExY9hDMhnnwE9ewL1\n6qltiXbhnB/DGIzffwfefx84ftwzurZwzo9hGKSlSSuwzZrlGcLnCix+MmHEvFE2RvYNMJZ/48ZJ\noW5g4D+fGck/OeF+fgxjEI4elRo44uPVtkQfcM6PYQzAo0dAkyZSyNu/v9rWuBfO+TGMBzN9OlCi\nhHEmKHUHLH4yYeS8ipF9A/Tv319/Ad9+C8ybBwgL9R+9+6cULH4Mo2OIgA8+AEaNAqpUUdsafcE5\nP4bRMfPnA3PnAnv3AkU8tPmS1+1lGA/j8mWgTh1gxw7gX/9S2xr14AYPlTFyXsXIvgH69I8I+Ogj\naRU2W8KnR//cgYdWlBlG3/zyC3DyJLBihdqW6BcOexlGZ9y8CdSuLQlg06ZqW6M+nPNjGA+hRw+g\nQgVg6lS1LdEGnPNTGSPnVYzsG6Av/9auBWJjgW++sf87evLPnXDOj2F0QnKy1MixahXw5JNqW6N/\nOOxlGJ3w7rvAM89IQ9mYf3A27OWaH8PogI0bpY7MPGOLfHDOTyaMnFcxsm+A9v0zmYCBA4EFC4Cn\nnnL8+1r3Ty1Y/BhG4wwdCnTuDLRoobYlxoJzfgyjYTZskMQvPp6npbcG9/NjGINx86Y0dG3lSqB5\nc7Wt0S7cz09ljJxXMbJvgDb9IwIGDQL++1/XhU+L/mkBbu1lGA2ycqW09OSSJWpbYlw47GUYjXHl\nClC3LhARATRsqLY12ofDXoYxANkzMwcFsfApDYufTBg5r2Jk3wBt+ffTT9IkpV9+Kd85teSfluCc\nH8NohPPngZEjgT/+AIoWVdsa48M5P4bRABkZwOuvA+3bSwLI2A/n/BhGx0ybJuX7hg9X2xLPgcVP\nJoycVzGyb4D6/sXHA5MmSd1aHntM/vOr7Z9WsSl+QoiFQohrQohjhZSZIYQ4K4SIF0K8Iq+JDGNc\n0tOBnj2lWZn9/NS2xrOwmfMTQjQHcBfAEiKqbeF4ewAfE1F7IURjAN8Tkb+FckREMKebsfvCbgRW\nD5TJBYbRL8OGARcvShOUCoezVgygYM6PiKIBmAop0hHA4qyy+wF4CyHKWSpoTjcjeEcwAnwDHLWT\nYQzHjh2S6M2Zw8KnBnLk/HwAXMy1/zeAipYKBu8IRmirUHh7ectwWW1h5LyKkX0D1PHPZAL69pXm\n6HvmGWWvZfTn5yxy9fPL/3fLYix9eellTD8ozcHt7e2NunXromXLlgD+eUB63Y+Li9OUPbyv3X0i\noGvXKNSvD7Rtq749etuPiorCokWLAAB+LiRK7ernJ4TwA7DRSs5vDoAoIgrP2j8FoAURXctXjipO\nq4j97+9HheIVnDaYYfTOkiXAxInAwYO8EJEcqLmGxwYAHwMIF0L4AzDnF75smlZqikZhjbD//f3w\nKeEjw6UZQ0AELFoEXL9e8FiLFoB/gfYz3ZKQAHz2mZTvY+FTF3u6uqwAsAfAi0KIi0KIfkKIICFE\nEAAQ0WYAiUKIBABzAQyydq7lXZajSukqeGPJG7hz/07O5+Z0Mzad2eSqL6qSXS03Ior6RgR8/DEw\ncyaQkpJ3u3ED6NgR2L5duevDfc/u4UNpfr6vvpImKXUXen43N53ZBHO6Oc9nsukFEbllky5FlJya\nTLVm1qJGYY3oVvotMqWZaFDEIDKlmUjPREZGqm2CYijmW2Ym0aBBRP7+RGaz5TLR0URlyxJt26aM\nDeS+ZzdqFNGbb0puuxM9v5u59SEzM9OiXmRpi8OapMrY3pS0FLz606so8n9FUL9CfUxtM9WQLcBM\nIWTX+A4fBrZsAUqWtF42Jgbo0gX4+WegdWv32SgjUVFSrS8uDnj2WbWt0RfmdDP6rOuDG/duoM5z\ndTC+1fg8eqGrsb2li5XGxh4bEX8tHjEXYpD2ME0NMxi1cET4AKBZM+DXX6WhEAqHwEqQkgK89x6w\ncCELnzNEn49G9Plo7Pl7Dz4P+Fy2ipIq4mdON2PKnilIHJwIby9vBCwMwDnzOTVMkQ0951VsIatv\njgpfNgoKoJLPLnty0n//G2jXTrHLFIqe383lx5aj/4b+eP3515E0JAmTd08ukAN0FreLX/Yoj9BW\noXi+1PPY2msrfEr4oMG8Bjh46WCBsnpvCGFy4azwZaPDGuDcuUBiIjBhgtqWaJ/8jRtzD83FsK3D\nULNMTYR1DIOftx9CW4UieEewPALoTKLQmQ1ZDR4RpyMKNG6Y0kw0YOMAejL0Sdp1bleBRCdjAOxp\n3LAXNzSCyEF8PFGZMkSnT6ttiT7I/ZufFDOJfKf5UsufWtI507kC5SJOR+TsQ08NHtZYErcEAzYN\nwOK3FyPqXJRhh8J5HK7W+Cyh8UaQe/eABg2A0aOBd99V2xr9YEoz4Y2lb+BW+i00822G6e2m29QA\nXTV4WOO9uu9hzltz0G1NN9R5ro6uhE/PeRVbuOSbEsIHyBoCK/HsPv4YaNxYG8Knl3fzUeYjDN82\nHI8yHuEv018IaRmiqAZoSvzM6Wbs/3s/1v9nPYZuHYrp+6arbRLjCkoJXzYazQEuXQrs3Sv122bs\nI/VhKjqv7Izzt86jccXGsjduWMSZWNmZDVk5P3vifSKi2MuxVGxcMfr0t08pM1ev0PzxPqNR5Mzx\n2SI6WkquaSAHeOqUZEp8vNqWaJf8ef+b925Sw3kNqfnC5jRg44CcY/bm/eFkzk8z4mepIeTA3weo\nxPgS1HdtX3qU8YgbQfSCO4UvGw00gqSlEdWpQ/Tjj6qZoAty/47Pm89TtRnV6JU5r1D4sXCLjaG2\nKju6Fz9rnDefJ5+pPtR6SWv6cMOHmhU+PQ8hsoVDvqkhfNk4KYByPbuPPiLq2tX9w9dsocV305Rm\nom6rutFzU56jgAUBLv2unRU/za/b61vSF3/0/gMvznwRdcrVwf1H99U2ibGG0jk+W2TnAFVoBV65\nUnL50CGeldke9l7ci98Tf0dKegr29t+rTuOmM4rpzAYna37ZVeTElERqOK8hVZpWiX488KNT1WNG\nQdSs8eXHzSFwdp7v8GG3XE5XWEpnTd0zlUpMKEFdwrtQkinJ5VQWjBj25s/xmdJM1GpxKyr1bSnq\nuLyjw4lRRiG0JHzZuEkA790jevllorlzFb2Mbsn928zIzKCPN31MJcaXoHdWvSPb79eQ4mdtNMj4\n6PFUZlIZem3Ra7L85ZADLeZV5KJQ37QofNnYKYDOPrvMTKLevYl69dJeni83ar+bpjQTfbjhQ2q3\ntB2Vn1Ke5h+aL2vk5qz4aTrnZ2l5S28vb4xqNgqda3RGm6Vt8Pz3zyPhkwRddYg2DGrn+GyhcA5w\n4UJpKvoDBzjPVxj3HtzDvr/34ej1ozj10Sm8WObFAmW8vbzdvpytpoa3OYI53YxhW4fhzxt/4urd\nq9jVdxd8S/rKdn7GBloXvtwoMBQuPh544w1g506gVi1ZTmlIDlw6gLfD34ZvSV+s6LoCU/ZMkX3Y\nqrPD2zQd9lojd47gwaMH1HdtX3o69GnambSzQDluBFEALYe61pAxB2g2E1WtSrRsmQx2GQRLKap5\nh+ZR8fHF6c1lbyqan4cRc37WsHSjv476mop+XZTWn1pPRO5vBFE7r6IkeXzTo/BlY0UAHXl2GRlE\nnToRDRwos20K4o53M3/Dxqe/fUrFxxenSTGTFO+Z4az4aTrnZw1LuYGvWnyFOs/VwX/W/AdfBHyB\nq/euFpjumnERPYW6lsjOAXbuDCxf7lQIPHEicO0asGqVAvbpGG8vb4S2CsWIbSOQZE7CyRsncSTo\nCKqUrmKxrLvzexZxRjGd2SBjza8wopKiCCGgLiu70L0H99xyTY9AzzW+/DgZAm/bRlS+PNHffytk\nl845fu04PT/9eUII6PRN901iCCdrfpqa1cVVzOlmrDqxCicGnsCxa8fQOKwxFhxeoNzSd56C3mt8\n+XFiNpjz56XpqZYvB3w8fMlpS8tJLjy8EE0WNMELpV5A0pAkfL/ve2VnZJEDZxTTmQ0K1/zy5/hS\nUlOo2cJmVHpiaQr8OVDxDtGGzfllZlJkp07GqPHlJ6sGGDl5cqHF0tKI6tcnmjLFTXbJjNzvZv4G\nx4ERA+np0Kfp7fC3VRl4AE9q8LCEtQ7RE2Mm0nNTnqOG8xrSXyl/KfZADCl+WaFuZM2axhO+bKKj\nKbJkyUJD4P79ibp103ZH5sJQ4t00pZmoz9o+1CisEfl+50vzY+XtuOwIzoqfbvv5OcKVO1fQMbwj\nDl0+hIMfHESDCg1UsUNXGC3ULYxC+gGGhQHTpwP79wNPP62SfRrkj6Q/0OOXHrh+7zoSPkmw2LDh\nLgwxjb1SFHu8GBpUaICBDQai5aKWGBs5lvOAheFJwgdYzQHu2QMEBwNr13qu8OXP72VkZmDk9pHo\nsrILGvtIMy5P2ztN+/k9SzhTXXRmg5tae/OTP/ew5sQaKjauGL0y5xW6ee+mxTLOYJiw10KrrmF8\ns0KOf7lagS9eJKpQgWjTJlVNkwVXnl/u38al25coYEEAlZ9Snnqs6aGZiUXg6Tk/a1jKBf55/U+q\nNbMW+Uz1of0X98vy4AwhEFa6sxjCt0LI4190NGWWLUsfVd9GEyaoZpKsuPr8TGkm6rC8A5WZVIYa\nzmtIy48u19SUcs6Kn0fk/CyRkZmBYVuHYcaBGQh7Kwzv139fbZPUxdNCXSsQAaFvxuCTqC4oseFn\niDbaWxbTndx/dB+jd4zGz8d+xrV715A0JAl+3n5qm5UHzvk5yJ0Hd/Ao8xFW/XsVhm0bhr7r+uKX\nP3/xzFzCiNvsAAAdWElEQVQgC18O338PrLnaDEU3/grRS1urwimJpb57ey7uQY1ZNXAm5QwCqwW6\nZ0U1d+JMddGZDSqFvZbIn6M4ZzpH1X+oTn7f+dG/V/3bqVyGbkNDO0Zu6NY3O8n2b/t2onLliJKS\nsg5oaFU4V7Dn+eUfmzt+13jyGudFE2Mm0sCIgZrJ71kCnPOzH2t9Aj/f9jk9M/EZahzWmBKSExx6\nyLoUCDuHrOnSNweIjIykhASiZ58lKuCqBlaFcxV7n58pzUS91/amVxe+Ss9OfpZiL8da/a1oabYk\nFj+ZuGC+QE3mNyGEgLYnbFfbHOUw0lhdFzGbiWrUIJo1y0oBAwigLTIzM2lJ3BIq/W1pQgjozM0z\naptkN86Kn8fm/KxR/IniqPtcXXzT8ht0CO+Anr/0RHJqcp4yus8Dco4vh0ePgO7dpYlJBw2yUsiJ\nscBaJn9+7/Kdy3jz5zcR/EcwWr3QCklDkjB933Tj5Pas4YxiOrNBBzW//PmM+KvxVH5KeSo7qSzt\nu7jPYplsdBMaOlHj041vTjB4MFGDBpH08KEdhXVaA8z//LLf4ZTUFFoct5iemfgM1Z1dl/qu66vp\n3J41wGGv61jKbySnJlO/df3oiW+eoBFbR9CAjQMsvhC6EAgnQ11d+OYEs2cTvfgi0caNkfZ/SYcC\naOn5Hb16lCp/V5lqzKxB76x6h8KPh2s+t2cNFj+FiT4fTQgBVf+hOu25sEcXieA8cI4vDzt2SA0c\nZ8868WWdCKCld/TmvZv0/ob36ZmJz9Cnv31KCAElmZLUMVAmnBU/zvnZgTndjBXHViBxcCL8Svqh\n88rO+PXkrxi+bXhOXsScbkbwjmAE+AaobK0FOMeXh7NngR49gPBwoGpVJ06gkxxggG8AgncE57yj\nO8/tRM1ZNXHqxils6bUFDzIeGK/vniM4o5jObNBpzc/Swun91/en3r/2puemPEdtlrahxJRE6jS+\nkzbzIzLU+IwU9t64QVStGtG8ef985rR/OqgBmtJMFBgaSP3W9aNi44rRrAOzKCU1pcA7rZf8niXA\nYa8yFBbexpyPoWozqhFCQJOXTdZeKCxTqGsU8UtLIwoIIBo5Mu/nLvmnEQG09O6lpKbQ8K3DqdSA\nUoQQ0KFLh6yW1XTKxgYsfipgSjPRgI0DaMwfY8hrnBf1XduX3l//vjb+onKOLw8ZGUT/+Y80KWlG\nhswn14AA5n/Xdp3bReWnlKdaM2tR5/DOlGRK0nXtrjBY/NxM/pftzM0zVHNmTSozqQy9tug1RWeN\ntgkLXwFGjyZq0oQoNVWhC2hgKJwpzUT91vWjnr/0pGLjitH4XeM1PzRNDlj83Ez+0CEyMpJMaSaa\ntmca1Z1dlxACWn50ufsNU0D49B72hoURValCdP265eOy+adiDTD9YTpN3TM1Z4TGkStHct7R3P7p\nOby1hrPix629ThJYPbDAmsDeXt7o+0pfNKnUBDPazUBQRBDqza2HvRf35imn2AgR4lbd/GzfDnz5\nJbB5M1C2rMIXc0MrcP7RGZmUifmx8+E73Rfb/9qONlXbIGlIEsJiwxDgG2DxHdXEmrlawBnFdGaD\nwWp+lsgfVly9c5UahzWmJ755gt795V26dPuScqEHh7oFiIuTKmK7drn5wgrWAHO/P5FJkVR3Tl0q\nO6ksrTi2wlAtuI4ADnvVx1or2oLYBfTKnFeo5ISSVHdOXUpITpC3xY2FrwDnzhFVrEi0cqVKBsgg\ngNbekXE7x1GlaZWo4tSK1HpJa0pOTTZcC64jsPipjK28UZIpiRAC6vlLTyo9sTQN2zKM+q3v5/pf\najcIn95yfjdvSrO0TJ9uX3nF/HNRAPO/EzuTdlLl7ypThakVaNzOcXaPztDb83MUZ8WPc35uwJxu\nxuTdk5E0JAklnyiJyN6RuH3/NtaeXIvWS1vj6LWjCN4RjNBWoQVyNIXCOb4CpKUBHTsCHToAQ4ao\nbIyLOUBvL2+EtgrFgIgBaLu0Ld5c/iYGNRyE2A9jcfnOZc8enSEHziimMxsMXvOzhqURItn7CckJ\n1HVlV0IIKGhDEC2NX2p/6MKhbgEePiTq1ImoZ08F+vK5go0aoLWQdVLMJHpr+VtUZlIZQgjo5PWT\nhb5Pngo47NUmheVisl/c6PPRVPvH2lRifAmqNasWHb58OKecxRebha8AmZlEQUFErVsT3b+vtjUW\nKEQAcz/nzMxMWnl8JZWfUp58p/nS1D1T6cMNH+Z0Utbz7CtKweKnMo7mVSz9Be+7ri99tOkjeuKb\nJ+itn9+iriu75hFKIsoRvoeNGtKWQ+7J5ushZxQSQlSvHtHt245/123+RUdTeumSdCfi1zwfm9JM\ntCxuGbVa3IpqzqxJz0x8hsIOhdGNezdkqeXp4fm5Aoufyjj6ghVWIzx69SghBOQz1YcahTWieYfm\nUdDGIDKlpuQI37CV/d0W6mj9xzNjBlHVqkRXrzr3fXf6d/v3zXS7ZLEcATx5/SQ1mNuAyk4qS80X\nNieEgBJTEolIvjG4Wn9+rsLiZxCy/7onmZJoYMRAWha/jF5b9BqVm/QsrW75LN2p9zINW9mfw58s\nli2TurTkrLimEQoTrlu/byZz8aI0ZlQTKvpNUeq3rh/tvbg357l7eg7PUVj8DIDVZHZqCt3s051O\nV3uGSnwBarOkDa08vtIjxm0WxsaN0lKTJ06obUlBLD3Lfuv70YToCVT7x9rU5eNn6dqToITw2dyI\n4SIsfiojR2hhsbaQmkJJ/w3MCXWPXT1GLX9qSfXm1qMKUytQg7kNKPp8dM6PRYnOrloMm3bulNoP\n9u93/Vyu+mftnocfC6eBEQNpxbEVVG1GNSoxvgT1WNODNpzaQAMjBtLlzavodsliFDlvtKK1eC0+\nPzlh8VMZRV6wXI0buXN8uVuJe6/tTQgBNZzXkOYcnEN/Jf8ley1Caz+ew4cl4fv9d3nO56p/lmpu\n3Vd3p09/+5R8pvoQQkBj/hhDN+/dLFA2fw5QCbT2/OTGWfET0neVRwhB7rqWIaB/OjBvnTkUjV9q\nk6cDtDndjK0JW7Hr/C4M8R+CoVuGwutxL+xI3AH/iv7IpEx81/Y7/HjwR7zq9yraVmlb4Pu7L+zW\n3SD3U6eA118HZs4EunRx77U3ndlUYLKA7PsY4BuATzZ/gvJPl8eCuAV48vEn0bVmV1y9exXfvvEt\nJu+ejNBWoTllc5/jzo7fULRbDzwRvhpo3dq9ThkAIQSISDj8RWcU05kNBq/5yYod/fis5YkumC/Q\n0vil1OKnFoQQ0GuLXqMZ+2bodlnC3CQkSI0bixerc/389y0lNYW6repGI7ePpDqz6+RMJ7Xy2EpK\nTk12rAaugQlR9Qo47FUX2UILOzsw29N5Ov5qPLVe0preDn+bSowvQeWnlKdRv4+i7qu7U0pqit35\nQS2ETefPE1WuTDRnjvzntuSftXuzJG4JdVjegfqu60vFxxcn32m+NGzLMNp0ehMNjBjoWmdkhQRQ\nC89PSVj8VEaWF0yGkRvWaoRX7lyhhYcXEkJAFadWpPJTylOPNT3ojSVv0NmbZ3PKtl/Wns6ZzuU5\n58atG1XtQnPpktSP77vv5D93xOkI2rh1Y57Pshsrshf2PnbtWM5MKk+FPkX15tYjhIA2n9lMmZmZ\n8rbWKiCALH5Oih+AdgBOATgLYKSF4y0B3AJwJGv70sp53HAbdIxMQ9YKa3nM3Y8s9nIszTowi9ov\na09Fvy5K1X+oTi/Neom+2/MdvfvruwV+yGr1K7x2jahmTaLx410/l6V7c850jgJ/Dsz5/Ma9G/Tv\nVf+m8bvGU4flHajYuGLkM9WHas2qRUvjltI507kC/fFkb2HnENghFBE/AI8BSADgB+BxAHEAalJB\n8dtg80IsftZReKyurZrJmZtnCCGgr/74irqu7EplJpWhp0Ofprd+fouazG9C606uo8SUxALnsFRL\nNKWZaOwfY2URg+Rkojp1iL780noZS8ITfjycwo+FF7h+9h+A3Dm7nr/0pFn7Z1G9OfWo6fym9PjX\nj1PNmTVpwMYBtCx+Wc5i9UmmJPf2x2MBtBulxK8JgC259r8A8EW+Mi0BbLR5IYOLn9OhhRsmKbAn\nP5i7JpOZmUl/JP5BCAH1XdeXan9eOye/9cL3L9DgzYOp9ZLWFHE6gj7Y8EEBMciuHTkrlOHHw2nB\nvnCqX5/os8+kW2StbP6amynNRP3X9af+6/vnEbk+a/vQxtMbadqeafSvH/9F/mH+VPSbouQz1YcC\nvgqgT3/7NGfti9z25b43bq/9yiSAHPY6J37/BhCWa78XgB/ylWkBIBlAPIDNAGpZOZcbboN6OPWC\nqTw7i7WaTP7QbuPWjZSRmUGnb56mGftmEEJA7Ze1p9o/1iavb7zo6fFPU8CCAKoxswaN3D6SFh5e\nSOtOrqP//vJfOnr1KA2MGGhREK0JZa9V/emZvv1p4FBTjvBZKzsoYhAlJCdQ33V9adOZTdR+WXua\nFDOJBm8eTFVnVKW6s+tS0W+Kkve33hSwIID6r+9Po34fRQgBHbx0kIiknGb+PwCaGXUhw6pwLH6W\nt0L7+QkhugJoR0QfZO33AtCYiD7JVaY4gAwiShVCvAngeyKqbuFcVNi1PA5SfyJSS/3WzpvP46PN\nH2FZl2Xw9vKGOd2cM9EqAATvCMaIgBE5/daKFy2O3Rd3o8WiFhj/+nikPkzFhdsXcOHWBfyV8hcu\n3r6IIqIIyj5VFqW8SuHW/Vt4qexLSDInoUnFJihVrBT+T/wfdl/YjVcrv4rIxGhc2tUalSsDqLId\n/hX9sev8LjQo3wAZlAHzfTNiL8eidLHSSDQlAgDuPriL4kWLw3zfjJZ+LVHjmRqo7F0ZxYoUw6db\nP0Xsh7GoV74eAOT4k+3D5wGfY9LuSTkTyWYf11TfyJgYqVPjzz9zP0ALONvPz5b4+QMIIaJ2Wfuj\nAGQS0cRCvpMEoD4RpeT7nHr37g0/Pz8AgLe3N+rWrYuWLVsCAKKiogDAM/aJENW5M3DmDFru3QuU\nLKkZ++5VuIcA3wDE7YvLOW5ON2Paimk4evUoFg1dBG8vb0Rsi8D82Pn4fuD3mLR7EppTc4QfCy9w\nfPrA6fg25lv4P/THI3qEYlWLodfaXhj4zEB4Pe6FSv+qhNSHqTi4+yDWn16PMhiI58t7o2qZ87j3\n4C42PNiAz5t+jrSzaXiiyBOo16Qe7jy4g6AfgjClzRT0frs3BAT6f98f/6n9H0SLaIS2CkXMrpic\n60/ePRlvPv4mAOC3h78htFUo4vbF4e6Du5h9YzZ+DPwRSXFJefydvXo2mlRqovrzyNn/4Qfgq6/Q\ncrXUEVp1e1Tcj4qKwqJFiwAAfn5++N///ueU+NkKe4sA+AtSg0dRWG7wKId/RLQRgHNWzqVQpVcb\n2B1a6HAi0sjISLtaSm2Fp5ZaSrOP9/9lEL3cLInqjBpEKammAvm2/OfK/tzStfLn/Gy1WI9fIkNT\nsjtwMgfIYa8TOT/pvHgTwGlIrb6jsj4LAhCU9f+PABzPEsY9APytnMctN0It7HrBdCh8RNZ9s9aQ\nYm/DRLZ49ftlENVpbKKhQ4lSUq2Ll70NKdZae601TOhKHJwQQF355wSKiZ9cm9HFzyY6FT65sCaU\nn20cS7XqmWj4cOkWEVkXL7m60Oge7gaTB2fFjyc2cAcaaNzQIn//DbzxBtCjBzBmDCAcz9p4LtwI\nkoOzDR68dKVMZCdkC2AA4bPqmwucOwe0aAH06weMHauu8Cnhn+I4sCymLv1zAyx+SmIA4VOCM2eA\nV18Fhg4FPv9cbWt0jIvrAns6HPYqBQufRY4fB9q2Bb7+GujfX21rDIKHh8Ac9moJFj6L7N0LtGoF\nTJnCwicrXAN0ChY/mcjJqxhQ+OTIGW3ZAnTsCCxaJDVwaAlD5MQKEUBD+KcALH5yYkDhk4MVK4De\nvYH164E331TbGgPDNUCH4JyfXLDwWWTWLGDCBOC334DatdW2xkPwsBygImN75cTQ4sfCVwAiqVFj\n6VJg2zbghRfUtsjD8CAB5AYPtcgSvqjISMMKn6M5o0ePgKAgKcyNidG+8BkyJ5YrBI6aMkVtazRJ\nEbUN0DW5a3yTJhlS+Bzl7l2ge3cgIwPYuRMoXlxtizyYbAF86y2gTh3D1wAdhcNeZ+FQtwDXrkm/\ns9q1gblzgccfV9siBoDhQ2AOe90JC18BzpwBmjYFAgOBBQtY+DQFtwJbhMXPUawInyHzRlnY8m33\nbmm42ujRQEiI/iYoMPKzA7L8YwEsAIufI3CNrwBLlwKdOwOLF/OoDc3DApgHzvnZCwtfHjIzpWmo\nli8HNm4EXnpJbYsYuzFYDpD7+SkJC18eUlOlERtXrgBr1wJly6ptEeMwBhJAbvBQCjuFz8h5o9y+\nXbkizcPn5QXs2GEM4TPyswOs+MchMItfoXCNLw8HDgCNGgGdOgFLlgBPPKG2RYxLeLgActhrDRa+\nPCxaBIwYAcyfL4kfYyB0HgJzzk9OWPhyePgQ+Owz6TasWwfUqqW2RYwi6FgAOecnF04KnxHzRjdu\nSL+DAweicOCAcYXPiM8uN3b554EhMItfbrjGl8OhQ0DDhtJvIjQU8PZW2yJGcTxMADnszYaFD4B0\nG+bMkVZUmzNHioQYD0NnITDn/FyBhQ+ANCNLUJC0yNCaNUC1ampbxKiGjgSQc37OIpPw6T1v9Oef\nUjcWLy9g3768wqd332zB/lnAA0JgzxY/rvEBkP64t2ghdWVZsAAoVkxtixhNYHAB9Nywl4UPd+4A\nn3wi1fRWrpTmu2SYAmg8BOaw1xFY+BAbC9SrBxQpIv2fhY+xikFrgJ4nfgoJn17yRpmZwNSp0hKS\n48ZJIzaeeqrw7+jFN2dh/+zAgALoWWt4eHiN78oVoG9f4PZtaZyun5/aFjG6IlsANRwCO4Ln5Pw8\nXPhWr5bcDwqS5uEr4ll/9hg50VgOkPv5FYYHC5/JJLkeGyvNxNKokdoWMYYgJkaawnv5ctUFkBs8\nrOEm4dNi3mjrVuBf/wLKlJHcd1b4tOibnLB/TtCsmTSTrY5zgMYWPw+t8d2+DQwYAHz4oTQV1fff\nA08+qbZVjOHQeSOIccNeDxW+iAhg0CCgXTtg8mSPcZtRE5VzgJzzy40HCt+NG8CQIVIrblgY8Npr\nalvEeBQqCiDn/LJRSfjUyhsRSTnn2rUBHx/g6FH5hY9zYvrGLf7pMAQ2VocHD6vxnTkjuXvtmrR8\nZMOGalvEeDQ66wdonLDXg4QvNRWYMAGYPRsIDpbG53K/PUYzuDkE9uyw14OEL3uB8LNngfh4YOhQ\nFj5GY+gkBNa/+GlE+JTOq5w5A3ToIC0mNG8eEB4u5fjcAefE9I0q/ulAAPUtfhoRPiUxmaTaXdOm\nwKuvAseOaT6VwjASGhdA/eb8DC58Dx8Cc+cCX38tpU++/hp49lm1rWIYJ1B4KJxn9fMzsPARAevX\nA6NHS2HttGlSNxaG0TUKNoJ4ToOHRoVPjrxKZCTQpIm0ctrkycC2bdoQPs6J6RtN+KfBEFhf4qdR\n4XOV2FigbVvg/feBwYOBI0eAwEBAOPy3jGE0jMYEUD9hrwGFLy5Omk15717gyy+B/v2BokXVtoph\nFEbmENjYYa/BhO/QIaBTJ6B9e6kV98wZYOBAFj7GQ9BIDVD74qcT4bMnr7J3ryR4nTtLf/D++gsY\nNsz2Ghpqo4mckYKwfyqgAQHUtvjpRPgKIzNTar199VXgv/+VanwJCZJbvD4u49GoLIDazfnpXPjS\n0oDFi4HvvpNMHz5cSnPwUDSGyYeLOUBj9fPTsfBduCANPwsLA/z9peFozZtzyy3DFIoLAmicBg8d\nCl9mJjBpUhQ6dQJeeQW4cwfYteufcFfvwqfJnJGMsH8aQIUQWFtBWLbwxcZKq+9oXPiuXgWWLpWG\noQkBjBwpjeDRegMGw2iSbAHs3FmqAbZpo+jltBP26qTGl54uTSu1aBGwZ4/0nIKCpJXR9F7DYxhN\nkB0CL1tmlwDqO+enceHLzJS6qSxbJi3+Xbcu0Lu39Hy4lscwCuCAAOo356dR4cvMBHbvlhYF8vWV\nloKsWFEy8/ffgXffzSt8usirOImRfQPYP02SHQL36iUNclcAdXN+GhO++/elhoqNG6X7XqoU8M47\nUv61Zk1VTWMYzyP3miB2hsCOoF7YqxHhu3IF2LwZ2LQJ2LEDqFVLmlSga1cWPIbRBDZCYH3l/FQU\nPrMZ2LlTmj4qMlLql9emDfDWW9JC32XLus0UhmHsJSYG+N//pBA4X8uiYjk/IUQ7IcQpIcRZIcRI\nK2VmZB2PF0K8UugJ3Sh8RMDFi8CaNdIIiwYNgEqVgFmzgHLlpC4q168DK1dKOTxXhE+XeRU7MbJv\nAPunC5o1syh8rlCo+AkhHgMwE0A7ALUA9BBC1MxXpj2AqkRUDcCHAGZbPaGCwkckhbDbtgETJ0q1\n5IoVgfr1pWFm3t7SULObN6Uyo0ZJIzAef1ye68fFxclzIg1iZN8A9k83yNyXzFaDRyMACUR0Trq2\nCAfQCcDJXGU6AlgMAES0XwjhLYQoR0TXCpxNBuHLyAD+/luaESUhAThxQlrU5+hR6Xjt2lJXlG7d\ngKlTAT8/9/S/M5vNyl9EJYzsG8D+eSq2xM8HwMVc+38DaGxHmYoACopfIcKXmQncvg3cuiWFoleu\n5N2yBe/cOaBMGaBKFWnLbqCoXRt47jnuaMwwjH3YEj97W0PyS47F73UsugXmDiXx8KG0Otn9+5LY\n3boF3L0r9ZsrWVJapax8+X+2unUlgatSBXjhBW1OBXXu3Dm1TVAMI/sGsH+eSqGtvUIIfwAhRNQu\na38UgEwimpirzBwAUUQUnrV/CkCL/GGvEMI9zcoMw3gczrT22qr5HQJQTQjhB+AygO4AeuQrswHA\nxwDCs8TSbCnf54xxDMMwSlGo+BHRIyHExwC2AngMwAIiOimECMo6PpeINgsh2gshEgDcA9BXcasZ\nhmFcxG2dnBmGYbSE7BMbyN4pWkPY8k0I0TPLp6NCiN1CiH+pYaez2PPssso1FEI8EkJ0cad9rmLn\nu9lSCHFECHFcCBHlZhNdwo73s4wQYosQIi7Lvz4qmOkUQoiFQohrQohjhZRxTFeISLYNUmicAMAP\nwOMA4gDUzFemPYDNWf9vDGCfnDYotdnpWxMAJbP+304vvtnrX65yfwCIANBVbbtlfn7eAE4AqJi1\nX0Ztu2X2LwTAhGzfACQDKKK27Xb61xzAKwCOWTnusK7IXfPL6RRNRA8BZHeKzk2eTtEAvIUQ5WS2\nQwls+kZEe4noVtbufkj9HfWCPc8OAD4BsAbADXcaJwP2+PdfAL8Q0d8AQEQ33WyjK9jj3xUAJbL+\nXwJAMhE9cqONTkNE0QBMhRRxWFfkFj9LHZ597CijB5Gwx7fc9AewWVGL5MWmf0IIH0g/qOwhjHpK\nGNvz/KoBKC2EiBRCHBJCvOs261zHHv/CALwkhLgMIB7AEDfZ5g4c1hW55/OTtVO0xrDbRiHEawD6\nAQhQzhzZsce/6QC+ICISQggUfI5axh7/HgdQD0ArAE8C2CuE2EdEZxW1TB7s8W80gDgiaimEqAJg\nuxCiDhHdUdg2d+GQrsgtfpcAVMq1XwmSAhdWpmLWZ1rHHt+Q1cgRBqAdERVWTdca9vhXH1J/TkDK\nGb0phHhIRBvcY6JL2OPfRQA3iSgNQJoQYheAOgD0IH72+NcUQCgAENFfQogkAC9C6s+rdxzXFZmT\nkkUA/AUp6VoUths8/KGTRgE7ffOFlHT2V9teJfzLV/4nAF3Utlvm51cDwO+QGg+eBHAMQC21bZfR\nv2kAxmb9vxwkcSyttu0O+OgH+xo87NIVWWt+ZOBO0fb4BmAMgFIAZmfVjh4SUSO1bHYEO/3TLXa+\nm6eEEFsAHAWQCSCMiP5Uz2r7sfP5jQfwkxAiHlK+/3MiSlHNaAcQQqwA0AJAGSHERQBjIaUpnNYV\n7uTMMIxHov7qbQzDMCrA4scwjEfC4scwjEfC4scwjEfC4scwjEfC4scwjEfC4scwjEfC4scwjEfy\n/7MSuP4OUIeAAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x8414510>"
       ]
      }
     ],
     "prompt_number": 15
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The figure shows that the $\\mathbb{E}(\\xi|\\eta)$ is continuous over the entire domain. The code above solves for the conditional expectation using optimization assuming that $h$ is a quadratic function of $\\eta$, but we can also do it by using the inner product. Thus,\n",
      "\n",
      "$$ \\mathbb{E}\\left((2 x^2 - h(\\eta_1(x)) )\\eta_1(x)\\right)= \\int_0^{\\frac{1}{2}} (2 x^2 - h(\\eta_1(x)) )\\eta_1(x)  dx = 0$$\n",
      "\n",
      "where $\\eta_1 = 2x $ for $x\\in [0,\\frac{1}{2}]$. We can re-write this in terms of $\\eta_1$ as\n",
      "\n",
      "$$ \\int_0^1 \\left(\\frac{\\eta_1^2}{2}-h(\\eta_1)\\right)\\eta_1 d\\eta_1$$\n",
      "\n",
      "and the solution jumps right out as $h(\\eta_1)=\\frac{\\eta_1^2}{2}$. Note that $\\eta_1\\in[0,1]$. Doing the same thing for the other piece,\n",
      "\n",
      "$$ \\eta_2  = 2 - 2 x, \\hspace{1em} \\forall x\\in[\\frac{1}{2},1]$$ \n",
      "\n",
      "gives,\n",
      "\n",
      "$$ \\int_0^1 \\left(\\frac{(2-\\eta_2)^2}{2}-h(\\eta_2)\\right)\\eta_2 d\\eta_2$$\n",
      "\n",
      "and again, the optimal $h(\\eta_2)$ jumps right out as\n",
      "\n",
      "$$  h(\\eta_2) = \\frac{(2-\\eta_2)^2}{2} , \\hspace{1em} \\forall \\eta_2\\in[0,1]$$ \n",
      "\n",
      "and since $\\eta_2$ and $\\eta_2$ represent the same variable over the same domain we can just add these up to get the full solution:\n",
      "\n",
      "$$ h(\\eta) = \\frac{1}{2} \\left( 2 - 2 \\eta + \\eta^2\\right) $$\n",
      "\n",
      "and then back-substituting each piece for $x$ produces the same solution as `sympy`."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Example\n",
      "-----------\n",
      "\n",
      "This is Exercise 2.14\n",
      "\n",
      "![alt text](files/ex214.jpg)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x,a=S.symbols('x,a')\n",
      "half = S.Rational(1,2)\n",
      "\n",
      "xi = 2*x**2\n",
      "eta=S.Piecewise((2*x,    S.And(S.Ge(x,0),S.Lt(x,half))), \n",
      "                ((2*x-1),S.Ge(x,half)),\n",
      "               )\n",
      "\n",
      "v=S.var('b:3')\n",
      "\n",
      "h = (eta**np.arange(len(v))*v).sum()\n",
      "\n",
      "J=S.integrate((xi - h)**2 ,(x,0,1))\n",
      "sol=S.solve([J.diff(i) for i in v],v)\n",
      "hsol = h.subs(sol)\n",
      "\n",
      "print S.piecewise_fold(h.subs(sol))\n",
      "\n",
      "t = np.linspace(0,1,51,endpoint=False)\n",
      "\n",
      "fig,ax = subplots()\n",
      "fig.set_size_inches(5,5)\n",
      "\n",
      "ax.plot(t, 2*t**2,label=r'$\\xi=2 x^2$')\n",
      "ax.plot(t,[hsol.subs(x,i) for i in t],'-x',label=r'$\\mathbb{E}(\\xi|\\eta)$')\n",
      "ax.plot(t,map(S.lambdify(x,eta),t),label=r'$\\eta(x)$')\n",
      "ax.legend(loc=0)\n",
      "ax.grid()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Piecewise((2*x**2 + x + 1/4, And(x < 1/2, x >= 0)), (x + (2*x - 1)**2/2 - 1/4, x >= 1/2))\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAT8AAAE4CAYAAAAto/QTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXd4FFXbh++RnlAChqJUBVQQFOSVqlRfQRBDERAFwfah\nFBEEFVGxIdWGCNKrGDpiolSJAtKl9xI6UlKQkpD2fH9MNu8mpOzOzu7sbs59XXPB7s6eec7u5jfn\n+Z2miQgKhUKR27jD6gAUCoXCCpT4KRSKXIkSP4VCkStR4qdQKHIlSvwUCkWuRImfQqHIlWQrfpqm\nldc0bZ2mafs1TdunadqbWZw3TtO0o5qm7dY0rbZ7QlUoFArzyJvD64nAABHZpWlaYWCHpmmrReSg\n7QRN01oDVUSkqqZp9YCJQH33haxQKBSuk23LT0T+EZFdqf+/DhwE7s5w2jPArNRztgBBmqaVdkOs\nCoVCYRoOe36aplUCagNbMrxUFjhj9/gsUM7VwBQKhcKdOCR+qSnvIqB/agvwtlMyPFZz5hQKhVeT\nk+eHpmn5gMXAXBFZlskp54Dydo/LpT6XsRwliAqFwi2ISMYGWI7k1NurAdOAAyLyTRanLQdeTD2/\nPhArIhezCNBvj2HDhlkeg6qbql9urJ9Rcmr5NQK6AXs0TduZ+tz7QIVUMZskIr9qmtZa07RjwA3g\nJcPR+DAnT560OgS34c91A1W/3Eq24iciG3DAFxSRvqZFpFAoFB5AzfAwiZ49e1odgtvw57qBql9u\nRXMlZ3bqQpomnrqWQqHIPWiahpjd4eEJNE1ThwuHJ4iIiPDIdaxC1S93kuNQF0+gWoTG8JT4KRT+\niOVpb2qT1SMx+Bvqs1MofDjtVSgUCkcw+z6vxE+RI/7uGan6eT+HDkGTJuYKoBI/hULh9Xz8MbRu\nDWba3Mrz80KOHj3Kvn372LNnD23btuWRRx7J9Dz12SlyA3v2QMuWcOwYBAbe/rry/PyIsLAwypYt\ny8CBAxk7dqzV4SgUljJsGLz7bubC5wpK/DzAvHnz+PLLL+nSpQuhoaE5nj9gwADq1q3LmTNnuOee\nezwQYfb4g2eUHap+3sv27bBtG7z+uvlle8U4P2/n4MGDDBs2jPbt21O/fn2nBOnYsWNERUXx9ttv\nc+XKFapWrUq9evUcKmPp0qUMHTrUldAVCp/mww9h6FAoWND8slXLLwc2b97Mf//7X6pVq8a1a9co\nUKCAU+/fv38/o0ePBiA4OJgqVaqwY8eOHN+3fPly3nzzTc6du21pRI/TtGlTq0NwK6p+3smGDXov\n7yuvuKd81eGRA507d+aZZ56hW7du6Z4/ceIEU6ZMyfJ99evXJyQkhMTERA4fPkyNGjUQEcqXL09Y\nWBi1atXK8r1Lly7liy++ICgoiKZNm2bZ+vP2z06hMIoINGsGPXrASzkskme0w0OJXw4MGjSIAwcO\n8PLLL1OxYkUeffRRw2WFhYUxdepUli3TF8Revnw5efLkYf369dSsWZMVK1YwdOhQHnjgAYfK89Rn\nFxER4bOtB0dQ9fM+1q6FN96AAwcgbw7mnFHxU55fDlSvXp3k5GRSUlLIkyeP4XJiY2OZOXMmc+fO\nBeD06dNUr16dKlWq8NFHH/Hee+9RrFgxKlSoYFboCoVPIgIffACffJKz8LmC17f8zBrUaKSaCxYs\noFy5cjRs2PC21xxNe/VrC0OGDOG9994jKCiIU6dOUbFiRQAuXrxIly5dDPXIeXurWaEwQng4vPce\n7N4NdzjQK6HSXjfQrVs3pk2b5nQnR0bGjRtHo0aNKFu2LKdPnyYuLo7SpUsTHx/Pzp07iYyM5NNP\nPyUsLIynn37a4XK9+bNTKIyQkgJ16ui9vB06OPYeNcjZDXTv3p3HH3+cESNGcOrUKUNlbNiwgQED\nBvDoo49y991306BBA6pUqcKqVasIDw9HRIiPj2fp0qWUKlXK5BqYgy+PE3MEVT/vYeFCPdVt3979\n11KeXza0bNmSxx57jCVLltCjRw969uzp9JLgjz32GMnJybc9/+abb5oUpULhHyQl6S2+7783dw5v\nVqi01wF++eUXzpw5Q7Vq1WjWrJnV4aThC5+dQuEoU6fCvHl6T68z4qc8v1yI+uwU/kJ8PFStqqe9\n9es7917l+Snchi95RkZQ9bOeiRPhkUecFz5XUJ6fQqGwlGvXYORIWLPGs9dVaa8Poz47hT/w6adw\n5Aikjv93GuX55ULUZ6fwda5cgQcegC1boHJlY2Uoz0/hNnzBM3IFVT/rGDkSOnc2LnyuoDw/hUJh\nCWfPwvTpsG+fNddXaa8Poz47hS/z6qsQHKy3/lxBreqiUCh8hoMH4eef9Y4Oq1CenyJHvNkzMgNV\nP8/z/vvwzjtQvLh1MaiWn0Kh8CibNukbE82bZ20cyvMzgcjISJd2Wbtw4QLFihUjICDAqff5w2en\nyF2IQJMm0LMnvPyyOWWqoS4WceLECTZv3uxSGSVLlkzb5Eih8Gd+/RWiouDFF62ORImfy0yaNImu\nXbtme86oUaOyfT1v3ry0adOG2bNnmxmaaXijZ2Qmqn6eITkZhgyBL75w7/L0juKV4hd+JJzY+Nh0\nz8XGxxJ+JNyjZdh44YUXePzxx+nevXva0bdvX/bs2UO5cuVuO//gwYN07tyZn376icjISOLj43O8\nxqOPPsoaT09uVCg8yI8/QpEi8MwzVkeSioh45NAvdTuZPR8TFyO9w3pLTFxMpo8dwYwyRESioqJk\n/Pjx8sMPP9z22ueffy779u1L99ymTZukbNmy8tFHH8nkyZPl3Llz8vHHHzt0rYEDB8rRo0cdji2r\nz1Sh8Dbi4kQqVhRZv978slP/DpzXJCNvMnQhJ8RP5H9iFRkTaUi0zCpjwYIFcvjwYZk1a9Ztr4WE\nhEhKSkq65zp16iRz5sxJ95yj4jdr1iwJDQ11ODYlfgpf4auvRNq2dU/ZRsXPCzLvzAkqGMTgRoO5\n51u9F3XC9gmGy5qwfQKR/SMJKhjk9HtPnDjBfffdx9y5c1m9ejUA7dq1o2PHjty8eRMtw5KzFSpU\nYN68eRQsWPC2fX7//vtvNm/ezLlz56hfvz6JiYksWrSIeal9/sWLF+eIlaM+s8AX9311BlU/9xIb\nCyNGwLp1loWQKV4rfrHxsYzZOIbI/pGM2TiG4S2GOy1esfGxDF07lMGNBhsuI2/evGzZsoVVq1bd\n9lpme3Nkt8/vpUuXqFGjBitXrmT48OGkpKTw1ltvpb1eqFAhEhISnIpPofB2RoyAkBB48EGrI8mA\nkeaikQMf9PxOnDghy5cvl++++y7T15988sl0j+fPny8bN2687Tz7tPf999+XBQsWiIjuD3bo0CHt\ntUWLFmXqLWZFVp+pQuEtnDolUqKEyLlz7rsGBtNer+zt3Xh6Y7pWWlDBIIa3GM7G0xs9Wsb69eup\nW7cuhQoVyvT1MmXKcP369bTHy5cvp06dOtmWuXbtWlq0aAHAihUr0u3Te+HCBapUqeJwfAqFt/PB\nB9CnD9x9t9WR3I5Xil+b+9rclp4GFQyizX1tPFrGxo0beemll1i4cGG6YS7rUs2LJk2asHXr1rTz\nc9rnNy4ujiJFilCiRAkAVq9eTevWrdNe37VrF40aNXI4Pk/hLePE3IWqn3vYuRNWrYLBgy25fI54\nrefnDUyaNCnb1zt06MDYsWNp3rw5kPM+v4UKFUrrNImOjiYpKYnSpUsDEB8fT9GiRSlYsKB7KqNQ\neBARfeGCjz7Sx/Z5I17Z8vMVgoKCCA4O5sqVK2nPBQYGEhQUROfOnalYsWKW7928eTPt7balDw0N\npVevXm6N1yj+3BMKqn7uYNUqOH0aXnvN45d2GNXyc5H+/fszdepUXrP7ltu2bZvj+1q3bp2W8p45\nc4bixYtz//33uy1OhcJTJCfrrb6RIyFfPqujyRrV8nMRTdPSCV9m5HTnLV++PCEhISZGZS7KE/Nt\nPF2/OXP0VLddO49e1mlUy88DNGnSxOoQFAqPcPMmfPghLFgAmtOLTHkWtZ6fD6M+O4W3MXw47NoF\nCxd67ppq395ciPrsFN7EP/9AjRqu7cFrBLWYqcJtKE/Mt/FU/T76SF+h2Yo9eI2gPD+FQuEye/fq\nu7EdPmx1JI6j0l4fRn12Cm9ABFq2hLZtoV8/z19fpb0KhcISVqyAU6fg9detjsQ5lPgpckR5Yr6N\nO+uXlASDBsGYMd49oDkzlPgpFArDTJ0KpUvrKa+voTw/H0Z9dgor+fdfuO8++O03qF3bujiU52cB\nkZGROZ5z4cIFbt686YFoFArPMmIEtGplrfC5ghI/gzi6Wbk/bEiuPDHfxh31O3ECpkzR9+D1VZT4\nGcSRzcrB+zckVyiMMHgwDBjgnSs0O4oSPwPs3r07083Ks8LXNyRX6935NmbXLyICduyAgQNNLdbj\n5Ch+mqZN1zTtoqZpe7N4vammaVc1TduZenxgfpjeRVhYWNrqzY5SsmRJjh075qaIFArPkJwMb70F\no0dDFlvb+AyOtPxmAK1yOOcPEamdenxuQlxezbZt26hevbpT73n44YfZsWOHmyJyL8oT823MrN/0\n6VC0KHTqZFqRlpHj3F4RWa9pWqUcTnPfyl1mLQpmYEjIzp072bRpU7pNxhcuXJjpZuXLly8nT548\nrF+/npo1a7JixQqGDh3KAw88AHjvhuQKhaNcvaqv1Rce7v1r9TmCGZ6fAA01TdutadqvmqY51yTK\nsXQx5zDAxYsXqVGjBvv27aNt27a0a9eODRs2kJKSku6806dPU716ddq0acPq1atp06YNXbp0oUKF\nCmnn+PKG5MoT823Mqt/w4dC6NeSwO6vPYMaqLn8D5UXkpqZpTwHLgPtMKNdyWrVqxdChQ+nWrRsA\nW7dupW7duty4cSPdeTaRu3jxIkWKFCEoKCjdfrwAV69eTduyUqHwNY4d01PeffusjsQ8XBY/Eblm\n9//fNE2boGlaCRGJznhuz549qVSpEqDvfFarVi1XL+921q5dy9tvvw3AypUradOmDX/88QfXr1+n\ncOHCABw6dIhbt27x999/07hxY0DvFMm4IXm1atXcFqfN17Hd5c18bO8ZuaN8qx+r+uX8uGfPCDp0\ngDJlvKM+M2fOBEjTEyM4NL0t1fP7RURqZvJaaeCSiIimaXWBBSJyW0S+OL0tLi6OZ555Jm2v3UaN\nGrFkyRLCw8OpVKlSWo/vuHHjuHbtGnfddReHDh2iQYMGlC1blrp166aV9eqrrzJ+/HhT9+X11GcX\nERHh16mhql/2rF4NvXrBgQPgjdtKG53ehohkewA/AeeBBOAM8DLQC+iV+nofYB+wC/gLqJ9FOZIZ\nWT3vbURFRUndunVFRCQmJkaGDh3q8Hvj4uJkwIABpsfkK5+dwndJSBCpVk1k2TKrI8ma1L+DHLUs\n4+FIb2+20xhE5Hvge6dV18ew32TcfrPy4ODgHN/rzRuSKxTZ8f33UL48PPOM1ZGYj5rh4SCtW7fm\nvffeS3vcv39/li5dmuP7/GFDcjUOzrcxWr9Ll/Qe3m++8Y+hLRlRe3gYxJHNykHfkLx8+fIeiEih\nMJchQ+DFF8GN/XSWotbz82HUZ6dwF1u3QkgIHDoExYpZHU32qPX8FAqFKaSkwJtv6stVebvwuYIS\nP0WOKE/Mt3G2fnPm6JOievRwTzzegvL8FApFGv/+q3t9y5bBHX7eNFKenw+jPjuF2QwaBFFRMGOG\n1ZE4jlHPzytafhlXSFEoFJ5n/36YNUv/NzdgecPWyMhsbzzWrVtnyXU9gfLEfBtH6icCffvCsGFQ\nqpT7Y/IGLBc/hUJhPaGhEBsLr79udSSew3LPT6FQWMu1a/pA5gULoGFDq6NxHqOenxI/hSKX44ud\nHPaoQc4W48++kT/XDXJ3/fbvh9mzYdQoz8XjLSjxUyhyKbmxk8MelfYqFLmU0FC9xbd9O+TJY3U0\nxlGen0KhcJh//4Xq1X23k8Me5flZjD/7Rv5cN8id9fvoI2jVyveFzxW8YoaHQqHwHDt36ilvbpnJ\nkRUq7VUochHJyXprr1cvePllq6MxB5X2KhSKHJkyBfLnh549rY7EepT4mYQ/+0b+XDfIPfW7eFH3\n+iZO9P/lqhxBfQQKRS5h8GC9xVejhtWReAfK81MocgHr1unCt38/FC5sdTTmojw/hUKRKQkJ0Ls3\njBvnf8LnCkr8TMKffSN/rhv4f/3eeCOCqlX13dgU/0ON81Mo/JgjR2DRIti71+pIvA/l+SkUfooI\ntGgBbdvCgAFWR+M+lOenUCjSMXs2XL0K/fpZHYl3osTPJPzZN/LnuoF/1u/yZXj3XX1Q84YNEVaH\n45Uo8VMo/JC334YXXoBHHrE6Eu9FeX4KhZ+xZg28+irs25c7hrYoz0+hUBAXp+/A9v33uUP4XEGJ\nn0n4o29kw511Cz8STmx8bLrnYuNjCT8S7rZrZsSfvrvPP9dT3TZt/vecP9XPTJT4KSylUYVGDF07\nNE0AY+NjGbp2KI0qNLI4Mt9jzx69g+Pbb62OxDdQnp/Ccg5dPkSjGY1Y12Mdk7ZPYniL4QQVDLI6\nLJ8iKQkaNNBT3ldesToaz2LU81MzPBSWIiIMWj2I6LhoHv7hYSL7RyrhM8A330DRov6zQKknUGmv\nSfizr+LOuk3bOY1TV09RrEAxVrywgjEbx9zmAbobX//ujh+HkSNh8mTQMmn/+HL93OkJK/FTWEZk\nTCTvrH6HB0s+yL3F76V4oeIMbzE8nQeoyB4ReO01GDIEKle2OhrzsfeE4xLjTPWEleensIQUSaH5\nrObcU/wevm75Ne1C2zGsyTCa3dOM2PhYNp7eSJv72uRcUC5n6lSYNAk2bYK8fmpixcbH0nNZT45F\nH6NxxcZ80eKLdNaI8vwUXk34kXAaVWiU9qP9dvO3xCfF0/6B9gQVDCIwfyA3E28CEFQwSAmfA5w/\nr7f41q71X+ED0NDYfn47566dI+z5MNM8YZX2moQv+yo5YUbd7NOXA5cP8Nmfn1GlRBUaV2wMQEC+\nAG4k3nD5Okbwxe9OBPr00Xdhe+ih7M/1xfrZ80b4GwQVDCKyf6SpnrASP4VHCCoYxPAWwxmyZghd\nFnahZumajG89Pu0uHpAvIK3lp8iZxYvh4EH44AOrI3Evs3fN5tejv7Ky20oqBVUy1RNW4mcSTZs2\ntToEt2FW3YIKBlEgbwH2Xd7HzJCZ6dKXwHyB3EiwpuXna9/dlSv6MlXTp0PBgjmf72v1s3Hh2gX6\nr+jPwk4LKVu0LPC/m+jG0xtdLl+Jn8JjrDq2isk7JrP5lc2M/Wtsuru3avk5Tr9+8Pzz+ubj/oT9\nsBYR4dVfXuW1Oq+RkJyQ7jyzPGElfibh675KdphRtzNXz9BpUScmPT2JeuXq3Za+BOYLtEz8fOm7\nW7oUduyAzz5z/D2+Uj97X3jyjsmc+/cc125dc9tURyV+Co/wyvJXaFu1Ld0f7g7cnr5Y2eHhK0RF\n6Z0c06dDQIDV0ZiP7TfR99e+DFk7hOolqzPiiRFum/Gjxvkp3M7Sg0sZtHoQu3rtokiBIpmeM37r\neA5dOcT41uM9HJ3v0L073HmnPpXNX0lMTqTulLrsuriLyP6RVAqqlON71Hp+Cq/B3ru5cO0Cb4S/\nwcTWE/nz1J9Zvke1/LLnl1/0gczDh1sdiXt5f+37RMdHc+LNE26f6qjEzyR8xVcxgrN1s3k3MXEx\nvPTzS/Ss1ZOfD/+crXejPL+siYmBN96AadMgMND593t7/Wz8euRXJmyfwKpuq7in+D1un+qoxE9h\nOjbvJiQ0hIvXL3I1/mqOy1QF5AuwbKiLtzNgALRvD02aWB2J+4iNj6Xnzz2Z9sw07g++HzB3WEtm\nKM9P4Rb2XdpH4xmNiYmPcci7+T3ydz7/83N+7/G7ZwL0EZYv18Vv927/XZZeRHh+yfPcWehOQ56v\n8vwUXkNcYhydFnSidpnaDk9JUp7f7Vy5oi9OOnOmfwlfxmWq5u6Zy84LO2l+T3OPxqHEzyR8xVcx\ngrN16/trX5IkiUWdFzk8JcnKQc7e+N2JQO/e+mDmxx93rSxvq5/9eL4TMSd4a+Vb1CpTy+Pi58dr\nQSisYNmhZYQdDWPrq1spXqg4kN67yWpkvpUdHt7I/Pn61pOzZ1sdifnYz/PefG4z1YKr8cPTP3h8\nBW/l+SlM4+y/Z6kzuQ7LuiyjQfkGTr33wrUL1J5Um38G/eOm6HyHCxegVi0IC4NHH7U6GvfRJ7wP\nE7ZP4Pibx7m3+L2Gy1Gen8Lj2Hs3ySnJdFvSjddqv0Z0XLTTZdmv55ebsa3M3KuXfwvfz4d+Zubu\nmWx7bRtf/vWlJSt3K/EzCW/zVcwkq7rZezcjN4wkKSWJ6PhoQ3MxbR0eVmQH3vTdzZihL1Jq5lJV\n3lQ/gKNRR3l+yfP82OFH/nP3fyzbukCJn8IwNu/mlZ9f4evNX1O5ROXblhh3lLx35CXvHXlvW8Ej\nN3HqFLz7LsyaBfnzWx2Ne0iRFJ5b/Byv13mddg+0A9w/ni8rlOencInouGhqTqzJ+WvnHZ6LmRXF\nRxXn+JvHKVGohHkB+gjJydC8ObRurQugvzL2r7EsPbSUiB4R5MuTz5Qyleen8DgiQrcl3QguFGzK\nEuO5ucf3q690v2/QIKsjMRd7X3jrua2M3jiaiW0msur4KosjU+JnGt7mq5hJVnUbuWEk289vZ/WL\nq01ZYtyqKW5Wf3e7d8Po0fqwljx5zC/fyvrZfOFTsafourgrY/87lknbJ7ltjT5nyFH8NE2brmna\nRU3T9mZzzjhN045qmrZb07Ta5oao8Ea2ndvGyI0jWdV9FaUCSwGueze5cTXn+Hh44QX48kuoVMnq\naMwnqGAQnzf/nBazW1C/bH22nNuS4zxvT+FIy28G0CqrFzVNaw1UEZGqwP8BE02Kzafw1X0SHKFp\n06bp0pfY+Fi6LOrCuFbjOPfvuXTnurLEuFXDXaz87t5/H6pV09fqcxdW/zbn7Z1H/jz5mbdvHoMb\nDfYK4QMHxE9E1gMx2ZzyDDAr9dwtQJCmaaXNCU/hLdgvU/XaL6/RvFJztp7bamr6ktvm965dCwsW\nwA8/gOa0Xe8bbD+/nWERw3jkrkdM33rSVczw/MoCZ+wenwXKmVCuT2G1b+ROIiIi0lLadqHt2H9p\nP3nz5DU9fbGqw8OK7y4mBl56SV+j78473Xstq36bsfGxdJzfkf/c/R/Gtx5v+taTrmLW3N6M961M\nx7T07NmTSqnGRlBQELVq1Uprktu+IF99vGvXLq+Kxx2PD185zN5Le4mJj+Gjih+xa/MuU8u/dvga\nN6rf8Jr6uuuxCHTsGEGdOtCypfXxuOPxunXr+HDdhzxU7SHmtJ/Drs3/+/sY3mI4ExdOpEH5BobK\nj4iIYObMmQBpemIIEcnxACoBe7N47QfgObvHh4DSmZwnCt8l6maUlP+qvLSc01IiYyKld1hviYmL\nMfUar/78qkzePtnUMr2RWbNEqlcXuXHD6kjMI+xwWLrfw1d/fSW1f6gtSw8sdfu1U7XFIS2zP8xI\ne5cDLwJomlYfiBWRiyaUq/ASUiSF5xc/T4lCJQh9NtT89EUEli0jME8hv/f8jh2Dt9+Gn37yrx3Y\n7Kc6bj67mS82fMGDJR+k6T1NrQ4tSxwZ6vIT8Bdwv6ZpZzRNe1nTtF6apvUCEJFfgROaph0DJgG9\n3Rqxl2KVr+IJen3Xi1NXT7G6++o0j8+0KUki0LcvtG9PmWspfu35JSbq6/N9+CE89JBHLgl4pn62\n38PAlQPpuKAjj979KN+1/s5renYzI0fPT0S6OnBOX3PCUVhN+JFwGlVolPajXRe5jkUHF/FNr28o\nGVgy3bmuDGsB/id8f/8NZctSNDkv5/x4H49hwyA4GPr1szoS91AkfxGORR/j/LXzbHx5o1cLH6gZ\nHqZhM2Z9Hfv05fy183Rd3JXHGz9OyAMh5l7IXvhWrIDgYIom5vHbcX4REfpy9DNnen5Yi6d+m++t\neY/jMcc52u+o+UNa9u2D557TfzcmocRPkQ77VXbbhbajQrEKzG4/29y7eEbhK1YMAgIonHyHX3p+\n0dHw4oswfTqUKmV1NO4hdF8oE7dPZN2L66hSooq5nvC+ffDf/0JIiKl3DiV+JuFPnl9QwSASkxPZ\ndn4boc+Gpg1TMIXMhA8gMJDCCZrfeX62xUmffRZaZTlPyr24+7d5PPo4vcJ6sbjzYu4Lvg8w0RO2\nCd9XX0HXHB04p1Dip7iNqTumsuDAAnb22smXf33J9YTr5hSclfABBAQQmGSN+LmTSZPgxAkYMcLq\nSMzDfqrjzcSbdFzQkfcfe5+klKR057nsCbtR+ECJn2n4i+e34fQG+q3oR/jz4dQqU4vhLYbzW+Jv\nrqcv2QkfQGAgAYlYkva667vbs0fv2Z0/HwoUcMslHMLs+tlPdewd3psqJapw6uopc1dqcbPwgRI/\nhR0xcTF0WtiJb1t+y+MV9f0STUlfchI+gIAAAhLEb1p+N25Aly763+5991kdjbmkTXWc345NZzdR\nolAJwyt4Z4oHhA+U+JmGL3p+9ulLiqTwwpIXaHd/O8oWLZvuvF2bdxlPXxwRPoDAQAompPjNen59\n+0K9eu5drcVR3FG//Zf2s+/iPo5EHeH9x9/3OeEDJX65GvthLR9HfMzVW1fTnjcFR4UPICCAgreS\n/aLlN2cObNoE48dbHYl7OPfvOZ5d+CwNyjcwd6UWDwof4NjcXjMO1NxeryQmLkaemvuUlBlbRnou\n7WnefN2UFJHevUXq1xeJjc35/M8+k9hBfaXcV+XMub5FHDokEhwssnu31ZG4h7jEOHlk0iNSb0q9\ntN9KTFyM63O99+4VKVNGZN48p9+KhXN7FT7MhWsX2HJ2C/9c/4dhTYeZk7440+KzERBA/rhEn275\nxcfrPt+nn3p2+pqnEBH6hPehUN5C/PbCb+ZNdfR0iy8VJX4m4YueX0xcDE/Pe5paZWplm744VTcj\nwgcQEEC+W4k+7fkNGgRVqsDrr5tSnGm4Uj97X3jCtglpYz//OvNXuvMMD2uxSPhAiV+uJSkliY4L\nOlI4f2G3VlGMAAAgAElEQVQWd1lszkotRoUPIDCQPHG3SEhOIDkl2dj1LWT+fL3KU6f616rMNl84\n/Eg4n/75KbPazWLE+hHm+MIWCh+gPL/cysAVA6XWD7Xk8o3L6Z6PiYuRsMNhzhforMeXkUWLRNq3\nl8DhgXLt1jXn328hNp/v77+tjsQ97P5ntwQMD5BZu2aZt46jCx5fRlCenyIr7FMXgNm7Z7P00FLe\nbfQuwQHB6c41lL640uKzERgIN29atn2lUW7e1KeuDR8Otf1w38LrCdfpvrQ7A+sPpMeyHuZsQGR1\niy8VJX4m4c2en/2Qli1ntzBw5UDqlatHqyqOTTbNtm5mCB/oK3veuGHJDm5GvzsR6N0batXS5+96\nK0brlyIpdF/anYdKPURUXJQ5w1q8RPhAiV+uwNYb99aKtwgJDaFe2XpMbDPR9Tu4WcIH6Vt+PrKy\ny/TpsG2b/+6+9tG6j7h4/SKB+QP5osUXrvvCXiR8gPL8cgvXb12XB79/UPgYiYyJdL1AVz2+jBw4\nIHL//fLo5Edly9ktrpfnZnbt0n2+/futjsQcMu7B8eOeH6X8V+Vl0MpBt3l8hnxhEz2+jKA8P0VW\npEgKzy1+jqSUJE68ecL11MXMFp8NH/L8rl7Vfb5vvoHq1a2OxhzsrZGt57bS77d+NCzfkKGNh96W\nITjtC3tbi8+GEcU0cuDnLb9169ZZHUKWDFgxQO4ae5f8c+0fEXF+RH66upnd4rNx+bJIiRLS+sfW\nxnqbXcCZ7y45WSQkROSNN9wXj9k4Wr+YuBh5cemLUnpMaXlq7lNe16ubFaiWnwJu79mdtWsWc/fM\n5ZOmn1C6cGnAhRH57mjx2fARz2/UKLh4UW/1+Rt578jLtnPbuHjjIhPaTPCbXt0sMaKYRg78vOXn\nLdi36v48+afcOepO6bqoq+t3cXe1+OzL1zTpuai7zNg5w/zyTWDVKpG77hI5e9bqSMwnKTlJWs5p\nKdXGV5MT0ScsnavrLBhs+Snx80Ni4mLkhcUvSPCoYGk7r633C5+NgAB5a+Gr8v3W7913DYOcPClS\nurSIF7sbLtFreS8p+2VZuXT9koi4uFiBB4VPRKW9luNN4/xSJIVNZzdxJe4K454a51r6IkJE+/bu\nSXUzEhhIsZR8Hu/wyOm7i4+Hjh1h8GDwxQW7c6rfd1u+I/xYOBtf3pi2Palha8TbU107lPj5GbeS\nbtH2p7YUyV/E9UGpNo/vyBH3Cx9AQADFkvJ63coufftC5cowcKDVkZiDvS8cdiSMERtGENY1jH2X\n9qU7z296dbPCSHPRyIFKe91OckqydJzfUSp/W1mibkaJiAvpi6dSXXuqV5epswfI4FWDPXM9B5g8\nWaR6dZFrvjXdOFtsv4k/Iv+Q4NHBsurYKp/y+DKCSnsVH/z+Afsv7Wf9S+spUagEYDB9cWevbnYE\nBFAk6Q6vafn99RcMHQpLl0LhwlZHYx5BBYN44z9v0Hpeaz5p+gnLDi1jeIvhxu0RX2vxpaLEzyQ8\n7fllHNIyecdkQveF8mGTD7mryF3pznUqfclE+DxWt8BACid5fuPyzOp39ix06gQzZ/r+BkQZ6xcb\nH8tzi5+jf73+9Pm1j2uLFfio8IESP5/FfkT+b0d/48N1H9KwfENaV21tvFCrWnw2AgIonGj93r3x\n8dChA/TrB61d+Di9kVtJt2gX2o7HKzxOTHyMa76wDwsfgKanzB64kKaJp66VW4iNj+X/fvk/1p5Y\nS9NKTZkWMs34Hdxq4QPo1Intj93LxyX3E/Z8mOevj/4x9OwJt27BTz/514IFKZJC18VdiU+Kp2yR\nsmnbTcbGxzJ07VDnUl8vEj5N0xARp78p1fLzYaJuRhFxMoLo+Gi+bPmlbwsfeMXevd9+C7t3w7Rp\nvi98Ga2RwasGc+bqGWqUrJFun12nfWEvEj5XUOJnEp72/C5ev8gTs5+gesnqrqUuDgifx+oWEECh\nBLHM81uzBkaOhGXL9Nl2vo7NGglbFcbXm74m7EgY1UpWy9Tjc9gX9hPhAyV+Psm1W9d4cu6TBAcG\ns+y5ZcbXWfOWFp+N1I3LrWj5HT8OL7wAoaFQqZLHL+8WbC26TyI+YdTGUdQrV48vn3QhQ/Aj4QMl\nfqbR1I1D/+3Tl4TkBDos6EDJQiV5u/7bxlMXJ4TPnXVLR0AABW4leVz8atduytNPw7BhvjmDIzu2\nntvKiWInuHjjIp82+1QJnx1K/HwAW/oSHRdNz2U9yZ8nP/cF30erqumXoXc4dfG2Fp+NwEDy30ry\n6PS2pCR9r90nntCXpPcntpzdQtfFXWlcqXGu7tXNCiV+JuFOXyyoYBCfN/+cZjObcSz6GOWKlEtn\nWDuFAeHzpOeXP96zG5e//TZERUXw9dceu6RHOHD5AG1/akv9cvV5pfgrxq0RPxU+UOLnM4zbMo7E\nlES2nd/GkMeHeEz4PEpgIHnjE7iReANPDIv64QdYuVJPd/Pmdfvl3Iq9NXL66mlazW1Fp+qdePGh\nFymcX5+eklt7dbPEyJw4Iwdqbq9hvtn0jdz7zb3Sc2lPiYyJ9J25us7y008inTtL/s/yS3xivFsv\ntXatSKlSIkePuvUyHsM2X/folaNy/3f3yxd/fuHafF0L5+o6C2o9P/9kxs4ZUvbLstJ9Sfe0H7LT\nixX4gvCJiCxfLvL00xI0MihtYQZ3cOSILny//+62S1jCqdhTUnJ0SekT3ifXCJ+IEj/LMWMPj4w7\naC0+sFhKjSmV6Q/Z4R20TBA+j+1PsmaNSLNmUvbLsnLm6hm3XOLyZZGqVfXVWmx48/4rjnL91nV5\nbPpj0m1xt9t26HOqfj4mfCLGxU95fl6E/Xzd1cdX83+//B9NKjbJdNqRQz273u7xZcTNO7jFx0O7\ndvq8XW/eZNxZbiXdov389pQvWp4iBVxYx9HfPb6MGFFMIwd+3vIzi5i4GGkf2l6Kjywu7ULbGU9d\nfCXVtWf3bpGaNaXWD7Xk7/N/m1p0crLIc8+JdO6s/99XyZgdJCQlSOu5raX+lPryetjrxq0RH2zx\n2UC1/PyDw1cO88fJP4iJj+Hrll/7Z69uVgQEwI0bBOQLMH24y4cfwqlT+hJVd/jwr94+O0hOSabr\n4q7sv7yfvnX7MqLFCGOD3nNbi8+GEcU0cuDnLT8zfKMd53dI8OhgafNjG6/q1fWYJ3b+vEiZMvLE\n7Cdk5bGVphU7ZYpI5coily5l/rqveX4xcTHyRtgb0nlBZyn7ZVk5/+/5bM/Ptn4+3OKzgWr5+RYZ\nV9zYe3EvT855kgrFKjC3w1xjg1J9tcVnI7XlF5gv0LSW3+rV8MEH8OuvULKkKUVaTrECxbiRcIMF\nBxawstvK2xavdZjc2uKzYUQxjRz4ecvPWew9mQOXDkjpMaWl9g+15WTMydvO81SvruUkJIjkySNd\nFz4nc3fPdbm4XbtESpYU+fNPE2LzElJSUqTX8l5Sakwp2f3Pbp/ZXtKdoIa6+B4xcTHy/OLnpfSY\n0tJiVovc1bmRFfnyyeuLXpLJ2yfnfG42nDwpUq6cyPz5JsXlBaSkpMjrv7wupUaXSrtJGtqgyo+E\nT0SlvZZjZP7r5RuX+T3ydy7euMjUZ6Z6beeGR9cqDAykWHI+l9b0i4qCVq1g0CDo3Dnn871pz2V7\n7K0REWHQqkGsPL6SEU+MoGJQRcCxjo109cvtqa4dSvws4vCVwzSZ2YRqwdWMj8vydY8vMwICCErJ\nZ9jzi4uDZ56Btm2hf3+TY/Mwtp7dmLgY3ln9Dmsi19D8nuZ0qNYh3Xm5cSFSM1B7eFjAoSuHaDar\nGdWCq7GkyxJj+yj4o/AB3HcfE4a15lyZQIa3GO7UW5OS4Nln9W0mZ8/27SEtNmLiYmgxuwXxSfE0\nKN/A+GKkfix8ag8PLyVjr+7BywdpNqsZ9e6ulyZ84OS4LH8VPoCAAIom5XW65Wf7SG7ehOnT/UP4\nRITP/vyMW0m3OHjlIB82/lAJn4n4wU/EO8jKN7IflLr/0n6azWrGgyUfZGb7mT4zZc3Tnl/hJOe3\nr/z0U9i2DRYvhvz5nbukN3p+KZJCn1/78MepP2hQvoFLU9YiGjdWwpcJSvzcjK1F1+uXXjSb1Yya\npWqyqPMir+3csJyAAKc3Lv/uO5g7Vx/LV6SIG2NzExmzg+SUZHos7cHaE2t55K5HGPvkWGPjPm0t\nvj59lPBlhpEuYiMHuXioy+Yzm+XOUXfettqGU/jTcJbseOYZ+WvcO9I+tL1Dp8+dqw9piYx0b1ju\nxH64SmJyonRe2Fnu/vJumblzpvHVfPxsOEt2oIa6eCcRJyNoM68N9cvVV726jhAYSEAiDrX8wsL0\nZehXrvTtHdds2cGQNUNoF9qOTWc2se21bfSo1cOYNaI8PodQ4mcSERERt6Uvvx39jQ7zO1ClRBWf\nnrLmUU/MwY3L//wTXn4Zli+H6tVdu6Q3eH4F8hTgcNRhwo+Gs7r7au4ucrexgjIRPm+onzeixM9E\n7Ds3Fh9YzIvLXuSB4AeY/+x81avrKKl792a3nt/OnfqQlp9+grp1PRibm7gaf5Un5jzB2X/PcqTv\nEcZtGad2WfMERnJlIwe5xPOLiYuR5jObS/DoYOm8oLOasuYs770nl95/S+7/7v5MXz54UOSuu0QW\nL/ZwXCaRcT2+S9cvyYPfPygVvq6QtnS/mrLmHCjPzzuYsmMKh6IOceXmFUb9d5Tq1XWWwEDyJyRl\n6vkdP643bEaO1Fdj9kXss4MzV8/QcFpDbiXd4o8ef1CiUAlA7bLmKZT4mcS6desYvGow03ZO44l7\nn/Crzg1Pe37545Nu8/xOn4YWLfTlqV580dxLerJ+NmHr82sfGkxrQKnCpdj2f9uoVLzSbeeZNWVN\neX6Zo8TPIPadG4nJiYzaOIq1kWupUKwC37b61mc7NywnMJB88QnpxO/8eV343noLevWyMDaTOBJ1\nhFXHV3Hu2jl+7PCjsewAVIvPVXLKi4FWwCHgKPBuJq83Ba4CO1OPD7IoxwPZv+ew+TLn/j0nbX5s\nI0/MekJazmmZu9fjM4PZsyWlWze545M7JCk5SS5eFKlWTeSLL6wOzBzCDofJnaPulNY/tja+WrdI\nrvb4MoI71vMD8gDHgEpAPmAXUE1uF7/lOV7Iz8RPROR41HEpM7aMtA9tL6//8rrq3DCDRYtE2reX\nwOGBcurCNXn4YZEPPrA6KOfJ2LEhIjJu8zgpOqKodJzf0fhGQyJK+DJgVPxySnvrAsdE5KSIJAKh\nQEgm5zm9ooKvczz6OE/Ne4qQ+0JYemgpTWjit50bnp7by82bBOQNpG2HGzzxhD5v1524o372HRsi\nwvtr3mfo70P5pOkn6dZu9ETnhvL8Micn8SsLnLF7fDb1OXsEaKhp2m5N037VNM3FIafez5azW3hs\nxmP0qtOLPHfkIbJ/JKF7Q/2ic8NyAgJIvHqDq1EB/KfBTcaMAc0Hb61pszbWDqHbkm5M2TmFba9t\n4636bxmbtQHK4zOb7JqFQEdgit3jbsB3Gc4pAgSk/v8p4EgWZbm58eseMqYvSw4skTtH3Snvrno3\nXbridPqiUt1Mufr7djkQ8IiU+LC67Plnr9XhuMTV+KvSZEYT4WNkr6t1UalulmAw7c2bgzaeA8rb\nPS6P3vqzF89rdv//TdO0CZqmlRCR6IyF9ezZk0qpkzCDgoKoVasWTZs2Bf7XNPe2x43q6+nLU/me\nYsWxFSyOW0zze5qT/2x+nir7VNpdfNfmXTyV7yk2nt5Im/vaZF++CBHt28ORIzTdtAmKFfOa+lr5\n+No1mPVeab4vcIOSN2HDn+up2amG18TnzOP5YfN5d/W7BN4fyJG+Rxg8eTCv1nmVp5982vnybMtS\n9elD0wxT1rylvp58HBERwcyZMwHS9MQQ2SkjkBc4jt7hkZ/MOzxK878VoesCJ7Moy+13AHdx+cZl\nqTmhplT+trJ0X9I909adw3u/+mCLzxP72kZHi9SpI/LJK6ckpXx5aTKjifx+4ne3X1fE9fplzA62\nnt0qpcaUkmrjq0n0zWgRMdixIWJKi8/X9iV2FtzR4SEiSUBfYCVwAJgvIgc1TeulaZptxNWzwF5N\n03YB3wDPGZdi7yM2PpYXlrxAUMEgjscc59Nmnxofl6U8vkyJitKtrMaN4cMRAWg3bxKY37y9e92N\nfefGkoNLaPVjK8oXLc9vL/xG8ULFAQMdG6A8PndjRDGNHPhAyy/jHfzIlSNS5dsq0mpOK3k97HXX\nxmX5YIvPE/zzj0jNmiKDBukfkdy8KVKwoDy74FmZv8939p2Mvhkt9afUl9JjSkunBZ2MD3uyoTw+\nh0HN7XUd+zv475G/03B6Q0oXLk3ZomUZ0WKEsVkboFp8WXD2LDRpAh07wujRqb26BQtCQgKF8wT4\nTMsvLjGOvr/15VrCNS7euMjo/442nh2AavF5CCV+dthSk5DQEDov7Eyj8o3oV68fY58cm+O4rCzH\nUvmB8LljnNjJk7rwvfwyDBtmN5xF06BQIYJS8ntM/BytX8b1GgH2XdxHzYk1SUhOoFH5RsbndKcV\naL7wqXF+maPEz45bSbd4Z/U7nP33LFFxUXzT6hu6PNjF+LgsPxA+d3DkiO7vDRgA77yTyQkBARRL\nyZftmn5WYJ8ZAKw4uoJ60+rRqXonSgaUZNR/RxnPDkC1+DyNkVzZyIGXe35nr56VelPqydPznpZX\nl7/qmr8nojy+LNi7V+Tuu0WmTs3mpIoV5asf+8nH6z72WFyOYuu1HbF+hBT8vKAs2Lcg06lsDs/p\ntqE8PsOgPD/HyZi+rD+1njqT63BP0D2UL1qeMf8d49odXLX4MmXTJn11lrFj4ZVXsjkxMJCiyXmd\n2sHNUxTMW5DouGiGrB1CWNcwOj3YiTb3tTGeHYBq8VlErhQ/W/oSExfDd1u+o8OCDjx696O0q9aO\nL1p8YWjeZZqv4ofCZ4ZntGIFPPMMzJzpwN93QABFku7wOs/vZOxJ6k+tz7bz29jz+h6WHFxi3Nuz\n4QHhU55f5uRK8QsqGMSQx4dQb2o9JmyfQKvKrZjTYY5r/h74pfCZwU8/QY8e8PPP8NRTDrwhMNDp\nvXvNJrPNqOpMrkNCcgLbXttGzdI1jWcGNlSLz1qM5MpGDrzI89t7ca88MP4BeXb+s67tpWuP8vgy\nZfx4kbJlRfbsceJNTz0lf3z/rnRe2NltceWEzdu7cuOKDFs3TO4ae5fUn1rf+HqNGVEen2mgPL/M\nyXgHn7FzBk1mNqFpxaaUCizl+tAEUC2+TBCBTz6Br7/Wt5msWdOJNwcGEpAolvb2BhUMon+9/tSe\nVJsVx1bQsnJLfnvhNyoGVbztPIczAxuqxecV+L342fy989fO89LPLzFyw0iaVmxKYkoiw1sMd61j\nA9KEL2LdOr8VPmc9o6Qkfbn5n3+GDRvg3nudvGBAAAGJWOr5rTmxhmazmxFyfwhbzm1hWNNhrg1c\ntmGB8CnPL3P8XvyCCgbxfM3nqTGhBlfjr9KkUhM61+js0MDlHLFv8Y0e7ZfC5yzXr0NIiL7h0B9/\nQJkyBgoJDKTQrRSPiF/4kXCuJ1xPe5yYnMjAlQPpsqgLE1pPIEVSzMkOQLX4vA0jubKRAws8v5SU\nFPl287cSPDpYxmwcY56/pxeuPL4M/POPyH/+I/LSSyIJCS4U9Pbbcv6jgVJjQg3TYssK+9VWTsac\nlLpT6kqFryvI9nPbXVuvMSPK43Mb5HbPL6O3d/H6RZ6c8yTjt45nZbeVRMZEmncHVx7fbRw5Ag0b\nQps2MG0a5MvnQmGBgRS4dfv2le7A1urvsrALdSbXIY+Wh529dvLP9X8Y3mK469kBqBaft2JEMY0c\nuLnlZ39nDj8SLqXGlJJHJj0iR68cNfcOnkWLz5/XTMupbhs2iJQuncOsDWcYOVKuvdVHSo8pbVKB\nOpnNxIiMjpSH331Y7v3mXnMzAxte0OLz59+miGr5pY3de2z6Y7y6/FUalmvI2hfXcjjqsHl3cNXi\nu405c6B9e5g1K4dZG84QEED+DHv3mkHGubnLDy+n5sSalA4o7dpG81mhWnzejRHFNHJgYssvszt4\n+OFwKTO2jHSY38E9d3Dl8aUjOVlk6FCRe+4R2bfP5MKnTZPknj3kjk/ukJSUFFOLjomLkV6/9JKX\nl70sgcMDZdbOWeZmBja8oMWXWyA3tfzs7+BxiXH0Du9Np0Wd+LTpp5QJLGP+HVy1+NJx8yZ06QIR\nEbBlCzz4oMkXCAjgjrh48t6Rl4TkBFOLPnTlEGtOrGH6rumsf2k9dwbcaV5mYEO1+HwDI4pp5MBk\nzy8mLkY6zu8o935zr1QZV0V2nt/pnju4gy0+f/ZV7Ot2/rzeo9utm0h8vJsuuHy5yNNPS9DIIIm6\nGeX02zPLDM79e07a/dROSo0pJU/OeVJORJ9I+32Y+t15YYvPn3+bIrms5Xc94TofR3zMHyf/4ETs\nCVZ3X825a+fMv4OrFl86tm6FunX1cXyzZ0OBAm66UEAA3LxJYD5j+3hk9PbCjoRR7ftqJEsybaq2\nYf6z87mn+D1pg9vtx/m5hGrx+RZGFNPIgYGWX2Z38Pn75kvJ0SWl66Ku8vKyl11fdy8rlMeXjhkz\nRIKDRZYt88DFNm0SqVdPqo6rKoevHDZURExcjLy6/FXptribBAwPkLm755qz7l5WeGGLL7eAwZaf\nV4uffep68fpF6TC/gxQdUVTm7JrjnhTXhhK+NBISRPr1E6laVWT/fg9ddPdukZo1pdYPteTv839n\ne2pmghZ9M1reXf2ulBpTSvgY2XVhlzujVcJnMUbFz6vT3qCCQXzW/DPah7an+vfVORFzgoN9DlK8\nUHHzU1wbBlNdf5w/efmynsVt3RrB1q1QvbqHLhwYCDduEJAv502MMqa4Oy/spNakWiw9uDRtT43J\nOyZn2/nl0nfnA6muP/42zcBrxC+zzWF+P/E79afW5+qtq0TFRbG0y1LuLnK36yvnZoXy+NLYvh0e\nfRQeewyGD4cgE+b0O4yd55fTmn62G997a95j0MpBNJjWgOdrPE/ze5ozPWS66wtXZIcPCJ8iG4w0\nF40c5JD22qeuV25ckR5LekjA8AAZs3GMvBH2hvu8PRsq1RUR/WOYMEGkZEmRxYstCuLqVZEiRSTk\npxBZenBptqempKTIL4d/kUpfVxI+RtafWu9eb8+GSnW9BvzB87t847I0ntFY7hx1p9ScUFN2/7Pb\nvd6eDSV8IiJy7ZrI88+LPPSQyJEjFgaSkCCSJ490XficzN09N+3pjKJ24NIBaT6zudz95d3y9I9P\nu/8GaUMJn1dhVPw8nvZmlt7GxMXwScQnNJnZhFtJt4iKi2J51+WcuXrGfd6eDZNSXV/3VQ4c0Iex\nFCwImzdD1ar/e83jdcuXD+64g6JaoXSen83fOxlzkrdWvMVjMx7jRuINWlZuyZwOcwynuE7VzwdT\nXV//bboLj4tfRoN6w+kNPPTDQ8zZM4ehjw/lkbseSZuh0ahCI/d4ezaUxwfAjz/qG4gPHqyvyFKo\nkNURAYGBBCXnS+f5FcpbiHJFy1F9QnUu37hMyP0hDGgwgK9afuXeG6QNHxQ+RTYYaS4aObBLe2Pi\nYuTFpS9K5wWdpdDnhWTUhlFy6folz6S4NlSqK//+K9Kjh8j994vscvNoEGcIOxwmyXfdJSNC+8rw\nP4dLckqyTNo2SUqNKSVPz3tafjvym3vmb2eHSnW9FnzF8zt65aiE/BQiRb8oKnyM7P5nt4hkPl7L\ndJPahhI+2b5dpEoVkVdeEbl+3epo0hMTFyP/3F1MvpzZS1rNbSU1vq8hpceUTvuN9A7r7Tl/T0QJ\nn5fjE+L39sq3pcBnBaT74u7SY2kPz/6AbbhJ+Hxl/mRyssjYsXpvbmioY++xom5JD9WQVz/5j1T+\ntrK0mtNKom9G35YNmJUdZFs/PxA+X/ltGsWo+HnU8ws/Gs6aF9dQpEARvmn1jXvHYGVGLvf4LlyA\n1q1h8WJ9nm6XLlZHlDV5Chfl07rvcjzmOBOfnkjxQsXZeHqj+zvA7FEen39jRDGNHKB7NB5Nb+3J\n5anuggUipUqJfPihSGKi1dHkTELzpjLus7bWZAciftHiyy3gC2mvJT9ikVwtfNHR+ti9++8X2bLF\n6mgcIyYuRnbVv0euz5+b9tijvx0lfD6FUfHzaNrr0RTXhodSXW8cS7VyJTz0EAQH69WvW9dYOZ6u\n28bTG6lW4RECkzTA/eltuvr5Yarrjb9Nb8Cj4ud2jyYjudTj+/dfeP11+L//g5kz4dtv9emyvkKb\n+9qQv0iQvmR0KqaO78wKPxQ+RdZoeqvRAxfSNPHUtYBcK3xhYdC7N7RqBWPG+HC1+/eHe+/V//UE\nSvh8Fk3TEBHN2ffldUcwlpMLhe/yZV0ntm7Vd1Jr1szqiFwkdWUXj6CEL1fiNUtamYZFwmeVryIC\n8+ZBzZpQtizs2WO+8FlSt9Q1/dzOvn1ENG7s18KnPL/M8a+WXy5r8R05olf34kX45Rd9/T2/ISAA\noqPdew1bi69PH78VPkXW+I/nl4uE7+ZNGDECJk6EoUOhXz/I61+3MZg0Sf8uJ01yT/kq1fUbcrfn\nl4uE75df4M03oV492L1bT3X9End6fkr4FPiD5+clwuduX+XIEWjbFt5+GyZPhtBQzwmfX3l+mQif\nv3ti/l4/o/i2+HmJ8LmTmBgYMAAaNoTGjWHvXv1v1+9xR8tPtfgUdviu5+fnwpeYqNtdn34KHTro\n/5YqZXVUHmT9et3Q/PNPc8pTwue35C7Pz4+FTwR+/hnef19Pa9eu1Yex5DoCAsxLe5XwKTLB99Je\nLxU+M3yVdeugQQMYNkyfnbFqlXcInyWekVlprwPC5++emL/Xzyi+1fLzUuFzlR079JbesWPw2Wfw\n3HNwh+/dlszFjA4P1eJTZIPveH5+KHy7dsHnn8OmTfDBB/DKK5A/v9VReQlXrsADD+j/GkEJX67B\nqHVcfVYAAAduSURBVOfnG+0LPxO+7dshJERfVblhQ30YyxtvKOFLhystPyV8CgfwfvHzEeFzxFfZ\ntEkXvPbt9b/N48dh4ED979ybscQzKlgQEhIgJcW59xkQPn/3xPy9fkbxbs/PR4QvO1JS9FkZX34J\nZ87Ae+/B0qVQoIDVkXk5mqZvIHzzJhQu7Nh7VItP4QTe6/n5uPDFxelLS339tR76oEH6eD2/m4Pr\nTkqV0gXNkQGOSvhyLf41zs+Hhe/0aX362ZQpUL++/u/jj+sNGYWTBAY6NtxFCZ/CAN7n+fmg8KWk\nwOjREYSEQO3acO2aPjHh55/1KWm+LnyWeUaODHQ2Qfj83RPz9/oZxbtafjbh27FD333Hy4Xvn39g\nzhx9Gpqmwbvv6guLensHhs+QU8tPtfgULuA9np+PtPji4/UOjJkz4a+/9J7bXr30ndF8vYXndTRt\nCh9/rP+bESV8ilR82/PzcuFLSdGHqcydCwsXQq1a0KMHLFigWnluJaspbkr4FCZgvefnpcKXkgIb\nN+qbAlWooG8FWa6cHuaaNdC9e3rh82dfxbK6ZTbQ2Q3C58/fHfh//YxibcvPy4Tv1i29o+KXX2DJ\nEiheHDp1gtWroVo1S0PLnWRs+akWn8JErPP8vET4LlyAX3+F8HB9+ajq1aFNG+jYUQme5fTuDTVq\n6P8q4VNkgW95fhYKX2ws/PGHvnzUunX6uLwnn9Q7LiZNgpIlPRaKIidsLT8lfAo3kKPnp2laK03T\nDmmadlTTtHezOGdc6uu7NU2rnW2BHhQ+EX1K2aJF+gyL//wHypeH77+H0qV1sbt0CebP1z08V4TP\nn30VSz2/bdvcLnz+/N2B/9fPKNmKn6ZpeYDxQCugOtBV07RqGc5pDVQRkarA/wETsyzQjcInoqew\nq1bBqFH6VLJy5aBOHX2aWVCQPtXsyhX9nCFD9BkY+fKZc/1du3aZU5AXYlndAgL0O5ebW3z+/N2B\n/9fPKDmlvXWBYyJyEkDTtFAgBDhod84zwCwAEdmiaVqQpmmlReTibaWZIHzJyXD2rL4iyrFjsH+/\nvqnPnj366zVr6kNROnfWFxOoVMkz4+9iY2PdfxGLsKxubdvCww9Dq1ZuvYw/f3fg//UzSk7iVxY4\nY/f4LFDPgXPKAbeLXzbCl5IC//4LV6/qqeiFC+kPm+CdPAnBwVC5sn7YOihq1oQyZdRAY7+ienX9\nUCjcQE7i52hXcEbJyfR9z+RfQWzbYiQm6ruT3bqli93Vq3D9um7xFCumL+Jx113/O2rV0gWucmW4\n9159pSNv4+TJk1aH4Db8uW6g6pdbyXaoi6Zp9YGPRaRV6uMhQIqIjLI75wcgQkRCUx8fAppkTHs1\nTfPMmBqFQpHrcMdQl+1AVU3TKgHngS5ARud5OdAXCE0Vy9jM/D4jwSkUCoW7yFb8RCRJ07S+wEog\nDzBNRA5qmtYr9fVJIvKrpmmtNU07BtwAXnJ71AqFQuEiHpvhoVAoFN6E6QsbmD4o2ovIqW6apr2Q\nWqc9mqZt1DTtISviNIoj313qeY9qmpakaVoHT8bnKg7+NptqmrZT07R9mqZFeDhEl3Dg9xmsadoK\nTdN2pdavpwVhGkLTtOmapl3UNG1vNuc4pysiYtqBnhofAyoB+YBdQLUM57QGfk39fz1gs5kxuOtw\nsG4NgGKp/2/lK3VztH525/0OhAEdrY7b5O8vCNgPlEt9HGx13CbX72NghK1uQBSQ1+rYHazf40Bt\nYG8WrzutK2a3/NIGRYtIImAbFG1PukHRQJCmaaVNjsMd5Fg3EdkkIldTH25BH+/oKzjy3QH0AxYB\nlz0ZnAk4Ur/ngcUichZARAzumG4JjtTvAlA09f9FgSgRSfJgjIYRkfVATDanOK0rZotfZgOeyzpw\nji+IhCN1s+cV4Fe3RmQuOdZP07Sy6H9QtimMvmQYO/L9VQVKaJq2TtO07ZqmdfdYdK7jSP2mAA9q\nmnYe2A3091BsnsBpXTF7VRdTB0V7GQ7HqGlaM+BloJH7wjEdR+r3DfCeiIimaRq3f4/ejCP1ywc8\nArQAAoBNmqZtFpGjbo3MHByp3/vALhFpqmlaZWC1pmkPi8g1N8fmKZzSFbPF7xxQ3u5xeXQFzu6c\ncqnPeTuO1I3UTo4pQCsRya6Z7m04Ur866OM5QfeMntI0LVFElnsmRJdwpH5ngCsiEgfEaZr2J/Aw\n4Avi50j9GgLDAUTkuKZpkcD96ON5fR3ndcVkUzIvcBzddM1Pzh0e9fGRTgEH61YB3XSub3W87qhf\nhvNnAB2sjtvk7+8BYA1650EAsBeobnXsJtbvK2BY6v9Lo4tjCatjd6KOlXCsw8MhXTG15Sd+PCja\nkboBHwHFgYmpraNEEalrVczO4GD9fBYHf5uHNE1bAewBUoApInLAuqgdx8Hv7wtghqZpu9H9/ndE\nJNqyoJ1A07SfgCZAsKZpZ4Bh6DaFYV1Rg5wVCkWuxPrd2xQKhcIClPgpFIpciRI/hUKRK1Hip1Ao\nciVK/BQKRa5EiZ9CociVKPFTKBS5EiV+CoUiV/L/msea8Ip5xD4AAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x81753b0>"
       ]
      }
     ],
     "prompt_number": 16
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "As before, using the inner product for this problem, gives:\n",
      "\n",
      "$$ \\int_0^1 \\left(\\frac{\\eta_1^2}{2}-h(\\eta_1)\\right)\\eta_1 d\\eta_1=0$$\n",
      "\n",
      "and the solution jumps right out as \n",
      "$$h(\\eta_1)=\\frac{\\eta_1^2}{2}  , \\hspace{1em} \\forall \\eta_1\\in[0,1$$\n",
      "\n",
      "where $\\eta_1(x)=2x$. Doing the same thing for $\\eta_2=2x-1$ gives,\n",
      "\n",
      "$$ \\int_0^1 \\left(\\frac{(1+\\eta_2)^2}{2}-h(\\eta_2)\\right)\\eta_1 d\\eta_2=0$$\n",
      "\n",
      "with  \n",
      "\n",
      "$$h(\\eta_2)=\\frac{(1+\\eta_2)^2}{2} , \\hspace{1em} \\forall \\eta_2\\in[0,1]$$ \n",
      "\n",
      "and then adding these up as before gives the full solution:\n",
      "\n",
      "$$ h(\\eta)= \\frac{1}{2} +\\eta + \\eta^2$$ \n",
      "\n",
      "Back-substituting each piece for $x$ produces the same solution as `sympy`."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "xs = np.random.rand(100)\n",
      "print np.mean([(2*i**2-hsol.subs(x,i))**2 for i in xs])\n",
      "print S.integrate((2*x**2-hsol)**2,(x,0,1)).evalf()\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0.282679822464177\n",
        "0.270833333333333\n"
       ]
      }
     ],
     "prompt_number": 17
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Summary"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We worked out some of the great examples in Brzezniak's book using our methods as a way to show multiple ways to solve the same problem. In particular, comparing Brzezniak's more measure-theoretic methods to our less abstract techniques is a great way to get a handle on those concepts which you will need for more advanced study in stochastic process. \n",
      "            \n",
      "As usual, the corresponding [IPython Notebook](www.ipython.org) notebook for this post  is available for download [here](https://github.com/unpingco/Python-for-Signal-Processing/blob/master/Conditional_expectation_MSE_Ex.ipynb)  \n",
      "\n",
      "Comments and corrections welcome!"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "References\n",
      "---------------\n",
      "\n",
      "* Brzezniak, Zdzislaw, and Tomasz Zastawniak. Basic stochastic processes: a course through exercises. Springer, 2000."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}