Xiu Paper.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# (1) Compute $k_i$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ k=\\sum_{i=0}^{P} k_i \\Phi_i(\\xi)$$\n",
    "\n",
    "$k_i$ are deterministic coefficients. We use the Galerkin projection to compute them:\n",
    "$$ k_i=\\frac{<k\\Phi_i(\\xi)>}{<\\Phi_i^2(\\xi)>}=\\frac{1}{<\\Phi_i^2(\\xi)>}\\int_{-\\infty}^{\\infty} \\frac{1}{\\sqrt{2\\pi}\n",
    "}k\\Phi_i(\\xi)e^{-\\frac{\\xi^2}{2}}\\ d\\xi$$\n",
    "\n",
    "<font color = red> To be more general </font>\n",
    "    \n",
    "It is kind of meaningless becuase we don't know what $k$ is. We are gonna do something called an inverse sample transform, the reason why we can do that is because $k$ and $\\xi$ are fully dependent\n",
    "\n",
    "<font color = red> Numerator: </font> Using the CDF of both random variables, transform them to the same **uniform** random variables supported on $[0,1]$\n",
    "$$u=G(\\xi)=F(k)$$\n",
    "$$k=F^{-1}(u)=h(u)$$\n",
    "$$\\xi=G^{-1}(u)=l(u)$$\n",
    "Using this transform, the integral is transformed:\n",
    "$$\\int_{-\\infty}^{\\infty} \\frac{1}{\\sqrt{2\\pi}}k\\Phi_i(\\xi)e^{-\\frac{\\xi^2}{2}}\\ d\\xi \\Longrightarrow \\int_0^1 h(u)\\Phi_i(l(u))\\ du$$\n",
    "\n",
    "<font color = red> To be more general </font>\n",
    "\n",
    "For the inner product term, three options:\n",
    "- use the closed form representation from the orthogonality condition: $<\\Phi_i^2(\\xi)>=i!$ <font color = red> Here since the weight function is the real pdf of $\\xi$, but weighting function in the package is just $e^{-\\xi^2/2}$ </font>\n",
    "- Choose a scaling of $\\xi$ such that the orthogonality condition is just the Kronecker delta:\n",
    "- Just compute the integral\n",
    "\n",
    "<font color = red> Denominator: </font> Option 2 is best, but tedious. So let's go with option three using Gaussian-Hermite quadrature\n",
    "$$ \\int_{-\\infty}^{\\infty} f(x)e^{-\\frac{x^2}{2}}\\approx\\sum_{i=0}^{N}w_i f(x_i)$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.polynomial.hermite_e as H\n",
    "from scipy.stats import norm\n",
    "import scipy.integrate as integrate\n",
    "import scipy.integrate.quadrature as G_qua\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The reason why it uses hermite_2 is because they are probabilists' Hermite polynomials, usually use notation $H_{e_0}, H_{e_1},\\ldots$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 158,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5.0"
      ]
     },
     "execution_count": 158,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H.hermeval(2, [0,1,1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 159,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Herm(n):\n",
    "    #define H_n\n",
    "    coeffs = [0]*(n+1)\n",
    "    coeffs[n] = 1\n",
    "    return coeffs\n",
    "\n",
    "def inner_product(h1, h2):\n",
    "    return lambda x: H.hermeval(x, H.hermemul(h1, h2))  \n",
    "\n",
    "\n",
    "def trapezoid_int(f, a, b, n=100):\n",
    "    P = [a + i*(b-a)/n for i in range(0, n+1)]\n",
    "    F = [1/2*np.abs(P[i+1] - P[i])*(f(P[i+1]) + f(P[i])) for i in range(0, n)]  \n",
    "    return sum(F)\n",
    "\n",
    "def unif_icdf(params):\n",
    "    a = params[0]\n",
    "    b = params[1]\n",
    "    return lambda u: u*(b-a)+a\n",
    "\n",
    "def expo_icdf(params):\n",
    "    return lambda u:-np.log(1-u)\n",
    "\n",
    "def norm_icdf(params):                          #k=h(u)\n",
    "    return lambda u:norm.ppf(u, loc=0, scale=1)\n",
    "\n",
    "def factorial(n):\n",
    "    k=1\n",
    "    if n==0:\n",
    "        k=1\n",
    "    else:\n",
    "        for i in range(1,n+1):\n",
    "            k=k*i\n",
    "    return k\n",
    "\n",
    "# def h_phi_l(h):\n",
    "#     return lambda x: h(x)*H.hermeval(norm.ppf(x,loc=0,scale=1),Herm(i))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 160,
   "metadata": {},
   "outputs": [],
   "source": [
    "def approximate_rv_coeffs(P,h):\n",
    "    #P is the number of Hermite-Chaos\n",
    "    #h is the inverse pdf of k\n",
    "    #initialize lists for output to make syntax more canonical with the math\n",
    "    ki=[0]*P\n",
    "    \n",
    "    #Set up Gauss-Hermite quadrature\n",
    "    m=P**2\n",
    "    x, w=H.hermegauss(m)        #x is point, w is the corresponding weight\n",
    "    \n",
    "    #Compute the coefficients, and also build out k in the same pass\n",
    "    for i in range(0, P):\n",
    "        #compute the inner product with Gaussian-Hermite quadrature, denominator\n",
    "        ip=sum([inner_product(Herm(i), Herm(i))(x[idx]) * w[idx] for idx in range(m)])         \n",
    "        #compute the integral, numerator\n",
    "        integrand=lambda u:h(u)*H.hermeval(norm.ppf(u, loc=0, scale=1), Herm(i))     #norm.ppf, since expansion of \\xi\n",
    "        #ki[i]=np.sqrt(2*np.pi)/G_qua(h_phi_l(h), 0, 1)\n",
    "        ki[i]=trapezoid_int(integrand, 0.001, 1-0.001, 1000)/factorial(i)\n",
    "    return ki     #ki is the coefficient"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 161,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_rv(ki, S):\n",
    "    #different P will give different ki\n",
    "    #build out k termwise\n",
    "    k = [0]*len(S)\n",
    "    for i in range(len(ki)):\n",
    "        k=np.add(k, ki[i]*H.hermeval(S, Herm(i)))\n",
    "    return k"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 162,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Generate a bunch of Gaussian random variables to use\n",
    "N=5000\n",
    "S=np.random.normal(loc=0, scale=1, size=N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 163,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAEa9JREFUeJzt3WuQZGddx/Hvjyyg3ExgJxh2s06wNkhIYZEaYpASgSCEhMryAqhNCSwY3RIiIqBkkRextKgKoHIpEFxITGJhQoxItriIGIJRiwQ3XEIuRNYQkyELOwjEC2Ug8PdFn8CwzM709Omey7PfT9XW9HnO033+z87Mb04/fS6pKiRJ7brfahcgSZosg16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUuA2rXQDAxo0ba3p6erXLkKR15frrr/96VU0t1W9NBP309DR79+5d7TIkaV1J8h/D9HPqRpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGrcmzoyVpNZN7/rwgu23n3/GxLftHr0kNW7JoE9yYZIDSW48qP0VSW5NclOSN81rf12Sfd26Z02iaEnS8IaZurkIeAdwyX0NSZ4GbAMeX1X3JDm6az8B2A48DngU8A9Jjq+q7427cEnScJbco6+qa4BvHNT8MuD8qrqn63Oga98GXFZV91TVl4F9wMljrFeStEyjztEfD/xSkuuS/GOSJ3btm4A75/Wb7dokSatk1KNuNgBHAacATwQuT/JoIAv0rYVeIMlOYCfAli1bRixDkrSUUffoZ4EP1MCnge8DG7v2Y+f12wzctdALVNXuqpqpqpmpqSVvkCJJGtGoQf9B4OkASY4HHgB8HdgDbE/ywCTHAVuBT4+jUEnSaJacuklyKfBUYGOSWeA84ELgwu6Qy+8AO6qqgJuSXA7cDNwLnOMRN5K0upYM+qo66xCrXniI/m8A3tCnKEnS+HhmrCQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcUsGfZILkxzo7iZ18LrfTVJJNnbLSfL2JPuS3JDkpEkULUka3jB79BcBpx3cmORY4FeAO+Y1P5vBfWK3AjuBd/UvUZLUx5JBX1XXAN9YYNVbgNcCNa9tG3BJDVwLHJnkmLFUKkkayZL3jF1IkjOBr1TV55PMX7UJuHPe8mzXtn/kClfQ9K4PL9h++/lnrHAlkjQ+yw76JA8CXg88c6HVC7TVAm0k2clgeoctW7YstwxJ0pBGOermZ4HjgM8nuR3YDHwmyU8z2IM/dl7fzcBdC71IVe2uqpmqmpmamhqhDEnSMJYd9FX1hao6uqqmq2qaQbifVFVfBfYAL+6OvjkFuLuq1sW0jSS1apjDKy8FPgU8JslskrMX6f4R4DZgH/Ae4OVjqVKSNLIl5+ir6qwl1k/Pe1zAOf3LkiSNi2fGSlLjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuOGuZXghUkOJLlxXtubk3wxyQ1J/jbJkfPWvS7JviS3JnnWpAqXJA1nmD36i4DTDmr7OHBiVT0e+DfgdQBJTgC2A4/rnvNnSY4YW7WSpGVbMuir6hrgGwe1/X1V3dstXgts7h5vAy6rqnuq6ssMbhJ+8hjrlSQt0zjm6H8N+Gj3eBNw57x1s13bj0myM8neJHvn5ubGUIYkaSEb+jw5yeuBe4H33de0QLda6LlVtRvYDTAzM7NgH+lwNr3rwwu2337+GStcida7kYM+yQ7gOcCpVXVfUM8Cx87rthm4a/TyJEl9jTR1k+Q04FzgzKr69rxVe4DtSR6Y5DhgK/Dp/mVKkka15B59kkuBpwIbk8wC5zE4yuaBwMeTAFxbVb9ZVTcluRy4mcGUzjlV9b1JFS9JWtqSQV9VZy3QfMEi/d8AvKFPUZN2qLnP9aLludv1NLb1VKsOb54ZK0mNM+glqXEGvSQ1zqCXpMb1OmFKklq22IEb6+lDd4NeY+ERKNLa5dSNJDXOPXpJQ1vuOSi+o1sbDPrD3Gr94q73k9akQ1mLP9sG/Qqa9Dy28+Rrg98HrTXO0UtS49yjH4J7aNLa5O/mcAx6aZ0x3LRcBv1hYJwfDq3FD5okLc45eklqnEEvSY0b5g5TFzK4N+yBqjqxa3s48H5gGrgdeEFVfTOD2029DTgd+Dbwkqr6zGRKb5/TJG1xbl2rZZg5+ouAdwCXzGvbBVxVVecn2dUtnws8m8F9YrcCvwC8q/t6WDGgJa0lw9xK8Jok0wc1b2NwH1mAi4FPMgj6bcAlVVXAtUmOTHJMVe0fV8FSH61cjXA5fCehUY+6eeR94V1V+5Mc3bVvAu6c12+2a/uxoE+yE9gJsGXLlhHL0OFunCFmIB6+Wn8XPu4PY7NAWy3Usap2V9VMVc1MTU2NuQxJ0n1GDfqvJTkGoPt6oGufBY6d128zcNfo5UmS+ho16PcAO7rHO4Ar57W/OAOnAHc7Py9Jq2uYwysvZfDB68Yks8B5wPnA5UnOBu4Ant91/wiDQyv3MTi88qUTqFnSBLU+X304Guaom7MOserUBfoWcE7foiRJ4+OZsZLUOC9qtgb4Vnlt8/uj9c6g10QZkivH/2sdikGvJhl60g8Z9JJWnGchryw/jJWkxrlHL2linEJbGwx6aZUZhj80rikdp4Z+lFM3ktQ49+ilw9R6eiexnmpdi9yjl6TGuUcvSSNYT+8y3KOXpMa5R9/DevqLLunwZdBrXfCPqjQ6p24kqXG99uiTvAr4dQY3AP8CgztKHQNcBjwc+Azwoqr6Ts86Jam3w/Wd4ch79Ek2Ab8NzFTVicARwHbgjcBbqmor8E3g7HEUKkkaTd+pmw3ATybZADwI2A88HbiiW38x8Nye25Ak9TBy0FfVV4A/ZnBz8P3A3cD1wLeq6t6u2yywaaHnJ9mZZG+SvXNzc6OWIUlaQp+pm6OAbcBxwKOABwPPXqBrLfT8qtpdVTNVNTM1NTVqGZKkJfSZunkG8OWqmquq7wIfAH4ROLKbygHYDNzVs0ZJUg99gv4O4JQkD0oS4FTgZuBq4Hldnx3Alf1KlCT10WeO/joGH7p+hsGhlfcDdgPnAq9Osg94BHDBGOqUJI2o13H0VXUecN5BzbcBJ/d5XUnS+HhmrCQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcb2CPsmRSa5I8sUktyR5UpKHJ/l4ki91X48aV7GSpOXru0f/NuDvqurngJ8HbgF2AVdV1Vbgqm5ZkrRKRg76JA8DnkJ3T9iq+k5VfQvYBlzcdbsYeG7fIiVJo+uzR/9oYA74iySfTfLeJA8GHllV+wG6r0cv9OQkO5PsTbJ3bm6uRxmSpMX0CfoNwEnAu6rqCcD/soxpmqraXVUzVTUzNTXVowxJ0mL6BP0sMFtV13XLVzAI/q8lOQag+3qgX4mSpD5GDvqq+ipwZ5LHdE2nAjcDe4AdXdsO4MpeFUqSetnQ8/mvAN6X5AHAbcBLGfzxuDzJ2cAdwPN7bkOS1EOvoK+qzwEzC6w6tc/rSpLGxzNjJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mN6x30SY5I8tkkH+qWj0tyXZIvJXl/d/cpSdIqGcce/SuBW+YtvxF4S1VtBb4JnD2GbUiSRtQr6JNsBs4A3tstB3g6cEXX5WLguX22IUnqp+8e/VuB1wLf75YfAXyrqu7tlmeBTT23IUnqYeSgT/Ic4EBVXT+/eYGudYjn70yyN8neubm5UcuQJC2hzx79k4Ezk9wOXMZgyuatwJFJNnR9NgN3LfTkqtpdVTNVNTM1NdWjDEnSYkYO+qp6XVVtrqppYDvwiar6VeBq4Hldtx3Alb2rlCSNbBLH0Z8LvDrJPgZz9hdMYBuSpCFtWLrL0qrqk8Anu8e3ASeP43UlSf15ZqwkNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXF9bg5+bJKrk9yS5KYkr+zaH57k40m+1H09anzlSpKWq88e/b3Aa6rqscApwDlJTgB2AVdV1Vbgqm5ZkrRKRr6VYFXtB/Z3j/87yS3AJmAb8NSu28UMbjF4bq8qRzS968OrsVlJWlPGMkefZBp4AnAd8Mjuj8B9fwyOHsc2JEmj6R30SR4C/A3wO1X1X8t43s4ke5PsnZub61uGJOkQegV9kvszCPn3VdUHuuavJTmmW38McGCh51bV7qqaqaqZqampPmVIkhbR56ibABcAt1TVn85btQfY0T3eAVw5enmSpL5G/jAWeDLwIuALST7Xtf0+cD5weZKzgTuA5/crUZLUR5+jbv4ZyCFWnzrq60qSxsszYyWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1Lj+hxHvyZ44TJJWpx79JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1LiJBX2S05LcmmRfkl2T2o4kaXETCfokRwDvBJ4NnACcleSESWxLkrS4Se3Rnwzsq6rbquo7wGXAtgltS5K0iEkF/SbgznnLs12bJGmFTerqlQvdNLx+pEOyE9jZLf5PklsnVMtG4OsTeu21yPG2zfE2Jm/8kcXljvdnhuk0qaCfBY6dt7wZuGt+h6raDeye0PZ/IMneqpqZ9HbWCsfbNsfbtkmNd1JTN/8KbE1yXJIHANuBPRPaliRpERPZo6+qe5P8FvAx4Ajgwqq6aRLbkiQtbmJ3mKqqjwAfmdTrL8PEp4fWGMfbNsfbtomMN1W1dC9J0rrlJRAkqXHNBP1Sl1xI8sAk7+/WX5dkeuWrHJ8hxvvqJDcnuSHJVUmGOgxrrRr2khpJnpekkqzrIzWGGW+SF3Tf45uS/NVK1zhOQ/w8b0lydZLPdj/Tp69GneOQ5MIkB5LceIj1SfL27v/ihiQn9d5oVa37fww+8P134NHAA4DPAycc1OflwLu7x9uB96923RMe79OAB3WPX9b6eLt+DwWuAa4FZla77gl/f7cCnwWO6paPXu26Jzze3cDLuscnALevdt09xvsU4CTgxkOsPx34KIPzkU4Bruu7zVb26Ie55MI24OLu8RXAqUkWOrFrPVhyvFV1dVV9u1u8lsG5DOvVsJfU+CPgTcD/rWRxEzDMeH8DeGdVfROgqg6scI3jNMx4C3hY9/inOOi8nPWkqq4BvrFIl23AJTVwLXBkkmP6bLOVoB/mkgs/6FNV9wJ3A49YkerGb7mXmDibwR7CerXkeJM8ATi2qj60koVNyDDf3+OB45P8S5Jrk5y2YtWN3zDj/QPghUlmGRzN94qVKW1VjP0SMhM7vHKFLXnJhSH7rBdDjyXJC4EZ4JcnWtFkLTreJPcD3gK8ZKUKmrBhvr8bGEzfPJXBu7V/SnJiVX1rwrVNwjDjPQu4qKr+JMmTgL/sxvv9yZe34saeVa3s0S95yYX5fZJsYPD2b7G3T2vZMOMlyTOA1wNnVtU9K1TbJCw13ocCJwKfTHI7g3nNPev4A9lhf56vrKrvVtWXgVsZBP96NMx4zwYuB6iqTwE/weC6MC0a6vd7OVoJ+mEuubAH2NE9fh7wieo++ViHlhxvN5Xx5wxCfj3P38IS462qu6tqY1VNV9U0g88kzqyqvatTbm/D/Dx/kMEH7iTZyGAq57YVrXJ8hhnvHcCpAEkeyyDo51a0ypWzB3hxd/TNKcDdVbW/zws2MXVTh7jkQpI/BPZW1R7gAgZv9/Yx2JPfvnoV9zPkeN8MPAT46+4z5zuq6sxVK7qHIcfbjCHH+zHgmUluBr4H/F5V/efqVT26Icf7GuA9SV7FYBrjJet1Ry3JpQym3DZ2nzmcB9wfoKrezeAziNOBfcC3gZf23uY6/b+SJA2plakbSdIhGPSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXu/wHUSMpdx8LpdwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1a23feba58>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "h = unif_icdf([0,1])                          #index k follows uniform distr\n",
    "ki_uniform = approximate_rv_coeffs(10, h)     #change the order\n",
    "k = generate_rv(ki_uniform, S)\n",
    "out = plt.hist(k, bins=50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAADsxJREFUeJzt3X+MZWV9x/H3p/gjaTUBu6MlwHTQrKZq7FomxMZoqNh2RQNiooU0llbS0UQSTfxD0KSa9h9aRZPGFrPGDZjQVVpEScXqlhhJE7HO6na7dEGBrjKw2V2hERsMzeK3f8zZeF1md+7ec+/cmWfer+RmznnuOfd8z/747LPPPec5qSokSe36lWkXIEmaLINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1LhnTbsAgC1bttTc3Ny0y5CkDWXPnj0/rqqZ1bZbF0E/NzfH4uLitMuQpA0lyQ+H2c6hG0lqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJaty6uDNWasnctV9Zsf3g9W9e40qkZfboJalxBr0kNc6gl6TGOUYvrRHH7jUt9uglqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWrcqkGfZGeSI0n2D7R9Icne7nUwyd6ufS7Jzwbe+/Qki5ckrW6YG6ZuAj4FfO54Q1X90fHlJDcAPxnY/sGq2jauAiVJ/awa9FV1d5K5ld5LEuAdwBvGW5YkaVz6jtG/DjhcVT8YaDs/yfeSfDPJ63p+viSpp75z3VwJ7BpYPwTMVtVjSS4AvpTkFVX1xIk7JlkAFgBmZ2d7lqFWTHM+mGkd2zlwNGkjB32SZwFvAy443lZVTwFPdct7kjwIvBRYPHH/qtoB7ACYn5+vUevQ5mZISqvrM3TzRuC+qlo63pBkJskZ3fKLga3AQ/1KlCT1MczllbuAbwEvS7KU5OrurSv45WEbgNcD+5L8B/BPwHuq6vFxFixJOj3DXHVz5Una/3SFttuA2/qXJUkaF++MlaTGGfSS1DgfJagN4WRX10hanUGvTcd/NLTZOHQjSY2zRy+twpuytNHZo5ekxhn0ktQ4h240FX4hKq0de/SS1DiDXpIa59CNmuTQkPQL9uglqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWrcMA8H35nkSJL9A20fTfJIkr3d65KB965L8kCS+5P84aQKlyQNZ5ge/U3A9hXaP1lV27rXnQBJXg5cAbyi2+fvk5wxrmIlSadv1Ttjq+ruJHNDft5lwOer6ingv5M8AFwIfGvkCqV1yrtvtVH0GaO/Jsm+bmjnrK7tHODhgW2WurZnSLKQZDHJ4tGjR3uUIUk6lVGD/kbgJcA24BBwQ9eeFbatlT6gqnZU1XxVzc/MzIxYhiRpNSNNalZVh48vJ/kM8M/d6hJw3sCm5wKPjlydpGfw0YY6XSP16JOcPbB6OXD8ipw7gCuSPDfJ+cBW4N/7lShJ6mPVHn2SXcBFwJYkS8BHgIuSbGN5WOYg8G6Aqro3ya3AfwHHgPdW1dOTKV2SNIxUrTiEvqbm5+drcXFx2mWoB69AmT6HbjafJHuqan617XzwiLSJOd6/OTgFgiQ1zqCXpMYZ9JLUOINekhrnl7FSI/xiVSdjj16SGmePXqfF6+U3Hn/PZI9ekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuNWDfokO5McSbJ/oO1jSe5Lsi/J7UnO7Nrnkvwsyd7u9elJFi9JWt0wPfqbgO0ntO0GXllVrwK+D1w38N6DVbWte71nPGVKkka1atBX1d3A4ye0fb2qjnWr9wDnTqA2SdIYjGOM/l3AVwfWz0/yvSTfTPK6MXy+JKmHXtMUJ/kwcAy4pWs6BMxW1WNJLgC+lOQVVfXECvsuAAsAs7OzfcqQJJ3CyD36JFcBbwH+uKoKoKqeqqrHuuU9wIPAS1fav6p2VNV8Vc3PzMyMWoYkaRUjBX2S7cAHgUur6smB9pkkZ3TLLwa2Ag+No1BJ0mhWHbpJsgu4CNiSZAn4CMtX2TwX2J0E4J7uCpvXA3+Z5BjwNPCeqnp8xQ+WJK2JVYO+qq5cofmzJ9n2NuC2vkVJksbHO2MlqXE+HHyTONkDog9e/+Y1rkTSWrNHL0mNM+glqXEO3WxyDulI7bNHL0mNM+glqXEO3Uh6Bof02mKPXpIaZ9BLUuMMeklqnGP0knpzTH99s0cvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGjdU0CfZmeRIkv0DbS9IsjvJD7qfZ3XtSfK3SR5Isi/J70yqeEnS6obt0d8EbD+h7VrgrqraCtzVrQO8CdjavRaAG/uXKUka1VBBX1V3A4+f0HwZcHO3fDPw1oH2z9Wye4Azk5w9jmIlSaevzxj9i6rqEED384Vd+znAwwPbLXVtkqQpmMQUCFmhrZ6xUbLA8tAOs7OzEyhDfZzslnZJG0+foD+c5OyqOtQNzRzp2peA8wa2Oxd49MSdq2oHsANgfn7+Gf8QSFp/7ABsTH2Gbu4AruqWrwK+PND+J93VN68BfnJ8iEeStPaG6tEn2QVcBGxJsgR8BLgeuDXJ1cCPgLd3m98JXAI8ADwJ/NmYa5YknYahgr6qrjzJWxevsG0B7+1TlEbnf60lncg7YyWpcQa9JDXOJ0xJmhifPLU+2KOXpMYZ9JLUOINekhpn0EtS4/wyVtKa80vatWWPXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGuedsZLWPe+k7WfkoE/yMuALA00vBv4COBP4c+Bo1/6hqrpz5AolSb2MHPRVdT+wDSDJGcAjwO0sPwz8k1X18bFUKEnqZVxj9BcDD1bVD8f0eZKkMRlX0F8B7BpYvybJviQ7k5w1pmNIkkbQO+iTPAe4FPjHrulG4CUsD+scAm44yX4LSRaTLB49enSlTSRJYzCOHv2bgO9W1WGAqjpcVU9X1c+BzwAXrrRTVe2oqvmqmp+ZmRlDGZKklYzj8sorGRi2SXJ2VR3qVi8H9o/hGBpwskvNJGklvYI+ya8Cvw+8e6D5b5JsAwo4eMJ7kqQ11ivoq+pJ4NdPaHtnr4okSWPlFAiS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOB4+sA97pKmmS7NFLUuMMeklqnEEvSY1zjF7SuuH3VZNhj16SGmfQS1LjDHpJapxj9JI2rJON6R+8/s1rXMn6Zo9ekhpnj15Sc+zp/zKDfg156Zikaegd9EkOAj8FngaOVdV8khcAXwDmWH5A+Duq6n/6HmujMNAlrSfjGqP/varaVlXz3fq1wF1VtRW4q1uXJE3BpL6MvQy4uVu+GXjrhI4jSVrFOIK+gK8n2ZNkoWt7UVUdAuh+vvDEnZIsJFlMsnj06NExlCFJWsk4vox9bVU9muSFwO4k9w2zU1XtAHYAzM/P1xjqkCStoHePvqoe7X4eAW4HLgQOJzkboPt5pO9xJEmj6RX0SX4tyfOPLwN/AOwH7gCu6ja7Cvhyn+NIkkbXd+jmRcDtSY5/1j9U1b8k+Q5wa5KrgR8Bb+95HEnSiHoFfVU9BPz2Cu2PARf3+WxJ0ng4140kNc6gl6TGGfSS1DgnNevBOW0kbQT26CWpcQa9JDXOoZshOEQjte1Uf8dbeFiJQS9p09isnTaHbiSpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1buSgT3Jekm8kOZDk3iTv69o/muSRJHu71yXjK1eSdLr6zHVzDPhAVX03yfOBPUl2d+99sqo+3r88SVJfIwd9VR0CDnXLP01yADhnXIVJksZjLGP0SeaAVwPf7pquSbIvyc4kZ43jGJKk0fSepjjJ84DbgPdX1RNJbgT+Cqju5w3Au1bYbwFYAJidne1bhiRNxMmmNt5I89T36tEneTbLIX9LVX0RoKoOV9XTVfVz4DPAhSvtW1U7qmq+quZnZmb6lCFJOoU+V90E+CxwoKo+MdB+9sBmlwP7Ry9PktRXn6Gb1wLvBP4zyd6u7UPAlUm2sTx0cxB4d68KJ2CzPmVG0ubU56qbfwOywlt3jl6OJGncvDNWkhpn0EtS4wx6SWpc7+voJUm/sB6vuzfoJWkEG+nqvSaCfiP9gkvanKbZ03eMXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNm1jQJ9me5P4kDyS5dlLHkSSd2kSCPskZwN8BbwJezvIDw18+iWNJkk5tUj36C4EHquqhqvo/4PPAZRM6liTpFCYV9OcADw+sL3VtkqQ1NqkHj2SFtvqlDZIFYKFb/d8k90+olnHYAvx42kVMkefv+W/W85/4ueeve+3+m8NsNKmgXwLOG1g/F3h0cIOq2gHsmNDxxyrJYlXNT7uOafH8Pf/Nev6tnPukhm6+A2xNcn6S5wBXAHdM6FiSpFOYSI++qo4luQb4GnAGsLOq7p3EsSRJpzaxh4NX1Z3AnZP6/DW2IYaYJsjz39w28/k3ce6pqtW3kiRtWE6BIEmNM+iHlORjSe5Lsi/J7UnOnHZNaynJ25Pcm+TnSTb8VQjD2MzTeCTZmeRIkv3TrmUakpyX5BtJDnR/7t837Zr6MOiHtxt4ZVW9Cvg+cN2U61lr+4G3AXdPu5C14DQe3ARsn3YRU3QM+EBV/RbwGuC9G/n336AfUlV9vaqOdav3sHxvwKZRVQeqaj3f1DZum3oaj6q6G3h82nVMS1Udqqrvdss/BQ6wge/uN+hH8y7gq9MuQhPlNB4CIMkc8Grg29OtZHQTu7xyI0ryr8BvrPDWh6vqy902H2b5v3W3rGVta2GY899EVp3GQ+1L8jzgNuD9VfXEtOsZlUE/oKreeKr3k1wFvAW4uBq8LnW1899kVp3GQ21L8myWQ/6WqvritOvpw6GbISXZDnwQuLSqnpx2PZo4p/HYxJIE+CxwoKo+Me16+jLoh/cp4PnA7iR7k3x62gWtpSSXJ1kCfhf4SpKvTbumSeq+eD8+jccB4NbNNI1Hkl3At4CXJVlKcvW0a1pjrwXeCbyh+/u+N8kl0y5qVN4ZK0mNs0cvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJatz/A06twthOVjm4AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1a1a458630>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "h = norm_icdf([])                          #index k follows norm distr\n",
    "ki_norm = approximate_rv_coeffs(5, h)     #change the order\n",
    "k = generate_rv(ki_norm, S)\n",
    "out = plt.hist(k, bins=50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[4.770489558936219e-18,\n",
       " 0.9774860948298655,\n",
       " 2.8796409701214993e-16,\n",
       " -0.03233346505595982,\n",
       " 1.214306433183764e-17]"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ki_norm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.0133265865433016"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAEuBJREFUeJzt3W2MneV95/HvLzgkgTYxDwOybLOmipV2FSmEHVF3kaIuTiucRLFfhBXRbnGRJfcF7SbLSq3bN1GlfUGkVUkiVUgWTtfspiSEBGGlKC0yRN2+gMYmFEKcLA4leGoXT8pDSlE2S/vfF+fyZjI+9pzxnPGZufh+pKP7vq/7Ouf8j0G/c8117odUFZKkfr1l0gVIkpaXQS9JnTPoJalzBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUucMeknq3JpJFwBw+eWX16ZNmyZdhiStKocPH/5hVU0t1G9FBP2mTZs4dOjQpMuQpFUlyQ9G6efUjSR1zqCXpM4Z9JLUOYNekjpn0EtS50YK+iT/OckzSb6d5N4kb09ydZLHkzyb5EtJLmx939a2j7b9m5bzA0iSzm7BoE+yHvhPwHRVvRe4ALgZ+DRwZ1VtBl4GdrWn7AJerqp3A3e2fpKkCRl16mYN8I4ka4CLgBPADcD9bf9+YEdb3962afu3Jsl4ypUkLdaCQV9Vfwf8N+AFBgH/KnAYeKWq3mjdZoD1bX09cKw9943W/7L5r5tkd5JDSQ7Nzs4u9XNIks5gwTNjk1zCYJR+NfAK8GVg25Cup+4yPmz0ftodyKtqL7AXYHp6+pzvUL5pz58NbX/+jg+f60tKUldGmbr5IPC3VTVbVf8X+Crwb4G1bSoHYANwvK3PABsB2v53AS+NtWpJ0shGCfoXgC1JLmpz7VuB7wCPAh9rfXYCD7b1A22btv+RqjrnEbskaWlGmaN/nMGPqk8AT7fn7AV+D7g9yVEGc/D72lP2AZe19tuBPctQtyRpRCNdvbKqPgV8al7zc8B1Q/r+GLhp6aVJksbBM2MlqXMGvSR1zqCXpM4Z9JLUOYNekjpn0EtS5wx6SeqcQS9JnTPoJalzBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUucMeknqnEEvSZ1bMOiTvCfJk3MeP0ryySSXJnk4ybNteUnrnySfS3I0yVNJrl3+jyFJOpNRbiX4vaq6pqquAf4N8DrwAINbBB6sqs3AQX56y8BtwOb22A3ctRyFS5JGs9ipm63A96vqB8B2YH9r3w/saOvbgXtq4DFgbZJ1Y6lWkrRoiw36m4F72/qVVXUCoC2vaO3rgWNznjPT2iRJEzBy0Ce5EPgo8OWFug5pqyGvtzvJoSSHZmdnRy1DkrRIixnRbwOeqKoX2/aLp6Zk2vJka58BNs553gbg+PwXq6q9VTVdVdNTU1OLr1ySNJLFBP3H+em0DcABYGdb3wk8OKf9lnb0zRbg1VNTPJKk82/NKJ2SXAT8GvBbc5rvAO5Lsgt4AbiptT8EfAg4yuAInVvHVq0kadFGCvqqeh24bF7bPzA4Cmd+3wJuG0t1kqQl88xYSeqcQS9JnTPoJalzBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUucMeknqnEEvSZ0z6CWpcwa9JHXOoJekzhn0ktQ5g16SOmfQS1LnDHpJ6txIQZ9kbZL7k3w3yZEkv5Lk0iQPJ3m2LS9pfZPkc0mOJnkqybXL+xEkSWcz6oj+s8DXq+oXgfcBR4A9wMGq2gwcbNsA24DN7bEbuGusFUuSFmXBoE/yTuADwD6AqvpJVb0CbAf2t277gR1tfTtwTw08BqxNsm7slUuSRjLKiP4XgFngT5J8K8ndSS4GrqyqEwBteUXrvx44Nuf5M61NkjQBowT9GuBa4K6qej/wT/x0mmaYDGmr0zolu5McSnJodnZ2pGIlSYs3StDPADNV9Xjbvp9B8L94akqmLU/O6b9xzvM3AMfnv2hV7a2q6aqanpqaOtf6JUkLWDDoq+rvgWNJ3tOatgLfAQ4AO1vbTuDBtn4AuKUdfbMFePXUFI8k6fxbM2K/3wG+kORC4DngVgZfEvcl2QW8ANzU+j4EfAg4Crze+kqSJmSkoK+qJ4HpIbu2DulbwG1LrEuSNCaeGStJnTPoJalzBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUucMeknqnEEvSZ0z6CWpcwa9JHXOoJekzhn0ktQ5g16SOmfQS1LnDHpJ6pxBL0mdGynokzyf5OkkTyY51NouTfJwkmfb8pLWniSfS3I0yVNJrl3ODyBJOrvFjOj/XVVdU1Wnbim4BzhYVZuBg20bYBuwuT12A3eNq1hJ0uItZepmO7C/re8Hdsxpv6cGHgPWJlm3hPeRJC3BqEFfwF8kOZxkd2u7sqpOALTlFa19PXBsznNnWpskaQLWjNjv+qo6nuQK4OEk3z1L3wxpq9M6Db4wdgNcddVVI5YhSVqskUb0VXW8LU8CDwDXAS+empJpy5Ot+wywcc7TNwDHh7zm3qqarqrpqampc/8EkqSzWjDok1yc5OdPrQO/DnwbOADsbN12Ag+29QPALe3omy3Aq6emeCRJ598oUzdXAg8kOdX/T6vq60m+CdyXZBfwAnBT6/8Q8CHgKPA6cOvYq5YkjWzBoK+q54D3DWn/B2DrkPYCbhtLdZKkJfPMWEnqnEEvSZ0z6CWpcwa9JHVu1BOmVp1Ne/5saPvzd3z4PFciSZPliF6SOmfQS1LnDHpJ6pxBL0mdM+glqXPdHnVzJmc6Ggc8IkdSnxzRS1LnDHpJ6pxBL0mdM+glqXMGvSR1zqCXpM6NHPRJLkjyrSRfa9tXJ3k8ybNJvpTkwtb+trZ9tO3ftDylS5JGsZgR/SeAI3O2Pw3cWVWbgZeBXa19F/ByVb0buLP1kyRNyEhBn2QD8GHg7rYd4Abg/tZlP7CjrW9v27T9W1t/SdIEjDqi/wzwu8C/tO3LgFeq6o22PQOsb+vrgWMAbf+rrb8kaQIWDPokHwFOVtXhuc1DutYI++a+7u4kh5Icmp2dHalYSdLijTKivx74aJLngS8ymLL5DLA2yalr5WwAjrf1GWAjQNv/LuCl+S9aVXurarqqpqemppb0ISRJZ7Zg0FfV71fVhqraBNwMPFJV/wF4FPhY67YTeLCtH2jbtP2PVNVpI3pJ0vmxlOPofw+4PclRBnPw+1r7PuCy1n47sGdpJUqSlmJRlymuqm8A32jrzwHXDenzY+CmMdQmSRoDz4yVpM696W48cjZnuimJNySRtJo5opekzhn0ktQ5g16SOmfQS1LnDHpJ6pxH3YzAo3EkrWaO6CWpcwa9JHXOoJekzhn0ktQ5g16SOmfQS1LnDHpJ6pxBL0mdM+glqXMLBn2Styf56yR/k+SZJH/Y2q9O8niSZ5N8KcmFrf1tbfto279peT+CJOlsRhnR/x/ghqp6H3ANcGOSLcCngTurajPwMrCr9d8FvFxV7wbubP0kSROyYNDXwGtt863tUcANwP2tfT+wo61vb9u0/VuTZGwVS5IWZaSLmiW5ADgMvBv4Y+D7wCtV9UbrMgOsb+vrgWMAVfVGkleBy4AfjrHuFcGLnUlaDUb6Mbaq/rmqrgE2ANcBvzSsW1sOG73X/IYku5McSnJodnZ21HolSYu0qKNuquoV4BvAFmBtklN/EWwAjrf1GWAjQNv/LuClIa+1t6qmq2p6amrq3KqXJC1olKNuppKsbevvAD4IHAEeBT7Wuu0EHmzrB9o2bf8jVXXaiF6SdH6MMke/Dtjf5unfAtxXVV9L8h3gi0n+K/AtYF/rvw/4H0mOMhjJ37wMdUuSRrRg0FfVU8D7h7Q/x2C+fn77j4GbxlKdJGnJPDNWkjrnPWOXgYddSlpJHNFLUucMeknqnEEvSZ0z6CWpcwa9JHXOoJekzhn0ktQ5j6NfATzuXtJyMujPozMFuiQtJ6duJKlzBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUudGuWfsxiSPJjmS5Jkkn2jtlyZ5OMmzbXlJa0+SzyU5muSpJNcu94eQJJ3ZKCP6N4D/UlW/BGwBbkvyr4E9wMGq2gwcbNsA24DN7bEbuGvsVUuSRrZg0FfViap6oq3/I3AEWA9sB/a3bvuBHW19O3BPDTwGrE2ybuyVS5JGsqg5+iSbGNwo/HHgyqo6AYMvA+CK1m09cGzO02ZamyRpAkYO+iQ/B3wF+GRV/ehsXYe01ZDX253kUJJDs7Ozo5YhSVqkkYI+yVsZhPwXquqrrfnFU1MybXmytc8AG+c8fQNwfP5rVtXeqpququmpqalzrV+StIBRjroJsA84UlV/NGfXAWBnW98JPDin/ZZ29M0W4NVTUzySpPNvlKtXXg/8BvB0kidb2x8AdwD3JdkFvADc1PY9BHwIOAq8Dtw61orfRLx8saRxWDDoq+qvGD7vDrB1SP8CbltiXToLvwAkLYZnxkpS57zxSEcc6UsaxhG9JHXOoJekzhn0ktQ5g16SOmfQS1LnDHpJ6pyHV74JeNil9ObmiF6SOmfQS1LnDHpJ6pxBL0mdM+glqXMGvSR1zqCXpM4Z9JLUuQVPmEryeeAjwMmqem9ruxT4ErAJeB7491X1cru/7GcZ3ErwdeA3q+qJ5SldS3WmE6nAk6mknowyov/vwI3z2vYAB6tqM3CwbQNsAza3x27grvGUKUk6VwsGfVX9JfDSvObtwP62vh/YMaf9nhp4DFibZN24ipUkLd65ztFfWVUnANryita+Hjg2p99Ma5MkTci4L2qWIW01tGOym8H0DlddddWYy9BSeSE0qR/nGvQvJllXVSfa1MzJ1j4DbJzTbwNwfNgLVNVeYC/A9PT00C8DrTx+AUirz7lO3RwAdrb1ncCDc9pvycAW4NVTUzySpMkY5fDKe4FfBS5PMgN8CrgDuC/JLuAF4KbW/SEGh1YeZXB45a3LULMkaREWDPqq+vgZdm0d0reA25ZalCRpfLzDlMbCuXtp5fISCJLUOYNekjpn0EtS5wx6SeqcP8ZqWfkjrTR5juglqXMGvSR1zqkbrSiLnepxakhamEGviTjb3a0kjZdTN5LUOYNekjrn1I1WhcVO9Th3L/2UI3pJ6pxBL0mdc+pGbyqLnQJyqkc9MOilszjbF4NfAlotliXok9wIfBa4ALi7qu5YjveRJskffLVajD3ok1wA/DHwa8AM8M0kB6rqO+N+L2kl8gtAK81yjOivA45W1XMASb4IbAcMer2pjetsYL8wtFjLEfTrgWNztmeAX16G95HelJb7C2Ml/kXij+hLsxxBnyFtdVqnZDewu22+luR7S3jPy4EfLuH554M1Lt1Krw9WUY359OKetNj+S7Skf8fzVOtK+G/9r0bptBxBPwNsnLO9ATg+v1NV7QX2juMNkxyqqulxvNZyscalW+n1gTWOizWO13KcMPVNYHOSq5NcCNwMHFiG95EkjWDsI/qqeiPJbwN/zuDwys9X1TPjfh9J0miW5Tj6qnoIeGg5XvsMxjIFtMyscelWen1gjeNijWOUqtN+J5UkdcSLmklS51Z10Ce5Mcn3khxNsmfS9QyT5PNJTib59qRrGSbJxiSPJjmS5Jkkn5h0TfMleXuSv07yN63GP5x0TWeS5IIk30rytUnXMkyS55M8neTJJIcmXc8wSdYmuT/Jd9v/l78y6ZrmSvKe9u936vGjJJ+cdF1ns2qnbtqlFv43cy61AHx8pV1qIckHgNeAe6rqvZOuZ74k64B1VfVEkp8HDgM7VtK/Y5IAF1fVa0neCvwV8ImqemzCpZ0mye3ANPDOqvrIpOuZL8nzwHRVTfr47zNKsh/4X1V1dzty76KqemXSdQ3TcujvgF+uqh9Mup4zWc0j+v9/qYWq+glw6lILK0pV/SXw0qTrOJOqOlFVT7T1fwSOMDi7ecWogdfa5lvbY8WNUJJsAD4M3D3pWlarJO8EPgDsA6iqn6zUkG+2At9fySEPqzvoh11qYUUF1GqTZBPwfuDxyVZyujYl8iRwEni4qlZcjcBngN8F/mXShZxFAX+R5HA7O32l+QVgFviTNgV2d5KLJ13UWdwM3DvpIhaymoN+pEstaDRJfg74CvDJqvrRpOuZr6r+uaquYXCm9XVJVtQ0WJKPACer6vCka1nA9VV1LbANuK1NLa4ka4Brgbuq6v3APwEr9fe3C4GPAl+edC0LWc1BP9KlFrSwNu/9FeALVfXVSddzNu3P+G8AN064lPmuBz7a5sC/CNyQ5H9OtqTTVdXxtjwJPMBgCnQlmQFm5vzFdj+D4F+JtgFPVNWLky5kIas56L3Uwhi0Hzr3AUeq6o8mXc8wSaaSrG3r7wA+CHx3slX9rKr6/araUFWbGPy/+EhV/ccJl/UzklzcfnCnTYf8OrCijgarqr8HjiV5T2vaysq9xPnHWQXTNrCKbyW4Wi61kORe4FeBy5PMAJ+qqn2TrepnXA/8BvB0mwMH+IN2dvNKsQ7Y345weAtwX1WtyMMXV7grgQcG3+2sAf60qr4+2ZKG+h3gC20A9xxw64TrOU2Sixgc8fdbk65lFKv28EpJ0mhW89SNJGkEBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUucMeknqnEEvSZ37f4b+TVj19CkFAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1a1a68b1d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "h = expo_icdf([])                          #index k follows expo distr\n",
    "ki_expo = approximate_rv_coeffs(3, h)     #change the order\n",
    "k = generate_rv(ki_expo, S)\n",
    "out = plt.hist(k, bins=50)\n",
    "k.min()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# (2) Compute $y_i$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "def approx_inner2(P,i,j):\n",
    "    #P is the number of Hermite-Chaos\n",
    "\n",
    "    #Set up Gauss-Hermite quadrature\n",
    "    m=P**2\n",
    "    x, w=H.hermegauss(m)        #x is point, w is the corresponding weight\n",
    "    inner=sum([inner_product(Herm(i), Herm(j))(x[idx]) * w[idx] for idx in range(m)])         \n",
    "    \n",
    "    return np.sqrt(2*np.pi)*inner"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "-----\n",
    "Doesn't work"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "def approx_inner3(P,i,j,l):\n",
    "    #P is the number of Hermite-Chaos\n",
    "\n",
    "    #Set up Gauss-Hermite quadrature\n",
    "    m=P**2\n",
    "    x, w=H.hermegauss(m)        #x is point, w is the corresponding weight\n",
    "    inner=sum([inner_product(inner_product(Herm(i), Herm(j)),Herm(l))(x[idx]) * w[idx] for idx in range(m)])         \n",
    "    \n",
    "    return np.sqrt(2*np.pi)*inner"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "def inner3_product(i,j,l):\n",
    "    return lambda x: H.hermeval(x, H.hermemul(H.hermemul(Herm(i),Herm(j)),Herm(l))) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [],
   "source": [
    "def approx_inner3(P,i,j,l):\n",
    "    #P is the number of Hermite-Chaos\n",
    "\n",
    "    #Set up Gauss-Hermite quadrature\n",
    "    m=P**2\n",
    "    x, w=H.hermegauss(m)        #x is point, w is the corresponding weight\n",
    "    inner=sum([inner3_product(i,j,l)(x[idx]) * w[idx] for idx in range(m)])         \n",
    "    \n",
    "#     return np.sqrt(2*np.pi)*inner\n",
    "    return inner/np.sqrt(2*np.pi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4.602162686824381e-17"
      ]
     },
     "execution_count": 94,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "approx_inner3(3,1,1,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ode_system(y, t):\n",
    "    y0, y1, y2 = y\n",
    "    dydt = [(-approx_inner3(3,1,0,0)*y0 - approx_inner3(3,1,1,0)*y1 - approx_inner3(3,1,2,0)*y2),\n",
    "            (-approx_inner3(3,1,0,1)*y0 - approx_inner3(3,1,1,1)*y1 - approx_inner3(3,1,2,1)*y2),\n",
    "            (-approx_inner3(3,1,0,2)*y0 - approx_inner3(3,1,1,2)*y1 - approx_inner3(3,1,2,2)*y2)/2]\n",
    "    return dydt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ode_system(y, t):\n",
    "    y0, y1, y2 = y\n",
    "    dy0dt = (-inner3(3,1,0,0)*y0 - inner3(3,1,1,0)*y1 - inner3(3,1,2,0)*y2)\n",
    "    dy1dt = (-inner3(3,1,0,1)*y0 - inner3(3,1,1,1)*y1 - inner3(3,1,2,1)*y2)\n",
    "    dy2dt = (-inner3(3,1,0,2)*y0 - inner3(3,1,1,2)*y1 - inner3(3,1,2,2)*y2)/2\n",
    "    return dy0dt, dy1dt, dy2dt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.integrate import odeint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol = odeint(ode_system, [1.0, 0.0, 0.0], np.linspace(0,1,101))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1a1dcf3320>]"
      ]
     },
     "execution_count": 99,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3XmYHGdh5/Hv29f0TM99SJpDM6NjdFvWMZJ84PvgMjaYy4YcQBItEEM22WyWXXafZMk+uyTAE9iYAAYcQpLFsIQFm8uAsZGxLVmjw7pl3ZpLmvu++nj3j7dnNJJ1zGh6umemfx+eeqqqq9T1FpLfX9f7vlVlrLWIiEj68aS6ACIikhoKABGRNKUAEBFJUwoAEZE0pQAQEUlTCgARkTSlABARSVMKABGRNKUAEBFJU75UF+BqiouLbXV1daqLISIya+zatavNWlsykX1ndABUV1dTV1eX6mKIiMwaxpgzE91XTUAiImlKASAikqYUACIiaUoBICKSphQAIiJpSgEgIpKmFAAiImlKASAiMkP0DoV55rUmvvLCiaQcb0bfCCYiMte19A7x3OEWfnHwHC8db2ckGqMsL8gf3rYIv3d6f6MrAEREkshay4nWfn556Dy/PHSOPfVdWAsVBZn87s1VvHn1AjZWFeD1mGkviwJARGSahaMxdp7u4NeHW3juSAun2voBWFOey7+/Zxn3r57PigU5GDP9lf54CgARkWnQ0T/CC0ddhb/t9VZ6hyIEvB5uXlLER960iHtWzKMsPzOlZVQAiIgkQCxmOdDUzQtHW3n+aAt74007xdkZvGX1Au5ZOZ/baooJZcycajchJTHGPAk8ALRYa9dcZvudwI+AU/GPfmCt/Uwiji0ikiod/SO8eKyV3xxtZduxVtr6RjAGbqzI50/uqeHuFfNYU5aHJwnt+dcjUVH0LeBx4NtX2edFa+0DCTqeiEjShaMx9tZ3se31Vra93sq+xm6shYIsP7cvK+GO+FSUnZHqok5IQgLAWrvNGFOdiO8SEZkprLWcbOvnpeNtvHisjVdOtNM3HMFjYH1lAf/+nmXcvqyYtRX5SRm1k2jJbIy62RjzGtAE/Lm19mASjy0iMiEtPUO8dKKNl46389LxNpq7hwBYWJjJg+vKuL2mmJuXFJOX6U9xSacuWQGwG6iy1vYZY94G/BCoudyOxpitwFaAysrKJBVPRNJV18AI20928MqJNl460c7xlj4A8jL93LKkiMfuLuZNS4upKgqluKSJl5QAsNb2jFv+qTHmH4wxxdbatsvs+wTwBEBtba1NRvlEJH10D4R59XQH20+288qJdg6f68FayPR72bSokPdsrODWJcWsKsudlc06k5GUADDGLADOW2utMWYz7hlE7ck4toikt47+EV491cGrp1ylP1rhB3weNsbb8W9dWsTainwCvvR6PFqihoF+B7gTKDbGNAB/CfgBrLVfBd4DfMwYEwEGgUestfp1LyIJ19g1yM5THbx6uoO60x28ft416WT4PGyIV/hbFheybmE+Qb83xaVNrUSNAnr0Gtsfxw0TFRFJmGjMcvRcL7vOdLDzdCd1pztoinfa5mT42FBVwDvXl7NlUSFryvPI8KV3hX+pmXNLmojINXQPhtlb38XuM53sPtvJnrNd9A1HAJiXk8GmRYX8UVUBm6oLWVk699vwp0oBICIzUjRmOd7Sx976Tnaf6WJPfSfHWvqwFjwGls3P4Z3ry9hYVUBtVSEVBZlJf5jabKcAEJEZoaVniD31XbxW38Xe+i72NXSP/brPy/SzvjKft9/gKvwbF+aRE5z94/BTTQEgIknXPRBmf2M3rzW4Cn9fQzfnelzbvc9jWFHqft2vX1jAusp8FheH9Ot+GigARGRadQ+EOdDUzf7G+NTQzdmOgbHti4pDbFlcyNqKfNYtzGd1WW7aj85JFgWAiCRMS88QB5t7ONTUw4HGbg40dVPfMTi2vTw/kxvK83hk80LWludzQ3keeVlqykkVBYCITFo0ZjnV1s/h5h4OxSv8Q809tPYOj+1TWZjF2vJ8Ht1cyZqyPNaU51EYCqSw1HIpBYCIXFVn/whHzvVy5FwPR5rd/Oj5XobCMcC12dfMz+H2mhJWleWyuiyXVWW55KqTdsZTAIgIAIMjUY639HH0fC9Hz/Vw9HwfR8/1cL7nwq/6giw/K0tz+cDmKlaV5bKyNIel87J1g9UspQAQSTODI1FOtPZxrKWX18/3cey8m9d3DjD6gJaAz8PSkmxuXVrMigU5LF+Qy8oFOZTkZGg0zhyiABCZozr7RzjR2sfxlvgUX27sGhyr6H0ew6LiEDdU5PHuDRUsm5/N8gU5VBWFdBdtGlAAiMxiw5Eo9R0DnGzt52RbPydb+8aWO/pHxvbL8HlYXJLN+soC3le7kKXzsqmZl011cQi/N72egCkXKABEZrhINEZj1yCn2wc43dbPqfh0ur2f+o4BYuOeq1ucncHikhBvXr2AJSUhlpRks3ReNmX5mfpFL2+gABCZAQZHotR3DnCmfYAz7f2c7biw3NA5SGRcLR8KeKkuDrGmPI+HbixjUUmIRcXZLCoOzYnXFEryKABEkiAaszR3D9LQOUh9xwD1o/OOAc52DNAybvw8uEcZVxVnsbo8j7fdUEp1cYjqohDVRVnqiJWEUQCIJMBIJMa57iEaugZo7ByksWuQxk5X4Td0DdDcNXTRr3hjoDQ3SEVhFncsK6GyMIvKoiwqC7OoKgpRkOVXJS/TTgEgcg3RmKW9b5im7iGauwbH5s3dQzR2DdLUNUhr3zCXvuNuXk4GCwuz2FBZQPnaTBYWZrGwIIuKgkzK8jPT7vWDMvMoACStDYWjnO8Z4lz3EOd6huLLw24e//x8z8W/3gGCfg9l+ZmU5WVy5/ISyvOzKMsPUp6fSXlBJgvygro5SmY8BYDMOdZaeocjtPYO09o7TEvvMC09Q7T2DdPSM0xL7xAtPa6S7xmKvOHPZwW8LMgNMj83yJZFhZTmB1mQl0lpbpDS/CBleZnkq4lG5gAFgMwK0Zilc2CEjv4R2vqGae8bob1vmLY+t97WN0xr3whtvW55OBJ7w3cEvB5KcjKYl+uGSt6ypIh5uUHm5WSwIC/oKv28IDkZPlXukhYUAJJ01lp6hiJ0DYzQORCms39krHLvGgjT3j9CZ79b74h/3jkw8oY2dgCvx1AYClCcnUFxdoAlxSGKczKYl5NBSU4GxdlueV5OkNxMVewi4ykA5LpYaxkMR+kZjNAzFKZ7MEz3QHw+buoZDNM1GKZrYISu+D5dg2GiscvU5rgKvSDLT0FWgMJQgGXzsykMBSgMuQq+MBSgaNxyQVYAj25wErkuCQkAY8yTwANAi7V2zWW2G+BLwNuAAeBD1trdiTi2TM5oxd03HKF/OEr/cIS+4Qh9QxH6RyL0Dl1Y7x0K0zu2HKF3OEzvUISeQTe/tGP0UtkZPvKz/ORl+snP8lOan0l+fLkgK0B+VsBV9vGKvCDLT27QrwpdJEkSdQXwLeBx4NtX2P5WoCY+bQG+Ep/LONGYZTgSZSgce8N8KByNTxeWB0enETcNxJcHRiIMjH4WX++Pr/ePRC7blHIpj3EVeE7QT07QR07Qx7ycIEtKfOTGP8vN9I8t52W6ij53dB704dMzZkRmtIQEgLV2mzGm+iq7PAR821prge3GmHxjTKm1tjkRx7/ULw+dH9fEYONldEtubrEWYvbCPDa2bonGIGotsZglOjrZC8uRmCUaixGJWSJRSyQaIxyzRKOWcCxGePSzaIyRqCUcGV2OMRIZNx+3PByJXbFZ5FqMgUy/l6yAl8yAN77sIyvgJT/LTyjDN7YeyvARCnjJyvCRneElFPCRneEjO+gjlOGLV/o+Mv1etZeLzHHJ6gMoB+rHrTfEP3tDABhjtgJbASorK6/rYJ/8zh4Gw9Hr+rMT5fUYfPHJ6zH4vR58XoPP4yHg8+Aftxzwunl20Ic/vpwRn49uz/B7CHi9BP0eMnweMvyjy14yfB4y/d6xzzL9rqIP+tw8w+dRZS0ik5asALhc7XTZn7vW2ieAJwBqa2uv6yfxDz5+y0XNHKN1ozFgMPE5eDwGjzEYXIVuDHiM+8zjAa9xlfvo5DGuovcYVOGKyKyXrABoABaOW68AmqbrYCtLc6frq0VE5oxk9dI9DfyecW4Cuqer/V9ERCYmUcNAvwPcCRQbYxqAvwT8ANbarwI/xQ0BPY4bBvrhRBxXRESuX6JGAT16je0W+ONEHEtERBJDA7VFRNKUAkBEJE0pAERE0pQCQEQkTSkARETSlAJARCRNKQBERNKUAkBEJE0pAERE0pQCQEQkTSkARETSlAJARCRNKQBERNKUAkBEJE0pAERE0pQCQERkBonGotT31iflWMl6J7CIiFzBcHSYHc07eO7sc7xQ/wI+j49fvueXeMz0/kZXAIiIpED3cDfbGrbxfP3z/LbxtwxGBsn2Z3NbxW3cU3kPMRtTAIiIzBVne87yfP3zvFD/Anta9hC1UeZlzuMdi9/BXZV3sWXBFvxef9LKowAQEZkm4ViYvS172dawjRfqX+B0z2kAagpq+Miaj3DXwrtYXbx62n/pX4kCQEQkgdoH23mp6SW2NWzj5caX6Q334vP42DR/E4+seIQ7Ku6gIqci1cUEEhQAxpi3AF8CvMA3rLWfvWT7h4DPAY3xjx631n4jEccWEUmlaCzK/rb9vNT0Ei82vMjB9oMAFGcWc2/VvdxecTs3ld5EdiA7xSV9oykHgDHGC3wZuA9oAHYaY5621h66ZNfvWmsfm+rxRERS7Xz/eV5uenls6hnpwWM8rClew2PrHuO2ittYUbgiZU07E5WIK4DNwHFr7UkAY8xTwEPApQEgIjIrDUYG2X1+91iFf7zrOOB+5d+58E5uK7+Nm0pvIj+Yn+KSTk4iAqAcGH/XQgOw5TL7vdsYczvwOvCn1trk3OkgIjJJ0ViUwx2H2d68nVeaXmFPyx7CsTB+j58N8zfwjiXv4NayW1lWsAxjTKqLe90SEQCXO3t7yfozwHestcPGmI8C/wTcfdkvM2YrsBWgsrIyAcUTEbk6ay0nu0+yo3kHO5p3sPP8TnpHegFYXrCcR1c8yi1lt7Bh/gYyfZkpLm3iJCIAGoCF49YrgKbxO1hr28etfh34myt9mbX2CeAJgNra2kuDRERkyqy1nO09y85zO3n13Ku82vwq7UOumioLlXFf1X1sWbCFzaWbKc4sTnFpp08iAmAnUGOMWYQb5fMI8IHxOxhjSq21zfHVB4HDCTiuiMiEWGs503OGuvN17Dy3k7rzdbQMtABQklnCltItbCndwqYFm6jIrpjVzTqTMeUAsNZGjDGPAc/ihoE+aa09aIz5DFBnrX0a+KQx5kEgAnQAH5rqcUVEriRmYxzrPMbult3Unatj1/ldY7/wizOLqZ1fy6YFm9i0YBPVudVpU+Ffylg7c1tZamtrbV1dXaqLISIz3HB0mINtB9ndsps9LXvY07JnrA1/ftZ8ahfUsnH+RjbN30RVbtWcrvCNMbustbUT2Vd3AovIrNM22MZrra+xt2Uve1v2crD9IOFYGIBFeYu4v+p+Ns7fyPp56ynPLp/TFf5UKABEZEaLxCK83vk6+1r3jVX6DX0NAPg9flYVreKDKz/I+nnrWTdvHYXBwhSXePZQAIjIjGGt5fzAefa37Wd/6372te3jUPshBiODABQFi1g3bx3vX/5+bpx3I6uKVpHhzUhxqWcvBYCIpEz3cDcH2w9ysO0g+9v2c6DtAK2DrYD7db+ycCUP1zzMjSU3srZkLWWhMjXnJJACQESSom+kj8MdhznUfoiDbQc52H6Qs71nx7ZX5VaxpXQLa4rXsLZ4LcsLlxPwBlJY4rlPASAiCdc93M2RjiMcbj/MoY5DHG4/PPYsfIDSUCmri1bzzqXvZE3xGlYVrSIvIy91BU5TCgARuW6jbfZHOo5cNDX2NY7tUxoqZWXhSh5Y/ACrilaxqmgVRZlFKSy1jFIAiMiEDEeHOdF1gtc7X+dox1Fe73ydIx1H6BnpAcBgqMqt4obiG3jvsveysnAlK4tWUhAsSHHJ5UoUACJykZiN0djbyLGuYxzrPDY2P9NzhqiNAhD0BqkpqOG+qvtYXriclYUrWVawjCx/VopLL5OhABBJU6PNN8e7jnOi6wTHOo9xousEJ7pPjA27BCjPLqemoIZ7q+5lecFylhUsY2HOQrwebwpLL4mgABCZ46KxKE19TZzsPsnJ7pOc6Doxttwf7h/bryhYxNL8pby75t0szV/Kkvwl1BTUEPKHUlh6mU4KAJE5oj/cz+nu05zqOeXm3ac41XOKM91nGImNjO1XnFnM4rzFPLjkQZbmL2VR3iKW5i9VW30aUgCIzCLhaJj6vnrO9pzlTM8ZTvec5kzPGc50n6FlsGVsP4/xUJFdQXVeNbeW3cqivEUszlvMorxFGm4pYxQAIjPMUGSIxr5G6ntdRX+29+zYvLm/mZiNje2bn5FPVW4VN5XdxKK8RVTnVlOdW01lbqVuopJrUgCIJJm1lvahdhp6G2joa3Dz+HJ9b/3Yi0pG5fhzqMytZG3JWh5Y/ABVuVVjk37Ny1QoAEQSzFpLx1AHzf3NNPY10tzn5o19jTT1NdHU33TRKBuAeZnzqMip4KbSm1iYs3BsqsypJC8jT8+/kWmhABCZpOHoMOf7z3Ou/xznBs7R3NdMc7+bmvqaONd/jqHo0EV/JsefQ3lOuWuTL7+V8uxyKnIqKM8upzy7nKAvmKKzkXSmABAZZzAySMtACy0DLZwfOD9W0Z8fuDDvGOp4w58rChZRGiqlpqCG2ytupyy7jNJQKeXZ5ZRml5IbyE3B2YhcnQJA0sJwdJi2wTZaB1ppHWwdm7cMtNA64OYtgy1jrxEcLyeQw/ys+cwPzWdV0SoWhBawILSA0lApC0ILmJ81X7/gZVZSAMisFY6F6RjsoH2onfbB9ovmbYNttA+6eetg62Urdp/xUZxVTElmCdV51WxasIn5ofnMy5rH/KwLcz3eQOYqBYDMGCPRETqHOuka7qJjqIPOoU46hzvpHOocW+8Y6hibRh9CdqlMXybFmcVjNzxtXrCZkqwSSjJLKM4spiTLzQuDhXiMJ8lnKTJzKAAk4cKxMD3DPfSM9NA93P2GeddwF13DXReWh9z6QGTgst9nMORn5FMYLKQgWMCygmUUBgspDBZSlFlEUbDoorl+sYtMTEICwBjzFuBLgBf4hrX2s5dszwC+DWwE2oH3W2tPJ+LYkljWWgYjg/SF++gb6Rub94Z7x9Z7R3rpC/fRM9xD70gvPSM99IZ7xyr9S4c4jmcwZAeyyc/IH6vUl+QtIS8jj/yMfAqCBRQGC8nLyBur8PMCeXrwmMg0mHIAGGO8wJeB+4AGYKcx5mlr7aFxu/0B0GmtXWqMeQT4G+D9Uz12ugvHwgxFhsamgcgAg5HBsWl0fSA8cPFyeID+cD/9kf4Ly+H4cqT/ojtNL8dgyPZnk5uRS04gh2x/NguzF5Jb5NZzA7luysglPyN/bD0/I5+cQI4qc5EZIhFXAJuB49bakwDGmKeAh4DxAfAQ8Ffx5e8DjxtjjLXWJuD4b2QtJODGmZiNEbVRorHoRcsRGyEaixK1USKxCJFYhHAsPLYejoXH5uFo2M1jYUaiIxctj0RHGImNEI6GGY4OMxJznw1Fhtw8emE+HBl28+jwWIUfsZFJnU/AEyDkD5HlzyLTl0nIHyInkMOC0AJC/pDb5ssiJ5Aztu2iuT+H7EA2IX9Ibecic0AiAqAcqB+33gBsudI+1tqIMaYbKALaEnD8i1nLW/9xLUMeD3j9WI8P6/GB8WABO/o/G5+wxGxsbD46jb74YroZDBneDALeAAFvYGw5w5tBhjeDoDdIyB8i6AsS9AbJ8LnPgr4gGd4MMn2ZZPoyCfqCF5a9wbFKfrSiz/Rl4vOoy0dELkhEjXC5n9qX/rKfyD5uR2O2AlsBKisrJ1+aaJjbsqsJ97dCTzvEwhjAeHyYrCIIlWBCxZjQPDyZ+RjjxRiDNz734MFj3OQ1Xjf3eMeWfR4fPuPD4/HgMz63PjrF1/1e/9hywBvA7/FftBzwuMp+dD/d5i8iqZCIAGgAFo5brwCarrBPgzHGB+QBb7ydErDWPgE8AVBbWzv5JiJfgP/yvmfcciwKbcegeS807YWmPXBsF4Tjo038ISi9EcrWQdl6NxUuAY+aN0Rk7ktEAOwEaowxi4BG4BHgA5fs8zTw+8ArwHuAX09b+/94Hi/MW+GmGx9xn8Wi0Pa6C4PRUKh7EiLxZ7cEcuKBMC4UChYlpE9BRGQmmXIAxNv0HwOexQ0DfdJae9AY8xmgzlr7NPBN4J+NMcdxv/wfmepxr5vHC/NWumldPKeiEWg9Eg+F3S4YdnwNovG3KAXzLoTB6JS3UKEgIrOaScYP8etVW1tr6+rqUnPwyAi0Ho6Hwh5o3A0thyAWH3mTVRQPgw0XQiG3NDVlFRGJM8bsstbWTmRfDQu5El/A9Q+U3ggbP+Q+Cw9By0EXBqP9Ci9+AUZHDOWUjrtKiAdDqChlpyAicjUKgMnwB6F8o5tGjQzAuX0X+hOadsPRnzE2yCm/8kIglG+A0nUQ1KOBRST1FABTFciCypvcNGqoB5pfc2HQuNvND/0ovtFAcc2FQCjbAAtucOEiIpJECoDpEMyFRbe5aVR/+4UrhMbdcPJ52PeU2+bxwbxV8auLDW5essJ1WIuITBMFQLKEiqDmXjeBe1xFT9OFQGjcBQd+ALv+0W0fvUdhNBDKN7rmJI08EpEEUQCkijGQV+6mle9wn8Vi0HHShUHjLhcOr34doo+77VnFF8KgYqNrPsoqTN05iMispgCYSTweKF7qphvjD0uNjMD5AxeuFBrq4NgvGOtkLlwSD4RaN19wA/gyUnYKIjJ7KABmOl8g3gy0ATbFPxvqccNQG+rclcKpbbD/e26bN+BCoLwWKja5YCioVtORiLyBAmA2CubCotvdBPH+hEYXBqOhsOef4dWvue1ZRRcHQvlGDUUVEQXAnGAM5FW4adVD7rNoxN253FgHDbugYScce3b0D7hRRgs3xUNhExQv10PwRNKMHgWRTga7XF9CQx3Uv+pCYajLbcvIdVcGCzdDxWbXyZxZkNryisik6VEQcnmZ+bDkbjeBazpqPwEN8TCo3wnbPgejr4QsXu6uEhZucVNRja4SROYQBUA6M+bCqKPRJ6MO98ZHG73qrhKO/AT2/IvbFsx3VwgLN7tAKN8IgVDqyi8iU6IAkItl5MDiO9wE8auE41C/wwVC/Y74MFTAeN2Io8qbXCBU3gS5Zakru4hMivoAZPIGO10/wtntLhAax71lLa8y/mykLVB5M5SsVLORSBKpD0CmV2YB1NznJoBo2D0R9ewOqN8Op35z4b6EYN6Fq4PKm93dy3rwnciMoACQqfP6Lzyi4uaPu2ajztPxK4Ttbj7abOQNuBCouhkqb3FXCsG8lBZfJF2pCUiSo789HgavwJlX3J3MsQhgXD9C1S1uqrwFsktSXVqRWWsyTUAKAEmNkX439PTMK3D2ZTcENTLothUvd2FQ/SaoulWv2hSZBPUByMwXCMHiO90E7qF3TXtcGJx5GQ7824VHYxcudkFQfRtU3+rueBaRKdMVgMxMsajrWD79Epz+rQuGoW63rWCRuzqojr90R0NPRcaoCUjmnlgUzh90YXD6t3DmtxcCoXCJC4Lq29wD8rLnpbasIimUtAAwxhQC3wWqgdPA+6y1nZfZLwrsj6+etdY+OJHvVwDIFcWi7j0Jp3/rHod95mUY7nHbSlZceFpq9Zv0TCNJK8kMgL8FOqy1nzXGfAoosNb+p8vs12etzZ7s9ysAZMKiETj3Gpx60d2HcHZ7/OY0416tufgOWHSHuxchkJXq0opMm2QGwFHgTmttszGmFHjBWrv8MvspACS5IiPuUdgnf+OuEBp2Qizs7kNYuCXeAX0XlK0DjzfVpRVJmGQGQJe1Nn/ceqe19g3X28aYCLAXiACftdb+cCLfrwCQhBnuc/cgnHzBhcL5eItkMM81FS2+C5bc5UYcicxiCR0Gaoz5FbDgMps+PYkyVVprm4wxi4FfG2P2W2tPXOF4W4GtAJWVlZM4hMhVZGRf/PiKvlbXVHTyBTcdfsZ9XlDtwmDpPa5TOTP/Cl8oMvslpQnokj/zLeDH1trvX+v7dQUgSTH6XoSTz8OJX7t+hJFe97TTilpYco8LhLL1ai6SGS+ZTUCfA9rHdQIXWmv/4pJ9CoABa+2wMaYYeAV4yFp76FrfrwCQlIiGXZ/B8efgxHPQtBewbjTR4rtg6b0uEHIud2EsklrJDIAi4HtAJXAWeK+1tsMYUwt81Fr7h8aYW4CvATHAA3zRWvvNiXy/AkBmhP52d3Vw/Dk4/ivob3Gfz78Bau6FmvvdazS9urFeUk83golMl1jM3X9w/Fduqt/hHmqXkec6kWvud1cIOfNTXVJJU3oWkMh08XigdK2bbvszdzfyyd+4x10f+yUcig9wK10Hy97sptL1eimOzEi6AhBJFGvh3H449iy8/gt3H4KNQWieuzJY9mZ3lZCRk+qSyhymJiCRmaC/3TUTHXvWzYe63Y1o1W+C5W+DZW+B/IWpLqXMMQoAkZkmGnaPp3j953D0Z9ARvw1m/g2w/K1uKlsPxqS2nDLrKQBEZrq2Yy4Ijv7MvSnNxiCn1F0ZrHi7uwnNF0h1KWUWUgCIzCb97a4T+ciP3Y1o4QHIyHX9Bive7u5eVr+BTJBGAYnMJqEiWPeom8KDblTRkR+7q4MD33f9BovvgpUPuCuEUHGqSyxzhAJAZCbxZ8Lyt7gpFnX3GRx+Bg7/2HUmmz9xr8dc+aALBL0NTaZATUAis4G17hWZh5+BQ09D21H3+cItsOohFwgaUSSoD0Bk7ms96oLg8I/cvQcA5Rth1TtdIBRUpbZ8kjIKAJF00n4CDj8NB38IzXvdZ+UbYfW7XCDoyiCtKABE0lXHKfc4ioP/D5pfc59VbIY174bV79QTTNOAAkBE3JXBwf/npvMHAOPuQl7zbtdMlFWY6hLKNFAAiMjFWo/WMiZdAAAMJElEQVTCgR+4YaXtx8Hjcy+6ueG97i7kjEm/sltmKAWAiFze6Gii/d93gdDTAP4sd3/B2vfBkrvB6091KWUKFAAicm2xmHsMxb7vuX6DwU7IKoLVD8ONj7iOZD2baNZRAIjI5ERG3BNL933X3YEcHYbCJS4I1r4PCqpTXUKZIAWAiFy/oW53j8G+78LpF91nlbe4MFj9LgjmprZ8clUKABFJjK6zronotaeg/Rj4Mt0jKNZ9ABbdAR5vqksol1AAiEhiWQuNu2Dv/4ED/wZDXZBbDjc+6sKgaEmqSyhxCgARmT6RYTj6U9jzr3DiOfcug6pbYf3vuPsLAqFUlzCtKQBEJDl6mlzz0J5/cW85C+TADe+G9b8H5Rs0iigFJhMAnike6L3GmIPGmJgx5ooHNMa8xRhz1Bhz3BjzqakcU0RmkNwyuO3P4BO74MM/g5XvgNe+C9+4G75yK+z4mhteKjPSlAIAOAA8DGy70g7GGC/wZeCtwCrgUWPMqikeV0RmEmOg6hZ411fgz4/CA3/nbij72V/AF1bAD7bCmVdcX4LMGFN6IYy19jCAufpl3mbguLX2ZHzfp4CHgENTObaIzFDBPKj9iJuaX4Pd33YjifZ9F0pWwMYPuSGlmQWpLmnam+oVwESUA/Xj1hvin12WMWarMabOGFPX2to67YUTkWlUeiO8/QvwH47Ag4+7DuKff8pdFfzw49BQp6uCFLrmFYAx5lfA5Z4h+2lr7Y8mcIzLXR5c8W/cWvsE8AS4TuAJfL+IzHSBEGz4XTc174Nd/+iuCvb+qwuJTX8Ia94DgaxUlzStXDMArLX3TvEYDcD4N1JUAE1T/E4Rma1K17o+gvs+45qFdj4JT38Cnv2vsP6DLgx0X0FSJKMJaCdQY4xZZIwJAI8ATyfhuCIyk2XkuMr+Yy/Bh38ONffBq1+Hv98A//wwHP25e2CdTJupDgN9lzGmAbgZ+Ikx5tn452XGmJ8CWGsjwGPAs8Bh4HvW2oNTK7aIzBnGQNXN8J5vwp8ehLs+DS2H4Dvvh79fDy8/DoNdqS7lnKQbwURk5omG4fAz8OoTcPYV8Idg3aOw+d9BybJUl25Gm8yNYFMaBioiMi28fljzsJuaX3M3lO3+Nuz8Biy9F276uHt5je40nhJdAYjI7NDX6kYPvfp16G9x9xTc9HFY+37wB1NduhkjaY+CEBFJmuwSuOMv4E8PwDu/Ch4/PPNJ+LvV8Pz/cgEhk6IAEJHZxZfh+gM++iL8/jNQUQu/+Sx8cQ088yfQdizVJZw11AcgIrOTMbDodje1vg7bvwx7vwO7/sm95P7WP4HKLaku5YymKwARmf1KlsE7vuSGkd7+H+Hsy/Dk/fDN++HIT3U/wRUoAERk7sgugbs/7YLgrZ+D3mZ46lH4ys3ubWbRcKpLOKMoAERk7gmEYMtW+MQeePgb4PHBDz8G/3u9G1I6MpDqEs4ICgARmbu8Plj7Xvjob+GD34e8he4dBV+8AbZ9Hoa6U13ClFIAiMjcZ4x71tBHfuaeO1S2Hn791/B3N8Bzfw397akuYUooAEQkvVTdDL/zfdj6Aiy+HV78vLsi+MV/g76WVJcuqRQAIpKeytbD+/8FPr4dlr8VXnkcvrgWfv5foPd8qkuXFAoAEUlv81a6J5H+8U5Y/S7Y8VX40lr42aeg91yqSzetFAAiIgDFS91L7R/bCWve7Z5E+qUb4ef/ec5eESgARETGK1oC7/wH+ESdC4IdX3VB8Oynob8t1aVLKAWAiMjlFC52QfBYHax6CLb/g+sjeO4zMNiZ6tIlhAJARORqipbAw19zncXL7ocXvwBfvBG2fQ6G+1JduilRAIiITETJcnjvt+CjL0H1rfDr/+Gahl75BwgPpbp010UBICIyGQvWwKPfgT/4FcxfDc/+Z3i8Fvb8C0QjqS7dpCgARESux8JN8PtPw+/+EELF8KM/hq/cAkd+AjP4TYvjKQBERKZiyV3wR8/D+74NNgpPfQCefDOc3Z7qkl3TlALAGPNeY8xBY0zMGHPFd1AaY04bY/YbY/YaY/SSXxGZW4xxI4U+vgMe+CJ0nnEh8NQHZ/QbyqZ6BXAAeBjYNoF977LWrpvoy4pFRGYdrw9qPwyf3A13/1c4+Rv48hb48Z/NyHcWTykArLWHrbVHE1UYEZE5IRBybyb75B6o/Qjs/if3LoJtn4fwYKpLNyZZfQAW+IUxZpcxZmuSjikiklrZJfD2z7t7CBbf4R5B/fcbYd/3ZsRrKq8ZAMaYXxljDlxmemgSx7nVWrsBeCvwx8aY269yvK3GmDpjTF1r68y7ZBIRmbTiGnjkX+FDP4VQCfzgj+Cb98LZHSktlrEJGK5kjHkB+HNr7TU7eI0xfwX0WWs/f619a2trbV2d+oxFZA6JxWDfU+6REr3N7nlD9/53yF+YkK83xuyaaF/rtDcBGWNCxpic0WXgflznsYhI+vF4YN0H4BO74Pa/cPcNPL4Jnv+fSX9X8VSHgb7LGNMA3Az8xBjzbPzzMmPMT+O7zQd+a4x5DXgV+Im19udTOa6IyKwXCMHdn3YPm1vxNvjN37ggOPBvSbuRLCFNQNNFTUAikjbOvOxeWH9uP1Td6l5iH8ia9NdMpgnIN+lvFxGRxKu6Bbb+xg0ZbdpzXZX/ZCkARERmCo/X3TeQrMMl7UgiIjKjKABERNKUAkBEJE0pAERE0pQCQEQkTSkARETSlAJARCRNKQBERNLUjH4UhDGmFThznX+8GGhLYHFmA53z3Jdu5ws658mqstaWTGTHGR0AU2GMqUu310/qnOe+dDtf0DlPJzUBiYikKQWAiEiamssB8ESqC5ACOue5L93OF3TO02bO9gGIiMjVzeUrABERuYpZHwDGmLcYY44aY44bYz51me0ZxpjvxrfvMMZUJ7+UiTOB8/0zY8whY8w+Y8xzxpiqVJQzka51zuP2e48xxhpjZv2IkYmcszHmffG/64PGmP+T7DIm2gT+bVcaY543xuyJ//t+WyrKmSjGmCeNMS3GmMu+I904/zv+/8c+Y8yGhBfCWjtrJ8ALnAAWAwHgNWDVJft8HPhqfPkR4LupLvc0n+9dQFZ8+WOz+Xwnes7x/XKAbcB2oDbV5U7C33MNsAcoiK/PS3W5k3DOTwAfiy+vAk6nutxTPOfbgQ3AgStsfxvwM8AANwE7El2G2X4FsBk4bq09aa0dAZ4CHrpkn4eAf4ovfx+4xxhjkljGRLrm+Vprn7fWDsRXtwMVSS5jok3k7xjgr4G/BYaSWbhpMpFz/iPgy9baTgBrbUuSy5hoEzlnC+TGl/OApiSWL+GstduAjqvs8hDwbetsB/KNMaWJLMNsD4ByoH7cekP8s8vuY62NAN1AUVJKl3gTOd/x/gD3C2I2u+Y5G2PWAwuttT9OZsGm0UT+npcBy4wxLxljthtj3pK00k2PiZzzXwG/Y4xpAH4KfCI5RUuZyf73Pmmz/Z3Al/slf+mwponsM1tM+FyMMb8D1AJ3TGuJpt9Vz9kY4wH+DvhQsgqUBBP5e/bhmoHuxF3lvWiMWWOt7Zrmsk2XiZzzo8C3rLVfMMbcDPxz/Jxj01+8lJj2umu2XwE0AAvHrVfwxsvCsX2MMT7cpePVLrtmsomcL8aYe4FPAw9aa4eTVLbpcq1zzgHWAC8YY07j2kqfnuUdwRP9d/0ja23YWnsKOIoLhNlqIuf8B8D3AKy1rwBB3DNz5qoJ/fc+FbM9AHYCNcaYRcaYAK6T9+lL9nka+P348nuAX9t4D8ssdM3zjTeHfA1X+c/2dmG4xjlba7uttcXW2mprbTWu3+NBa21daoqbEBP5d/1DXIc/xphiXJPQyaSWMrEmcs5ngXsAjDErcQHQmtRSJtfTwO/FRwPdBHRba5sTeYBZ3QRkrY0YYx4DnsWNInjSWnvQGPMZoM5a+zTwTdyl4nHcL/9HUlfiqZng+X4OyAb+b7yv+6y19sGUFXqKJnjOc8oEz/lZ4H5jzCEgCvxHa2176ko9NRM85/8AfN0Y86e4ppAPzeIfcxhjvoNrwiuO92v8JeAHsNZ+FdfP8TbgODAAfDjhZZjF//+JiMgUzPYmIBERuU4KABGRNKUAEBFJUwoAEZE0pQAQEUlTCgARkTSlABARSVMKABGRNPX/AYPj4ORseBBkAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1a1dbf7c50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.plot(np.linspace(0,1,101),sol[:,0])\n",
    "plt.plot(np.linspace(0,1,101),sol[:,1])\n",
    "plt.plot(np.linspace(0,1,101),sol[:,2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def inner2_product(h1, h2):\n",
    "    return lambda x: H.hermeval(x, H.hermemul(h1, h2))     \n",
    "\n",
    "def inner3_product(h1,h2,h3):\n",
    "    return lambda x: H.hermeval(x, H.hermemul(H.hermemul(h1,h2),h3)) \n",
    "\n",
    "def integ1_func(h1,h2):\n",
    "    return lambda x:inner2_product(h1,h2)(x)*(norm.pdf(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def exact_mean(t):\n",
    "    f=np.exp(t**2/2)\n",
    "    return f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1a18232da0>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEKCAYAAAAFJbKyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xl4XXW97/H3N3MztGmThs50oC10CAXSiXKgyIxaBvGIolAUe1A4cK/CEY9HUa8i9xzUR0CGipzqARRlrMgVZKjM0rSWQuncBhpa2jRt0zRzdr73j70JaZph7yZ7r53m83qe9ew1/PZe35Wk+9M1/Za5OyIiItFKCboAERHpWxQcIiISEwWHiIjERMEhIiIxUXCIiEhMFBwiIhKTwILDzEab2YtmttbM1pjZ9R20MTO73cw2mdlqMzsxiFpFRORjaQGuuxn4pruvNLM8YIWZ/dXd323T5jxgYmSYDdwdeRURkYAEtsfh7jvcfWVkvBpYC4xs1+wC4Lce9gaQb2bDE1yqiIi0EeQeRyszGwucAPy93aKRwLY20+WReTs6+IxFwCKAnJyck4499th4lCoickRasWLFbncfGk3bwIPDzHKBR4H/5e772y/u4C0d9pHi7ouBxQAlJSVeWlraq3WKiBzJzOy9aNsGelWVmaUTDo0H3f2xDpqUA6PbTI8CtieiNhER6ViQV1UZ8Gtgrbv/rJNmS4HLI1dXzQGq3P2Qw1QiIpI4QR6qmgd8CXjbzFZF5v07MAbA3e8BngbOBzYBtcCVAdQpIiJtBBYc7v4KHZ/DaNvGgWsSU5GIiERDd46LiEhMFBwiIhITBYeIiMREwSEiIjFRcIiISEwUHCIiEhMFh4iIxETBISIiMVFwiIhITBQcIiISEwWHiIjERMEhIiIxUXCIiBwB3qus4cX1uxKyrsCfACgiIrFraA6xfOteXly/ixfX72JLRQ15WWn847tnkZYa330CBYeISB/xfmUtf9tYwd/W7+K1zZXUNobISE1hzoQCLp9zNKcfWxT30AAFh4hI0qppaOaNLZW8vHE3f9tQwdbdNQCMGjyAi08cyfxJRZx8TAHZGYn9KldwiIgkiVCL8/YHVbyysYKXN+5m5ft7aQo5WekpzBlfwOVzj+bUSUMZX5hD+OnbwVBwiIgExN3ZXFHDa5t38+qm3by+uZL99c0ATB0xkC+fMo7TJg7lpLGDyUxLDbjajwUaHGZ2P/ApYJe7T+tg+XzgSWBrZNZj7v7DxFUoItK7tu2p5fUtlby+uZLXNu9m5/4GAEbmD+C8acOZN7GQeRMKKMjNDLjSzgW9x7EEuBP4bRdtXnb3TyWmHBGR3rVtTy1/37qHNyJh8cG+OgAKcjKYO6GAuRMKOOWYQsYMyQ708FMsAg0Od3/JzMYGWYOISG9xd8oqa3lzayV/37qHv2/Z0xoU+dnpzB43hK/+0zjmTihk0lG5fSYo2gt6jyMac83sLWA7cIO7rwm6IBERgOZQC+s+rGZ52R6Wl+3hza172X0gfOipICeDmWPDQTFnQgGTivJISembQdFesgfHSuBodz9gZucDTwATO2poZouARQBjxoxJXIUi0m8caGhm1fv7KH1vDyve28vK9/ZS0xgCwucoTjmmgJnjhjB73BAmDO27exTdSergcPf9bcafNrO7zKzQ3Xd30HYxsBigpKTEE1imiByB3J33KmtZ+f7ecEi8v4/1H+6nxcEMJh+Vx8UnjqJk7GBKxg5hZP6AoEtOmKQODjMbBux0dzezWYT71qoMuCwROQJV1TXx1rZ9rIoM/3h/L3trmwDIy0xjxph8zvrEREqOHswJY/LJy0oPuOLgBH057u+A+UChmZUDNwPpAO5+D3AJ8DUzawbqgEvdXXsTItIj9U0h1n1YzVvb9oXDonwfWyrCd2WbwTFDczlrylGcMCYcEhOL8kg9Qs5P9Iagr6r6fDfL7yR8ua6IyGFpCrWwYWc1b5dXsfqDKlaX72P9h9U0hcL/By3MzWTG6HwuPmEkM0YPpnj0IAb2472JaCT1oSoRkVg0NIfY8OEB3tlexdsfVLHmgyrWflhNY3MLAHlZaRSPGsRV/zSe4pGDKB6dz4hBWUfsSex4UXCISJ+0v76Jtdv3syYyvLtjPxt3VtPcEt6TGJiVxrSRg1h48limjxxE8ahBfeomu2Sm4BCRpNbS4mzbW8vaHftZu6OatTvCIVG+t661TWFuJlNHDOT0yUOZNnIQU0cMVEjEkYJDRJLG3ppG1n1YzfoP97N+Z3VkvJrayL0SKQZjC3OYMTqfz88aw5ThA5k6YiBFA7MCrrx/UXCISMLtr29i485qNu48wIadB9iws5r1O6upqG5obZOfnc6xw/L455LRHDssj+OGD2TSUXkMyEieXmL7KwWHiMTNnppGNu06wKZdB9i4qzr8uvMAH+6vb22TlZ7CpKPyOG3SUCYflcfkYXkcOyyPoXmZOtSUpBQcItIjzaEWtu2tY0vFATZXHGBLRQ2bKw6wuaKGPTWNre0GpKdyTFEuJ08o4JijcplUlMeko/IYNXjAEdOHU3+h4BCRbrk7FdUNbNldQ9nuGrburmFzRQ1bdx/g/T21rfdEQLhzvwlDczln6lFMGJrLMUXhYcQgBcSRQsEhIkD46qVd1Q28V1nDe5W1lFXWUFZZw9bdtbxXWdN6ghogIzWFsYXZHFOUy1lThjF+aA4ThuYyYWgO+dkZAW6FJIKCQ6QfqWsMUb63lm17a3m/spb39tSybU8t71XW8v6eWhoiN8oBpKUYo4dkM7YgmznjhzC2IIdxheFhRP4AdcHRjyk4RI4gDc0hduyrp3xvHeV7aynfW8e2veFw2La37qCrliB83uHogmzGFuZw2qShHF2Yw9FDshlXmMPwQVmkpaYEtCWSzBQcIn2Eu7O/rpkP9tWxfV8dO6rqKN9Xxwd7w9Mf7KtjV3UDbbsBTU0xhg/KYvTgbE6fPJTRg7MZU5DNqMHZjB4ygKG5unJJYqfgEEkC7s7++mZ27q9nR1U9H1bVsX1fPTuq6thRVR8JivqDzjMApKcaI/IHMDJ/AKdOHMqowdmMHByeHj1kAMMGaq9Bep+CQyTO6ptCVFQ3sKu6np37G9i5v54P99eza38DH1bVt063DwWzcFcaIwZlMbEoj1MnDWXEoAGMyB/AiPwsRuYPoDA3U1cqScIpOEQOQ0uLU1XXRMWBBnZXN1BxoCESDuHXtkFRVdd0yPszUlMoGpjJUQOzOG74QOZPLmLYoEyGDxrA8EFZHDUwPGSkaW9Bko+CQySirjFEZU0De2oaqaxppPJAI5UHGqisaWT3gQZ2R6Z3H2ig8kBjay+sbWWkpVCUl8nQvEzGFeYwe1wBRw3MpCgvqzUoivIyGZKToXML0mcpOOSIVNcYYl9dI3trmthX28je2ib21jayt6aRPbWN7KttorImMh0Z6ppCHX5WZloKhbmZFORmcNTALKaOGEhhbiaFueGAKMzNpGhgeDwvM02BIEc8BYckrfqmEPvrm9hf18z++iaq6prYXxd+raqNvLYb9kZCoe39CO3lZqaRn51OQU4GBbkZTCzKZUhOBkNyMyjIyWBITjgkwuMZ5CoMRA6i4JBe5e40NLdQ09BMTUOIAw3N1DQ2h18bmjlQHx6vrv9oaAq/NjQdNG9/XTONoc6//CF8D8KgAemtw5gh2RSPGkR+dgb52enkD8hgcHY6+dkZDM5JZ3BkfmaaelcV6YlAg8PM7gc+Bexy92kdLDfgF8D5QC2w0N1XJrbKI4u709zi1DeFqG9qob4pRENziLrGFuqaQtQ3hVpfaxvDQ3i8mdrGEHWNH80PB0NtU4iahmZqG5qpaQyPd3TsvyO5mWnkZqaRlxUe8rMzGDMkm7ysdAZmpTFwQHp4+Gg8K51BA9IYNCCDgQPSFAAiAQl6j2MJcCfw206WnwdMjAyzgbsjr3FR3+4Ytzs4HnkNf+lCZLwFWtxxwq8t7rRE5oVavN0rhFrC0yF3Qi0tNIfC000t4emmkNMccppbx1toagm/NoecxlALTa2D09jcQkNzeLqxuSUyHaIx1EJDU0vra0NziIbmjwIi/Brl9/pBUlOM7IzUyJDWOj5oQDoj87PIzkgjJyOVnMy08BAZz41M52alkRcZz8tKIycjTZeRivRRgQaHu79kZmO7aHIB8FsPf2O/YWb5Zjbc3Xd0+cHr18P8+V2v/FOfghtuCI/Pnw8LF3LC5mFkVe3h7id+0m3tj0w/k0emn8ng2irufuIn/GrWRTx/zGzGV5ZzyzN3dvv+9u3/89QrWDnqOE4sX8u/vfSbTt9nZhhwx5lXsnbcdE76YC1fe/bX3PHPN7B7xFhOfvc1Lnjh96SYYQYpZqS0vobHzYxXb/wRjRMmMu71Fzj2oft45xf3kVZUxMgnf0/RY79vbZ+a8vH7Dv4BPAKFhbBkCfxqCSxbFp5/223w1FPdbv9B7V9/HR59NDz97W+Hp7tSUHBw+8pKWLw4PL1oEWzY0PX7J006uH1BAfwk8jv/zGfCn9eVuXMPbj937sF/S93p4G+PhQth92645JLu39++/Te/CZ/+dPjv/l/+pfv3t29/yy1w8snw2mvw7//e/fvbt7/3Xpg8Gf70J/jpT7t/f/v2bf+Wlizp/v3t2+tvLzzd07+9GAS9x9GdkcC2NtPlkXmHBIeZLQIWARRnZh7Wyr5x1iRS9+5m9EvZH3/uQSPWOr3g+BFM/eQUsvbtYewrOXxx9tGccep08rdlM35lLkb4Bi4If4FHRiPzjRvOmsxXz5pDbtkmxr4ziP/6bDHNs+eSsyKTwq1PhgPCwmsMf9l/HBoA//OV2ZF/vNmw9lF+vXBm5B/jXlg7sNtt/dzMMTB5LOweCrkZzJ9cFP7HWJgLGcn+ZyEiQTL3wzhu0ZsFhPc4nurkHMefgZ+4+yuR6eeBf3P3FV19ZklJiZeWlsahWhGRI5OZrXD3kmjaJvttqeXA6DbTo4DtAdUiIiIkf3AsBS63sDlAVbfnN0REJK6Cvhz3d8B8oNDMyoGbgXQAd78HeJrwpbibCF+Oe2UwlYqIyEeCvqrq890sd+CaBJUjIiJRSPZDVSIikmQUHCIiEhMFh4iIxETBISIiMVFwiIhITBQcIiISEwWHiIjERMEhIiIxUXCIiEhMFBwiIhITBYeIiMREwSEiIjFRcIiISEwUHCIiEhMFh4iIxETBISIiMVFwiIhITBQcIiISk0CDw8zONbP1ZrbJzG7qYPlCM6sws1WR4aog6hQRkY8F9sxxM0sFfgmcBZQDy81sqbu/267pw+5+bcILFBGRDgW5xzEL2OTuW9y9Efg9cEGA9YiISBSCDI6RwLY20+WRee19xsxWm9kjZja6sw8zs0VmVmpmpRUVFb1dq4iIRAQZHNbBPG83/SdgrLsXA88Bv+nsw9x9sbuXuHvJ0KFDe7FMERFpK8jgKAfa7kGMAra3beDule7eEJn8FXBSgmoTEZFOBBkcy4GJZjbOzDKAS4GlbRuY2fA2kwuAtQmsT0REOhDYVVXu3mxm1wLPAKnA/e6+xsx+CJS6+1LgOjNbADQDe4CFQdUrIiJh5t7+tELfV1JS4qWlpUGXISLSZ5jZCncviaat7hwXEZGYKDhERCQmCg4REYmJgkNERGIS1VVVZjYJuBE4uu173P0TcapLRESSVLSX4/4RuIfwTXih+JUjIiLJLtrgaHb3u+NaiYiI9AnRnuP4k5l93cyGm9mQj4a4ViYiIkkp2j2OKyKvN7aZ58D43i1HRESSXVTB4e7j4l2IiIj0DdFeVZUOfA04NTJrGXCvuzfFqS4REUlS0R6quhtIB+6KTH8pMk/PABcR6WeiDY6Z7n58m+kXzOyteBQkIiLJLdqrqkJmNuGjCTMbj+7nEBHpl6Ld47gReNHMthB+5OvRwJVxq0pERJJWtFdVPW9mE4HJhINjXZtHuoqISD/SZXCY2Sfc/QUzu7jdoglmhrs/FsfaREQkCXW3x3Ea8ALw6Q6WOaDgEBHpZ7oMDne/OTL6Q3ff2naZmfX4pkAzOxf4BeFnjt/n7re2W54J/BY4CagEPufuZT1dr4iIHL5or6p6tIN5j/RkxWaWCvwSOA+YAnzezKa0a/YVYK+7HwP8HPi/PVmniIj0XHfnOI4FpgKD2p3nGAhk9XDds4BN7r4lsq7fAxcA77ZpcwHw/cj4I8CdZmbu7j1ct4iIHKbuznFMBj4F5HPweY5q4Ks9XPdIYFub6XJgdmdt3L3ZzKqAAmB3+w8zs0XAIoAxY8b0sDQREelMd+c4ngSeNLO57v56L6/bOlrlYbQJz3RfDCwGKCkp0R6JiEicRHsD4CIzO2QPw92/3IN1lwOj20yPArZ30qbczNKAQcCeHqxTRER6KNrgeKrNeBZwEYd+ycdqOTAxcnXWB8ClwBfatVlK+FkgrwOXAC/o/IaISLCivXP8oKuqzOx3wHM9WXHknMW1wDOEL8e9393XmNkPgVJ3Xwr8GvgfM9tEeE/j0p6sU0REei7aPY72JgI9PgPt7k8DT7eb97024/XAZ3u6HhER6T3RPsipmvBJaYu8fgh8K451iYhIkor2UFVevAsREZG+obsbAE/sarm7r+zdckREJNl1t8fx0y6WOfCJXqxFRET6gO5uADw9UYWIiEjfEO3J8XTga8CpkVnLgHvdvSlOdYmISJKK9nLcu4F04K7I9Jci866KR1EiIpK8og2Ome5+fJvpF8zsrXgUJCIiyS3a53GEzGzCRxNmNh4IxackERFJZtHucdwIvGhmWyLTY4Er41KRiIgktS73OMxsppkNc/fnCXcz8hiwH3gW0KEqEZF+qLtDVfcCjZHx2cBNwG+AnUSefSEiIv1Ld4eqUt39o+dffA5YHOkp91EzWxXf0kREJBl1t8eRGnmAEsAZwAttlh1uz7oiItKHdffl/zvgb2a2G6gDXgYws2OAqjjXJiIiSai7Lkd+bGbPA8OBZ9s8fS8F+Nd4FyciIsmn28NN7v5GB/M2xKccERFJdtHeACgiIgIEFBxmNsTM/mpmGyOvgztpFzKzVZFhaaLrFBGRQwW1x3ET8Ly7TwSej0x3pM7dZ0SGBYkrT0REOhNUcFxA+EZCIq8XBlSHiIjEKKjgOMrddwBEXos6aZdlZqVm9oaZdRkuZrYo0ra0oqKit+sVEZGIuN3EZ2bPAcM6WPSdGD5mjLtvj/TG+4KZve3umztq6O6LiXSDUlJS4h21ERGRnotbcLj7mZ0tM7OdZjbc3XeY2XBgVyefsT3yusXMlgEnAB0Gh4iIJEZQh6qWAldExq8AnmzfwMwGm1lmZLwQmAe8m7AKRUSkQ0EFx63AWWa2ETgrMo2ZlZjZfZE2xwGlkScNvgjc6u4KDhGRgAXSUaG7VxLuNLH9/FIizzF399eA6QkuTUREuqE7x0VEJCYKDhERiYmCQ0REYqLgEBGRmCg4REQkJgoOERGJiYJDRERiouAQEZGYKDhERCQmCg4REYmJgkNERGKi4BARkZgoOEREJCYKDhERiYmCQ0REYqLgEBGRmCg4REQkJgoOERGJSSDBYWafNbM1ZtZiZiVdtDvXzNab2SYzuymRNYqISMeC2uN4B7gYeKmzBmaWCvwSOA+YAnzezKYkpjwREelMWhArdfe1AGbWVbNZwCZ33xJp+3vgAuDduBcoIiKdSuZzHCOBbW2myyPzOmRmi8ys1MxKKyoq4l6ciEh/Fbc9DjN7DhjWwaLvuPuT0XxEB/O8s8buvhhYDFBSUtJpOxER6Zm4BYe7n9nDjygHRreZHgVs7+FniohIDyXzoarlwEQzG2dmGcClwNKAaxIR6feCuhz3IjMrB+YCfzazZyLzR5jZ0wDu3gxcCzwDrAX+4O5rgqhXREQ+FtRVVY8Dj3cwfztwfpvpp4GnE1iaiIh0I5kPVYmISBJScIiISEwUHCIiEhMFh4iIxETBISIiMVFwiIhITBQcIiISk0Du4whCU1MT5eXl1NfXB12KtJGVlcWoUaNIT08PuhQRiVK/CY7y8nLy8vIYO3Zsd925S4K4O5WVlZSXlzNu3LigyxGRKPWbQ1X19fUUFBQoNJKImVFQUKC9QJE+pt8EB3T74CgJgH4nIn1PvwoOERHpOQWHiIjERMEhIiIxUXAkqbq6Ok477TRCoVAg629sbOTUU0+lubk5kPWLSPJScCSp+++/n4svvpjU1NRA1p+RkcEZZ5zBww8/HMj6RSR59Zv7ONr6wZ/W8O72/b36mVNGDOTmT0/tss3bb7/N1VdfzauvvgrAypUrueGGG3jhhRcOafvggw/y0EMPAVBWVsa5557LKaecwhtvvMHxxx/PlVdeyc0338yuXbt48MEHmTVrFgAPPPAAt99+O42NjcyePZu77rqL1NRULrzwQrZt20Z9fT3XX389ixYtoqysjPPOO49TTjmF1157jZEjR/Lkk08yYMAAAC688EK+/e1vc9lll/Xmj0pE+jjtcSTQ1KlT2bx5c+vhp29+85vcdttth7RrbGxky5YtjB07tnXepk2buP7661m9ejXr1q3joYce4pVXXuG2227jlltuAWDt2rU8/PDDvPrqq6xatYrU1FQefPBBILwHs2LFCkpLS7n99tuprKwEYOPGjVxzzTWsWbOG/Px8Hn300dZ1Tps2jeXLl8frxyEifVQgexxm9lng+8BxwCx3L+2kXRlQDYSAZncv6Y31d7dnEC8pKSlMnTqVNWvWsHHjRsaMGcOJJ55ITU0NX//618nIyGD+/Pmcfvrp5OfnH/TecePGMX36dCAcQGeccQZmxvTp0ykrKwPg+eefZ8WKFcycORMInycpKioC4Pbbb+fxx8NP6922bRsbN25k2LBhjBs3jhkzZgBw0kkntX4WQGpqKhkZGVRXV5OXlxfPH42I9CFBHap6B7gYuDeKtqe7++4415Mwc+bM4dVXX+Wuu+7iL3/5CwCPPfYYl1xyCZ/+9Kf53Oc+x/nnn3/I3dSZmZmt4ykpKa3TKSkprSew3Z0rrriCn/zkJwe9d9myZTz33HO8/vrrZGdnM3/+/NbPb/u5qamp1NXVHfTehoYGsrKyemnrReRIEMihKndf6+7rg1h30ObMmcN//Md/cNFFFzFy5Egg3I/W6NGjgfCX9+DBgwmFQjF3xXHGGWfwyCOPsGvXLgD27NnDe++9R1VVFYMHDyY7O5t169bxxhtvRPV5lZWVDB06VB0QishBkv0chwPPmtkKM1sUdDG94dhjjyUzM5NvfetbrfNGjRpFeXk5AC0tLQCcffbZvPLKKzF99pQpU/jRj37E2WefTXFxMWeddRY7duzg3HPPpbm5meLiYr773e8yZ86cqD7vxRdf5Pzzz4+pBhHpB9w9LgPwHOFDUu2HC9q0WQaUdPEZIyKvRcBbwKldtF0ElAKlY8aM8fbefffdQ+YF4ZprrvElS5YcNO/AgQO+cOFCv/rqq/2BBx5wd/eVK1f6F7/4xSBKbHXRRRf5unXr4r6eZPndiPRnQKlH+f0et3Mc7n5mL3zG9sjrLjN7HJgFvNRJ28XAYoCSkhLv6bp72+bNm/nkJz/JvHnzuOKKKw5alpOTw3//938fNO+EE07g9NNPJxQKBXIvR2NjIxdeeCGTJ09O+LpFJLkl7X0cZpYDpLh7dWT8bOCHAZd12CZMmMC6deties+Xv/zlOFXTvYyMDC6//PLA1i8iySuQcxxmdpGZlQNzgT+b2TOR+SPM7OlIs6OAV8zsLeBN4M/u/pcg6hURkY8Fssfh7o8Dj3cwfztwfmR8C3B8gksTEZFuJPtVVSIikmQUHCIiEhMFh4iIxETBISIiMVFwJFBqaiozZsxg6tSpHH/88fzsZz9rvVO8M2VlZa3dq8fD0qVLufXWW7tss337di655JJOl+/bt4+77ror6vYi0rcpOBJowIABrFq1ijVr1vDXv/6Vp59+mh/84AddvudwgiOWpwYuWLCAm266qcs2I0aM4JFHHul0efvg6K69iPRt/Tc45s/vfmj7rIz582HJkvD47t2Hto1RUVERixcv5s4778TdCYVC3HjjjcycOZPi4mLuvTfccfBNN93Eyy+/zIwZM/j5z3/eabtly5Zx+umn84UvfKG1q/Vjjz2Wq666imnTpnHZZZfx3HPPMW/ePCZOnMibb74JwJIlS7j22msBWLhwIddddx0nn3wy48ePb/3yLysrY9q0aQCsWbOGWbNmMWPGDIqLi9m4cSM33XQTmzdvZsaMGdx4440HtQ+FQtxwww1Mnz6d4uJi7rjjjph/ViKSXJL2zvH+YPz48bS0tLBr1y6efPJJBg0axPLly2loaGDevHmcffbZ3Hrrrdx222089dRTACxevLjDdgBvvvkm77zzDuPGjaOsrIxNmzbxxz/+kcWLFzNz5szWhz8tXbqUW265hSeeeOKQmnbs2MErr7zCunXrWLBgwSGHnO655x6uv/56LrvsMhobGwmFQtx666288847rFq1CuCgZ3osXryYrVu38o9//IO0tDT27NkTp5+miCRK/w2OZcsOv31hYezv70S4bzF49tlnWb16dev/8quqqti4cSMZGRkHte+q3axZsxg3blxr22ge/tTehRdeSEpKClOmTGHnzp2HLJ87dy4//vGPKS8v5+KLL2bixIldbt9zzz3H1VdfTVpa+E9tyJAhUfxURCSZ9d/gSAJbtmwhNTWVoqIi3J077riDc84556A2y9oFVFftcnJyDpoXzcOf2mv7no9Cra0vfOELzJ49mz//+c+cc8453HfffYwfP77TbXR3zKzT5SLS9/TfcxwBq6io4Oqrr+baa6/FzDjnnHO4++67aWpqAmDDhg3U1NSQl5dHdXV16/s6a5coW7ZsYfz48Vx33XUsWLCA1atXH1JjW2effTb33HNPa1DpUJVI36c9jgSqq6tjxowZNDU1kZaWxpe+9CW+8Y1vAHDVVVdRVlbGiSeeiLszdOhQnnjiCYqLi0lLS+P4449n4cKFXH/99R22S5SHH36YBx54gPT0dIYNG8b3vvc9hgwZwrx585g2bRrnnXce11xzTWv7q666ig0bNlBcXEx6ejpf/epXW09p8diEAAAGZ0lEQVTGi0jfZB0djujrSkpKvLS09KB5a9eu5bjjjguoIumKfjciwTOzFe5eEk1bHaoSEZGYKDhERCQm/So4jsTDcn2dficifU+/CY6srCwqKyv1RZVE3J3KykqysrKCLkVEYtBvrqoaNWoU5eXlVFRUBF2KtJGVlcWoUaOCLkNEYtBvgiM9Pf2gu6pFROTwBHKoysz+y8zWmdlqM3vczPI7aXeuma03s01m1nUXriIikhBBneP4KzDN3YuBDcC32zcws1Tgl8B5wBTg82Y2JaFViojIIQIJDnd/1t0/6izpDaCjg9yzgE3uvsXdG4HfAxckqkYREelYMpzj+DLwcAfzRwLb2kyXA7M7+xAzWwQsikweMLP1h1lPIbD7MN/bV2mbj3z9bXtB2xyro6NtGLfgMLPngGEdLPqOuz8ZafMdoBl4sKOP6GBep9fSuvtiYPFhlHrwSs1Ko73t/kihbT7y9bftBW1zPMUtONz9zK6Wm9kVwKeAM7zjmyvKgdFtpkcB23uvQhERORxBXVV1LvAtYIG713bSbDkw0czGmVkGcCmwNFE1iohIx4K6qupOIA/4q5mtMrN7AMxshJk9DRA5eX4t8AywFviDu69JQG09PtzVB2mbj3z9bXtB2xw3R2S36iIiEj/9pq8qERHpHQoOERGJSb8Nju66MzGzTDN7OLL872Y2NvFV9p4otvcbZvZupBuY580s6mu6k1W0XdaY2SVm5mbW5y/djGabzeyfI7/rNWb2UKJr7G1R/G2PMbMXzewfkb/v84Oos7eY2f1mtsvM3ulkuZnZ7ZGfx2ozO7HXi3D3fjcAqcBmYDyQAbwFTGnX5uvAPZHxS4GHg647ztt7OpAdGf9aX97eaLc50i4PeIlwDwYlQdedgN/zROAfwODIdFHQdSdgmxcDX4uMTwHKgq67h9t8KnAi8E4ny88H/h/he+HmAH/v7Rr66x5HNN2ZXAD8JjL+CHCGmXV0U2Jf0O32uvuL/vGl0Z11A9OXRNtlzf8B/hOoT2RxcRLNNn8V+KW77wVw910JrrG3RbPNDgyMjA+ij98P5u4vAXu6aHIB8FsPewPIN7PhvVlDfw2OjrozGdlZGw9fGlwFFCSkut4Xzfa29RXC/2Ppy7rdZjM7ARjt7k8lsrA4iub3PAmYZGavmtkbkXuq+rJotvn7wBfNrBx4GvjXxJQWmFj/vccsGfqqCkI03ZnE1OVJkot6W8zsi0AJcFpcK4q/LrfZzFKAnwMLE1VQAkTze04jfLhqPuG9ypfNbJq774tzbfESzTZ/Hlji7j81s7nA/0S2uSX+5QUi7t9d/XWPI5ruTFrbmFka4V3crnYPk1lU3beY2ZnAdwjf0d+QoNripbttzgOmAcvMrIzwseClffwEebR/10+6e5O7bwXWEw6Sviqabf4K8AcAd38dyCLcGeCRKu7dNfXX4IimO5OlwBWR8UuAFzxy5qkP6nZ7I4dt7iUcGn39uDd0s83uXuXuhe4+1t3HEj6vs8DdS4Mpt1dE83f9BOELITCzQsKHrrYktMreFc02vw+cAWBmxxEOjiP5GdJLgcsjV1fNAarcfUdvrqBfHqpy92Yz+6g7k1TgfndfY2Y/BErdfSnwa8K7tJsI72lcGlzFPRPl9v4XkAv8MXINwPvuviCwonsoym0+okS5zc8AZ5vZu0AIuNHdK4Orumei3OZvAr8ys/9N+JDNwj78n0DM7HeEDzUWRs7b3AykA7j7PYTP45wPbAJqgSt7vYY+/PMTEZEA9NdDVSIicpgUHCIiEhMFh4iIxETBISIiMVFwiIhITPrl5bgivcXMCoDnI5PDCF/i+tE9ArXufnIghYnEkS7HFeklZvZ94IC73xZ0LSLxpENVInFiZgcir/PN7G9m9gcz22Bmt5rZZWb2ppm9bWYTIu2GmtmjZrY8MswLdgtEOqbgEEmM44HrgenAl4BJ7j4LuI+Pe2v9BfBzd58JfCayTCTp6ByHSGIs/6i/IDPbDDwbmf82kb6jgDOBKW0e+zLQzPLcvTqhlYp0Q8EhkhhtextuaTPdwsf/DlOAue5el8jCRGKlQ1UiyeNZ4NqPJsxsRoC1iHRKwSGSPK4DSsxsdaT32quDLkikI7ocV0REYqI9DhERiYmCQ0REYqLgEBGRmCg4REQkJgoOERGJiYJDRERiouAQEZGY/H/YuWtgFwq72QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1a161e55c0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "t=np.linspace(0,1,100)\n",
    "plt.ylim([-2,2])\n",
    "plt.plot(t, exact_mean(t),label=\"$y_0$(mean)\")\n",
    "\n",
    "plt.axhline(y=1.0, color='r', linestyle='-.',label='Deterministic')\n",
    "\n",
    "#plt.yticks(np.arange(-2, 2, step=0.5))\n",
    "plt.ylabel(\"Solution\")\n",
    "plt.xlabel(\"Time\")\n",
    "plt.legend(loc=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}