{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Betrachtet wird die stationäre Verteilung des Kapitalstocks bei Technologie-Schocks (AR(1)-Modell für $\\theta_t$, siehe Gleichung (6.2)).\n",
    "\n",
    "Gegebene Parameter: mit α = 0.3, β = 0.9 und $σ_ε^2$  = 0.01\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy import pi\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%load_ext ipydex.displaytools\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Normalverteilung** einer Zufallsgröße $x \\sim \\mathcal N(\\mu, \\sigma)$:\n",
    "$$\n",
    "f_{normal}(x)=\\frac{1}{\\sqrt{2\\pi\\sigma^2}}\\operatorname{exp}\\left(-\\frac{(x-\\mu)^2}{2\\sigma^2}\\right)\n",
    "$$\n",
    "\n",
    "**Lognormalverteilung**\n",
    "Verteilung einer Zufallsgröße $y$ mit $\\operatorname{ln}(y) \\sim \\mathcal N(\\mu, \\sigma)$:\n",
    "$$\n",
    "f_{lognorm}(x)=\\frac{1}{\\sqrt{2\\pi}\\sigma x }\\,\\exp\\Big( -\\frac{(\\ln(x)-\\mu)^2}{2\\sigma^2}\\Big) \\quad, x > 0\n",
    "$$\n",
    "\n",
    "Siehe auch:\n",
    "\n",
    "https://de.wikipedia.org/wiki/Logarithmische_Normalverteilung\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def norm_dist(x, mu, sigma):\n",
    "    return 1/np.sqrt(2*pi*sigma**2) * np.exp(-(x-mu)**2/(2*sigma**2))\n",
    "\n",
    "def lognorm_dist(x, mu, sigma):\n",
    "    return 1/(np.sqrt(2*pi)*sigma*x) * np.exp(-(np.log(x)-mu)**2/(2*sigma**2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "muk := -1.8704761714053748"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "---\n"
     ]
    }
   ],
   "source": [
    "xx = np.linspace(0.1, 0.22, 10000)\n",
    "\n",
    "sigma = np.sqrt(.01)\n",
    "alpha = 0.3\n",
    "beta = 0.9\n",
    "\n",
    "# epsilon wird als Mittelwertfrei angenommen\n",
    "mu_eps = 0\n",
    "\n",
    "muk = np.log(alpha*beta)/(1-alpha) + mu_eps/(1-alpha) ##:\n",
    "ff_ln = lognorm_dist(xx, muk, sigma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "maxval := 0.15251725172517253"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "---\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3Xnc1XP+//HHq81SUnKVVGRMsnwtJVka/CyRECEiZJjJOoxlMPgOMzcGjf1rLKE0jDCWZqihNEhGSBIVkq1CRWUpWt+/P17ncl1trnNdZ3l/zuc877fbuX0+17nOuc7rU9d5Xp/z/rwXCyEgIiKlr17sAkREJD8U6CIiKaFAFxFJCQW6iEhKKNBFRFJCgS4ikhIKdBGRlFCgi4ikhAJdRCQlGhTzxTbddNPQvn37Yr6klIuvvvJtixZx6xApgDfeeOPLEEJFTY8raqC3b9+eCRMmFPMlpVzcf79vTzklZhUiBWFmn2TzODW5iIikhAJdRCQlFOgiIimhQBcRSQkFuohISijQRURSQoEuIpISRe2HLpKTlSvhpZfgmWdgyhRYvNgHEu28M4QAbdrErlAkKgW6JF8I8MQTcNll8P77a37/0Ud9u+22fttjj+LWJ5IQCnRJtoUL4bTTPNAB2rWDvn2ha1do2hS++AJefBGGDYN334W99oLzz4drr4VGjeLWLlJkCnRJrk8+ge7dYfp0D+9rr4Vf/xoaNlz1cSefDLvuCiNGwLPPwk03wYQJ/kdAc7tIGdFFUUmmGTNg7709zHfeGSZOhLPOWjPMK224IfTpA//9L2y+OYwd68+fO7e4dYtEVGOgm1k7M3vezKaZ2RQzOy9z/1VmNtvMJmVuPQtfrpSFL7+EHj1g5kzo1s2bVLbeOrvndu0Kr70GO+wA06b5Gf78+YWtVyQhsjlDXw5cGELYDtgDONvMts987+YQwi6Z28iCVSnlY+lS6N0bPvgAOnXyHi0bb1y7n9GmDYwZAx07wuTJcOyxsHx5YeoVSZAaAz2E8HkIYWJm/1tgGqD+YVIYl18O48ZB27bw9NPQpEndfk6rVjB6NLRs6eF+4YX5rVMkgWrVhm5m7YFOwKuZu84xs8lmNtjMmue5Nik3o0bBDTdA/freFXHzzXP7ee3a+YXRhg3httvgscfyU6dIQmUd6GbWBHgc+G0I4RvgTmBrYBfgc+DGdTxvgJlNMLMJ8+bNy0PJkkpff121OMVVV8Gee+bn53br5r1eAAYMgNmz8/NzRRIoq0A3s4Z4mP89hPAEQAhhTghhRQhhJXAP0HVtzw0hDAohdAkhdKmoqHEFJSlXl18On3/ug4J+//v8/uyzz4ZDDoEFC/yPRgj5/fkiCZFNLxcD7gOmhRBuqnZ/62oP6w28k//ypCy89hrccYc3tdx9t2/zyQwGD4ZNN4XnnoO//S2/P18kIbI5Q+8GnATsv1oXxYFm9raZTQb2A84vZKGSUitX+hl0CHDBBbDTToV5nc02q2p6ueiiqkWlRVKkxpGiIYRxgK3lW+qmKLl79FEf1dm6NVx5ZWFf68QTYcgQeP55uPhiuO++wr6eSJFppKjEs3Spt50D/PGP0LhxYV/PDO680+d4GTzY/5CIpIgCXeK5+2748EOfIfGXvyzOa3bsCOed5/sXXaQLpJIqCnSJ4/vv4ZprfP/aa6FBEeeJu+wy2GQTn1Lg6aeL97oiBaZAlziGDIE5c6BzZzjiiOK+drNm8Ic/+P7FF8OyZcV9fZECUaBL8S1bBgMH+v5ll3nbdrGdeaZP+PXuu/Dgg8V/fZECUKBL8T30kM91vu22PhFXDI0aVfWqufpqnaVLKijQpbhWroTrrvP9Sy+FehF/BY8/Hjp08Auzf/97vDpE8kSBLsX17LPezNGuHZxwQtxaGjSAK67w/auv1hS7UvIU6FJct93m27PPXvfqQ8V0wgnw85/7CklqS5cSp0CX4nn/fV+wYv314Ve/il2Nq36WPnCgNwmJlCgFuhTPX//q2379krV48wkn+IIa06bBSM1oIaVLgS7F8e233vcc4De/iVvL6ho2hN/+1vdvuCFuLSI5UKBLcTz4oIf63nvDzjvHrmZNv/41NG3qo0dffz12NSJ1okCX4rj3Xt+edVbcOtalaVM4/XTf11m6lCgFuhTepEkwcSI0bw5HHhm7mnU791y/SPrYY/DRR7GrEak1BboU3uDBvu3Xz3u4JFXbtj7YaOVKX0FJpMQo0KWwfvihqn/3aafFrSUblRds77sPFi+OW4tILSnQpbCGD/fFmTt3hl12iV1NzXbbzW8LFsDDD8euRqRWFOhSWJXLvJXC2Xmlc87x7e23awEMKSkKdCmcmTNhzBhYbz1vmy4Vxx4Lm24Kb74J48fHrkYkawp0KZxhw/wM94gjvIdLqVh//apPFJWjW0VKgAJdCuehh3wbe1bFujjjDJ/a9x//8JWVREqAAl0KY8oUeOstX+6tR4/Y1dRe+/Zw2GGwdCncf3/sakSyokCXwhg2zLfHHONt6KVowADf3nuvLo5KSVCgS/6FUNXc0q9f3Fpy0aMHtGkDH3wAY8fGrkakRgp0yb/x433ofJs2PhlXqapfH0491ffvuSduLSJZUKBL/lWenfft66FYyk49Fcx8fpcFC2JXI/KTFOiSX8uXwyOP+H4p9m5ZXfv20L07LFmihaQl8RTokl8vvgjz5kGHDtCpU+xq8qNyubx77tHFUUm0GgPdzNqZ2fNmNs3MppjZeZn7NzGz0WY2PbMtoZEjUjCPP+7bPn28qSINevXykaOTJ8OECbGrEVmnbM7QlwMXhhC2A/YAzjaz7YFLgTEhhA7AmMzXUs5WrIAnnvD9o4+OW0s+rbcenHyy71cu1CGSQDUGegjh8xDCxMz+t8A0oA1wBDA087ChQIJXLpCi+O9/fVTlVlulp7mlUmWzy7BhmlZXEqtWbehm1h7oBLwKtAohfA4e+kDLfBcnJaayueXoo9PT3FJpu+1g9919XdThw2NXI7JWWQe6mTUBHgd+G0L4phbPG2BmE8xswrx58+pSo5SClStXDfQ0qmx2GTr0px8nEklWgW5mDfEw/3sIIdNIyhwza535fmtg7tqeG0IYFELoEkLoUlFRkY+aJYlefx1mzfJl3Lp2jV1NYfTtC40awXPPwezZsasRWUM2vVwMuA+YFkK4qdq3/gX0z+z3B/6Z//KkZFSenR91lM9SmEabbAKHH+6fRtQnXRIom3deN+AkYH8zm5S59QSuA7qb2XSge+ZrKUch+EhK8Mm40qx/5hxm6FD1SZfEaVDTA0II44B1XeE6IL/lSEmaNMnnbmnVCvbaK3Y1hdWjB1RUwNSp8MYb0KVL7IpEfpTSz8ZSVJXNLb17l/7cLTVp2LBqSgNdHJWEUaBL7p56yre9e8eto1gqm12GDfMFMEQSQoEuufn0Ux8S37gx7Ltv7GqKY5ddYMcd4auvYMSI2NWI/EiBLrmpDLSDDirdlYlqy2zVi6MiCaFAl9w8/bRvDzssbh3F1q+fXy8YMQK+/DJ2NSKAAl1ysWgRjBnj+z17xq2l2DbbDA4+2Od/f/jh2NWIAAp0ycV//uMLP+y2mwdcuTnxRN9qkJEkhAJd6q5cm1sq9erlF4PHj4cZM2JXI6JAlzoKQYHeuLFPdQA6S5dEUKBL3UyaBJ99Bptvnr65z2ujXz/fPvigpgKQ6BToUjeVZ+eHHpq+uc9r44ADfMqD6dO1PJ1Ep0CXuin35pZKDRr4tLrgZ+kiESnQpfbmzIHXXvOBRAdofrYfe7s8/LB3YxSJRIEutTdypG/3398vDJa7XXeFbbaBuXN98QuRSBToUntqblmVmfqkSyIo0KV2liyBUaN8/9BD49aSJJVT6j75pI+gFYlAgS61M3YsfPedzza45Zaxq0mOrbeGPff0MP+nVmOUOBToUjtqblm36n3SRSJQoEv2QqhazEKBvqZjj/VujKNG+QVSkSJToEv23n3X1w5t0QJ23z12NclTUeEzMK5YAY88ErsaKUMKdMleZXNLz57pXzu0rtTbRSJSoEv21H5es169oEkTePVVnw5ApIgU6JKd+fPh5Ze9jfigg2JXk1wbblg1A+NDD8WtRcqOAl2y8+yz3ja8997QrFnsapJNMzBKJAp0yY6aW7K3//6+gtMHH8Drr8euRsqIAl1qtnw5/Pvfvq9Ar5lmYJRIFOhSs1degQULoEMHn4RKalZ9BsZly+LWImVDgS41U3NL7XXuDNtuC/PmwejRsauRMqFAl5op0Guv+gyMDzwQtxYpGzUGupkNNrO5ZvZOtfuuMrPZZjYpc+tZ2DIlmg8/hKlToWlT+MUvYldTWip7uwwfDt98E7cWKQvZnKHfD/RYy/03hxB2ydxG5rcsSYwRI3x78MHQqFHcWkpN+/bezfOHH+CJJ2JXI2WgxkAPIYwF5hehFkkiNbfk5qSTfKtmFymCXNrQzzGzyZkmmeZ5q0iS49tv4YUXvD34kENiV1Oa+vTxtVeffx5mzYpdjaRcXQP9TmBrYBfgc+DGdT3QzAaY2QQzmzBv3rw6vpxE8dxzsHQp7LGHzyQotdesGRx+uI8Y1YRdUmB1CvQQwpwQwooQwkrgHqDrTzx2UAihSwihS4VCobSouSU/qje7aCoAKaA6BbqZta72ZW/gnXU9VkrUypVVF0QV6Lnp0cPnkJ8yBd56K3Y1kmLZdFscBrwCdDSzWWZ2GjDQzN42s8nAfsD5Ba5Tiu2NN2DOHGjXztcPlbpr1KhqKgBdHJUCyqaXy/EhhNYhhIYhhLYhhPtCCCeFEHYMIewUQugVQvi8GMVKEVVvbjGLW0saVA4yeughnxtHpAA0UlTWTu3n+bX77j4XzhdfwJgxsauRlFKgy5pmz4aJE2GDDWC//WJXkw6aCkCKQIEuaxqZGfh74IEe6pIflYH+5JPw3Xdxa5FUUqDLmtTcUhg/+xl06waLF3uoi+SZAl1W9f33PqAI4NBD49aSRpoKQApIgS6reuEFP4Ps1AnatIldTfoce6x3YxwzBj77LHY1kjIKdFmVmlsKq3lz/+SzcqV3YRTJIwW6VAlBgV4ManaRAlGgS5W334ZPP4WWLaFLl9jVpFfPnn6mPnmypgKQvFKgS5XqZ+f19KtRMOutByec4PtDhsStRVJF71qp8tRTvlVzS+GdeqpvH3zQpygWyQMFuri5c+HVV70HRvfusatJv06dYKed4Kuvqv6QiuRIgS5uxAi/KLr//tCkSexq0s+s6ix98OC4tUhqKNDFVbafH3543DrKSb9+0LAhPPOMz58jkiMFusCSJTBqlO+r/bx4Nt0UevXyPunqwih5oEAXHx363XfeprvFFrGrKS/Vm120PJ3kSIEuVRfl1NxSfAcdBJtvDtOnw8svx65GSpwCvdxVHx2qQC++Bg2gf3/fV590yZECvdy98w588gm0agW77Ra7mvJ0yim+feQRzZMuOVGgl7vK5pZDD9Xo0Fi22QZ+8QtYtAj+8Y/Y1UgJ0zu43Kn9PBnUJ13yQIFezqqPDj3wwNjVlLc+fXxA17hxMHVq7GqkRCnQy9nIkRodmhRNmvhAI4BBg+LWIiVLgV7O1NySLKef7tuhQ30pQJFaUqCXq++/9yHn4KMVJb5OnaBrV1i4EB59NHY1UoIU6OVq1ChfO3S33aBt29jVSKUzzvDtXXfFrUNKkgK9XD35pG+PPDJuHbKq446DjTeG8eO1mpHUmgK9HC1fXtV+3rt33FpkVRtuCCef7Pt33x23Fik5CvRyNHYszJ8PHTvCdtvFrkZWV3lx9MEHNXJUaqXGQDezwWY218zeqXbfJmY22symZ7bNC1um5NXw4b7V2Xky7bCDjxz99lsYNix2NVJCsjlDvx/osdp9lwJjQggdgDGZr6UUhFAV6Go/T67qF0c1ra5kqcZADyGMBeavdvcRwNDM/lBAyVAq3ngDZs70KVs1GVdyHX00tGgBEyf6BVKRLNS1Db1VCOFzgMy2Zf5KkoKq3rtFk3El1/rrw4ABvn/bbXFrkZJR8He0mQ0wswlmNmHevHmFfjmpSWWgq/08+c48E+rXh8ce05qjkpW6BvocM2sNkNnOXdcDQwiDQghdQghdKioq6vhykhfvvQfTpkGzZrDvvrGrkZq0a+dNL8uXw513xq5GSkBdA/1fQGaZFfoD/8xPOVJQlWfnhx3mq81L8p17rm/vvht++CFuLZJ42XRbHAa8AnQ0s1lmdhpwHdDdzKYD3TNfS9JVzg9yzDFx65Ds7bUXdO4MX36pLoxSo2x6uRwfQmgdQmgYQmgbQrgvhPBVCOGAEEKHzHb1XjCSNNOnw5tvwkYbwcEHx65GsmUG553n+7fdpi6M8pPUzaFcVC5tdsQR3oNCSsdxx0HLljBpErz0UuxqJMEU6OWisrnl2GPj1iG1t956VQONbrklbi2SaAr0cvDeez5zX9OmcNBBsauRujjjDF8qcPhweP/92NVIQinQy0Flc8uRR/rZnpSe1q19FsYQ4MYbY1cjCaVALwdqbkmHiy7yi6RDh8IXX8SuRhJIgZ5206bB22/7ogndu8euRnLRsaN/ylqyRNMByFop0NOusrmld29vg5XSdvHFvr3jDp9eV6QaBXqahQAPP+z7ffrErUXyY489YO+94euvYdCg2NVIwijQ0+zNN73JZdNN1dySJpdc4tubb4alS+PWIomiQE+zBx7wbd++mrslTQ45xFc1mj0b/va32NVIgijQ02r58qq5P046KW4tkl/16sFll/n+NdfAsmVx65HEUKCn1ZgxMGcOdOiglYnS6LjjvNfLxx9XfRKTsqdAT6vKN/mJJ3rfZUmX+vXhf//X96++WmfpAijQ0+m776rmPj/xxLi1SOH07QvbbAMffaSzdAEU6Ok0fDgsXuxzaf/sZ7GrkUKpfpautnRBgZ5O1ZtbJN0qz9I//FBn6aJAT51PP4XRo31UqOZuSb8GDarO0v/4Ry1TV+YU6GkzZIiPED3qKGjRInY1UgzHHw//8z/+x1yLSZc1BXqarFgBgwf7/q9+FbcWKZ769eG6zLK+V1/t0wJIWVKgp8lzz/lZ2lZbwX77xa5GiqlnT9hnH5g/HwYOjF2NRKJAT5P77vPtqaf6aEIpH2Zw/fW+f/PN8NlnceuRKPSuT4t587y7Yr16cMopsauRGPbYw6+dfP+9XyCVsqNAT4sHHvB+yIccAm3bxq5GYvnzn71N/d57YfLk2NVIkSnQ02DlSrjrLt/XxdDy1rEjnHOO/06ce673eJKyoUBPg1GjYPp0aNcODjssdjUS21VX+Rz4L75YtWKVlAUFehrcfrtvzzrLB5pIeWvWzJtewBeWXrQobj1SNAr0UjdjBowcCeutB6edFrsaSYpTT4XOnWHmzKreL5J6CvRSd8cd3k7aty9UVMSuRpKifn247TbfHzgQ3n8/bj1SFAr0UrZoUdXI0N/8Jm4tkjzdukH//rBkCZx+ui6QloGcAt3MPjazt81skplNyFdRkqWhQ2HhQth9d9h119jVSBLdeKN/cnvhhao//pJa+ThD3y+EsEsIoUsefpZka/lyf7MCXHhh3FokuVq0gFtv9f2LLoIvvohbjxSUmlxK1eOP+xzYW2/towNF1qVvXx9wtnCh902X1Mo10AMwyszeMLMB+ShIshBC1QRMF13kF8BE1sXMp9Vt3Nj7patvemrlGujdQgidgUOAs81sn9UfYGYDzGyCmU2YN29eji8nAPznPzBxIrRs6Re9RGqy5Zbwl7/4/umnw+zZceuRgsgp0EMIn2W2c4Enga5recygEEKXEEKXCnWry4/KfsXnnQcbbBC3FikdZ5zhTS8LFng/9ZUrY1ckeVbnQDezxma2UeU+cBDwTr4Kk3V4+WVfYm6jjeDMM2NXI6XEzKdYbtHCp4v4619jVyR5lssZeitgnJm9BbwGjAghPJOfsmSd/vAH355/PjRvHrcWKT2tW8OgQb5/8cXw9ttx65G8qnOghxA+DCHsnLntEEK4Jp+FyVq88IK3n2+8sQe6SF0cdZQ3ufzwAxx9NHzzTeyKJE/UbbFUhABXXun7F17oEzCJ1NX//R/suKPP0nnaaRpFmhIK9FIxejSMHQubbOIXQ0VyseGG8Nhjfi3msceq5n2RkqZALwUrVnh/c4BLLoGmTePWI+mwzTYwZIjvX3SRN+lJSVOgl4IhQ/zi1ZZbaqSf5NfRR3uYL1/u+9Onx65IcqBAT7pvv4UrrvD966+H9dePW4+kz3XX+UpX8+f7dsGC2BVJHSnQk+6662DOHF/R/dhjY1cjaVS/Pjz0EOy0k8+bfswxsHRp7KqkDhToSfbuu3DDDb5/000+MESkEDbaCJ56Clq18q6xJ57o126kpCjQkyoEH6q9dCn86lew556xK5K022ILX86waVOfwOuMM9SdscQo0JNq6FBftb2iQmtCSvF07uxn6uuvD/fe66NJFeolQ4GeRF98UdVN8aabvO+5SLHss4/Pt9+ggTf5XXKJQr1EKNCTJgRvYvnqK+jeHfr1i12RlKOePWHYMA/1v/zF16zV7IyJp0BPmkGDYMQIH9o/eLAuhEo8xxwDTzwBjRr5zIy//rX3V5fEUqAnyfvvwwUX+P6dd0LbtnHrETn8cHj6aZ93f/BgOOIIHxshiaRAT4pFi3yk3uLFcPzxvg6kSBJ07+5zCbVo4b1g9tlHKx4llAI9CUKAAQPgnXegY0e4667YFYmsqls3eOUV6NABJk2C3XeH8eNjVyWrUaAnwa23+ki9xo29zVKTb0kSdejgob733n6Gvs8+Pg2vesAkhgI9tieeqGo3HzIEtt8+bj0iP6VFC3juOZ/Cedkynyyub19YuDB2ZYICPa6XX/ZuiSHANddAnz6xKxKpWaNGcMst8OijPmXAo4/6YhmjRsWurOwp0GOZMMFntvvhBzj9dPj972NXJFI7ffrAG2/4xHGzZsHBB/t0AV9/HbuysqVAj+H11+HAA/1jau/ecPvt6m8upalDB3jpJbj2WmjYEO6+2y/sP/CA2tYjUKAX2wsveDewr7/2boqPPOKj8URKVYMGcOmlfra+554+3fPJJ/vF0wkTYldXVhToxTR0KBx0kId5nz4+tLphw9hVieTHjjvCuHFw//3QsqVfI9ptNzjqKO+SKwWnQC+GZcvgd7+DU07x/QsuUJhLOtWrB/37w3vv+e/8BhvAk0/64hl9++qMvcAU6IX2ySfeX/eGG3xlmNtvhxtv9H2RtGrWDAYOhBkz4JxzvFnmkUf8jH3ffeGf/9QCGgWgQC+UFSt8QqMdd/QRde3awdixcPbZsSsTKZ7WrX3w0YwZPiV006b+PjjySF/0/PLLtTB1HinQC+GVV3yo9Dnn+ERGRx3lw6X32it2ZSJxtGvn0/DOnOlz/G+9tY82/fOfYZtt4Be/8L7tn3wSu9KSpkDPp0mToFcvD+5XX4U2bbz98PHHtUiFCPgZ+vnn+1n5iy96e/uGG/oF1PPPh/btYddd4U9/8hMjTddbKwr0XK1YAcOHw377QadOvnxX48ZwxRUwdap/tBSRVZn5taX77/cVuoYN855fjRvDxIlw5ZV+YrTJJn6SdMst3nT5ww+xK080dYCui5Ur/Zdr2DBfTHfOHL+/SRNfbejSS331dBGp2UYbeQ+Yvn3h++99qt5nnoExY3yNgKee8ht4z7Cdd4auXX390+2399vGG8c9hoTIKdDNrAdwK1AfuDeEcF1eqkqaFSv8I+K4cf7LNmaMLxFXqUMHOOss+OUv9YslkosNNvAz8l69/OuZM/399uKLPsJ66lTv+rh698e2bWGHHbw9vn172Gor37ZvD82bF/kg4qlzoJtZfeCvQHdgFvC6mf0rhDA1X8UV3ZIl8PHHfkV+xgwP8Tff9NuiRas+dsstfYmu44/3MwUN3RfJv3btfPzGKaf419984yNSX3sN3n4bpkyBadN8LplZs+DZZ9f8GU2bwuab+6fmVq1gs82q9isqvItls2Z+MtasmX9iqFeardG5nKF3BT4IIXwIYGYPA0cA+Q/0BQtg6VK/QLJiRdWt+tdr+97SpR7Ea7t9/TXMmwdz51Zt589f9/wT7dp5H9oDDvCh+z//uUJcpNiaNvXrVfvtV3XfihXw4Yce7jNm+ElZ5e2jj/yPwDffwLvvZvcaZv46zZp5M+oGG/iF28pt9f0NNvBbo0Z+a9iwalt9v1Ej/yS/7bYF+EepkkugtwFmVvt6FrB7buWsw157Zf+fkYt69fwj2tZbV9122snPwCsqCv/6IlJ79et7WHbosOb3QvATtS++8NucOVXbOXP8ZO7rr/22cKHfvvuu6r58uvhiuP76/P7M1eQS6Gs7PV3j9NbMBgADALbYYou6vVKLFj43RP36fmvQILv9hg39qvnabhtt5D+zosK3LVv662gEp0h6mPn7ukULb2PPxvLlfka/cKF/ml+82C/W/tR22TJvEVi2bNX96vcVYfGaXAJ9FtCu2tdtgc9Wf1AIYRAwCKBLly51m09z3Lg6PU1EpNYaNPDukiU4diSXlv/XgQ5mtpWZNQL6Av/KT1kiIlJbdT5DDyEsN7NzgGfxbouDQwhT8laZiIjUSk790EMII4GReapFRERyUJqdLUVEZA0KdBGRlFCgi4ikhAJdRCQlFOgiIilhYV1zlxTixczmAXVdkmRT4Ms8lhOTjiV50nIcoGNJqlyOZcsQQo3zjxQ10HNhZhNCCF1i15EPOpbkSctxgI4lqYpxLGpyERFJCQW6iEhKlFKgD4pdQB7pWJInLccBOpakKvixlEwbuoiI/LRSOkMXEZGfkIhAN7MeZvaemX1gZpeu5fv7mNlEM1tuZses9r3+ZjY9c+tfvKrXVNfjMLNdzOwVM5tiZpPN7LjiVr6mXP5PMt9vamazzez24lS8bjn+fm1hZqPMbJqZTTWz9sWqe21yPJaBmd+xaWZ2m1m8NRSzOI4LMv/ek81sjJltWe17iXnPZ+qp07EU5H0fQoh6w6fenQH8DGgEvAVsv9pj2gM7AX8Djql2/ybAh5lt88x+8xI8jm2ADpn9zYHPgWal+H9S7fu3Ag8Bt5fq71fmey8A3TP7TYANS/FYgL2bEdXeAAAC5ElEQVSAlzM/oz7wCvD/Enwc+1X+WwNnAo9k9hPzns/DseT9fZ+EM/QfF5sOISwFKheb/lEI4eMQwmRg5WrPPRgYHUKYH0JYAIwGehSj6LWo83GEEN4PIUzP7H8GzAViLmKay/8JZrYr0AoYVYxia1DnYzGz7YEGIYTRmcd9F0JYXKS61yaX/5cArI+HznpAQ2BO4Uteq2yO4/lq/9bj8RXRIFnvecjhWArxvk9CoK9tsek2RXhuvuWlFjPrir/pZuSprrqo87GYWT3gRuB3BairLnL5f9kGWGhmT5jZm2b2FzOLuehsnY8lhPAK8Dx+Fvg58GwIYVreK8xObY/jNODfdXxuoeVyLD/K1/s+pwUu8iSrxaYL8Nx8y7kWM2sNPAD0DyGsceZbRLkcy1nAyBDCzIhNtNXlciwNgL2BTsCnwCPAKcB9eams9up8LGb2c2A7qs50R5vZPiGEsfkqrhayPg4zOxHoAuxb2+cWSS7HUnl/3t73SThDz2qx6QI8N99yqsXMmgIjgCtCCOPzXFtt5XIsewLnmNnHwA3AyWZ2XX7Lq5Vcf7/ezHycXg4MBzrnub7ayOVYegPjM81G3+FniXvkub5sZXUcZnYgcDnQK4SwpDbPLaJcjiX/7/tYFxOqXTBogF/Y2Iqqiwo7rOOx97PmRdGP8IsjzTP7m5TgcTQCxgC/jf3/keuxrPa9U4h/UTSX/5f6mcdXZL4eApxdosdyHPBc5mc0zPy+HZ7U48A/Fc0gc9Gw2v2Jec/n4Vjy/r6P8o+wln+UnsD7mYO+PHPfn/C/ZgC74X8JFwFfAVOqPfdU4IPM7ZeleBzAicAyYFK12y6leCyr/YzogZ6H36/uwGTg7UxINirFY8H/ON0NTAOmAjcl/Diewy/aVr4f/lXtuYl5z+dyLIV432ukqIhISiShDV1ERPJAgS4ikhIKdBGRlFCgi4ikhAJdRCQlFOgiIimhQBcRSQkFuohISvx/KBOhDQqv6lgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f9c426de710>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plt.plot(xx, ff1)\n",
    "plt.plot(xx, ff_ln, \"r-\", lw=2)\n",
    "axis = plt.axis()\n",
    "plt.vlines([0.15405], -5, 30, \"r\", alpha=0.5)\n",
    "plt.axis(axis)\n",
    "\n",
    "maxval = xx[np.argmax(ff_ln)] ##:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f9c42f330f0>"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3Xd4lFXawOHfSYHQi3QRIkhvQQIiovQiRUASendBBURRV/1UFMvusquiIBZABGQhlAQFAVcEQQVBCIjSRIoBIgghIi0Eksz5/jiTSkLazLxTnvu65pqZM+/MPG+SeXLmVKW1RgghhOfzszoAIYQQjiEJXQghvIQkdCGE8BKS0IUQwktIQhdCCC8hCV0IIbyEJHQhhPASktCFEMJLSEIXQggvEeDKN6tQoYIODg525VsKIYTH27Vr1zmtdcXcjnNpQg8ODiY6OtqVbymEEB5PKXU8L8dJk4sQQngJSehCCOElJKELIYSXcGkbuhDeLikpidjYWBITE60ORXigoKAgqlevTmBgYIGeLwldCAeKjY2lVKlSBAcHo5SyOhzhQbTWxMfHExsby+23316g15AmFyEcKDExkVtuuUWSucg3pRS33HJLob7dSUIXwsEkmYuCKuzfjjS5CM9is5Hw4yEObYwl4WISt9xajDu61SagVg2rIxPCclJDFx5BX7jIqiFL6V7iWyqG1mDEs1X4+z/K0mv8bVSoXZoR5dfw0yufQVKS1aFarmTJklaH4BDBwcGcO3cuU1lMTAzVq1fHZrNlKg8JCWHHjh15fu09e/awbt26XI+Ljo5m0qRJACxYsICJEyfm+T2sIAlduL3Ds76kbcVfeCmiASMTZ3O2SjP2PvAi3w/7gCOdH+WX4i1oen4z3aa25qFKq7m4aZfVIfu8lJQUp7xucHAwt912G999911a2S+//MKlS5do1apVnl4jOTk5zwk9NDSUmTNnFjheV5OELtyXzUZkrwW0eexOBid9wo+txzN4w98o8fuvsGoVLFoEX31Flfj9PP1JMw7X6o7/X/E071iWgy8utjp6t3L8+HE6depE06ZN6dSpEydOnADg6NGjtG7dmpYtW/LSSy+l1e43b95M+/btCQsLo379+gwdOhStNQAbN26kefPmNGnShDFjxnDt2jXAJNtXX32Vtm3bsmLFCtq3b8/kyZO57777aNCgATt37uTBBx+kTp06vPjii2mx9e3blxYtWtCoUSPmzJmT67kMHjyYpUuXpt1funQpgwcPBiAuLo7+/fvTsmVLWrZsydatWwGYOnUq48aNo2vXrowYMYKXXnqJZcuWERISwrJly7hy5QpjxoyhZcuWNG/enFWrVqX9HHr16nVDDKNGjSIyMjLtfl5+buvWraN+/fq0bduWSZMmZfu6haa1dtmlRYsWWog8SUrS77Wcr6sRq38MCNX6nXe0ttlu/pyEBK0nT9YLGa4r8Yf+ftCM3J/jYAcOHEi/A8655KJEiRI3lPXq1UsvWLBAa631vHnzdJ8+fbTWWvfs2VMvWbJEa631Bx98kPbcTZs26dKlS+uTJ0/qlJQU3bp1a/3dd9/pq1ev6urVq+tDhw5prbUePny4fvvtt7XWWtesWVP/+9//TnvPdu3a6WeeeUZrrfU777yjq1atqk+dOqUTExP1rbfeqs+dO6e11jo+Pl5rrXVCQoJu1KhRWnnNmjV1XFzcDedy+vRpXaVKFZ2UlKS11rp+/fp67969WmutBw8erL/77juttdbHjx/X9evX11pr/fLLL+s777xTJyQkaK21nj9/vp4wYULaa/7f//2fXrRokdZa6/Pnz+s6deroy5cv602bNumePXve8JyRI0fqFStW3PAzz+3nduzYMa211oMGDUp73awy/Q3ZAdE6DzlWaujC/WjNx+0/4Y2d7dhSvBsh6/8Djz8OuY0AKFYMpk9nxCddWOD/N/ouHchPo98xadDHbdu2jSFDhgAwfPhwtmzZklYeHh4OkPZ4qlatWlG9enX8/PwICQkhJiaGQ4cOcfvtt1O3bl0ARo4cybfffpv2nIEDB2Z6jQceeACAJk2a0KhRI6pWrUrRokWpVasWJ0+eBGDmzJk0a9aM1q1bc/LkSQ4fPnzTc6lSpQqNGjVi48aN7Nmzh8DAQBo3bgzAhg0bmDhxIiEhITzwwANcvHiRS5cupcVSrFixbF9z/fr1TJs2jZCQENq3b09iYmLat5j8yu7n9ssvv1CrVq208eWp3ygcLddRLkqp24BPgCqADZijtZ6hlJoKjAXi7Ic+r7XOvVFKiFysGb6MF7bezzdB3bl940fQunX+XmD4cO4vW5ZZfR/n/oXT2V5tNjX++Yhzgr0ZN/5HkpfhcUWLFk277e/vT3JyclrzQU5KlCiR7Wv4+fllej0/Pz+Sk5PZvHkzGzZsYNu2bRQvXjwtmWb03nvvMXfuXMA0W1SrVi2t2aVy5cqZkqPNZmPbtm3ZJu6ssWWktSYqKop69eplKj9z5ky2xwcEBKR1zGqtuX79+g3nDHn/uTlKXmroycBTWusGQGtgglKqof2xt7XWIfaLJHNRaEfmbmL04k58yoPUXfZa/pN5qt69CY94kMm8Tf9/tSDxs/85NlAP06ZNm7R258WLF9O2bVsAWrduTVRUFECmdumc1K9fn5iYGI4cOQLAokWLaNeuXYHjunDhAuXKlaN48eL88ssvbN++/YZjJkyYwJ49e9izZw/VqlUDoH///qxbt45ly5YxaNCgtGO7du3KrFmz0u7v2bMn2/ctVapUWs0doFu3brz77rtpiffHH3+8adzBwcHs2mU631etWkVSLqOr6tevz7Fjx4iJiQFg2bJlNz2+oHJN6Frr01rr3fbbl4CDwK1OiUb4tMTfThP+aAWmMpXW/+oD9q/rBTZgAE9PLUVNjvPEgFNw9KhjAnVzCQkJVK9ePe0yffp0Zs6cyfz582natCmLFi1ixowZALzzzjtMnz6dVq1acfr0acqUKXPT1w4KCmL+/PmEh4fTpEkT/Pz8eOSRgn/76d69O8nJyTRt2pQpU6bQOo//wMuWLUvr1q2pXLlypmnyM2fOJDo6mqZNm9KwYUM+/PDDbJ/foUMHDhw4kNYpOmXKFJKSkmjatCmNGzdmypQpN33/sWPH8s0339CqVSt++OGHm9b+AYoVK8b7779P9+7dadu2LZUrV871Z10geWloT70AwcAJoDQwFYgBfgY+Bsrl8JxxQDQQXaNGjWw7AYTQNpt+okaUDmeZtnXqrHVKimNeNyVFX+g1RAdzTK+t/6TW9o40Z8muQ8udXblyRdvsHccRERH6gQcesDgi73Xp0iWttdY2m00/+uijevr06dke55JOUaVUSSAKeEJrfRH4AKgNhACngbdy+IcxR2sdqrUOrVgx1x2UhI/6/uUvWXaiNR+Wex616BPwc1B/vZ8fpRe+y8cVnmXcL5M5P2W6Y17XS+zatYuQkBCaNm3K+++/z1tvZfsxFg4wd+5cQkJCaNSoERcuXODhhx92+HsonYfGeqVUILAG+FJrfcMnQikVDKzRWje+2euEhoZq2YJOZJX4ezzNa8bzWsrzhM3vBaNGOf5NNm7ksc4HuKxKMX/PndC0qePfAzh48CANGjRwymsL35Dd35BSapfWOjS35+ZaDVKmO3wecDBjMldKVc1wWD9gX54jFiKDfz+wlfop++l/3zkYOdI5b9KpE/8Yd4L1ugvbhrwLWaaOC+EN8vK99h5gONBRKbXHfukB/EcptVcp9TPQAZjszECFd4rd8Aszd9/DO/5Po2Z/mPtY80Io/Z8X+XeZf/HY/odJmf2R095HCKvkOg5da70FyO5TJsMURaE9P+oUjxBFzUd7QP36zn2zMmUYOqcdswde5eOn9jM2/BxUqODc9xTChWSmqLDMztm72fB7fZ4rMQtyGSbmKCo8jLdbLWXq1WdIeGmaS95TCFeRhC6soTXPPWPjFV6m1N8fgUqVXPO+ShE692Fas51ZswMgl2nmnkgpxVNPPZV2/80332Tq1KkujSHr4lWFsWDBAk6dOpXrcS+99BIbNmwAoH379vjiAAxJ6MIS3874kZiL5RhV/nPIkHxcomlTXu+3mzdtT3Lh6ddc+94uULRoUVauXHnDWuJ5lZyc7OCICi4lJSXPCf3VV1+lc+fOLojKfUlCF5Z4+RU/XuR1AidPBAs2ZGjw7nh6+n/JW6vvgGymm3uygIAAxo0bx9tvv33DYzktoztq1CiefPJJOnTowLPPPsvUqVMZOXIkXbt2JTg4mJUrV/LMM8/QpEkTunfvnjbV/dVXX6Vly5Y0btyYcePG3bBmyRdffMGAAQPS7m/evJnevXsDZkGsu+++mzvvvJPw8HAuX74MZF6GNyIigujoaIYOHUpISAhXr15l165dtGvXjhYtWtCtWzdOnz6ddg7ZfSvIuOFHZGQko+zDYkeNGsWkSZNo06YNtWrVSnuuzWZj/PjxNGrUiF69etGjRw+HfdtwNknowuU2v3+Ak3+VZHjJz2DCBGuCuPVWpow9zfuM5+KUN6yJwYkmTJjA4sWLuXDhQqbyiRMnMmLECH7++WeGDh2athsPwK+//sqGDRvSJhcdPXqUtWvXsmrVKoYNG0aHDh3Yu3cvxYoVY+3atWmvt3PnTvbt28fVq1dZs2ZNpvfr0qUL27dv58qVK4BZw2TgwIGcO3eO119/nQ0bNrB7925CQ0OZPj19iktQUBBbtmxh2LBhhIaGsnjxYvbs2UNAQACPPfYYkZGR7Nq1izFjxvDCCy8U+Od0+vRptmzZwpo1a3juuecAWLlyJTExMezdu5ePPvqIbdu2Ffj1XU0SunC5qS+lMIXXCJj4CJQrZ1kctV4bQzf/jXyw4Q744QenvIdSjr/kRenSpRkxYsQNu+3ktIwuQHh4OP7+/mn377//fgIDA2nSpAkpKSl0794dMEvhpi4ytWnTJu666y6aNGnC119/zf79+zO9X0BAAN27d+fzzz8nOTmZtWvX0qdPH7Zv386BAwe45557CAkJYeHChRw/fjzteVmX4U116NAh9u3bR5cuXQgJCeH1118nNjY2bz+UbPTt2xc/Pz8aNmyYtrLili1bCA8Px8/PjypVqtChQ4cCv76rySbRwqV+WH6c3+JLM7RoFDxh8WJZFSrw3MjTdP34CSa99CjFvvzM4W9h5Qq6TzzxBHfeeSejR4/O8ZiMy+jebOnbwMDAtGNTl75NTExk/PjxREdHc9tttzF16tQblr4Fk5zfe+89ypcvT8uWLSlVqhRaa7p06UJERES2ceW02JXWmkaNGuWr1pzxHLPGl3Gp29TmorzMnndXUkMXLvXWC/FM5m0CRgyBypWtDocm/x5GqP+PzF9fzWm1dKuUL1+eAQMGMG/evLSynJbRLYjU5FihQgUuX76cYztz+/bt2b17N3Pnzk2rebdu3ZqtW7emLcObkJDAr7/+mu3zMy51W69ePeLi4tISelJS0g3fCrKqXLkyBw8exGaz8emnn+Z6Xm3btiUqKgqbzcaZM2fYvHlzrs9xF5LQhcv89vMlvj5Sk4eYB489ZnU4RoUK/N/g47zB30l+9Z9WR+NwTz31VKbRLjkto1sQZcuWZezYsTRp0oS+ffvSsmXLbI/z9/enV69efPHFF2n7aFasWJEFCxYwePBgmjZtSuvWrfnll1+yff6oUaN45JFHCAkJISUlhcjISJ599lmaNWtGSEgI33///U3jnDZtGr169aJjx45UrVr1pseCWWu9evXqNG7cmIcffpi77rrLOUvdOkGeFudyFFmcy7c93m4Pxb79H9PafwmbNlkdTrqzZ2lT5Rh/1/+m34F/QiEW15LFubzD5cuXKVmyJPHx8bRq1YqtW7dSpUoVl7y3UxfnEsIRzsfbWLQlmMd4131q56kqVWJSx33MZBJMl+V1BfTq1YuQkBDuvfdepkyZ4rJkXljSKSpcYu4zh+ll+4FbawQUficiJ+g/4z6ealyCnxfspulrf4CHfICFc3hSu3lGUkMXTmezwYdLyzCRWfDIIxDgfvWIwEZ1ebTBZt5NfgQy7ElZEJ48SkJYq7B/O5LQhdN9ufQ85RJO0dJvt3M2r3CQcW/UJZIw4mdFgH3WYn4FBQURHx8vSV3km9aa+Ph4goKCCvwa7ldVEl7ng3/E8yjvo3r2gDyMMrBKpZ4t6VNxLfPi+vPMJ5/A+PH5fo3q1asTGxtLXFycEyIU3i4oKIjq1asX+PmS0IVTnTiu2fpLBSKIgIeWWB1Orh6ZVJThU8by91l9UY8+mu8NNwIDAzPtQi+EK0mTi3CqOS/HMtT2CSUql4IePawOJ1d3PdOOoIBkNh+sBB7aMSZ8lyR04TTXr8O85aV4hA/NXqGBgVaHlCtVJJCxnWOYy9hCd44K4WqS0IXTrF5+lbrX9tGQgzBmjNXh5Nmwt5qzjh7Ef/YdnDxpdThC5JkkdOE0C976k4dsc6BNG6hXz+pw8qx8wyr0qrmXRbYhMHu21eEIkWeS0IVTnD4N3+8vTX+iYOhQq8PJt7FPlWEuY9Fz5sK1a1aHI0SeSEIXTvHfDy/zYPJySvhfg/Bwq8PJt/smNCG5SAm2xdWGqCirwxEiTyShC4fTGubPuc4oPR+6doWKFa0OKd+Un2L0/X+wkJHw0UdWhyNEnkhCFw63Ywck/3WFe9jqkc0tqYb+oyGRhJG46Xuwr9sthDuThC4cbsG7FxmV+AGqeHHo08fqcArstkalaV71D9bQCzJsEiGEu5KELhzq6lVYvjKQ4SwyyTzDjuueaMSYQD5hBCxYAPad7oVwV5LQhUOtWgWh/ru5jViwb0bsyR58tg7f+rXn7B8psG6d1eEIcVOS0IVDLXr/EiMuvw/ly5sOUQ9XspSid0gsSxkknaPC7UlCFw4TFwdbdwbSh1XQty8UKWJ1SA4x4rlqfKJGmhp6bKzV4QiRo1wTulLqNqXUJqXUQaXUfqXU4/by8kqpr5RSh+3X5ZwfrnBnkZHQo+jXlOQKhIVZHY7DdHywLKeL3s4BWz3Tli6Em8pLDT0ZeEpr3QBoDUxQSjUEngM2aq3rABvt94UPi/g4gSEX3ocyZaBTJ6vDcRh/fxj2wAUWMdyMdrHZrA5JiGzlmtC11qe11rvtty8BB4FbgT7AQvthC4G+zgpSuL8TJ+DAAejKerNnqJc0t6Qa8lxNlvoPQ8fEwLffWh2OENnKVxu6UioYaA78AFTWWp8Gk/SBSjk8Z5xSKlopFS27uHivpUuhf/EvKEIS9O9vdTgO1zTEj6CyRdlBK1i4MPcnCGGBPCd0pVRJIAp4Qmt9Ma/P01rP0VqHaq1DK3rgFHCRN0sWXGPIuZlm3LkXjG7JSikYNMSPCAabzoIrV6wOSYgb5CmhK6UCMcl8sdZ6pb34jFKqqv3xqsBZ54Qo3N3+/XDu1HXu5Tvo1QuKFbM6JKcYOL4CywOHknI5AVauzP0JQrhYXka5KGAecFBrPT3DQ6uBkfbbI4FVjg9PeIKICBhcfDV+aK9sbklVvz5UrqL4jnul2UW4pbzU0O8BhgMdlVJ77JcewDSgi1LqMNDFfl/4GK1hyaJkhpx+09TM77/f6pCcatCYEiz1Gwpffy27GQm3k5dRLlu01kpr3VRrHWK/rNNax2utO2mt69iv/3RFwMK9/PADFL12kRD2mGReooTVITnVwFHFiAocSJL2h0WLrA5HiExkpqgolIgIGFR8NQrgwQetDsfpgoOhzu0pbKAzfPKJ+YoihJuQhC4KzGaDqEgbA068ZWbfeHlzS6pB40qzNGg0HDpkFn8Xwk1IQhcFtn07lPW/RIOUfXDPPWZBLh8QPsif1boXiRSVzlHhViShiwJbsQLCy39t7vTqZW0wLlS1KjRvamMdPcyMKtlEWrgJSeiiQGw2iIzUhB1/yxT4UEIHGDCmJJFlx8L58/D551aHIwQgCV0U0I4dUCrgKo3+2gq1aplB2j6kXz9Yl9iRqwTBf/9rdThCAJLQRQGtWAHh1bebO716mbnxPqRyZWh+J6xX3c066X/KqF1hPUnoIt+0NsuZhJ+ZZQp8rLklVdiQokRWmWD2Gl2xwupwhJCELvJv504oFpBEo8OfmsW47rvP6pAs0a8frLlwL9coIs0uwi1IQhf5tmIFhNf9yUwm6tYNiha1OiRLVKsGjZsFsKFoT9iyBWJirA5J+DhJ6CJftLYn9EsfmwIfbW5JFT7InxXVHjd3liyxNhjh8yShi3yJjoYigTaa7PzYdIT6yOzQnDz4IKyOu5vrBJq1XWQpAGEhSegiXyIjITzkCOr6NWjVygz38GHVq0P9xgF8Xbof/PIL/Pij1SEJHyYJXeRZWnOLbakp8PHmllRh4X5EVrc3uyxebG0wwqdJQhd5tns3+Ptrmn3/oSmQhA5AWBh89ntLkggw7egpKVaHJHyUJHSRZ5GRENb2D9Qfp+HWW6FZM6tDcgs1akDtegF8U3Uw/PGH2fxCCAtIQhd5ktbcUmS1KfDB2aE3ExamWFF1krkjzS7CIpLQRZ7s2WOSevNdH5mC3r2tDcjN9O8Pn8aEkIw/REVBQoLVIQkfJAld5MmKFRDW/TJqV7TZO7RjR6tDciu1asFtwQF8V38cXL4Mq1dbHZLwQZLQRa7S1m4pv9EUdOpkkrrIJDwcIis8Yu5Is4uwgCR0kat9++D6dWixd4EpkNEt2QoLg6hDjUjxC4T//Q/i4qwOSfgYSegiV5GR0L9PMuqr9aagZ09rA3JTd9wBVar5833LxyE5GZYvtzok4WMkoYtcRUVB/5rRpqMvJMRMjxTZCguDyLJ/M3ek2UW4mCR0cVOHDpld1loftScnaW65qbAwiNpbB1vxkrBtGxw9anVIwodIQhc3FRUFD/bT+K2175spwxVvqn59KFvOjx/ufdoUSC1duJAkdHFTkZHQv0UMHD8OlSpBaKjVIbm9sDCILDHS3Fm8WFZgFC4jCV3k6NgxiI2Fe0/bO/d69gQ/+ZPJTVgYREbXRFesBL/+atYcFsIF5NMpcrRyJfTtC/7r7M0t0n6eJ40aQbFiiuiOz5iCRYusDUj4DEnoIkeRkRDW5YLp3AsMhC5drA7JIyhlr6UXHWIKli41G0kL4WS5JnSl1MdKqbNKqX0ZyqYqpX5XSu2xX3o4N0zharGxcPgwdLj8Odhs0L49lCpldVgeIywMIrdUQderbyYYrV9vdUjCB+Slhr4A6J5N+dta6xD7ZZ1jwxJWW7nSDGgJ/J80txREs2aglGJPl7+bgk8+sTYg4RNyTeha62+BP10Qi3AjUVEQ1jfZTGEHSej5lNbsovubglWr4MIFa4MSXq8wbegTlVI/25tkyuV0kFJqnFIqWikVHSdrW3iEP/6An3+GLsW2wMWL0LChWU5Q5EtYGKxYXwbdrj1cu2Y6JYRwooIm9A+A2kAIcBp4K6cDtdZztNahWuvQihUrFvDthCt99hncfz8UXS/NLYXRooVZ1GxfR/vGFzLaRThZgRK61vqM1jpFa20D5gKtHBuWsFJUlKldsmaNKZCEXiBpzS4JPSAoCL75BmJirA5LeLECJXSlVNUMd/sB+3I6VniW+HjYsQO61z5sJsWUKwd33211WB4rLAwiPy9qBvSDLAUgnCovwxYjgG1APaVUrFLqIeA/Sqm9SqmfgQ7AZCfHKVxk1Soz3Lz4Rntzy/33Q0CAtUF5sFatTDfEgXaPmoJFi2QpAOE0uX5StdaDsyme54RYhBuIioJhw4C50tziCH5+Zr/RqDP30LBSJbN8ZXQ0tGxpdWjCC8lMUZHmwgXYsgV63vMXfPcd+PtD9+ymIIj86N8fIlf6w2B73Ug6R4WTSEIXaT7/HNq1g9Lb15sdd9q2NW3oolDatIGzZ+Fw+7GmICJClgIQTiEJXaSJijK1SRnd4lj+/vDggxB1sKEZ03/uXPqELSEcSBK6AODyZfj6a3igZwqss6/kIAndYcLCIDJKwfDhpkCaXYQTSEIXgMnhd98N5Q5tN2MXa9eGevWsDstr3HsvnDgBv907wgxQX70a/vrL6rCEl5GELgD7zkRZm1uUsjQmbxIQAP36QdS2amblSlkKQDiBJHTBlSvw5ZemnVfaz50nLMyew0eMMAWyAqNwMEnogrVrTXPLLZdiYN8+s+75ffdZHZbXad8ejhyBE63CoFgxMzT0t9+sDkt4EUnoguXLYcAAzLhFgG7doEgRS2PyRoGB0KcPrFxf0rS/gHSOCoeShO7jLl2Cr76yLzWSmtB797Y0Jm+W1uwyerQpmD/f7AglhANIQvdxa9aY+UPlAy7C5s1mrnoP2VHQWTp1ggMH4Pd6HaFGDbP64ubNVoclvIQkdB+X1tyyfr2ZvXj33VChgtVhea0iRcwXoE9X+aXX0j/+2NqghNeQhO7DLl40k4n69EGaW1wordll1ChTEBUlY9KFQ0hC92GrV5vBLGVLZZgdKgnd6bp0gZ9+gjPFgk0bTGIiLF1qdVjCC0hC92HLl8PAgcAPP5j1RWrVggYNrA7L6wUFmW6KqChgzBhTKM0uwgEkofuov/4yO6I98ACZm1tkdqhLDBpkr5T36wdlysDOnbB3r9VhCQ8nCd1HrVoFHTpA6dJI+7kFunWD/fvhRFwxGDLEFM6fb21QwuNJQvdRaaNbfvvNZJbSpc0KUsIlihQxSy0sW0Z6s8uiRXD9uqVxCc8mCd0HnT9vdibq3Zv02nn37jI71MUGDzZ7XdCiBTRpYvoxUtfSEaIAJKH7oKgo6NzZLNkizS3WadcO/vgDDv2qpHNUOIQkdB+0eDEMHYoZiP7NN2Z26P33Wx2Wz/H3N6OMIiIwv5DAQPjiCzh50urQhIeShO5jYmPNGOgePTBr5iYlmU0vb7nF6tB80uDBsGQJ6AoVzYgXmw3mzbM6LOGhJKH7mKVLTWdcUBDS3OIGWrY0OXz3buCRR0zh3Llmk24h8kkSuo9ZvNg+Si4pKT2h9+ljaUy+TCkzJj0iArNget26cOqUdI6KApGE7kMOHICzZ01nHN9+a2YXNWgge4dabPBgM3zRplV6Lf3DD60NSngkSeg+ZPFikzz8/YFPPzWFqRstCMs0agTly5uhpIwcCUWLmv7B78+wAAAaD0lEQVSNY8esDk14GEnoPkJr0/k2dCim0fazz8wDktDdQmrnKOXL2xfYAebMsTQm4XkkofuI778321iGhADR0fD771C9upnUIiw3eLBZUvfaNdKbXT7+2F4gRN7kmtCVUh8rpc4qpfZlKCuvlPpKKXXYfl3OuWGKwkode64U6c0tffvKYlxuomZNaNzYbNhN69bQtCnExaX/roTIg7zU0BcA3bOUPQds1FrXATba7ws3lZQEK1akrwEl7efuaeRIWLgQ809WOkdFAeSa0LXW3wJ/ZinuAyy0314I9HVwXMKBvvjCjIa7/Xbg4EE4dMi01d53n9WhiQzCwszE3bg4zNepEiVMwf79VocmPERB29Ara61PA9ivKzkuJOFoCxakb1+ZVjvv3RsCAqwKSWSjVCnza1myBLP65YgR5oF337U0LuE5nN4pqpQap5SKVkpFx8XFOfvtRBZxcWbf0AED7AXS3OLW0ppdAB57zFx/8gn8mfVLshA3KmhCP6OUqgpgvz6b04Fa6zla61CtdWjFihUL+HaioJYsMbW+0qUxiz5FR0Px4tC1q9WhiWx06GD+Ce/di5n01bUrXL0KH31kdWjCAxQ0oa8GRtpvjwRWOSYc4Wjz56dvLp9WO+/e3YxhFG7H3x+GDzeVcgAef9xcz5ol67uIXOVl2GIEsA2op5SKVUo9BEwDuiilDgNd7PeFm9mzx2xm0aGDvWDFCnPdv79lMYncjRhhhpkmJ2P++data75dpU4GEyIHeRnlMlhrXVVrHai1rq61nqe1jtdad9Ja17FfSwOfG5o/37TJ+vlhJhJt2WKWWZTVFd1a/fpw223w1VeYX15qW/qMGZbGJdyfzBT1Utevm/bzkakNY5GR5rpHD/tWRcKdjRqVYc/okSNNJ8iWLfZ1doXIniR0L7VmDTRsCLVr2wuWLTPXacNdhDsbPBjWr7ePSS9VCh56yDwgtXRxE5LQvVSmztATJ2DbNtMR2rOnlWGJPCpb1owsTRvCOHGiaX6JiDDNZ0JkQxK6Fzp5ErZuzVAZT21u6dkTSpa0LC6RP2PHmgUXtQZq1TKd2UlJ8M47Vocm3JQkdC/00Udm3ZYSJewFy5eba2lu8Sh33w1FipjZ/wA8+6y5nj3bbE4iRBaS0L1McrLZY/jhh+0FMTHwww9mMpE0t3gUpWDcuAzLordoAZ06waVL8MEHlsYm3JMkdC+zdi3UqAFNmtgLUsee9+5tkrrwKMOGwbp1cO6cvSC1lj5jBiQmWhaXcE+S0L3M7NnpK68C6aNbwsMtiUcUTvny5n/xokX2gs6doXlzOHMmQ4+pEIYkdC8SEwM7dmTI3QcPwq5dZgxzjx5WhiYKIbXZRWtMO0xqLf3NNyElxdLYhHuRhO5F5s41X9HTlmn573/NdXi4rN3iwdq2NXl882Z7Qf/+ZtTLkSPpI5iEQBK617h+3WxBOW6cvcBmS0/ow4ZZFpcoPKXMMPS0ZdEDAuCZZ8zt114zv2shkITuNVasMDNDGza0F2zZYiYU3Xab7EzkBUaMMMMXjx+3F4webX63+/dLLV2kkYTuBbQ2gx5SV1oF0nvRhg61r84lPFnJkmZJl/fftxcUKQIvvGBuv/KK1NIFIAndK2zfDvHxGYaZJyamD1ccPtyyuIRjTZhgmtUSEuwFqbX0Awekli4ASeheYcYMs8Kqv7+9YM0auHDBDG9La4MRnq52bWjd2r7nKEgtXdxAErqHi401q/KlbQIN6c0t0hnqdR57DGbOtA9hBKmli0wkoXu49983ebtMGXvBH3+YqYX+/mYNVuFVunQx63Olre+SsZb+0kuyTZ2Pk4TuwRISzEJcqRvaAGYzyuRk6NULqla1LDbhHErB5MnwxhsZCkePNuPSDx3KsCuG8EWS0D3Y/PlmRb46dewFWqfvDv+3v1kWl3CuESPMxkV799oLihSBf/zD3H75ZbhyxbLYhLUkoXuopCRTS/u//8tQ+N13cPgwVKtmNhcWXikoCCZNylJLHzDArMZ4+rTsauTDJKF7qOXLITjYjHpIk1o7Hz3azCYUXuvRR83KmidO2Av8/OA//zG3//3vDMszCl8iCd0DaQ3TpsFzz2Uo/Ouv9LHnY8ZYEpdwnbJlza/57bczFHbsaL6ZXbyY3gQjfIokdA+0bp2pgHfrlqFwyRIzoahTJ9NBJrzeE0+YFXT//DND4bRppuf0vffg118ti01YQxK6B0qtnStlL9A6fQcb6Qz1GbfeajaSnjkzQ2GzZqbqnpRkMn7agHXhCyShe5ivvzZ7G/Tvn6Hwm29g3z6oUgUefNCy2ITrPf88zJoF589nKPznP80a+F98YRrahc+QhO5BtDZzR15+OUufZ+q6qg8/bIawCZ9Ruzb07ZulLb1SJbMUAJha+rVrlsQmXE8SugdZv94swjVoUIbCEyfgs89Mhk9bDF34khdfNDOGM7WlT5hg1vE5ehSmT7csNuFaktA9RGrtfOrUDItwAXz4oVmUKSzMjD8XPic42LS0vfVWhsLAwPTG9ddfzzC+UXgzSegeYt06M9U/017PiYlm3znIMv9f+JoXXjD/2zMNP+/UyfzBJCSYgevSQer1CpXQlVIxSqm9Sqk9SqloRwUlMrPZYMoUUzvPtFfFkiXmE9y8uVkDQPismjVh4ED417+yPDBjhlm5bd06WLbMktiE6ziiht5Bax2itQ51wGuJbCxZYvo6Mw1gsdnSZwY++WSGMYzCV738shmXfuxYhsKqVeHNN83tSZNMJ4zwWtLk4uauXjVD0958M0vOXr3arK5Xo4apmgmfV7my2Ybw+eezPPDQQ9C+PcTFmX/+wmsVNqFrYL1SapdSSoZYOMGMGdCyJbRtm6FQa7NeB5gPaGCgJbEJ9/Pkk2Z/8B07MhQqBXPmQNGiZnllGZvutQqb0O/RWt8J3A9MUErdsL28UmqcUipaKRUdFxdXyLfzLXFxpmY+bVqWB7ZsMRuJli8vM0NFJiVKwKuvwtNPZ+kDrVPHjHYBM5P07FlL4hPOVaiErrU+Zb8+C3wKtMrmmDla61CtdWjFihUL83Y+56WXYOjQDOudp0qtnU+caD7BQmQwcqTZUjYqKssDkyebppezZ2HsWBn14oUKnNCVUiWUUqVSbwNdgX2OCszXRUeb+UJTp2Z5YOdO85W5eHGT0IXIwt/fLAcweTJcupTlgYULzaiX1ath3jzLYhTOUZgaemVgi1LqJ2AHsFZr/T/HhOXbUlLMsOFp06BcuSwPvvyyuX7sMZBvPCIH995rVtN99dUsD9SoYaaVglkW4NAhl8cmnEdpF37tCg0N1dHRMlw9Nx98ABERZs2tTCNbtm2DNm2gZEn47TeoUMGyGIX7O3MGGjeGTZvMdSZDhpg/ssaNTX+MNN25NaXUrrwMDZdhi27m7FlTCX///WyGlr/0krl+/HFJ5iJXlSubJrvx4820hUxmz4Z69cwqnTKL1GtIQnczEyeaHeRuqFFt3gwbNphlUZ96yorQhAd65BGz2OLs2VkeKFXK9JoWLw6LFplhjcLjSUJ3I8uXm53cU1c+TWOzpSfxp5/OpmFdiOz5+8OCBWbpiN9+y/Jgo0bpiXzSJNP0IjyaJHQ3ceaM+UwtXGh2dc/kv/+F3bvNFjUy00/kU4MG8MwzZvj5DU0vQ4eapXavX4c+fSAmxooQhYNIQncDWpt2ztGjoVXWkfxXrqTP5f7Xv6TzShTIU0+ZZSRSB7hk8vbb0KWL6cDp3dtsMi08kiR0N/DRR3D4cPqIxEzefBN+/x1atDC1KSEKIHUI+iuvmH7QTAIDTXtfgwbmwYEDITnZkjhF4UhCt9jevaYCvnx5Nk0tR46kr4c6fXqWtXOFyJ969Uz9YMAA88Uvk7JlYc0aM3rqf/8zC3rd0D4j3J1kCAtduWIqQ2++CfXrZ3lQazOc7No1GDEC7rthmRwh8m3kSLPYW7b7odSqZZJ6iRJmEa/HH5fhjB5GErpFUvN1aKj5kN1gyRIzTLF8+fT1rIVwgPfeM3PUFizI5sG77oJVq8wC/LNmmeExwmNIQrfIW2+Z5soPPsjmwbg4sxAHmGQuU/yFA5UsaYagP/OMSew36NTJtAH6+8M//mEmtElN3SNIQrfAunWmSfyzz7IZtKK1WQkvLs6sjDdqlAURCm/XsKGpoffvDydPZnNAnz5mwpG/P7z2Wjbr8Qp3JAndxfbvNzk6MtKsk3SD+fPNV97Spc0nTraWE07So4eZ1vDAA3D5cjYHDB5s9iENDDQ1kPHjzcpxwm1JQnehmBjo3h3eecessXWDo0dNRxSY9suaNV0ZnvBBTz1l+nH69TP97zfo3998lQwKgg8/hLCwbIbICHchCd1Fzp6Frl1Nu+WQIdkckJBgPjyXL0N4OAwb5vIYhe9RyuTpMmXMNIdsK+A9esAXX5ihjZ99Bu3awenTLo9V5E4SugvEx0O3bjBoUA7DxbQ2qyj99BPccYdZX0OaWoSL+PvD4sVml6OxY3MYft6+velBrVULdu0yo2FkKWy3Iwndyc6cMZ+F7t2zWXQr1axZpgOqeHH49FNTExLChYoWNX96x46ZYbTZThStX98s4NWmjelJveceM0xLOkvdhiR0J/r9d/PtNCwM/vnPHCrdq1aZnWPAbAl2w7q5QrhGyZJmBNbZs+bb5PXr2RxUsSJ8/bXpIL1+3VwPHSrrv7gJSehO8vPPcPfdZoW7l1/OIZlv22Y+OTabqb4PGuTyOIXIqHhxs91oSoppOj9/PpuDihY1s5MiIsy424gIaNIENm50ebwiM0noTvDFF9C5M7zxhukEzdbu3dCzJyQmmoZLmZEn3ETRorBihfmy2KaNWVIoW4MGmfb00FA4ccL80Y8fn2VnauFKktAdyGYza2mNGWMGAwwcmMOB0dFmNt7582YCR7b7zQlhnYAAM7x20iRo2xa+/DKHA+vVM980X3vNjFf/4ANTtnixtK1bQBK6g8TFma+o69bBzp05jDMHs/Nz587w119m8O/y5ebTI4QbevRRWLrULL743HOQlJTNQQEB8OKL5g+/VSszpHHYMLOg3O7dLo/Zl0lCd4BVqyAkxFw2bYLq1XM4cNEis5HAhQump3TZMrMIkhBurH17+PFH0y90331w8GAOBzZrZmrrH39sOk+3bDHr+IeFmSnSwukkoRdCXJyZHf3006ZfaNq0HCrbSUmmMX3ECHN78mRT7QkMdHnMQhRExYpmZd1hw+Dee00LS7ajYPz8zNZbv/5qpqEGBZmVwJo0MR8WqbE7lST0Arh2zayW2LCh2ebzp59uslz58ePmwTfeMDM43n3XrIvh7+/SmIUoLD8/s/3ojz/Cjh2mQr5mTQ5N5WXLmpVCjx6FiRNNTWfpUlNjb9fOdDLJujCOp7V22aVFixbakyUna71kida1a2vds6fWBw/mcvDMmVqXKqU1aF29utZbtrgsViGcyWbT+vPPtW7QQOv27bXevj2XJxw/rvWTT2pdurT5PKR+Jp5/XutDh1wSsycDonUecqwk9Dy4dk3rjz7S+o47tL7nHq2/+iqXJ2zdqnWrVul/uP36aX3unEtiFcKVkpK0nj1b69tu07pDB62//NIk+xxduKD122+bWlHq5wO0vvtuU/7bb64K3aNIQneAw4e1fvZZrStX1rprV62/+SaXJ/z4o9a9eqX/kd56q9YrV7okViGsdP261gsXat2okdZNmpgvp/HxN3mCzWY+UKNHa12iRObk3ry51q+8ovX335sXFnlO6Eq7cKxoaGiojnbzBX3++MOsabF8uemYHzHCzPupVy+HJyQnw+efw4wZZkgimNlzkyfD3/9u1jUXwkfYbGak17x5Zgjv/febxUO7dctmM5dUly+bz9Cnn5pZeRkXZy9VyvRBdexoxgKHhGSzm7r3U0rt0lqH5nqcryf05GTT8b5xo/lb2rvXTODs39+MKy9aNJsnpaSY4VkREWZKXVycKS9VCv72N3j2Wahc2aXnIYS7iY83FaOVK00naseOJsF37Ai1a+cwly4x0XwY1641a8YcOpT58YAAaNrU7HTdogU0amRGJ3j5gnYuSehKqe7ADMAf+EhrPe1mx7tDQo+LM7OVo6PNH9l335lx4506mSHinTtnk8STk80f1pYtZuPmr7+GP/9Mf7xuXTPlefRoqZELkY0//zQjYr76ytTglTLj21u1Mnm5WbMcavCxseYJmzebD+z+/dkPq6la1ST2+vXh9tvN5jDBweZyyy0ePxPb6QldKeUP/Ap0AWKBncBgrfWBnJ7jqoR+7RqcOmXWoDh0yAyJ/fVXMyHiwgXzB9SihVmCol27DJXpq1fht9/MUKujR+HwYTNGa88e81hGt99uJkwMGgTNm3v8H4wQrqK1+Wx+842pWO3aZfJ0cDA0aGDqR3XrQp06piZfqVKGUb6XLpmv1Dt2mPHCBw6YD3ZiYs5vWLy4qbVVrmwulSql365c2dTuy5RJvy5d2u1mb7siod8NTNVad7Pf/z8ArfW/cnpOQRN6wolz/H48mb/Oa86fN7Pmz/+lzPUFP+L+9OPU2QBOxQXye1wRLl7xp2r5a9xRNYG6lS9Q95Z46pY9S71Sp6gVdAq/hMtmuc+4OLNWaOr1X3/lHERwsPma16mTqcbXrp3v8xBCZC8pyeTm1ArY4cPm+tgxs+RRpUpmzke1aqYyXr48lCtnLuXLplDu2hnKnTtMyT+OUOxMDMVPHaH474cJPH7E1OLyq0QJk9zLlIFixcylePH021kvQUFm1ndgYPol4/0iRUxHXMOGBfr55DWhF+bf0K1Axv3CY4G7CvF6OdrYZgpP/P405ThPOc5Tlr8yXdfmHLfyO9U4RTVOUYFz+J3VcDafbxQQYBJ37drpl6ZNTQ38lluccWpCCEzOa9bMXLK6ft0MVjh1yuwxcPq0acKJjTV9XufP+3P+fDXOn6/GlSvtuHrV7OiYkGC+DRQvpSlWJIXigUkUUUkE6CQCbdcJsF0jIOUagSnmOiAlkYDkawQmXyXgShL+V1JQpzQKc/HDlnY7b5ckFNfT7s98ageBb+ZY33WIwiT07NoYbqjuK6XGAeMAamS7zX3uetc+QG99n/ne5e9vEm+2t4uBf13wb2DKAwPNf9rsLqVKmfnMlSqlX5crZ6bDCSHcRpEiUKOGueRXUhJcvapISAggISGApKRiJCWZbrHkZLK/naRJvpxI8oUr6ISr6OvX0deS0NeuY0u8jr6eZC7XMty+noROSUEn29Kvk1PQKen3VeNGjv/hZFGYhB4L3JbhfnXgVNaDtNZzgDlgmlwK9E6pwwGFECIfUls88jdWQQHF7BfPUpjq6E6gjlLqdqVUEWAQsNoxYQkhhMivAtfQtdbJSqmJwJeYYYsfa61ljUwhhLBIocbmaK3XAescFIsQQohCkB5AIYTwEpLQhRDCS0hCF0IILyEJXQghvIQkdCGE8BIuXT5XKRUHHC/g0ysA5xwYjpXkXNyPt5wHyLm4q8KcS02tdcXcDnJpQi8MpVR0Xhan8QRyLu7HW84D5FzclSvORZpchBDCS0hCF0IIL+FJCX2O1QE4kJyL+/GW8wA5F3fl9HPxmDZ0IYQQN+dJNXQhhBA34RYJXSnVXSl1SCl1RCn1XDaP36eU2q2USlZKhWV5bKRS6rD9MtJ1Ud+ooOehlApRSm1TSu1XSv2slBro2shvVJjfif3x0kqp35VSs1wTcc4K+fdVQym1Xil1UCl1QCkV7Kq4s1PIc/mP/W/soFJqplLWbYSbh/N40v7z/lkptVEpVTPDY27zmbfHU6BzccrnXmtt6QWz9O5RoBZQBPgJaJjlmGCgKfAJEJahvDxwzH5dzn67nAeeR12gjv12NeA0UNYTfycZHp8BLAFmeerfl/2xzUAX++2SQHFPPBegDbDV/hr+wDagvRufR4fUnzXwKLDMftttPvMOOBeHf+7doYbeCjiitT6mtb4OLAX6ZDxAax2jtf4ZsGV5bjfgK631n1rr88BXQHdXBJ2NAp+H1vpXrfVh++1TmN1Qc51E4ESF+Z2glGoBVAbWuyLYXBT4XJRSDYEArfVX9uMua60TXBR3dgrze9FAECbpFAUCgTPODzlbeTmPTRl+1tsxO6KBe33moRDn4ozPvTsk9Ow2m77VBc91NIfEopRqhfnQHXVQXAVR4HNRSvkBbwF/d0JcBVGY30td4C+l1Eql1I9KqTeUUv4OjzDvCnwuWuttwCZMLfA08KXW+qDDI8yb/J7HQ8AXBXyusxXmXNI46nNfqA0uHCRPm0074bmOVuhYlFJVgUXASK31DTVfFyrMuYwH1mmtT1rYRJtRYc4lALgXaA6cAJYBo4B5Doks/wp8LkqpO4AGpNd0v1JK3ae1/tZRweVDns9DKTUMCAXa5fe5LlKYc0ktd9jn3h1q6HnabNoJz3W0QsWilCoNrAVe1Fpvd3Bs+VWYc7kbmKiUigHeBEYopaY5Nrx8Kezf14/2r9PJwGfAnQ6OLz8Kcy79gO32ZqPLmFpiawfHl1d5Og+lVGfgBeABrfW1/DzXhQpzLo7/3FvVmZChwyAA07FxO+mdCo1yOHYBN3aK/obpHClnv13eA8+jCLAReMLq30dhzyXLY6OwvlO0ML8Xf/vxFe335wMTPPRcBgIb7K8RaP976+2u54H5VnQUe6dhhnK3+cw74Fwc/rm35IeQzQ+lB/Cr/aRfsJe9ivlvBtAS85/wChAP7M/w3DHAEftltCeeBzAMSAL2ZLiEeOK5ZHkNyxO6A/6+ugA/A3vtSbKIJ54L5p/TbOAgcACY7ubnsQHTaZv6eVid4blu85kvzLk443MvM0WFEMJLuEMbuhBCCAeQhC6EEF5CEroQQngJSehCCOElJKELIYSXkIQuhBBeQhK6EEJ4CUnoQgjhJf4fkZu123DTezIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f9c42bb47f0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Zum Vergleich: Normalverteilung mit gleichem Mittelwert und gleichem Maximum.\n",
    "\n",
    "# Man sieht, dass die rote Kurve leicht schief, d.h. nicht symmetrisch ist\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# sigma manuell so angepasst, dass Maxima übereinstimmen \n",
    "ff_n = norm_dist(xx, maxval, 0.0153)\n",
    "\n",
    "plt.plot(xx, ff_ln, \"r-\", lw=2, label=\"Lognormal-Verteilung\")\n",
    "plt.plot(xx, ff_n, \"b-\", lw=1., label=\"Normalverteilung\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Tabelle aus dem Buch Seite 175\n",
    "\n",
    "Lebensnutzen: $\\Delta = \\frac{1}{2}\\frac{\\sigma_\\varepsilon^2}{1- \\alpha^2}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "def delta(sigma, alpha):\n",
    "    return - 1.0/2* sigma**2/(1- alpha**2)\n",
    "\n",
    "sigma_values = np.r_[0.1, 0.2, 0.3, 0.4]\n",
    "delta_values = delta(sigma=sigma_values, alpha=0.3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.1, 0.2, 0.3, 0.4])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "___\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([-0.00549451, -0.02197802, -0.04945055, -0.08791209])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "___\n"
     ]
    }
   ],
   "source": [
    "sigma_values ##\n",
    "delta_values ##"
   ]
  }
 ],
 "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.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}