MALA-HMC.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Metropolis-adjusted Langevin algorithm (MALA)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ Y \\sim \\mathcal{N}(X_n + h \\frac{\\nabla\\pi(X_n)}{\\pi(X_n)}, 2hI) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ \\alpha(X_n, Y) = \\frac{\\pi(Y)}{\\pi(X_n)} \\exp \\Big\\{ \\frac{1}{4h}\\Big( \\big|Y - X_n - h \\frac{\\nabla\\pi(X_n)}{\\pi(X_n)}\\big|^2 - \\big|X_n - Y - h \\frac{\\nabla\\pi(Y)}{\\pi(Y)}\\big|^2 \\Big) \\Big\\}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mala(pi, grad_pi, x_0, h, n_of_steps):\n",
    "    #x = x_0\n",
    "    x = np.zeros((n_of_steps, x_0.size))\n",
    "    x[0] = x_0\n",
    "    #sumx = np.zeros(x_0.shape)\n",
    "    scale = np.sqrt(2*h)\n",
    "    for i in range(1, n_of_steps):\n",
    "        y = np.random.normal(x[i-1]+(h/pi(x[i-1]))*grad_pi(x[i-1]), scale)\n",
    "        sqsum1 = np.sum((y-x[i-1]-(h/pi(x[i-1]))*grad_pi(x[i-1]))**2)\n",
    "        sqsum2 = np.sum((x[i-1]-y-(h/pi(y))*grad_pi(y))**2)\n",
    "        alpha_cand = (pi(y)/pi(x[i-1])) * np.exp( (sqsum1 - sqsum2) / (4*h) )\n",
    "        alpha = min(alpha_cand, 1)\n",
    "        if np.random.rand() < alpha:\n",
    "            x[i] = y\n",
    "        else:\n",
    "            x[i] = x[i-1]\n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [],
   "source": [
    "def norm(x):\n",
    "    return (1/np.sqrt(2*np.pi))*np.exp(-x**2/2)\n",
    "\n",
    "def grad_norm(x):\n",
    "    return -(x/np.sqrt(2*np.pi))*np.exp(-x**2/2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [],
   "source": [
    "mala_samples = mala(norm, grad_norm, np.array([0]), 0.01, 10**5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x168e747b588>]"
      ]
     },
     "execution_count": 80,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl8VOXZ//HPNZOEfdW4sRhQVhFZIoh7KyK4gK1tBbXV1pbaau3zWNvi3upjH6vW+vSnttqqtW4U10YDsoOCIiSsBggERAiLhn1PMnOu3x9n0DEEMoGZ3LNc79drXplz5pyZbyBzzZn73Oe+RVUxxhiTGQKuAxhjjGk4VvSNMSaDWNE3xpgMYkXfGGMyiBV9Y4zJIFb0jTEmg1jRN8aYDGJF3xhjMogVfWOMySBZrgPUdOyxx2peXp7rGMYYk1KKi4s3q2puXdslXdHPy8ujqKjIdQxjjEkpIvJZLNtZ844xxmQQK/rGGJNBrOgbY0wGsaJvjDEZJKaiLyJDRaRURMpEZMxhtvuOiKiI5EetuyOyX6mIXBKP0MYYY45Mnb13RCQIPAlcDJQD80SkQFWX1tiuBXAr8HHUup7ASOA04CRgioh0VdVw/H4FY4wxsYrlSH8AUKaqq1W1ChgLjKhluweAh4H9UetGAGNVtVJVPwXKIs9njDHGgViKfjtgXdRyeWTdl0SkL9BBVd+t777GOOfZF0+TOWK5OEtqWfflxLoiEgD+DNxQ332jnmM0MBqgY8eOMUQy5ijs3AhznoJV02DXRti7Bdp2hs4XQrdL4dTBILX96RqT+mIp+uVAh6jl9sCGqOUWQC9ghvhvlBOAAhEZHsO+AKjqM8AzAPn5+TZTu0mMyl0w6W5Y+Ap4Iej8DegwAJoeA5uWwOJxUPQcnHoxXPYnaHOy68TGxF0sRX8e0EVEOgHr8U/MXnPgQVXdARx7YFlEZgC3q2qRiOwDXhGRx/BP5HYB5sYvvjEx2rYGXh0FFaXQ/wY4+xfkPfxVX4Q1D90N4WqY+3eY9j/w1Flw+eNwxtXkjSmM2u6yhs9uTBzV2aavqiHgFmAisAwYp6olInJ/5Gj+cPuWAOOApcB7wM3Wc8c0uLVz2Pr4OWz//DOuqfwtXP4YtO108HbBbBj0c7j5YzipH96bP+U3d93e8HmNSSBRTa7WlPz8fLUB10zcVJTCPy5m9b4m/Kj616zRE2vd7KAj+Op9zLx/CBcEF3Nn9Y28Er6o9u2MSRIiUqyq+XVtZ1fkmvS1ZzO88j3IyuH7VXccsuDXKrsJo6tvY2q4L3/IfpYLAosSl9OYBpR0QysbExehKhh7LezaBDcUsv6Jz2PaLbr9HnL4efUveVvu4bHsp7is8g+JyWpMA7IjfZOe3n8E1s2BK5+C9nV+4yVvTGGNgu+rJIebq39JY6r4S84TEA4lIq0xDcaKvkk7l93xBKGZj/JG+DzoddVRP99qPYk7qn/MgEApzPxjHBIa444VfZNeQlU8mv00W2jJ76u/f8gj+Poq8M7hzfC5MOvPsHllHIIa44YVfZNeZv2ZHoG13Fl9IztpHten/kP1tZDdFMbfDknW682YWFnRN+ljx3qY9RjvhM9iqtc/7k+/mVbwzbth9QwoeSvuz29MQ7Cib9LH9D+AevwxNCpxr3HmjXBCb5h4J1TuTtzrGJMgVvRNevi8BBa+DANGU665CXuZvDvf49uffdsfqG3uMwl7HWMSxYq+SQ+T74PGLeG8XyX8peZrV+gyBD78C+zfmfDXMyaerOiblJY3ppCr73gUyibz4K7LyLv/o4Z54W/cCfu2wZy/NszrGRMnVvRNyrs5620qtCX/Cg9puBc9qS90vxw+etIv/sakCCv6JqX1ktWcH1zCs6FLqSSnYV/8wjugcodf+I1JEVb0TUr7eVYBO7UpL4UHN/yLn9ALelzhj8FftafhX9+YI2BF36SuihUMDczjX+GL2U1TNxnOvhX2b4cFL7t5fWPqyUbZNKlr9v9RRRbPh4Y2+Et/bTatU8/059w980YIBBs8izH1YUf6JjXt2QxLxvFa+AK20MptlkG3wLZPoXS82xzGxCCmoi8iQ0WkVETKRGRMLY/fJCJLRGShiMwSkZ6R9Xkisi+yfqGI/C3ev4DJUPNfgHAV/wxf4jqJ367f+mT48AnXSYypU51FX0SCwJPAMKAnMOpAUY/yiqqerqp9gIeBx6IeW6WqfSK3m+IV3GSuU8YUUD7lKWaHT2OVtnMdx2/SOevn/vj964tdpzHmsGI50h8AlKnqalWtAsYCI6I3UNXoyxKbATYEoUmYiwLzaS+bG7Zffl36XAPZzaDoeeCrSVniMayzMfEUS9FvB6yLWi6PrPsaEblZRFbhH+nfGvVQJxFZICIzReS8o0prDPCD4CTW6zFM8fq5jgJECvzvPoDTr4JP3oD9O1xHMuaQYin6Usu6g47kVfVJVT0F+C1wd2T1RqCjqvYFbgNeEZGWB72AyGgRKRKRooqKitjTm8xTsYJzgyW8HBpMmOTqKXPFR12hei93P3Cf6yjGHFIsRb8c6BC13B7YcJjtxwJXAqhqpapuidwvBlYBXWvuoKrPqGq+qubn5iZuhESTBua/QJUG+Xf4QtdJDrJEO7PEy+Pa4BSshdMkq1iK/jygi4h0EpEcYCRQEL2BiHSJWrwMWBlZnxs5EYyIdAa6AKvjEdxkoFAVLBrLFK+/+26ah/BK+CJ6BNbRT2xKRZOc6iz6qhoCbgEmAsuAcapaIiL3i8jwyGa3iEiJiCzEb8a5PrL+fGCxiCwCXgduUtWtcf8tTGZYORH2bmZc+ALXSQ6pIHw2u7Ux12RNcx3FmFrFdEWuqo4HxtdYd2/U/V8eYr83gDeOJqDJbNG9X/6R/Sd6BdrwgdfbYaLD20MTCsKDGBH8kHu5gb00dh3JmK+xK3JNSjiObXwjsJA3wucl3Qncmt4In08zqWRYYK7rKMYcxIq+SQlXBT8gKMprSdy0c0CxdmWNdzzfDn7gOooxB7Gib1KA8p3gTOZ63VijJ7oOEwPhzfB5DAos5SQ220VaJqlY0TdJr7es5pTARt4In+86Ssze9M4lIMq3grNcRzHma6zom6R3ZXA2lZrFhPAA11FiVq7H8bHXnauC72N99k0ysaJvklqQMFcEP2K615edNHMdp15eD59P58Am+kqZ6yjGfMmKvklqZwdKyJUdvB0+x3WUepsQHkClZjM8+KHrKMZ8yYq+SWpXBmezU5sy3evjOkq97aYp07w+XBb8mACe6zjGAFb0TTKr2sslgXn+ETM5rtMckXfCgzhOtjMwsMx1FGMAK/omma14j+ayn7e91GvaOWCq14/d2pjhAWviMcnBJkY3yavkTb7Q1nzs9XCd5IhVksMkL59hwbl0GfMfqiNvuTUPXeY4mclUdqRvklPlLlg5mcLwQLwU/zMtCA+itezhvMBi11GMSfF3k0lfKyZCaD+F4YGukxy12d7pbNPm1ovHJAUr+iY5lbwFLU6kWA+acyflVJPFhPCZDA7MpxFVruOYDGdF3ySfyl1ULpvIc9vOQNPkT3S8dxbNZT/nBZa4jmIyXHq8o0x6KX2PRlLNu+GzXCeJmzleD7Zpcy4Nfuw6islwVvRN8il5iw3algV6quskcRMii0nhfAYHismh2nUck8FiKvoiMlRESkWkTETG1PL4TSKyREQWisgsEekZ9dgdkf1KReSSeIY3aahyF5RNYUJ4YNo07RwwwRtAS9nHOYFPXEcxGazOd1VkYvMngWFAT2BUdFGPeEVVT1fVPsDDwGORfXviT6R+GjAUeOrAROnG1GrlJAhXMiF8puskcTfb68VObcqlAWviMe7Ecig1AChT1dWqWgWMBUZEb6CqO6MWm/HVWLIjgLGqWqmqnwJlkeczpnbL3oFmxzE/DXrt1FRNFpO9/gwJFkHIevEYN2Ip+u2AdVHL5ZF1XyMiN4vIKvwj/Vvrs68xAFTvgxWToMflKX9B1qGMDw+glezl+nsftdm0jBOxvLOklnUHzQqhqk+q6inAb4G767OviIwWkSIRKaqoqIghkklHP77vEajew3UfHu86SsLM8k5ntzbmksA811FMhoql6JcDHaKW2wMbDrP9WODK+uyrqs+oar6q5ufm5sYQyaSjocF57NCmzEnhsXbqUkkOM7wzuDhYbMMtGydiKfrzgC4i0klEcvBPzBZEbyAiXaIWLwNWRu4XACNFpJGIdAK6AHOPPrZJO+FqBgeKmeL1J5Tm4wBOCp9Jruygr6z8ctJ0a+oxDaXOd5eqhkTkFmAiEASeU9USEbkfKFLVAuAWERkMVAPbgOsj+5aIyDhgKRACblbVcIJ+F5PK1syitezhvTTstVPTdK8PVRrkkmARxaFuruOYDBPTIZWqjgfG11h3b9T9Xx5m3weBB480oMkQy99lrzbifa+36yQJt4umfOj14pLAPP7ANdR+6suYxEjPLhImtXgeLB/PTK93ys6QVV8TvXxODnxBd1lX98bGxJEVfePc8LuegF0bmBzu7zpKg5kczsdTYUigyHUUk2Gs6BvnLg4WE9IA07y+rqM0mM20oli7cEnwq66bdkLXNAQr+sa5IYEi5nrd2U4L11Ea1KRwPqcFPqMddm2KaThW9I1bW1bRLVDOZC9zmnYOOPA7Dw7Od5zEZBIr+satUr9T2GQv33GQhrdGT2Sl146LrV3fNCAr+sat5YWUeCdTrpl5JfZkrz8DA8tpyR7XUUyGsKJv3NmzGdbOycimnQOmhPuRLWEuDCx0HcVkCCv6xom8MYXc/oeHAWVyOPOadg5YoKdSoa24OFjsOorJEFb0jTODA/PZoG0p0ZNdR3FGCTAl3I8LAovIJuQ6jskAVvSNE42o4vzAYqaE+5PpwxBM9vrTUvZxVmCp6ygmA1jRN04MCpTQVCqZ4vVzHcW52V4v9mojBgesiccknhV948TFgfns1sbM8WpOt5x5KslhlteLi4ILqGWOIWPiyoq+aXiex0XB+bzv9aaKbNdpksJkrz/tZTM9ZK3rKCbNWdE3DW/jQk6QbUwJW9POAdPDffFUuChgV+eaxLKibxpe6QTCKkz3+rhOkjQ204qFegqDreumSTAr+qbhlU6gSLuxjZaukySVKeF+9Amshl2bXEcxaSymoi8iQ0WkVETKRGRMLY/fJiJLRWSxiEwVkZOjHguLyMLIraDmvibDbF8Hny9hajhzhlGO1ZQDVyaveM9tEJPW6iz6IhIEngSGAT2BUSJSs8vFAiBfVXsDrwMPRz22T1X7RG7D45TbpKi7H/kTAFOtq+ZBVmh71nm5UGpF3yROLEf6A4AyVV2tqlXAWGBE9AaqOl1V90YW5wDt4xvTpIvBgfl86h3PKj3JdZQkJP51C6unQ9Xeujc35gjEUvTbAdETeZZH1h3KjcCEqOXGIlIkInNE5MojyGjSReVuBgVKIkf5mX0V7qFM9fpBaD83/u5Rm0XLJEQsRb+2d2etV5CIyHVAPvBI1OqOqpoPXAM8LiKn1LLf6MgHQ1FFhc0ilLZWTaORhKxp5zA+9nqwS5tY102TMLEU/XKgQ9Rye2BDzY1EZDBwFzBcVSsPrFfVDZGfq4EZwEFn8FT1GVXNV9X83NzMHFc9I6x4j53alHleN9dJklY1Wbzvnc5FwQUInus4Jg3FUvTnAV1EpJOI5AAjga/1whGRvsDT+AX/i6j1bUSkUeT+scA5gI0qlYm8MKyYyAzvDEJkuU6T1KaG+3G8bKeXrHEdxaShOt99qhoSkVuAiUAQeE5VS0TkfqBIVQvwm3OaA6+JCMDaSE+dHsDTIuLhf8A8pKpW9DPR+mLYu5kp4atdJ0l6070+hFVs7lyTEDEdcqnqeGB8jXX3Rt0ffIj9PgROP5qAJk2UTgAJMsM7w3WSpLeNlhRrVxt10ySEXZFrGsaK9+Dks9lJc9dJUsLUcD9OC3wGO8pdRzFpxoq+Sbxtn8EXS6HrUNdJUsaX8wzY1bkmzqzom8Q7ULi6DXObI4Ws0pNY4x1vV+eauLOibxKvdAIc2xWOOegSDXNI4l/P8On7ULXHdRiTRqzom8Tav5OqVe/zt01d7QrTepri9YNwJaya7jqKSSNW9E1irZpGjoSZahOm1Ns8rxs0auV/UzImTqzom8Ra8R7btDnztYvrJCknRBYFe3tSsaAAPLs618SHFX2TOJGrcKd5fQgTdJ0mJU0J9yNXdvoXtxkTB1b0TeKsmwv7tlrTzlGY4Z1BSAOwwpp4THxY0TcJ87e/P0GVBnnf6+06SsraSXPmed2t66aJGyv6JmEuDhQzx+vJbpq6jpLSpnj94IsS2LbGdRSTBqzom8TYXMYpgY1fXVlqjtjkA3Pn2tG+iQMr+iYxIm3Q1p5/9Nbq8ZDbHUrH172xMXWwom8So3QCy7yOrMcmxYmLbsPgs9mwb7vrJCbFWdE38bd3K6ydw2Rr2omfbpeCF4KyKa6TmBRnRd/E38pJoGGmhPu7TpI+2vWHZrl2da45alb0TVzljSlk/OvPsknbsEQ7uY6TPgJB6HoJrJwM4WrXaUwKi6noi8hQESkVkTIRGVPL47eJyFIRWSwiU0Xk5KjHrheRlZHb9fEMb5JPI6q4ILCIKeF+qB1TxFe3S6Fyh9+2b8wRqvNdKSJB4ElgGNATGCUiPWtstgDIV9XewOvAw5F92wL3AQOBAcB9ItImfvFNshkUKKGZVDLZy3cdJa3kjSmkxz/3s1+zYbn14jFHLpZDsQFAmaquVtUqYCwwInoDVZ2uqnsji3OA9pH7lwCTVXWrqm4DJgM2fVIaGxIoZpc24SOv5nGBOVr7aMwH3ul+101V13FMioql6LcD1kUtl0fWHcqNwIGzTfXd16Qyz2NwcD4zvd5Uke06TVqa5OXDjnWwaYnrKCZFxVL0pZZ1tR5miMh1QD7wSH32FZHRIlIkIkUVFRUxRDJJacN8jpPtTLZeOwkzNdwPJADLbUIac2RiKfrlQIeo5fbAhpobichg4C5guKpW1mdfVX1GVfNVNT831y7mSVnLCwlpgOleH9dJ0tZWWkKHs6zomyMWS9GfB3QRkU4ikgOMBAqiNxCRvsDT+AX/i6iHJgJDRKRN5ATukMg6k46WF/Kx14OdNHedJL11vxQ+X2IDsJkjUmfRV9UQcAt+sV4GjFPVEhG5X0SGRzZ7BGgOvCYiC0WkILLvVuAB/A+OecD9kXUm3WxeCZtLmWi9dhKv26X+T7tQyxyBrFg2UtXxwPga6+6Nuj/4MPs+Bzx3pAFNilj+LgCTw1b0E+6YU+C4nn4Tz1k/c53GpBi7esbExYJJL7HI68xGjnEdJe3ljSnkLxu6+hdp7bUvzqZ+rOibo7dzA30DZUy0o/wGMymcD+pZE4+pNyv65uhFepJM9M50HCRzfKKdoFUHWPaO6ygmxVjRN0dv+bus8k5kldp1dw1HoMcVsGoaVO5yHcakECv65ujs2wZrZtlRvgs9roBwpT/ypjExsqJvjs6KieCFrD3fhQ4D/TH2rYnH1IMVfXN0lv4HWrZnkZ7iOknmCQSh+2X+pDXV+12nMSnCir45cpW7oGwq9BxO7cMsmYTrcQVU7YbVM1wnMSnCir45cism+m3KPYbXva2Ju7wxhXT5xx52alNYVlD3DsZgRd8cjWUF0Px4v23ZOFFNlj8B/fJCm0bRxMSKvjkyVXv9XiM9roCA/Rm5ND48EPZvh9UzXUcxKcDerebIlE2B6r2Mmn08eWNsmF+XZnmnQ04LWPq26ygmBVjRN0dm6X/Yoi2Y63V3nSTjVZID3Yb5g95ZE4+pgxV9U3/V+2DFe0wK5xMm6DqNAX5S3N6/UO7T911HMUnOir6pv5WToWo373pnuU5iIt73erNbG/vXTRhzGFb0Tf2VvAlNj2WO19N1EhNRSQ5TvX6RJp6Q6zgmiVnRN/VTtcfvn99zhDXtJJnx4YGwdwussSYec2gxFX0RGSoipSJSJiJjann8fBGZLyIhEflOjcfCkSkUv5xG0aSwFROhei+c9i3XSUwNM7wzoFFL+OQN11FMEquz6ItIEHgSGAb0BEaJSM3v9WuBG4BXanmKfaraJ3KzSzdT3IR/P8UX2prOT293HcXUUEkOdL8clr4DoUrXcUySiuVIfwBQpqqrVbUKGAuMiN5AVdeo6mLAS0BGkywqd/GNwEIKwwPxrGUwOfW6Cip3+GMiGVOLWN657YB1UcvlkXWxaiwiRSIyR0SurFc6k1xKJ9BYqikM27ALSavzBdCkLXzyuuskJknFUvRrGz5R6/EaHVU1H7gGeFxEDhqDV0RGRz4YiioqKurx1KZBLR5HuR5LsXZ1ncQcQt5dk3hpV1/2LnnXP+luTA2xFP1yoEPUcntgQ6wvoKobIj9XAzOAvrVs84yq5qtqfm5ubqxPbRrS7gpYNY2C8NmoNe0ktXfCZ9NUKmHFe66jmCQUy7t3HtBFRDqJSA4wEoipF46ItBGRRpH7xwLnAEuPNKxxqOQt0DD/CZ/tOompwzztxiZtA0usicccrM6ir6oh4BZgIrAMGKeqJSJyv4gMBxCRM0WkHPgu8LSIlER27wEUicgiYDrwkKpa0U9FS16D406jVDu6TmLq4BGgIHy2f+X03q2u45gkkxXLRqo6HhhfY929Uffn4Tf71NzvQ+D0o8xoXNv6KZTP5aHqka6TmBi9HT6H0V6h32d/wE9cxzFJxBpnTd0izQTvhAc5DmJitVRPhuN6wuJxrqOYJGNF3xyeKiwZx8ded9ZjJ9lTh0Dvq6F8LmxZ5TqMSSJW9M3hrZ8Pm1fwZvg810lMfZ3+XUDsaN98jRV9c3gLX4asxv5gXia1tGoHnc6Dxf/2v7EZgxV9czjV+/0TgT2uYBdNXacxR6L3SNj2Kayb6zqJSRJW9M2hrZjgT7jd5xrXScwRyBtTyGn/bsQebQQLX3IdxyQJK/rm0Ba+Ci1Ogk4XuE5ijtAemlAYPgs+edOGZTCAFX1zKLs+h7IpcMZICNhkKalsXPgCqNptUykawIq+OZTFY0HDfHPqSeSNKXSdxhyFIu0GbU+BBdbEY6zom9qoQvELzPW6sVpPcp3GHDWBvtfBZ7Otz76xom9qsWYWbF3Fq6Fvuk5i4uWMUSABvwuuyWhW9M3B5r8AjVsx3rO++Wmj5Ylw6sWw4GUIV7tOYxyyom++bu9W/4Rf76v9OVdN+sj/IezeBKUTXCcxDlnRN1+3aCyEq6Df9a6TmHjrMgRatoei51wnMQ5Z0TdfUfWbdtr1hxN6uU5j4i0QhP43wOrpdkI3g1nRN19Z8wFULIf8H7lOYuIsb0wheWMKOXPCiYQ0wNN/vsd1JOOIFX3zlbnPQJO20Osq10lMglTQhklePt8NzvTHVjIZJ6aiLyJDRaRURMpEZEwtj58vIvNFJCQi36nx2PUisjJys4biZLV9HSwvhH4/gOwmrtOYBHopPJi2stuf99hknDqLvogEgSeBYUBPYJSI9Kyx2VrgBuCVGvu2Be4DBgIDgPtEpM3RxzZxV/y8/9OadtLeh95prPTawcd/tSGXM1AsR/oDgDJVXa2qVcBYYET0Bqq6RlUXA16NfS8BJqvqVlXdBkwGhsYht4mn6v1Q/E/oOgzanOw6jUk44fnwUNi4CNbOcR3GNLBYin47YF3UcnlkXSyOZl/TUEregr1bbALtDPJm+Fy2azPG/+Ne11FMA8uKYRupZV2s3wlj2ldERgOjATp27BjjU5u4UIWPnoDc7tD5QgAbYC0D7KcRr4a/yejgu7B9LbS2912miOVIvxzoELXcHtgQ4/PHtK+qPqOq+aqan5trk283qNXT4fNP4OxfgNT2GW3S1b9CQ1DE77VlMkYsRX8e0EVEOolIDjASKIjx+ScCQ0SkTeQE7pDIOpMsPvx/0Pz4yCTaJpNs5Bje8wZA8Quwf4frOKaB1Fn0VTUE3IJfrJcB41S1RETuF5HhACJypoiUA98FnhaRksi+W4EH8D845gH3R9aZZLDpE1g1DQb+FLIauU5jHHg6dDlU7oSi511HMQ0kljZ9VHU8ML7Gunuj7s/Db7qpbd/nABvsIxl99ARkN7NumhlsiXb2z+XMeQoG3gTZjV1HMglmV+Rmqu1rYclr0O/70KTNl5fp20nczDNq+Tmw+3Pu+N2drqOYBmBFP1PNehwQ/wSuyWgfeT1Z6HXmp8F3wQu7jmMSzIp+JtqxHha8CH2vI+9/F9nRfcYT/hoaTl7gcxuaIQNY0c9Es/8P1INz/9t1EpMkJnn5lHrtYebDdrSf5qzoZ5pdm/wx888YaUMumC8pAR4PXQWbS+GTN13HMQlkRT/TzHrcnxnr3NtcJzFJ5j3vTDi+F8x8CMIh13FMgljRzyTb10LRs9DnWjjmFNdpTJJRAoxedzFsKfN7dpm0ZEU/k0z/X0DgwoOmRDAG8Nv2P/Hy/KP9UJXrOCYBrOhnis+XwqJXYeBoaFXrdXTGAMIjoath2xr/W6FJOzFdkWvSwLT/gUYt4NzbrIumOayZXm/o/A2Y+Uc4YxQ0ae06kokjO9LPBGtmQWkhnPNLaNrWdRqT9ASGPAD7tsMHf3IdxsSZFf10Fw7BhN9Cq44w6GbXaUyKyHt8La+FzqNy9lN+U49JG1b00938f/rj5Q95gLx7plnTjonZo6Hv4RGASXe7jmLiyIp+Otu7lW3v3sdH4Z7kvWinb0z9fE5b/l/oW7DsHVg52XUcEydW9NPZtP+hJXv4fegH1D5zpTGH9/fwZazyTmTNizdD9X7XcUwcWNFPV2vnQNGzvBC+hOVq85+aI1NNFveEfugPxjb7cddxTBxY0U9HoUoo+AW06sijoe+5TmNS3IdeLwrCg+CDx6Ci1HUcc5RiKvoiMlRESkWkTEQOupxTRBqJyL8jj38sInmR9Xkisk9EFkZuf4tvfFOrDx6DzSu4vmIUe7GZkMzRu7/6B5DTDN7+mY3Lk+LqLPoiEgSeBIYBPYFRItKzxmY3AttU9VTgz8Afox5bpap9Ireb4pTbHMrGxfDBn3gzfC4zvTNcpzFpYjOt+MWOa2B9MQ/da11/U1ksR/oDgDJVXa2qVcBYYESNbUYAL0Tuvw5cJCJ25rChVe+DN3462g8HAAAMY0lEQVQMTY/hgerrXKcxaeYdbxDjwwP476zX4YtlruOYIxRL0W8HrItaLo+sq3UbVQ0BO4BjIo91EpEFIjJTRM47yrzmcCbf64+HfuVTbKOl6zQm7Qj3VP+Q3TTxDy6q97kOZI5ALEW/tiN2jXGbjUBHVe0L3Aa8IiIHVSMRGS0iRSJSVFFREUMkc5CVk2HuMzDwZ3DqRa7TmDS1hVb8qvom/4K/iXe5jmOOQCxFvxzoELXcHthwqG1EJAtoBWxV1UpV3QKgqsXAKqBrzRdQ1WdUNV9V83Nzc+v/W2S67WvhzdFw3Gkw+Heu05g0N8Pry99Cl/ujcNqcuiknlqI/D+giIp1EJAcYCRTU2KYAuD5y/zvANFVVEcmNnAhGRDoDXYDV8YluAP8r9r+vAy8EV78I2dZbxyTeo6HvMd87FQpuhc0rXccx9VDntfmqGhKRW4CJQBB4TlVLROR+oEhVC4BngRdFpAzYiv/BAHA+cL+IhIAwcJOqbk3EL5KRVKHwdti4CEaNJe+R5cBy16lMBgiRxS+qfsHs5g/AK1fDT6ZCkzauY5kYiGrN5nm38vPztaioyHWMpJc3ppAfBwu5O/tlOP/X8M27bTA10+DW/KwtvHAF5J0L174OQRvjyRURKVbV/Lq2sytyU9RlgTncnf0yheEBdJp0hhV840TeX7fy68ofwurpMOHX/rdPk9Ss6Keizz7kseynmOt147bqn6P232gcei18IX8NXQFFz/GXe25wHcfUwapFqikvhleuplxz+UnVr6gkx3UiY/hjaCSvhr7BrVlvw0dPuY5jDsMa4FLJ+mJ48VvQ9Biu3XEbO2juOpExEcJdoRtpJXu4dOIdEAjCwJ+6DmVqYUf6qWLdPHY8czmf7WvEoI23senLC56NSQ4eAf6r+mYmhvNhwm946C4r+snIin4qWD4eXriCrdqCUVV3s9EKvklSVWRzc/WtFIQHMSZ7LE/e/X06jXnHdSwTxZp3kt28Z2H87XBiH65aPZqtNqaOSXIhsviv6pvZrU24OauAk+ULqB4M2U1cRzPYkX7yqt7vT4RSeBucOhhueNcKvkkZHgHuDN3Ig9XXcGngY/jnZbBjvetYBiv6yWnbGnh+KMz/F5z3Kzov+QF5985wncqYehL+Hr6cn1b/N7vLS9j62AB+dOcDrkNlPCv6yUQVil+Av54DW1bB1S/DRffi2X+TSWGTvXyuqHqQjXoMz+U8yvN3f48eY96wCwodsTb9ZLH1U5jwG1g5iQ/DPfl19U9Z/4IA9sYwqe9TPZFvV/2eMVmv8sOsiVwUmM8doR8Dl7mOlnHsENK16n0w4yF4ciCsmc3vq7/PtdV3sh4bYtqkl0py+H3oer5XeQ8hgryc87+8e88lnHvH866jZRQ70nclVAXzX4D3H4Xdm3gnfBYPVl9r/e9N2purPRhW9RA/yyrgp8F3uTinCCaWwdm3QovjXcdLezbKZkPbvwPmvwgf/w12rIOOg/jeysHM1R6ukxnT4I5nK7/OHsd3smZBMAf638C573enXI9jzUPW9FMfsY6yaUf6DUHVH/N+wUvsmvsSLWQfH3vdeSr0W2au6E3ts00ak/4+py23V9/EE6ER/DxYwLfm/J33czyme32gNOhP/RnMdh0zrdiRfgIc6JWw5lddYNl/oORt+PwTKjWb8d4AngsNY4l2dpzSmORzAlsYlTWNa4LTyJUd0PRY6PVt6DkCOpxl4/UfRqxH+lb042nfNlg7h3+++BznBxbTObAJgGKvC2+Fz6UgPIidNkiaMXXKJsSFgYWMCM7m4sB8Gkk127UZrXtfCp0ugE7nQ5uTXcdMKnEt+iIyFPg//OkS/6GqD9V4vBHwL6A/sAW4WlXXRB67A7gRf7rEW1V14uFeK2WKfuUu+GI5bFrkN92UF8EXSwHYq42Y4/VghncGE8Nn8jltHYc1JnU1Yx/nBZYwODifCwILyZWdAGzUtiz0TmHYJZfBCafD8adD8+NAMrO5NG5FPzKx+QrgYqAcf6L0Uaq6NGqbnwO9VfUmERkJfEtVrxaRnsCrwADgJGAK0FVVw4d6vaQp+tX7YPcXsGsj7Cj3b9vWwLZPYXMZ7Cz/ctOt2pwlXmfmet0p1q7M97pQhbVDGhN/ShdZz9mBEvoFVtJXVtIxUPHlo9u1Gav1RPr1PRNanwytO0Kr9tDiRGhxAjRqkbYfCvE8kTsAKFPV1ZEnHguMAJZGbTMC+F3k/uvAEyIikfVjVbUS+DQycfoA4KNYf5GYhUOwYy2EqyO3qq9uoUoI7ffHs6ne69+qdkPlbv+IvXKn36tm71bYu8X/WbnjoJfYqs1Zq8fzqZ7MSu8cVmo7lnons55jsZOxxjQEYaW2Z2W4PS+ELwGgNbvoHlhHD/mMU2QDnWUjGxZO5AS2EZAaB7VZjaFZLjRt60/k3qQNNGrpfxg0agE5zfxbdlN/2+wmfq+irMaRnzkQyPZPLgeyvn4T8ecRkCBIIHI/4N8Q//Ek+MCJpei3A9ZFLZcDAw+1jaqGRGQHcExk/Zwa+7Y74rSHs28b/KVv/fYJ5rAl1Iid2pRdNGW7NmcbJ7JVu1KhraigNRXamg16DBv1GHbRNCHRjTFHbjstmOP1ZA49v7Y+h2pOks2cKFs5jm0cL9s4JrSTYyt30nrrblrLelpTSgvZR3P20VQqGy70gQ8CiHwQRO63z4cfvZfQl46l6Nf20VSzTehQ28SyLyIyGhgdWdwtIqW17HcssPkwOV2xXLFLxkxgueorGXPVmmmlgyA11PPfaiLceMTfBmI6sx1L0S8HOkQttwc2HGKbchHJAloBW2PcF1V9BnjmcCFEpCiW9qqGZrlil4yZwHLVVzLmSsZMkJy5Yhl7Zx7QRUQ6iUgOMBIoqLFNAXB95P53gGnqnyEuAEaKSCMR6QR0AebGJ7oxxpj6qvNIP9JGfwswEb/L5nOqWiIi9wNFqloAPAu8GDlRuxX/g4HIduPwT/qGgJsP13PHGGNMYsV0eZuqjgfG11h3b9T9/cB3D7Hvg8CDR5HxgMM2/zhkuWKXjJnActVXMuZKxkyQhLmS7opcY4wxiWPj6RtjTAZJyaIvIreLiIrIsa6zAIjIAyKyWEQWisgkETkpCTI9IiLLI7neEpHWrjMBiMh3RaRERDwRcd6rQUSGikipiJSJyBjXeQBE5DkR+UJEPnGd5QAR6SAi00VkWeT/75euMwGISGMRmSsiiyK5fu86UzQRCYrIAhF513WWA1Ku6ItIB/whIda6zhLlEVXtrap9gHeBe+vaoQFMBnqpam/8YTTucJzngE+AbwPvuw4SGWLkSWAY0BMYFRk6xLV/AkNdh6ghBPxKVXsAZwE3J8m/VSXwTVU9A+gDDBWRsxxnivZLYJnrENFSrugDfwZ+Qy0XebmiqjujFpuRBNlUdZKqhiKLc/CvkXBOVZepam0X37nw5RAjqloFHBhixClVfR+/F1zSUNWNqjo/cn8XfiFLzNX19aC+3ZHF7MjN+fsPQETa408C/A/XWaKlVNEXkeHAelVd5DpLTSLyoIisA64lOY70o/0ImOA6RBKqbYgR54Us2YlIHtAX+NhtEl+kCWUh8AUwWVWTIhfwOP4Bquc6SLSkm5FARKYAJ9Ty0F3AncCQhk3kO1wuVf2Pqt4F3BUZSvoW4D7XmSLb3IX/1fzlROepT64kEdMwIeYrItIceAP4rxrfcJ2JXPvTJ3Le6i0R6aWqTs+HiMjlwBeqWiwiF7rMUlPSFX1VHVzbehE5HegELPIH8KQ9MF9EBqjqJle5avEKUEgDFP26MonI9cDlwEXagH1z6/Fv5VpMw4QYn4hk4xf8l1X1Tdd5alLV7SIyA/98iOuT4OcAw0XkUqAx0FJEXlLV6xznSp3mHVVdoqrHqWqequbhv2H7NUTBr4uIdIlaHA4sd5XlgMjEN78FhqvqXtd5klQsQ4wYIDJU+rPAMlV9zHWeA0Qk90DPNBFpAgwmCd5/qnqHqraP1KqR+EPTOC/4kEJFP8k9JCKfiMhi/OanZOjO9gTQApgc6Ur6N9eBAETkWyJSDgwCCkXksDOpJVLkRPeBIUaWAeNUtcRVngNE5FX8OSe6iUi5iNzoOhP+kev3gW9G/p4WRo5iXTsRmB55783Db9NPmu6RyciuyDXGmAxiR/rGGJNBrOgbY0wGsaJvjDEZxIq+McZkECv6xhiTQazoG2NMBrGib4wxGcSKvjHGZJD/D9baVlcx36lsAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "mala_dens, bins, _ = plt.hist(mala_samples, bins=120, density=True)\n",
    "dens = np.zeros((bins.size, 1))\n",
    "for i in range(dens.size):\n",
    "    dens[i] = norm(bins[i])\n",
    "plt.plot(bins, dens)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Hamiltonian Monte Carlo (HMC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Startujemy z rozkładu $$ \\pi(q) = e^{-V(q)}. $$\n",
    "Dodajemy sztuczną zmienną (pęd) $p$:\n",
    "$$ \\pi(q,p) = \\pi(q)\\pi(p|q) = e^{-V(q)-K(p,q)}.$$\n",
    "Zazwyczaj $K(p,q) = \\frac{p^Tp}{2}.$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ H(q,p) = V(q) + K(p,q) $$\n",
    "Chcemy\n",
    "$$ \\frac{\\partial}{\\partial t} H = 0 $$\n",
    "czyli\n",
    "$$ \\frac{\\partial q}{\\partial t} = \\frac{\\partial H}{\\partial p} $$\n",
    "$$ \\frac{\\partial p}{\\partial t} = - \\frac{\\partial H}{\\partial q} $$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "bo wtedy zachowamy łączny rozkład, więc w szczególności zachowamy brzegowy."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### LFI:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$p^0, q^0 \\rightarrow p^1, q^1$\n",
    "\n",
    "$\\varepsilon$ - krok dyskretyzacji"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$p^{\\frac 1 2} = p^0 - \\frac \\varepsilon 2 \\frac{\\partial V}{\\partial q}(q^0)$$\n",
    "$$ q^1 = q^0 + \\varepsilon p^{\\frac 1 2}$$\n",
    "$$ p^1 = p^{\\frac 1 2} - \\frac \\varepsilon 2 \\frac{\\partial V}{\\partial q}(q^1)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "def hamiltonian_lfi(p0, q0, grad_v, epsilon, l):\n",
    "    p = p0\n",
    "    q = q0\n",
    "    for i in range(l):\n",
    "        p_half = p - (epsilon/2)*grad_v(q)\n",
    "        q += epsilon * p_half\n",
    "        p = p_half - (epsilon/2) * grad_v(q)\n",
    "    return np.stack((q,p))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Algorytm:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$ q_n \\rightarrow q_{n+1} $"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. $q^0 = q_n, p_0 \\sim \\mathcal{N}(0, \\text{Id})$\n",
    "2. $q^0, p^0 \\rightarrow q^l, p^l$ - Leap Frog Integrator\n",
    "3. $\\alpha = \\exp\\{-V(q^l) - K(-p^l) + V(q^0) + K(p^0)\\} \\wedge 1$ (dajemy $-p^l$, ale to jest symetryczne)\n",
    "4. $q_{n+1} = q^l$ z prawdopodobieństwem $\\alpha$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "def hmc(q_0, epsilon, l, v, grad_v, n, k=lambda p: np.sum(p**2)):\n",
    "    q_vec = np.zeros((n,q_0.size))\n",
    "    q_vec[0] = q_0\n",
    "    for i in range(1, n):\n",
    "        q = q_vec[i-1]\n",
    "        p = np.random.randn(q.size)\n",
    "        alpha_2nd = v(q) + k(p)\n",
    "        q, p = hamiltonian_lfi(p, q, grad_v, epsilon, l)\n",
    "        alpha_1st = -v(q) - k(-p)\n",
    "        alpha_cand = np.exp(alpha_1st + alpha_2nd)\n",
    "        alpha = min(1, alpha_cand)\n",
    "        if np.random.rand() < alpha:\n",
    "            q_vec[i] = q\n",
    "        else:\n",
    "            q_vec[i] = q_vec[i-1]\n",
    "    return q_vec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "def v_for_normal(q):\n",
    "    n = q.size\n",
    "    return np.sum(q**2)/2 + (n/2)*np.log(2*np.pi)\n",
    "\n",
    "def grad_v(q):\n",
    "    return q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "qs = hmc(np.zeros((1,1)), 0.001, 100, v_for_normal, grad_v, 10**5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x168e6d2f908>]"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl8VOW9x/HPbyYJq2yC4hVsQHFBNjGAoOJSURQFrRvaWrfqtW61rdeL1Wpr21trN5eqFde6IgpoBBQXXAFZZV8sIAKigApElpDMOc/94wQbMJAhzOSZmXzfr1deZiaHma+SfH3ynOc8x5xziIhIbon5DiAiIqmnchcRyUEqdxGRHKRyFxHJQSp3EZEcpHIXEclBKncRkRykchcRyUEqdxGRHJTn641btmzpCgsLfb29iEhWmj59+pfOuVbVHeet3AsLC5k2bZqvtxcRyUpm9mkyx2laRkQkB6ncRURykMpdRCQHqdxFRHKQyl1EJAep3EVEcpDKXUQkB6ncRURykMpdRCQHebtCVSRTFQ4Z8+3ny+4c4DGJSM1p5C4ikoNU7iIiOUjlLiKSg1TuIiI5SOUuIpKDVO4iIjlISyFFUkDLJyXTaOQuIpKDNHKXOk0jbslVSZW7mfUH7gHiwCPOuTt3ctw5wAtAD+ecbpAqWaVy0Vf1nMpfskm10zJmFgfuB04FOgIXmFnHKo7bC7gemJzqkCIisnuSGbn3BBY755YCmNkwYBAwf4fjfgfcBdyY0oQiaVDVKL0mf0ajeclUyZT7/sCKSo9XAr0qH2BmRwBtnXOjzUzlLlmlHmUcGfuYvSmhjDw20YDpYQe2UN93NJEaS6bcrYrn3LdfNIsBfwcuqfaFzK4ErgQ44IADkksokiLbj7wdx8dmcVn8VXrFFlLPyrc7ttTl817YheeD43kr7E7VPwYimSuZcl8JtK30uA2wqtLjvYBOwDtmBtAaKDazgTueVHXODQWGAhQVFTlEPDjYVnBr3tP0jc9hRdiKp4KT+CDszArXinwCWtoGvh+bwSnxqZwcn874oBu3Jy5mhdvXd3SRpCVT7lOBDmbWDvgMGAxcuO2LzrkNQMttj83sHeBGrZaRTHRGbCJ/yX+ILRRwR/lFPBX0o3zHHwMHH4Sd+X3iR1wSH8fP817kjYKbuKn8SorDo/0EF9lN1Za7cy5hZtcC44iWQj7mnJtnZncA05xzxekOKbLHnONn8RH8PH8Ek8ND+WnZDXxNk13+kYA4jwanMTo4insL/sG9Bfezb/k6Hg4GsG2apiYnZkVqgznnZ3akqKjITZumwb3UAufgtZth8oO8GPTlV+WXU0b+br1EPcr4a/6DnB6fzD8Tp3Nn4gKqm4fXShpJBzOb7pwrqu44XaEquW/ivTD5QR5L9OeOxEXU5OToVgq4rvw6vnZNuCpvNGtccx4LTt3ln9EFUOKTyl1y2+zh8MZtcPhZ/G762ezJqhdHjNsTF9PK1nNr3tOscK14I6x2ACXihTYOk9z1+Wx4+RooPBbOegiXgm93R4yfl1/NbNeOe/Lvp6Mt2/OcImmgcpecdNiQESx58DxWJxrCuf+CvHope+1S6vGTsv+hhIbck38/9ShL2WuLpIrKXXLS7XlP0s6+4Ibya6DR3il//S9pyo3lV9Eh9hlD8p5L+euL7CmVu+Se+S8zOO8dHggGMik8PG1v80HYmccS/bk0bxx9Y7PS9j4iNaFyl9xSugHG3sScsJC7E2en/e3+lBjMorANd+UPpTGb0/5+IslSuUtuGf8H2LiaX5X/hETFYrDCIWPSdrHRVgq4qfxK9mE91+eNSst7iNSEyl1yx6qZMPVh6PET5rj2tfa2s9xBPB8cz6Xx1+hgK2vtfUV2ReUuuSEMYfTPoWFLOPHWWn/7uxLns4n63JH3BJU2Tf2Obb9FaNsCSTeVu+SGuSNg1Qw4+XfQoFmtv/06mvDnxPn0js/n9NiHtf7+IjvSFaqS/RJbYfwd0LozdD7PW4znghP5UfxNbswbzmtlPb6d8wdtMCa1TyN3yXq/uf1GWL+ci5YPgJi/b+mQGH9OnEdhbDXnxd/1lkMEVO6S7UpLuC5vFB8Eh/N+2Nl3GsaHRzAtPJjr80bqylXxSuUu2W3SP9jbvuFPSWzBWzuMP5UPprWt4+L4ON9hpA5TuUv2Kt0AH/6T14Ietbr0sTpT3aG8HXTl6rxiGrHFdxypo1Tukr2mPAxbN3Bf4kzfSb7j74lzaGabuDD+1k6P0ZJISSeVu2Snsk0w6X7ocDLzXDvfab5jtjuQ94NOXJE3VnPv4oXKXbLTtMdhy9fQ93+2ezqTRsMPBIPYx9ZzrlbOiAcqd8k+5aXRrfPa9YW2PX2n2alJYUdmhAfx3/HR5JHwHUfqGJW7ZJ85L8DG1XDML3wnqYZxf2IQbWNrGRib6DuM1DG6QlWyi3PRXPu+nSl8eBOQGVMwO/NW2J1FYRt+kvcqI8uOJTOWa0pdoJG7ZJcl42HtAuh9NdlRlMZjwal0jH3KUbEFvsNIHaJyl+wy6X5ovC90Sv+NOFLlpeBovnJ7cVn8Vd9RpA5RuUv2WLMAlrwFPa9I6Q2v020rBTwTfJ+TYjM4wFb7jiN1hMpdsseHD0JeAzjyMt9JdtvTiX4kiHFJFVsSaI93SQeVu2SHLeth9nDofA402tt3mt22huaMDntzXvwd3WtVaoXKXbLDrOcgsSWakslS/0qcTGMr5cz4BN9RpA5QuUvmcw6mPgJtesJ+XX2nqbFZ7kDmhIX8KP4mu7oVn0gqqNwl8y19B75aDD1+4jvJHjKeDvpxaGwFRbbIdxjJcbqISTLf1Eeg4d7QcVDWn3QsDnpzS94z/CjvTaaVH+o7juQwjdwls234DBaNhSMugvz6vtPssS3UZ0RwLKfGptCCEt9xJIep3CWzzXwGXAhHXuI7Sco8E3yfepbgvPg7vqNIDlO5S+YKQ5jxFLQ/Hlpk3p7tNbXYtWFyeCjnx99GJ1YlXVTukrmWvg0blkP3H/tOknLPJ46nXWw1vWyh7yiSo1TukrlmPAkNWsChp/tOknJjw16UuAacl/e27yiSo1Tukpk2fQkLx0DXCyCvXs5dnl9KPV4OjmZAbDJN2OQ7juQglbtkplnDICyH7hf5TpI2w4ITqG/lDIz/50YeufY/MfFH5S6Zxzn46Clo0wP2Ocx3mrSZ59oxL/xexYlVkdRKqtzNrL+ZLTKzxWY2pIqvX2Vmc8xsppl9YGYdUx9V6oxVM2DtQoZ80jXnR7LDghPoHFtGR1vmO4rkmGrL3cziwP3AqUBH4IIqyvtZ51xn51w34C7gbylPKnXHzGchrz5jgqN8J0m74qAPW10eZ8ff9x1FckwyI/eewGLn3FLnXBkwDBhU+QDnXOVL7RqhxbtSU+Wl0Q2wDzuDb2joO03abaAxb4bdGRSfQB4J33EkhyRT7vsDKyo9Xlnx3HbM7BozW0I0cr8+NfGkzlk0Fko3QLcLfSepNSOCvrS0Eo6PzfIdRXJIMuVe1V2IvzMyd87d75w7EPhf4NYqX8jsSjObZmbT1q5du3tJpW6Y+Sw02R/aHec7Sa15L+zCWteUc+Lv+Y4iOSSZcl8JtK30uA2wahfHDwPOrOoLzrmhzrki51xRq1atkk8pdUPJ59E9UrsOhljcd5pakyCPl4KjOTE2g+baTExSJJlynwp0MLN2ZlYADAaKKx9gZh0qPRwA/Dt1EaXOmDM82iSsa92ZktlmRNCXAgsYGJ/kO4rkiGrL3TmXAK4FxgELgOHOuXlmdoeZDaw47Fozm2dmM4FfABenLbHkrlnPw/5F0PIg30lq3UJ3AHPDQs7W1IykSFI363DOjQXG7vDcbZU+/1mKc0ld88UcWDMPTvuL7yTejAqO4df5T3OgfeY7iuQAXaEqGWHoP/5IuYtzxMjGOX/h0s4UB70JnOkG2pISKnfxLwwYFJ/AO2E31tHEdxpv1tKcD8LOnBX/INrLXmQPqNzFv0/eZV9bz8jgGN9JvBsZHEMb+xJWfOg7imQ53SBb/Jv1PCWuIePDI3wn8e71sIhNrh4vP3wXv0pcAcCyOwd4TiXZSCN38atsEyx4hTFBL7ZS4DuNd1uoz2thD06PT6YeZb7jSBZTuYtfi16F8k28pCmZb40KjqWJbeb42EzfUSSLqdzFr9nDoUkbprhDfCfJGJPCjqx1TbVqRvaIyl382fRVtN1A57Nx+lb8VkCcV4LenBj7SLfgkxrTT5T4M28khAnofJ7vJBnnpeBo6lmCU+JTfUeRLKVyF3/mvAD7dITWnXwnyTizXXuWhq05M6apGakZLYUUP9YtgxWT4fu3+06SoYzisA/Xx0fRa8hTrKYFoGWRkjyN3MWPOS9G/+x8jt8cGeyl4Ghi5jhDO0VKDajcpfY5F03JHNAbmh3gO03GWub2Y2bYXqtmpEZU7lL7Vs+DtQs1ak9CcXA0nWLLaG+7uj+OyHep3KX2zXkBYnnQ8SzfSTLe6OAoQmcMjE/0HUWyjMpdalcYwtwR0P4EaLS37zQZbw3NmRR2ZGBsIlXculhkp1TuUrtWToENK6Dzub6TZI2Xwz60j31BZ/vEdxTJIloKKbVrzgtscQUUPRdn03N174YcNfFa0IPf5z1WMTVzne84kiU0cpfaE5TDvFG8FXZnEw18p8kaJTTmnbAbZ8Qn0X7IK3XyLlWy+zRyl9qz9F3Y/BUvB7p/+u4qDvpwcnw6vWILmBQevl3B68ImqYpG7lJ75rwA9ZvybtjVd5Ks82bYnU2uHmfEtGpGkqNyl9pRvgUWjobDzqCMfN9psk4p9Xg9LOK0+BTySfiOI1lA5S614+NxULZRq2T2QHHQh2a2ib6xWb6jSBZQuUvtmPMCNN4XCo/1nSRrvR925mvXmIHaa0aSoHKX9CvdAP9+Aw4/C2Jx32myVoI8Xg160S82nQaU+o4jGU7lLum3YDQEW6GT9pLZU8VBHxraVvrFZviOIhlO5S7pN+cFaPY9aFPkO0nWm+IO4XPXgjO014xUQ+Uu6bVxDXzybnQi1cx3mqzniPFK0JvjYrNoykbfcSSDqdwlveaNAhdqe98UejnoQ4EFnBqf4juKZDCVu6TXnBdh306wz2G+k+SMea6QJeF+DNIFTbILKndJn3XLol0gO53tO0mOMYqDPvSKLWBfvvYdRjKUyl3SZ+6I6J8q95QrDvsQM8fp8UkUDhmjzcTkO1Tukj5zXoS2R0Hz7/lOknM+cfsxK2zPIK2akZ1QuUt6rJ4Ha+brRGoaFQd96BL7hHb2ue8okoG05a+kx+zhYPHoqlTQtEEavBL05pa8ZxgUn8DdCf1PVLankbuk3rb7pB54IjRq6TtNztp2f9VBsQno/qqyI5W7pN6KD6P7pHY5z3eSnPdy2Id2sdV0saW+o0iGUblL6s15AfIbwiGnaSVHmr0W9GSry+PM+ATfUSTDaM5dUitRFl2VeshpUK+x7zQ5r4RGjA+P4Iz4JA4cUkxAtOumbr0nGrlLai0ZD1vW6aYctejl4Gha2QZ6x+b7jiIZJKlyN7P+ZrbIzBab2ZAqvv4LM5tvZrPN7C0z08LmumrOcGjQPDqZKrXi7bAbJa6BpmZkO9WWu5nFgfuBU4GOwAVm1nGHwz4CipxzXYAXgbtSHVSyQGkJLBwTXZGaV+A7TZ2xlQJeDXrRPzaF+mz1HUcyRDIj957AYufcUudcGTAMGFT5AOfc2865zRUPPwTapDamZIWFoyFRCl3O952kznkpPJrGVkq/2HTfUSRDJFPu+wMrKj1eWfHczlwOvFrVF8zsSjObZmbT1q5dm3xKyQ6zn4fmhdCmh+8kdc6H4WGsci00NSPfSqbcq7rDQpVXTJjZj4Ai4M9Vfd05N9Q5V+ScK2rVqlXyKSXzlXwOS9+NRu26KUetc8QoDo7muNgsWlDiO45kgGTKfSXQttLjNsCqHQ8ys5OAW4CBzjlN/NU1c18EHHTWhUu+jAyOIc9CzohP8h1FMkAy69ynAh3MrB3wGTAYuLDyAWZ2BPAQ0N85tyblKSXzzX4e9j8SWh6ki5Y8+di1ZX74Pc6Kf+A7imSAakfuzrkEcC0wDlgADHfOzTOzO8xsYMVhfwYaAy+Y2UwzK05bYsk8X8yFL+boRGoGGBkcQ7fYEvjy376jiGdJXaHqnBsLjN3hudsqfX5SinNJNpk9DGJ50Ek7E/pWHPTh5rxnic9+Hk681Xcc8UhXqMqeCRLR9r4dTqHwdx9qSsazNTRnQtgJZg2LdueUOkvlLntm6TuwcTV0Hew7iVR4Megb7cr5qebe6zKVu+yZWc9B/WZw8Cm+k0iF18MivnENePGxKlckSx2hcpeaKy2JrkrtfA7k1fOdRiqUUo+xQS9OjU2Gsk2+44gnKnepufkvRdsNdL3AdxLZwYjgWBrZVljwiu8o4onKXWpu5rOwd4dofbtklKnuEJaHraK/I6mTVO5SM18tgeWT4IgfaruBDOSIMTI8Fj55D9Yv9x1HPFC5S83MfAYsBl20SiZTvRj0BRzMfM53FPFA5S67LwyiwjjoJGiyn+80shMr3T5MCA5n+fiHaTdEc+91je6hKrtv6dvwzSro/0ddtJThhgfHcU/BAxwVWwCc4TuO1CKN3GX3ffR0dCu9Q071nUSq8VrYkxLXkHPi7/qOIrVM5S67Z/PX0a30Op+nte1ZYCsFvBL05rTYFDoPeUG/adUhKnfZPbOHQ1AG3S/ynUSSNDw4jgZWxuna571OUblL8pyDGf+C/+oOrTv7TiNJmuUOZGHYlsHxt31HkVqkcpfkrZwGa+bDkRf7TiK7xRgWnEDX2FI62jLfYaSWqNwleTOegPxG0Ols30lkN40KjqHU5Wv0Xoeo3CU5pSUwdyR0+gHU28t3GtlNG2jM2LAXZ8Y/0GZidYTKXZIz90Uo3wxHXuI7idTQsMQJNLEtMO8l31GkFqjcpXrOwbTHYZ/DtUlYFpviDmVJuB9Mf8J3FKkFKnep3mfT4YvZ0OOybzcJKxwyRmums47xbHAirJwS3dRccprKXao39VEoaAxdzvedRPbQiKAv5NWHaY/6jiJppnKXXdv8NcwbCV3O04nUHLCevaLVTrOej06SS85SucuuzXw2uttS0eW+k0iKDJx8KJRv4te/+7XvKJJGKnfZuTCMfn1v2wtad/KdRlJktjuQWWF7Loq/EZ0sl5ykcpedWzoevl6qUXsOejo4iYNjn8GnE3xHkTRRucvOTX4IGu0Dh5/pO4mk2CtBb9a7RjBlqO8okiYqd6naV0vg369Dj8u1tW8OKqUew4ITYcFoWL/CdxxJA5W7VG3yQxDLhyMv9Z1E0uTJRD8SoeOff7nZdxRJA91mT76rtCS6AXanH8Be+377tC5ayi2raMlrYQ8uiL8V7TdT0Mh3JEkhjdzlu2Y+C2Ubodd/+04iafZ4oj9NbTPMGuY7iqSYyl22FwYw+UFo01P7yNQB093BzArbw+R/RktfJWeo3GV7C16Bdcugz3W+k0itMB5L9IcvP4bFb/oOIymkcpf/cA4m3gst2sOhA3ynkVoyJjwKmrSBCff4jiIppHKX/1j+YbQDZO9rIBb3nUZqSYI86H01fPoBrJzuO46kiMpd/mPivdCgBXS90HcSqW3dfwz1msJEjd5zhcpdIms/hkVjoecVUNDQdxqpZYW3v8f9m44jmFccXcAmWU/lLpEP/g55DaDnlb6TiCdPJE4hQRwm3uc7iqSALmKSaHXM7Oeh11XQqCWgC5bqorU058XgOH448xk47iZo8l++I8ke0Mhd4IO7oxOoWv5Y5z0YnBFd6zDhXt9RZA8lVe5m1t/MFpnZYjMbUsXX+5rZDDNLmNk5qY8paVOyKtpqoNsPocl+vtOIZyvdPtB1MEx/HDau8R1H9kC15W5mceB+4FSgI3CBmXXc4bDlwCXAs6kOKGk28b5opHbMDb6TSKY45hcQlGnuPcslM+feE1jsnFsKYGbDgEHA/G0HOOeWVXxN1y9nk5LPYdpj0Y2vmxcCmmsXoOVBcPgPohujH30DNNrbdyKpgWSmZfYHKm/4vLLiOcl27/8VwkR08kykQuGQMZw0vSeUb4YJd/uOIzWUTLlbFc/V6MaLZnalmU0zs2lr166tyUtIqqxfDtOfgCMughbtfKeRDLPYtYl+o5syNPoNT7JOMuW+Emhb6XEbYFVN3sw5N9Q5V+ScK2rVqlVNXkJS5b0/gxn0vdF3EslUxw+JfrN7/y++k0gNJFPuU4EOZtbOzAqAwUBxemNJWn21BD56Booug6ZtfKeRDFV413yeLjuO8imPR9dCSFapttydcwngWmAcsAAY7pybZ2Z3mNlAADPrYWYrgXOBh8xsXjpDyx5667eQVz9aFSGyC/clziIgBm//0XcU2U1JXaHqnBsLjN3hudsqfT6VaLpGMt2KKTD/ZTj+5u1uoSdSldW04IngFK6a/Xy0c+R+XX1HkiTpCtW6xDl4/VZovC/0vtZ3GskSDyQGQYPmMO6W6HtIsoLKvS5ZUAwrJsMJv4J6jX2nkSxRQqPoN71l78PHr/mOI0lSudcV5aXwxu3Q6lDo9iPfaSTbFF0Ke3eA138NQbnvNJIE7QpZV0y6D9Z9AheNgvj2f+26KlWqFc+Hk38Hzw2O1r73vsZ3IqmGRu51wfoV8N5f4bCBcOCJvtNItjq4PxzUL1o5880XvtNINVTudcHrt0T/POUPfnNIdjPjuHmnsnVrKS/ddanvNFINlXuuWzI+Wvp47C+h2QG+00iW+9S15p/BGZwZnwifvO87juyCyj2XlW2CV34GLQ7UjTgkZR5IDGR52ArG/CI6US8ZSSdUc9nb/xdtEHbJWMivv92XdBJVdkfl75etFHBr4jKe/PJP8N5d8P3bdvEnxReN3HPVZ9Phwwei/WMKj/adRnLMe2HX6O5dH9wNn8/yHUeqoHLPRYmt8PK10Lg1nPRb32kkV53yh+iG6i9fo7XvGUjlnoveugPWzIcz7oH6Tbb7UuGQMZqSkdRo0BwG/A2+mAPv3uU7jexA5Z5rlr4Lk/4BRZfDwSf7TiO57rDToeuF0Z7vn07ynUYq0QnVXLJlHbz00+gy8ZN//+3TGqlLWp12FyyfCCOvhJ9+APWb+k4kqNxzh3PRPPvG1XD5G1DQ0HciyXGVBw3LrnkEHjsFRv8czn40usuXeKVpmVwx8T5YOBr63QH7d/edRuqatj3ghJth7giY+ojvNIJG7rlh2QR48zfR3jFHXf3t05qOkdpSOGQMxqE8kn8Efcf8L+eP2sjIP/7cd6w6TSP3bLfhM3jxUmheCIPu16/D4o0jxs/Lf8oqtzcPFNwDG9f4jlSnmfN0Z5WioiI3bdo0L++dM7ZuhMf7w9fLOPmbW/nYtfWdSITD7FNGFtzOQncAg8tuZdGdZ/mOlFPMbLpzrqi64zRyz1ZhwBu/H0jw+Vwu2fhTFbtkjAXue9xQfg1dbQl/y38AwtB3pDpJ5Z6NnINxv6JffDp3JH7MO2E334lEtjMu7MH/JS5kQHwKvPUb33HqJJV7Nnr3Lpj8Tx5NnMq/Al2oJJnpkeA0nkqcBBPuiT6kVqncs83kh+Cd/4NuP+T3iR8COoEqmcq4PXEJxUFveOM2br3lBt+B6hSVezaZ8jC8ehPjgiIO/LA/Tn99kuFCYvyi/Ke8EXTn9/mPc9MtN2qJbi1RO2SLCffC2BvhkAFcV34dAXHfiUSSkiCPa8uv592gC3flP8yl8Vd9R6oTVO6ZzjkY/wd449dw+A/gvH9RRr7vVCK7ZSsFXFH+S14LenB7/lPReSNPy7DrCpV7JkuUwairorvdHHERnP0IxFXskp3KyOea8usZERwDb/8Biq+LvsclLVTumWrTl/D0D2D2MDjxVhh4H8Q0FSPZLSDOjeVXcW/iTPjoKSb+9jjY/LXvWDlJe8tkohVT4YWLo4L/wcPQ5TydhJKc4Yjxt8R5LAtbc2f+wzD0eDj3CW14l2IauWeSMIRJD1D2yCksX1/GgM23U/hsIxW75KSRYV/OL7sNwiDaLnjyUM3Dp5DKPVOsXw5PDoRxN/Nu2JXTy/7APFfoO5VIWn3kOtBtza8ZX9YRXv2faCpyw0rfsXKCyt23IBFdmPRAH1j1EQy8jyvKf0kJjX0nE6kV69mLy8tv5NbyS9m8eAIlfyuCaY9HI3qpMZW7T8snw8MnwKs3QZsi+OkE6P5jdNWp1DWOGE8H/Til7E7mhu1g9A3w8InR+SepEZW7D2sWwnMXwmMnRydNz30CLhoV7ckuUoetcPtyYfkt0a36Nq6GR0+C4T+GtYt8R8s6Wi1Tm1bNhAl3w/yXoaBxtMTxqKuhoBGgOyeJRAw6nwMHnxJdmf3hA7DgFUYlejM0cTqv/vHq6l9CdLOOtAvKYdFYmPoofPIu1GsCPS6HPtdDwxYqdJFdWHbngOi32wl3s3HCwzS2Ut4LOvNU0I/x4REsuXPgt8du+1laducAX3FrRbI361C5p4Nz8MUcmPMCzB4OG7+Apm2jUi+6DOo3/fZQlbtIcpqwkR/Gx3Nx3jha2zpWu2aMCPpSHPRhoWtLVeeqcrHoky13TcukShjAymnRKH3RWPjyY4jlwUH94MhLoEO/b68wVaGL7L4SGvNgMJChwQBOiM3kgvh4royP5uq8Yv4d7s/r4ZG8FXRnpjuIUKcTVe41FiRg9VxYMSWabln2PpRuiAq98BjodRUcfhY0bOE7qUhOCYjzZngkb4ZH0oISTotPZkBsMv8dH801ecWsc435MDyMSWFH+Lwt7NMR4nWv6pKaljGz/sA9QBx4xDl35w5frwc8CRwJfAWc75xbtqvXzKppmdKSaCS+dmE03fL5LPh8NpRvir7e9AA48Hhofzwc+H1o0GyX838auYukXhM2cnxsNsfE5tAnPo829iUAm1w9GhUWQevO0Uerw6BlB6jfxHPimknZnLuZxYGPgX7ASmAqcIFzbn6lY64GujjnrjKzwcBZzrnzd/W6GVPuQQK2fA3ffBEtvSpZFV0ht2ElrPsEvl4aPb9NfqPoG2S/rtC2Z/TRtC2YqbRFMoajja2luy3miNi/6RJbymG2nIa29T+H7LUfNG8XLUFu1hb1R0jwAAAF8ElEQVSa7B997LUvNG4d/dadgZv1pbLcewO/cc6dUvH4ZgDn3B8rHTOu4phJZpYHfAG0crt48RqX++avYdNaCMqilShBecXnWyGxFcq3QKI0+mf5FijbBGUbYes3FR8lfPTxMpqwiQMblcKW9cAOMS0W/eW2aFfx0R5aHQotD6H9XxZoPk8kC8UI+Z6tpoOt5CBbRfvY57S1NRxga9iXdcRs+x4InbGeRrRo2RrqN4MGzaLFEPX2ipYyFzSGgoaQ3xDyG0QfefUhrx7E60G8APIKIJYfbdUdy/vPR4Nm3y6B3l2pPKG6P7Ci0uOVQK+dHeOcS5jZBmBv4Mvk4u6GGU/Cm7fv1h/Z7OqxifqUuIZspAHfuEaspCUTv2nMVzThK9eEta4Za1wzVrvmrKY5iS15sHbHV1qErvsSyU4hMT5x+/GJ24/XASrtbpBHgn1ZR2v7mla2gX1sHXvbNzTnG1qs+YYmbKKpLWUvNrOXbaExW2hge7AX/YC/Qo+f7Om/0i4lU+5VXQu/44g8mWMwsyuBKysebjSzbZedtSQd/yNILWXcc5meD5QxFTI9H1SRcUltvvtvrwCuqO6onf13/F4yb5FMua8E2lZ63AZYtZNjVlZMyzQFvrMDv3NuKDB0x+fNbFoyv2b4pIx7LtPzgTKmQqbng7qRMZk5hqlABzNrZ2YFwGCgeIdjioGLKz4/Bxi/q/l2ERFJr2pH7hVz6NcC44iWQj7mnJtnZncA05xzxcCjwFNmtphoxD44naFFRGTXklrZ75wbC4zd4bnbKn1eCpy7Bzm+M1WTgZRxz2V6PlDGVMj0fFAHMnrbW0ZERNJH6/pERHJQxpW7md1oZs7MWvrOsiMz+52ZzTazmWb2upn9l+9MlZnZn81sYUXGUWbWzHemHZnZuWY2z8xCM8uY1Qpm1t/MFpnZYjMb4jtPVczsMTNbY2ZzfWepipm1NbO3zWxBxd/xz3xn2pGZ1TezKWY2qyLjb31nqoqZxc3sIzMbXdPXyKhyN7O2RNscLPedZSf+7Jzr4pzrBowGbqvuD9SyN4BOzrkuRFtG3Ow5T1XmAj8A3vMdZJuKLTbuB04FOgIXmFlHv6mq9ATQ33eIXUgAv3TOHQYcBVyTgf8dtwInOue6At2A/mZ2lOdMVfkZsGBPXiCjyh34O3ATVVwAlQmccyWVHjYiw3I65153ziUqHn5IdE1CRnHOLXDOZdo903oCi51zS51zZcAwYJDnTN/hnHuPKq4fyRTOuc+dczMqPv+GqJz295tqey6yseJhfsVHRv0cm1kbYADwyJ68TsaUu5kNBD5zzs3ynWVXzOwPZrYC+CGZN3Kv7DLgVd8hskRVW2xkVCllGzMrBI4AJvtN8l0VUx4zgTXAG865TMt4N9EgN9yTF6nVTY7N7E2gdRVfugX4FXBybeapyq4yOudeds7dAtxSsYHatcDubXST5nwVx9xC9CvyM7WZbZtkMmaYpLbPkOSYWWNgBHDDDr/tZgTnXAB0qzgnNcrMOjnnMuI8hpmdDqxxzk03s+P35LVqtdydcydV9byZdQbaAbPMDKLphBlm1tM590UtRtxpxio8C4yhlsu9unxmdjFwOvB9X1cJ78Z/w0yRzBYbkgQzyycq9meccyN959kV59x6M3uH6DxGRpQ7cDQw0MxOA+oDTczsaefcj3b3hTJiWsY5N8c5t49zrtA5V0j0w9a9tou9OmbWodLDgcBCX1mqUnFTlf8FBjrnNvvOk0WS2WJDqmHRyOxRYIFz7m++81TFzFptW0VmZg2Ak8ign2Pn3M3OuTYVPTiYaCuX3S52yJByzyJ3mtlcM5tNNIWUaUu9/gHsBbxRsVzzn74D7cjMzjKzlUBvYEzFvQC8qjgJvW2LjQXAcOfcPL+pvsvMngMmAYeY2Uozu9x3ph0cDVwEnFjx/TezYgSaSfYD3q74GZ5KNOde4+WGmUxXqIqI5CCN3EVEcpDKXUQkB6ncRURykMpdRCQHqdxFRHKQyl1EJAep3EVEcpDKXUQkB/0/WoQPegKIO6gAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n, bins, _ = plt.hist(qs[:,0], bins=120, density=True)\n",
    "dens = np.zeros((bins.size, 1))\n",
    "for i in range(dens.size):\n",
    "    dens[i] = np.exp(-v_for_normal(bins[i]))\n",
    "plt.plot(bins, dens)"
   ]
  },
  {
   "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}