Sum of non-identical uniform distributions

sum_uniformly_distributed_variables.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The distribution of the sum of N uniformly distributed variables came up in a recent RdTools PR.  The approach taken there was to just sample each distribution and sum the samples, but I wondered how it could be done analytically.  This notebook is inspired by [this StackExchange post](https://math.stackexchange.com/a/2966816/197200) and extends the derivation from just summing two uniform distributions on $[0,1]$ to two distributions with arbitrary bounds.\n",
    "\n",
    "Given two independent random variabiles $A$ and $B$ with probability densities $f_a(x)$ and $f_b(x)$, the probability density $f_c(x)$ of their sum $A + B = C$ is given by the convolution $f_a * f_b$:\n",
    "\n",
    "$$\n",
    "f_c(x) = \\int_{-\\infty}^{\\infty} f_a(x') f_b(x - x') \\mathrm{d}x'\n",
    "$$\n",
    "\n",
    "For the case where $A$ and $B$ are uniformly distributed, we have:\n",
    "\n",
    "$$\n",
    "f_a(x) = \\begin{cases} \n",
    "      \\frac{1}{a_2 - a_1} & a_1 \\le x \\le a_2 \\\\\n",
    "      0 & \\mathrm{otherwise} \n",
    "   \\end{cases}\n",
    "$$\n",
    "\n",
    "$$\n",
    "f_b(x) = \\begin{cases} \n",
    "      \\frac{1}{b_2 - b_1} & b_1 \\le x \\le b_2 \\\\\n",
    "      0 & \\mathrm{otherwise} \n",
    "   \\end{cases}\n",
    "$$\n",
    "\n",
    "Using these definitions, let's visualize in the $X'-X$ plane where the integrand is nonzero:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEKCAYAAAAMzhLIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAARrklEQVR4nO3df6zddX3H8ecLqAIDhh1sVn6Im+BUMn5YQcEtpZggHRlZwh+4DTNi0uBkwUQzN/8Q9iNxJotxilIaNYpzElIRmYEZCTIgXZGW0UJtZjp10rSTwASsOEbhvT/OYbnc3vu95557vufXfT6Sm54f3++570/anGe/55z7vakqJEmazyGjHkCSNN4MhSSpkaGQJDUyFJKkRoZCktTIUEiSGrUeiiSHJvm3JN+c474k+VSS3Ul2JDm77XkkSYszjCOKa4Bd89x3MXBq92s9cMMQ5pEkLUKroUhyIvC7wOfm2eRS4Kbq2AIcm2RVmzNJkhbnsJYf/5PAnwFHz3P/CcBjM67v6d62b+ZGSdbTOeLgUA59y5EcM/hJJWmK/YyfPlFVx/ezb2uhSHIJ8HhVbUuyZr7N5rjtoHOKVNVGYCPAMVlZ5+bCgc0pScvBXbXpP/vdt82Xns4Hfi/Jj4CbgbVJ/mHWNnuAk2ZcPxHY2+JMkqRFai0UVfUXVXViVZ0CXA7cXVV/NGuz24H3dD/99Dbg6araN/uxJEmj0/Z7FAdJchVAVW0A7gDWAbuBZ4Erhz2PJKnZUEJRVfcA93Qvb5hxewHvH8YMkqT++JPZkqRGhkKS1MhQSJIaGQpJUiNDIUlqZCgkSY0MhSSpkaGQJDUyFJKkRoZCktTIUEiSGhkKSVIjQyFJamQoJEmNDIUkqZGhkCQ1MhSSpEaGQpLUyFBIkhoZCklSI0MhSWpkKCRJjQyFJKmRoZAkNTIUkqRGhkKS1MhQSJIatRaKJIcn+W6S7Ul2JvnLObZZk+TpJA93vz7a1jySpP4c1uJjPwesrar9SVYA9ye5s6q2zNruvqq6pMU5JElL0FooqqqA/d2rK7pf1db3kyS1o80jCpIcCmwDXg98pqoemGOztyfZDuwFPlRVO5se89nTVrH1xusGPqskTbULNvW9a6uhqKoXgDOTHAt8PcnpVfXojE0eAl7bfXlqHXAbcOrsx0myHlgPkNPe2ObIkpaBM87bf9Bt2zcfNYJJJkM6rxAN4Rsl1wI/r6q/a9jmR8Dqqnpivm2Oyco6Nxe2MKGkafatvQ833n/Ra84c0iSjcVdt2lZVq/vZt7UjiiTHA89X1VNJjgDeCXx81javBn5SVZXkHDqfwnqyrZkkLR8LhUG9a/Olp1XAl7rvUxwC3FJV30xyFUBVbQAuA96X5ADwC+DyGtYhjqSpYxza0eannnYAZ81x+4YZl68Hrm9rBknTyygMT6tvZkvSIBmH0TAUksaWYRgPhkLS2DAM48lQSBoZwzAZDIWkoTEMk8lQSGqVcZh8hkLSwBiF6WQoJPXNMCwPhkLSohiH5cdQSGpkGGQoJL2MYdBshkJa5gyDFmIopGXIOGgxDIW0DBgGLYWhkKaQYdAgGQppChgGtclQSBPIMGiYDIU0IYyDRsVQSGPKMGhcGAppDBgFjTNDIY2IcdCkMBTSkBgGTSpDIbXIOGgaGAppQIyCppWhkPpkGLRcGAppEYyDliNDIc3DKEgdhkKawThIB2stFEkOB+4FXtn9Ppuq6tpZ2wT4e2Ad8Czwx1X1UFszSbMZBmlhbR5RPAesrar9SVYA9ye5s6q2zNjmYuDU7te5wA3dP6XWGAdpcVoLRVUVsL97dUX3q2ZtdilwU3fbLUmOTbKqqva1NZeWF6MgLV2r71EkORTYBrwe+ExVPTBrkxOAx2Zc39O97WWhSLIeWA9wOEe2Nq8mn2GQBq/VUFTVC8CZSY4Fvp7k9Kp6dMYmmWu3OR5nI7AR4JisPOh+LW/GQWrXUD71VFVPJbkHeBcwMxR7gJNmXD8R2DuMmTS5DIM0XG1+6ul44PluJI4A3gl8fNZmtwNXJ7mZzpvYT/v+hGYzDNJotXlEsQr4Uvd9ikOAW6rqm0muAqiqDcAddD4au5vOx2OvbHEeTQjDII2XNj/1tAM4a47bN8y4XMD725pBk8EwSOPNn8zWSBgHaXIYCg2FYZAml6FQKwyDND0MhQbCMEjTy1CoL4ZBWj4MhXpiGKTly1BoToZB0ksMhf6fcZA0F0OxjBkGSb0wFMuIYZDUD0MxxQyDpEEwFFPEMEhqg6GYcMZBUtsMxQQxCpJGwVCMOeMgadTS+ZUQk+OwN7y5jr7xllGP0Zozzts/6hGkZWn75qNGPUKrnrrg9G1VtbqffQ8Z9DCSpOkycUcUx2RlnZsLRz1Ga3ypSRqNi15z5qhHaNVdtckjCklSOwyFJKmRoZAkNTIUkqRGhkKS1MhQSJIaGQpJUiNDIUlqZCgkSY1aC0WSk5J8J8muJDuTXDPHNmuSPJ3k4e7XR9uaR5LUnzbPHnsA+GBVPZTkaGBbkm9X1fdmbXdfVV3S4hySpCVY8IgiyZvmuG3NQvtV1b6qeqh7+WfALuCEPmaUJI1QLy893ZLkw+k4IsmngY8t5pskOQU4C3hgjrvfnmR7kjuTvHme/dcn2Zpk6/M8t5hvLUlaol5CcS5wErAZeBDYC5zf6zdIchTwNeADVfXMrLsfAl5bVWcAnwZum+sxqmpjVa2uqtUreGWv31qSNAC9hOJ54BfAEcDhwA+r6sVeHjzJCjqR+EpV3Tr7/qp6pqr2dy/fAaxIclyvw0uS2tdLKB6kE4q3Au8A3p1k00I7JQnweWBXVX1inm1e3d2OJOd053myx9klSUPQy6ee3ltVW7uX/wu4NMkVPex3PnAF8EiSl34bz0eAkwGqagNwGfC+JAfoxOjymrTfpCRJU27BUMyIxMzbvtzDfvcDWWCb64HrF3osSdLo+JPZkqRGhkKS1MhQSJIaGQpJUiNDIUlqZCgkSY0MhSSpkaGQJDUyFJKkRoZCktTIUEiSGhkKSVIjQyFJamQoJEmNDIUkqZGhkCQ1MhSSpEaGQpLUyFBIkhoZCklSI0MhSWpkKCRJjQyFJKmRoZAkNTIUkqRGhkKS1MhQSJIatRaKJCcl+U6SXUl2Jrlmjm2S5FNJdifZkeTstuaRJPXnsBYf+wDwwap6KMnRwLYk366q783Y5mLg1O7XucAN3T8lSWOitVBU1T5gX/fyz5LsAk4AZobiUuCmqipgS5Jjk6zq7junZ09bxdYbr2tr7JFbu3L/qEeQlqXtdx816hHadcGmvncdynsUSU4BzgIemHXXCcBjM67v6d42e//1SbYm2foi1daYkpap7ZunPBJL1OZLTwAkOQr4GvCBqnpm9t1z7HJQCapqI7AR4JisrNVrrxv0mGPjW3sfHvUI0tS76DVnvuz66hHNMUx3LWHfVkORZAWdSHylqm6dY5M9wEkzrp8I7G1zJknL0+w4qHethSJJgM8Du6rqE/NsdjtwdZKb6byJ/XTT+xOS1CvDMDhtHlGcD1wBPJLkpddTPgKcDFBVG4A7gHXAbuBZ4MoW55E0xQxDe9r81NP9zP0exMxtCnh/WzNIml6GYXhafzNbkgbFOIyGoZA0tgzDeDAUksaCURhfhkLSyBiHyWAoJA2NYZhMhkJSq4zD5DMUkgbGKEwnQyGpb4ZheTAUkhbFOCw/hkLSvIyCwFBImsU4aDZDIS1zhkELMRTSMmQctBiGQppyRkFLZSikKWMYNGiGQpoCxkFtMhTSBDIMGiZDIU0Aw6BRMhTSGDIMGieGQhoDhkHjzFBII2IcNCkMhTQkhkGTylBILTEMmhaGQhoQw6BpZSikPhkGLReGQuqRYdByZSikeRgGqaO1UCT5AnAJ8HhVnT7H/WuAbwA/7N50a1X9VVvzSL0wDtLB2jyi+CJwPXBTwzb3VdUlLc4gNTIM0sJaC0VV3ZvklLYeX+qHYZAWb9TvUbw9yXZgL/Chqto54nk0ZQyDtHSjDMVDwGuran+SdcBtwKlzbZhkPbAe4HCOHN6EmjiGQRq8kYWiqp6ZcfmOJJ9NclxVPTHHthuBjQDHZGUNcUxNAOMgtWtkoUjyauAnVVVJzgEOAZ4c1TyaDEZBGr42Px77VWANcFySPcC1wAqAqtoAXAa8L8kB4BfA5VXl0YIOYhyk0WrzU0/vXuD+6+l8fFZ6GcMgjZdRf+pJMgzSmDMUGjrDIE0WQ6HWGQZpshkKtcI4SNPDUGggDIM0vQyFFs0oSMuLoVBPjIO0fBkKzckwSHqJoRBgGCTNz1AsU4ZBUq8MxTJiHCT1w1BMMcMgaRAMxRQxDJLaYCgmmGGQNAyGYoIYBkmjkEn7XUGHveHNdfSNt4x6jNaccd7+l13fvvmoEU0iaZo8dcHp26pqdT/7ekQxZgyDpHEzcaE48vv7WL32ulGPIUkT5a4l7HvIwKaQJE0lQyFJamQoJEmNDIUkqZGhkCQ1MhSSpEaGQpLUyFBIkhoZCklSI0MhSWrUWiiSfCHJ40kenef+JPlUkt1JdiQ5u61ZJEn9a/OI4ovAuxruvxg4tfu1HrihxVkkSX1qLRRVdS/w3w2bXArcVB1bgGOTrGprHklSf0Z59tgTgMdmXN/TvW3f7A2TrKdz1AHw3F21ac6Xs6bEccATox6iRa5vck3z2mD61/eGfnccZSgyx21z/halqtoIbARIsrXfX74xCVzfZJvm9U3z2mB5rK/ffUf5qac9wEkzrp8I7B3RLJKkeYwyFLcD7+l++ultwNNVddDLTpKk0WrtpackXwXWAMcl2QNcC6wAqKoNwB3AOmA38CxwZY8PvXHgw44X1zfZpnl907w2cH3zStWcbwtIkgT4k9mSpAUYCklSo7ENRZJ3Jfn37ik+/nyO+yf6FCA9rO8Pu+vakWRzkjNGMWc/FlrbjO3emuSFJJcNc76l6mV9SdYkeTjJziT/MuwZl6KHf5u/nOSfkmzvrq/X9xdHbtpPLdTD+vp7XqmqsfsCDgX+A/h14BXAduBNs7ZZB9xJ5+cx3gY8MOq5B7y+84BXdS9fPCnr62VtM7a7m86HGi4b9dwD/rs7FvgecHL3+q+Oeu4Br+8jwMe7l4+ncwaGV4x69h7X9zvA2cCj89w/sc8rPa6vr+eVcT2iOAfYXVU/qKr/BW6mc8qPmSb5FCALrq+qNlfVT7tXt9D5OZNJ0MvfHcCfAl8DHh/mcAPQy/r+ALi1qn4MUFWTtMZe1lfA0UkCHEUnFAeGO2Z/aspPLbTQ+vp9XhnXUMx3eo/FbjOuFjv7e+n8L2cSLLi2JCcAvw9sGOJcg9LL391pwKuS3JNkW5L3DG26petlfdcDb6TzA7KPANdU1YvDGa91k/y8slg9P6+M8hQeTXo5vUfPpwAZQz3PnuQCOn+h72h1osHpZW2fBD5cVS90/lM6UXpZ32HAW4ALgSOAf02ypaq+3/ZwA9DL+i4CHgbWAr8BfDvJfVX1TNvDDcEkP6/0bLHPK+Mail5O7zHJpwDpafYkvwV8Dri4qp4c0mxL1cvaVgM3dyNxHLAuyYGqum04Iy5Jr/82n6iqnwM/T3IvcAYwCaHoZX1XAn9bnRe6dyf5IfCbwHeHM2KrJvl5pSf9PK+M60tPDwKnJnldklcAl9M55cdMk3wKkAXXl+Rk4Fbgign5n+hLFlxbVb2uqk6pqlOATcCfTEgkoLd/m98AfjvJYUmOBM4Fdg15zn71sr4f0zlaIsmv0Tkr6Q+GOmV7Jvl5ZUH9Pq+M5RFFVR1IcjXwLTqfwvhCVe1MclX3/qWcAmTkelzfR4FfAT7b/Z/3gZqAM1v2uLaJ1cv6qmpXkn8GdgAvAp+rqok4NX6Pf39/DXwxySN0Xqr5cFVNxOm5Wzy10FjoYX19Pa94Cg9JUqNxfelJkjQmDIUkqZGhkCQ1MhSSpEaGQpLUyFBIA5bklCT3jHoOaVAMhSSpkaGQlqD7OzV2JDk8yS8l2UnnFBBNZyiVJoo/cCctUZK/AQ6ncwLAPVX1sRGPJA2UoZCWqHtOpAeB/wHOq6oXRjySNFC+9CQt3Uo6v8DnaDpHFtJU8YhCWqIkt9P5TXCvA1ZV1dUjHkkaqLE8e6w0Kbq/ve5AVf1jkkOBzUnWVtXdo55NGhSPKCRJjXyPQpLUyFBIkhoZCklSI0MhSWpkKCRJjQyFJKmRoZAkNfo/+khhg6R7EW0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "a1 = 0.2\n",
    "a2 = 1.0\n",
    "b1 = 1.0\n",
    "b2 = 2.5\n",
    "\n",
    "Xlim = [1, 4]\n",
    "Xplim = [0, 1.2]\n",
    "Xp, X = np.meshgrid(np.linspace(*Xplim, 1000), np.linspace(*Xlim, 1000))\n",
    "nonzero = (a1 <= Xp) & (Xp <= a2) & (b1 <= X - Xp) & (X - Xp <= b2)\n",
    "\n",
    "plt.imshow(nonzero, extent=Xplim + Xlim, aspect='auto', origin='lower')\n",
    "plt.xlabel(\"x'\")\n",
    "plt.ylabel(\"x\")\n",
    "\n",
    "plt.axhline(b1 + a1, c='cyan')\n",
    "plt.axhline(b1 + a2, c='cyan')\n",
    "plt.axhline(b2 + a1, c='cyan')\n",
    "plt.axhline(b2 + a2, c='cyan');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The cross-section of the parallelogram is the range of $x'$ where the integrand is nonzero, so the lower and upper bounds of each cross-section define the limits of integration as a function of $x$.  Notice that it can be divided into three domains with nonzero probability, drawn here separated by the cyan lines. Going up from the bottom, they are defined by:\n",
    "\n",
    "| Domain | ---- Domain bounds ---- | Integration bounds |\n",
    "| --- | --- | --- |\n",
    "| 1 | $b_1 + a_1 \\le x \\le b_1 + a_2$ | $a_1 \\le x' \\le x - b_1$ |\n",
    "| 2 | $b_1 + a_2 \\le x \\le b_2 + a_1$ | $a_1 \\le x' \\le a_2$ |\n",
    "| 3 | $b_2 + a_1 \\le x \\le b_2 + a_2$ | $x - b_2 \\le x' \\le a_2$ |\n",
    "\n",
    "Finally, note that with the restricted bounds of integration $x_1$ and $x_2$, the integrand becomes the constant $(a_2-a_1)^{-1}(b_2-b_1)^{-1}$ and so the integral simplifies:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "f_c(x) &= \\frac{1}{(a_2-a_1)(b_2-b_1)} \\int_{x'_1}^{x'_2} \\mathrm{d}x' \\\\\n",
    "       &= \\frac{x'_2 - x'_1}{(a_2-a_1)(b_2-b_1)}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "where $x'_1$ and $x'_2$ are given by the integration bounds in the above table as a function of $x$.\n",
    "\n",
    "Note that we have assumed $b_1 + a_2 \\le b_2 + a_1$.  If that is not the case, swap the $a$s and $b$s.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEJCAYAAACE39xMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd3xUdb7/8dcnCVWqEBUpBhVUkKJEQFek2IJKh128Xvfq3SuXVezsT7DsZVdEV6kCgiiIdBEEEXWxEHpN6BGQLhGBECAQSDsz398fEzCEhEzCTL5TPs/HI49HZuY7c945Gd6cnDnne8QYg1JKqeAXYTuAUkop39BCV0qpEKGFrpRSIUILXSmlQoQWulJKhQgtdKWUChFeFbqIxInIThHZLSIDCnj8byKyKfdrm4i4RORK38dVSilVGCnqOHQRiQR+Bu4HkoH1wKPGmJ8KGd8JeNEY08HHWZVSSl1ClBdjWgK7jTF7AURkFtAFKLDQgUeBmUW9aM2aNU1MTIyXMZVSSgEkJiYeM8ZEF/SYN4VeGziY53Yy0KqggSJSEYgD+hXyeB+gD0C9evVISEjwYvFKKaXOEZEDhT3mzT50KeC+wvbTdAJWGmOOF/SgMWaCMSbWGBMbHV3gfzBKKaVKyJtCTwbq5rldBzhUyNjeeLG7RSmllO95U+jrgQYiUl9EyuIp7QX5B4lIVaAt8KVvIyqllPJGkfvQjTGOiPQDFgGRwCRjTJKI9M19fHzu0G7Ad8aYM35Lq5QqsZycHJKTk8nMzLQdRXmhfPny1KlThzJlynj9nCIPW/SX2NhYox+KKlV69u3bR+XKlalRowYiBX00pgKFMYbU1FROnz5N/fr1L3hMRBKNMbEFPU/PFFUqTGRmZmqZBwkRoUaNGsX+a0oLXakwomUePEryu9JCV0qpEKGFrpQqVfPmzUNE2LFjx2W9zhNPPMGcOXMuOWbIkCEX3L7rrrtKtKxBgwYxdOjQAu+vXbs2zZs359Zbb2XBggUX3d+gQQO6d+/OTz/9fnJ9u3btuOmmm2jevDnNmzcv8ufwljdniipVfFvnwJfPgKNHVJQegTYvwb1/tx3kkmbOnMndd9/NrFmzGDRokF+XNWTIEF599dXzt1etWuXzZbz44ov079+f7du306ZNG44ePXrB/QCfffYZHTp0YOvWrZw7qXL69OnExhb42WaJaaEr38s+A/8eCNXqQePu5+8e+ePPpRbhhXsbWlnu5bqc3C/cmAIrRkDTP0H0Tb6O5hPp6emsXLmS+Ph4OnfufL7QlyxZwqBBg6hZsybbtm2jRYsWTJs2DRHhn//8J1999RUZGRncddddfPjhhxfsX/7xxx8ZM2YM8+bNA+D7779n3LhxNGzYkIyMDJo3b07jxo2ZPn06lSpVIj09HYB3332XqVOnEhERQceOHXnnnXf46KOPmDBhAtnZ2dx4441MnTqVihUrevWz3XLLLURFRXHs2LGLHvvTn/7E119/zYwZM3j++ecvcy0WTgtd+dzb/3iZgWWO0v3402xIbpjnkaallmHkory3Sm+5l+tycr/QqxWMbApL3oZeky89+NsBcHhrceNd2jVNoOM7lxwyf/584uLiaNiwIVdeeSUbNmzg9ttvB2Djxo0kJSVx7bXX8oc//IGVK1dy9913069fP/7+d89fHY8//jgLFy6kU6dO51+zQ4cOPPPMM6SkpBAdHc0nn3zCk08+SadOnRgzZgybNm26KMe3337L/PnzWbt2LRUrVuT4cc9sJd27d+epp54C4PXXX2fixIk8++yzXv34a9euJSIigsKmNbn99tsv2M302GOPUaFCBcDzn1KNGjW8Ws6l6D505VuZp/jfqK+IdzVjg2lY9HjlMzFvrmV0xv2QNA8Ob7Mdp0AzZ86kd+/eAPTu3ZuZM3+fKaRly5bUqVOHiIgImjdvzv79+wGIj4+nVatWNGnShMWLF5OUlHTBa4oIjz/+ONOmTePkyZOsXr2ajh07XjLHDz/8wJNPPnl+6/vKKz2Xb9i2bRtt2rShSZMmTJ8+/aJlFWTEiBE0b96c/v3789lnnxV6dEr+c36mT5/Opk2b2LRpk0/KHHQLXfna2vFcKekMd3rZThKWPnIe4r8iv6NK/BB4dEbhA4vYkvaH1NRUFi9ezLZt2xARXC4XIsK7774LQLly5c6PjYyMxHEcMjMzefrpp0lISKBu3boMGjSowGOzz22Rly9fnl69ehEVdelqM8YUWLxPPPEE8+fPp1mzZkyePJklS5YU+XPl3Vd+KRs3bvT5PvP8dAtd+U7GCVg1hkWuWLaa622nCUunqMRHzkOw82v4dYPtOBeYM2cOf/7znzlw4AD79+/n4MGD1K9fnxUrVhT6nHPlXbNmTdLT0ws9GuTaa6/l2muvZfDgwTzxxBPn7y9Tpgw5OTkXjX/ggQeYNGkSZ8+eBTi/y+X06dPUqlWLnJwcpk+fXtIf9SJz587lu+++49FHH/XZaxZEC135zqoxkJXGCKen7SRh7RNXHFSoDvFv2Y5ygZkzZ9KtW7cL7uvRowczZhT+l0S1atV46qmnaNKkCV27duWOO+4odOxjjz1G3bp1adSo0fn7+vTpQ9OmTXnssccuGBsXF0fnzp2JjY2lefPm5w9JfPPNN2nVqhX3338/N998c0l+zPPO7Ypp0KAB06ZNY/HixYXuX/cVnctF+caZYzCqGTR4gJhELXTb9j+yC374P/jvRVCvNQDbt2/nlltusZzMf/r168dtt93GX/7yF9tRfKag35nO5aL8b+VIyDkL7QbaTqIAWj4FV0TD4sG2k5SKFi1asGXLFv7zP//TdhSrtNDV5Tt9GNZ9nHv8sx7ZEghi/r6Ef5yMg/3LYe9S23H8LjExkWXLll3wwWo40kJXl2/5cJycLO5Z15KYAV/bTqNyzXDdy2/mSs++dEu7VlXp0kJXlyctGRI/YbarLb+Yq22nUXlkUZYxTlc4uBZ2/2g7jioFWujq8ix7D4AxTrciBiobZrvaeaZgiA+PfenhTgtdldzxfbBxGrR4gkPUtJ1GFSCHKGj7ChzaCDkZtuMoP9MzRVWJzRnxHI9ECG2WNbMdRV1K096wfDhkpnn2peeeIenrzzv2v/OwT1/PHxYsWMBPP/3EgAEDiv3cmJgYEhISqFmz5kX3V65cmYiICK6++mqmTJnCNddcc/5+AJfLRffu3XnjjTcoV64c+/fv55ZbbuGmm36fRG3dunWULVv2sn4+3UJXJZPyM90iljPF9QApVLedRl1CzGuLeO5wHLiyIfOk7TjWOI5D586dS1TmRYmPj2fz5s3ExsZeMAd7fHw8W7duZd26dezdu5c+ffqcf+yGG244P5fLpk2bLrvMQQtdldTSd8ikLOOdTkWPVdYtdN/p2f1y6jdrR7yc2yp96qmnaNy4MQ888AAZGZ7dQO3atePciYbHjh0jJiYGgMmTJ9O1a1c6depE/fr1GTNmDMOHD+e2226jdevW50/Z37NnD3FxcbRo0YI2bdqcn9XwiSee4KWXXqJ9+/a88sorTJ48mX79+gFw5MgRunXrRrNmzWjWrNn5udK7du1KixYtaNy4MRMmTCjWz3jPPfewe/fui+6vVKkS48ePZ/78+ecz+4MWuiq+I0mwbS6fuOI4ThXbaZQX3ERwylQEVxZk+K9QirJr1y6eeeYZkpKSqFatGnPnzi3yOdu2bWPGjBmsW7eO1157jYoVK7Jx40buvPNOpkyZAnhO8R89ejSJiYkMHTqUp59++vzzf/75Z3744QeGDRt2wes+99xztG3bls2bN7NhwwYaN24MwKRJk0hMTCQhIYH333+f1NRUr3++hQsX0qRJkwIfq1KlCvXr12fXrl2A5z+hc1cseuaZZ7xexqV4tQ9dROKAUUAk8LEx5qKp2kSkHTASKAMcM8a09UlCFXjih0C5qkzIDPx9pup3mZSFqAqeE8Eq2NlNVr9+fZo3bw54zu48N0XupbRv357KlStTuXJlqlaten4u9CZNmrBlyxbS09NZtWoVvXr9PsNnVlbW+e979epFZGTkRa+7ePHi8/8hREZGUrVqVQDef//98xfLOHjwILt27Spyetv27dsTGRlJ06ZNGTy48COK8k61cm6Xiy8VWegiEgmMBe4HkoH1IrLAGPNTnjHVgA+AOGPMLyJylU9TqsBxaCPsWAjtXuXUvyvZTqOKwSBQpRYc3wtn7Wyl558i99wul6ioKNxuN8BF0+PmfU5ERMT52xERETiOg9vtplq1aoWW4xVXXOF1viVLlvDDDz+wevVqKlasSLt27Qqcrje/+Pj4iz4sze/06dPs37+fhg0bkpaW5nWm4vBml0tLYLcxZq8xJhuYBXTJN+Y/gC+MMb8AGGOO+jamChjxQzxbd63/ajuJKoEtxwxnTTmyT/5mO8oFYmJiSExMBCj2BZPP7cr4/PPPAc9W8ObNm4t83r333su4ceMAz1Eop06dIi0tjerVq1OxYkV27NjBmjVrivmTFCw9PZ2nn36arl27Ur26//468maXS23gYJ7byUCrfGMaAmVEZAlQGRhljJmS/4VEpA/QB6BevXolyatsOrgOdn3HOzm9GT9oue00qoQOm+pcH3GY/a+3hEr+nc7VW/379+ePf/wjU6dOpUOHDsV+/vTp0/nrX//K4MGDycnJoXfv3jRrdunDaUeNGkWfPn2YOHEikZGRjBs3jri4OMaPH0/Tpk256aabaN26dUl/JMCzK8YYg9vtplu3brzxxhuX9XpFKXL6XBHpBTxojPmf3NuPAy2NMc/mGTMGiAXuBSoAq4GHjTGFXuVWp88NQp92JmXvRu7JGkkG5W2nUcX0UedaXF3Pc+GR6+U3KkU6cFUjiLh4/7IKDP6YPjcZqJvndh3gUAFj/m2MOWOMOQYsA/Rsk1CybznsW8o4p4uWeQg4YqqD24GzF1+hXgUvbwp9PdBAROqLSFmgN7Ag35gvgTYiEiUiFfHsktnu26jKGmM8M/ZVvpbprnttp1E+cIbyULYypB8Ft8t2HOUjRRa6McYB+gGL8JT0bGNMkoj0FZG+uWO2A/8GtgDr8BzaGJiXHVfFt+dH+GU13PMyWVz+2WzKDoO58MrzVWp5ttLPpNgLpQpVkqvJeXUcujHmG+CbfPeNz3f7PeC9YidQgc0YNk35GzWlJu3nBsYHaKpkDpzMoUaNU0RVrIKIsOVoDjFSkYqnjhB1RU2I0KmdAoUxhtTUVMqXL97uTf0Nqkvb+S3NI/byt5w+nlPHVdAavfYEzwLXVTuG4Jmg6wQOV8kJOJoB5avaDaguUL58eerUqVOs5+i/UFU4txvih7DPfTVfuNrYTqMu06ksN28tu/g09g/KjOShCtvhhS1Q8UoLyZSv6FwuqnDbv4QjWxnp9MCFHtoWqkY4PSE7HVaOsh1FXSYtdFUwtwvi34bom/nKfZftNMqPdpk60KQnrJvgOepFBS0tdFWwrXPg2E5oNxC3vk1CX9sB4GTBihG2k6jLoP9S1cVcOeyf+wZJ7uuoP1V3tYSDmKE7mZ1zN1mrP6LVgKm246gS0kJXF9s8k5iIIwxzemH0LRI23nd1R3DTL2q+7SiqhPRfq7qQkwVL32Wj+0YWu2+znUaVomQTzWxXO/4UGQ8nDtiOo0pAC11daMMUSDvIMKcX5B6rrMLHGKer56+yZe/ajqJKQAtd/S4nA5YNhXp3scJ9q+00yoLD1PDM17NpJqTusR1HFZMWuvpdwiRIPwwdXke3zsPXOKczRJWDJRddaVIFOC105ZGVzrF/v8Ny163EjD9pO42yKIVqjM+4F/eWz7l/4Ie246hi0EJXHusmUFNOMdzpVfRYFfLGO49whvK8EFW8y8Epu7TQFWSmwcpR/Oi6jY2mge00KgCcpDKTXHE8HLkOfttiO47ykha6gjXjIPMkw52etpOoADLReYg0U9FzYXAVFLTQw93Z47B6LNzSiSRT33YaFUBOcQUTnEfg528hWa//Gwy00MPdqtGQdRravWo7iQpAk10PQsUanksQqoCnhR7O0lM4u3wsX7ruJGbEPttpVAA6QwUGp8XBnsX0GjjMdhxVBC30cLZyJOXIZqTTw3YSFcCmue7jqKlG/zKzPRcMVwFLCz1cnToE6z/mC1cb9plattOoAJZJOcY4XWgVsQP2LrEdR12CFnq4Wj4M3A6jXN1tJ1FBYJarA7+a3H3pupUesLTQw9HJXyDxU7jtcZLNVbbTqCCQTRlGO90geT3s+s52HFUIrwpdROJEZKeI7BaRAQU83k5E0kRkU+7X330fVfnM0ndBIuCev9lOooLIHNc9UD0GFg/WrfQAVWShi0gkMBboCDQCHhWRRgUMXW6MaZ779U8f51S+kroHZ8N0PslqR8zbm2ynUUHEIYqXjsTB4S3872uDbMdRBfBmC70lsNsYs9cYkw3MArr4N5bym6X/IocoPnD0V6iKb777bva4a/FS1BzPhcRVQPGm0GsDB/PcTs69L787RWSziHwrIo0LeiER6SMiCSKSkJKSUoK46rIc3QFbZvOp60FSqGY7jQpCbiIY4fTkpohkSJpnO47Kx5tCL2hi7Pw70DYA1xljmgGjgQIvSmiMmWCMiTXGxEZHRxcvqbp8S96GslfwofOw7SQqiH3tbsV2d13P+8nl2I6j8vCm0JOBunlu1wEO5R1gjDlljEnP/f4boIyI1PRZSnX5Dm+Fn+ZD66c5QRXbaVQQM7lb6aTuhi2f2Y6j8vCm0NcDDUSkvoiUBXoDC/IOEJFrRERyv2+Z+7qpvg6rLkP8EChfFe58xnYSFQK+c8eyxV2fg/P+jwYDvrQdR+UqstCNMQ7QD1gEbAdmG2OSRKSviPTNHdYT2CYim4H3gd7G6HFNAePXRNj5De+dfpCYf6y0nUaFBGG404u6ESn0ilxqO4zKJbZ6NzY21iQk6JScpWJqd47vXkubrFGcoYLtNCpkGOaWHcS1kkqt17dDmfK2A4UFEUk0xsQW9JieKRrqDqyGPT8yzumsZa58TBjq/JFachwSJ9sOo9BCD33xb0Glq5nqut92EhWCVrsbs9rVyDM3UPZZ23HCnhZ6KNu7FPYvhzYvk0k522lUiBrm9IQzR2H9R7ajhD0t9FBljGfOjSq14fb/sp1GhbAEczPceB+sGOm5+pWyRgs9RD3x2hBIXsfA1Dhi3vjRdhwV4jontYWM4wz954u2o4Q1LfRQZAwvR83mF3c0n7va2k6jwsAWcwPfuVrQJ+pryDhhO07Y0kIPRTsW0iRiP6OcHjhE2U6jwsRwpxdV5CysHms7StjSQg81bjfED2GPuxbz3X+wnUaFkR2mHgtdrWHNODijJ4rboIUeapK+gKM/MdLpgYtI22lUmBnh9ICcs7BypO0oYUkLPZS4HFjyDlzViIXu1rbTqDC0x9RmrnMXGSvHc8eA6bbjhB0t9BDy8huvQeou/jf5QYz+apUl7zvdKYPD01E6aVdp03/1ocKVw/NRc9nqjmGRu8BpHpQqFQfMNXzuast/RP4Iacm244QVLfRQsXEa9SJSGOb0ouBrkihVesY4XREMLBtqO0pY0UIPBTmZsOw9Et0NWOJubjuNUvxKNDNdHWDjVDi+z3acsKGFHgo2fAqnftWtcxVQxjpdISIKlr1nO0rY0EIPdtlnPTPdxbRhlbvAa3MrZcVRqvNRZgdcG2fQYaBO3FUatNCD3Fv/eBnSj9BzZ3t061wFmvFOJzIpywtRc21HCQta6MEs6zR9o75iqaupZ8Y7pQJMKlWZ7HqQRyLWwJEk23FCnhZ6MFs7nhpymuFOT9tJlCrUBOcR0invuVC58ist9GCVcRJWjeZ7Vws2mxttp1GqUGlUYqLzEOxYCIc22o4T0rTQg9XqsZCZplvnKihMcnWECtV1K93PtNCD0ZlUWPMBNOrKdnOd7TRKFek0FeGu52DXd3Bwne04IcurQheROBHZKSK7RWTAJcbdISIuEdHNRj8a//bzuLPOcN/Gu2xHUcprt3x9HSmmCismvGA7SsgqstBFJBIYC3QEGgGPikijQsb9C1jk65Aqj9NH+K/IRcx3/4Hdpo7tNEp5LYPyjHc6c3dkEuxbbjtOSPJmC70lsNsYs9cYkw3MAroUMO5ZYC5w1If5VH4rRlAGh1FOd9tJlCq2aa77OGyqQ/xbnguZK5/yptBrAwfz3E7Ove88EakNdAPGX+qFRKSPiCSISEJKSkpxs6q0XyFhInNc93DAXGM7jVLFlkVZxjhd4ZfVsGex7Tghx5tCL+j0w/z/tY4EXjHGuC71QsaYCcaYWGNMbHR0tLcZ1TnLh4IxjHa62U6iVInNdrWDqvVg8WDdSvcxbwo9Gaib53Yd4FC+MbHALBHZD/QEPhCRrj5JqAC4e+An5Kz/lCnZ7fgV/c9QBa9syvD/jj0Ihzbwl9fetB0npHhT6OuBBiJSX0TKAr2BBXkHGGPqG2NijDExwBzgaWPMfJ+nDWPPRc7DTQRjnYI+vlAquHzhasM+99W8HDXHc2Fz5RNFFroxxgH64Tl6ZTsw2xiTJCJ9RaSvvwMq4NhuekQuY6rrPo5wpe00Sl02hyhGOT1oFHEAti8o+gnKK2Is7cOKjY01CQkJVpYddOb+D2e3LKBN1khSqWo7jVI+EYGbRWVfocFVleDp1RARaTtSUBCRRGNMgdeZ1DNFA92Rn2DrHCa7HtQyVyHFTQQjnB5wbCds0+l1fUELPdAteRvKVeZD5xHbSZTyuW/dLeHqJp73ucuxHSfoaaEHsIcHjoXtCxh55j7SqGQ7jlI+Z4iA9q/C8b2weabtOEFPCz2AvRj1OSfNFZ6pR5UKUTGT3WxyX0/yl4PAybIdJ6hpoQeqg+u5L3IjE5xHPDPVKRWyhGHOH6kjx2DDFNthgpoWeqCKf4tjpgqTXQ/aTqKU3y13N2Gd+ybPBc9zMmzHCVpa6IFo/wrYG884pxNnKW87jVKlQBiW80c4/RskTLIdJmhpoQcaY2DxW1DpGqa57redRqlSs9bcAvXbwooRkJVuO05Q0kIPNHvj4ZdVcE9/sihrO41SpavD63AmBdZNsJ0kKGmhBxJj2PRpf5JNTRp+oRNwqfATMzaFxa7mnPxhKGSm2Y4TdLTQA8nPi2gesYfRTjeyKWM7jVJWDHd6Uk3OwJpxtqMEHS30QOF2Q/xg9ruvZq6rje00SlmzzVzPv113wOqxcPa47ThBRQs9UOz4Cg5vZZTTHYco22mUsmqE0wOyTsOq0bajBBUt9EDgdkH8EKjZkC/df7CdRinrdpp6LHC15szysbQYMMN2nKChhR4Annv9DUjZwTOH4nDrr0QpAEY6PShPNn2jvrIdJWhoe9jmcnghai7b3fX4xt3SdhqlAsZecy3z3G14PPJ7OPWb7ThBQQvdti2zuD7iMMOdnp6Z55RS541yuhGJ2zMlgCqSNohNTjYs+Reb3dfzvbuF7TRKBZyD5mo+d7WFxMlw8hfbcQKeFrpNG6dC2i8Md3oBYjuNUgFptNMNRGDZe7ajBDwtdFtyMmHZUKjbmqXuprbTKBWwfqMGtHgSNk6H1D224wQ0LXRL/jGoP5w+xKN77kW3zpW6tDuWNSXDHcnckc/ZjhLQtNBtyD7D01FfstLVmNXuxrbTKBXwUqjOp64H6BqxElJ22o4TsLwqdBGJE5GdIrJbRAYU8HgXEdkiIptEJEFE7vZ91BCybgLRcophTi/bSZQKGh86j5BBOc8FpVWBiix0EYkExgIdgUbAoyLSKN+wH4FmxpjmwH8DH/s6aMjIPAUrRxHvasYG09B2GqWCxgmqMMkVB0nz4PBW23ECkjdb6C2B3caYvcaYbGAW0CXvAGNMujHG5N68AjCogq0ZBxkndOtcqRL42HkIylWFeN1KL4g3hV4bOJjndnLufRcQkW4isgP4Gs9W+kVEpE/uLpmElJSUkuQNbmePw+oxcPMjbDPX206jVNA5RSW461nY+TX8mmg7TsDxptALOgTjoi1wY8w8Y8zNQFfgzYJeyBgzwRgTa4yJjY4Ovws4jBnyAu7M0zy4WT9iUKqkGn9bn+OmEkvGv2g7SsDxptCTgbp5btcBDhU22BizDLhBRGpeZrbQcuYYT0b+m6/drdhp6tlOo1TQOkMFxjudaBe5GX5ZYztOQPGm0NcDDUSkvoiUBXoDC/IOEJEbRURyv78dKAuk+jpsUFsxgvJkM9LpYTuJUkFviusBUkxVWDzYdpSAUmShG2McoB+wCNgOzDbGJIlIXxHpmzusB7BNRDbhOSLmT3k+JFWnD8P6j5nvvps95qKPH5RSxZRJOcY6XWD/cti71HacgCG2ejc2NtYkJCRYWXap++ZvkDCJezLe4xdzte00SoWEcmSz86rXoGod+O9FnvlewoCIJBpjYgt6TM8U9beTBz0zxTV/TMtcKR/Koizc0x8OroXdP9iOExC00P1sxtBnyXLc3LVKp8dVytcazKnBQXc0W6b+DXQvrxa6Xx3fS6/IZcx0deAQetCPUr6WQxSjXN1pGrEPdnxtO451Wuj+tPRdXER4PrxRSvnFPNfd7HHX8lxo3e22HccqLXR/SfkZtnzmObyK6rbTKBWyXEQyyukBR5Pgp3m241ilhe4vS96GKM8JEEop//rK3RquauSZ48Xl2I5jjRa6PxzeBklfQOu+HKeK7TRKhTxDBLQbCKm7YOvntuNYo4XuB4vGPs8pU5Gm399kO4pSYSNmSgTb3DEc+OLv4MqxHccKLXRf+3UDD0Ym8LHzkGdmOKVUKRGGOb24LuIobJpuO4wVWui+Fj+EE6aSZyJ+pVSpinc3Z4P7Rlj6HjhZtuOUOi10X/plLez+nvFOJ9KpaDuNUmHIs5XOqWRI/NR2mFKnhe5L8YPhimimuO63nUSpsLXSfStcdzcsHwrZZ23HKVVa6D7y6Kvvwr5l/ONkHBmUtx1HqTAm0OE1SD8CCRNthylVWui+YAwvRX3Ob+ZKZrjutZ1GqbAXM+4Ey1xNSF30L8g6bTtOqdFC94XdP3JHxM+Mcbp6ZoBTSlk3zOlFDTkNaz+0HaXUaKFfLmMgfjAH3dHMdrWznUYplWuzuZHvXbfDqvch46TtOKVCC/1y7acoirAAAA9RSURBVPwGDm3kfVc3coiynUYplccIpydkpsGaD2xHKRVa6JfD7fbM8HblDXzhamM7jVIqn59MDDTqAqs/gLPHbcfxOy30y/HTfDiyDdoNxEWk7TRKqYK0GwjZ6bBylO0kfqeFXlJuF7tnv8bP7tpcP0MPU1QqUMUM38s8112cXTEO0o/ajuNXWugltfVzbow4xHCnF25djUoFtFFOd8qSAytG2I7iV9pEJeHKgSXvkOS+jkXuAi++rZQKIPtNLea67oH1EyHtV9tx/MarQheROBHZKSK7RWRAAY8/JiJbcr9WiUgz30cNIJtmwIl9DHN6eeZhVkoFvNGubmDcsHyY7Sh+U2QbiUgkMBboCDQCHhWRRvmG7QPaGmOaAm8CE3wdNGA4WbD0Xagdy2L3bbbTKKW8lGyi4fY/w4YpcOKA7Th+4c3mZUtgtzFmrzEmG5gFXHDVY2PMKmPMidyba4A6vo0ZON74v/8Hp5L5z333A2I7jlKqOO7pDxIBy961ncQvvCn02sDBPLeTc+8rzF+Abwt6QET6iEiCiCSkpKR4nzJQ5GTQL2o+a903s8J9q+00SqliihmykYlZHXA2zIDUPbbj+Jw3hV7QZqgpcKBIezyF/kpBjxtjJhhjYo0xsdHR0d6nDBTrJ3K1nGRYTi9061yp4DTO6Uw2ZWDJO7aj+Jw3hZ4M1M1zuw5wKP8gEWkKfAx0Mcak+iZeAMlKhxXDWe66lXXmFttplFIldIyqfOp6wHMx6aPbbcfxKW8KfT3QQETqi0hZoDewIO8AEakHfAE8boz52fcxA8C6D+FsKsOdXraTKKUu04fOI1C2Eix523YUnyqy0I0xDtAPWARsB2YbY5JEpK+I9M0d9negBvCBiGwSkQS/JbYhMw1Wvg8N49hoGthOo5S6TCepDHc+DT99Cb9tth3HZ8SYAneH+11sbKxJSAiO3h/5+pO8EPUFD2cNIcnE2I6jlPKBKpxhWbkXSHA35L43l9iO4zURSTTGFHhGo54VU5Szx/lL5Ld842qpZa5UCDnFFUxwHua+yI2QHBwbl0XRQi/Kqve5gkzPvMpKqZAy2RVHqqkMiwfbjuITWuiXkn4U1n7IAved7DIhe66UUmHrLOUZ53SGvfGwf6XtOJdNC/1SVowEJ5NRTg/bSZRSfjLNdR9Uugbi3/JcUjKIaaEX5tQhWP8xNPsP9plattMopfwkk3LQ5mU4sBL2LrEd57JooRdm+TAwLmj7N9tJlFL+1uK/oEodz770IN5K10IvyMlfyF73CdOy2xLzryTbaZRSfhbz+g8MSI2DXxNg13e245SYFnpBlr6LIYIxTlfbSZRSpWSO6x4OuK8K6q10LfT8UvfAphlMd93LYWrYTqOUKiUOUYx0esDhLbD9K9txSkQLPb8l70BkWc+hTEqpsPKl+w9QsyHEDwG3y3acYtNCz+voDs8MbK36kEI122mUUqXMTQS0Gwgp2yFpnu04xaaFnsfXo58l3ZTjth91elylwlX9aWXY7q7H3s9fA5djO06xaKGf89sWHo5cx0RXR05QxXYapZQlhgiGOz25PuIwbPnMdpxi0UI/J34IaaYiE52HbCdRSln2vbsFm93Xw9J3wMm2HcdrWugAyYnw87dMcB7hFFfYTqOUsk48F7M5+QtsmmY7jNe00AHiB0OFK5nsetB2EqVUgFjqbgp1W8HS9yAn03Ycr2ihH1gFexbD3S9yhgq20yilAoZAh9fh9CFInGw7jFfCu9CNYc3ElzhqqnHzVzo9rlLqQjEfnmaVqxEp3w6B7DO24xQpvAt931JaR2xnrNPFM+OaUkrlM8zpRbSkwbqPbEcpUvgWujGweDC/mhrMdHWwnUYpFaASzU0scTWDlaMg85TtOJcUvoW+63tIXs8YpyvZlLGdRikVwIY5vSDjOKwdbzvKJYVnoRvjObKlegyfu9raTqOUCnBbzfVw08OwagxknLAdp1BeFbqIxInIThHZLSIDCnj8ZhFZLSJZItLf9zF9bMdC+G0ztH0FhyjbaZRSwaD9q5CV5in1AFVkoYtIJDAW6Ag0Ah4VkUb5hh0HngOG+jyhr7ld7Jw5gD3uWtwwq5LtNEqpIBEz8gBfuVqTvmwMtw+YaTtOgbzZQm8J7DbG7DXGZAOzgC55Bxhjjhpj1gM5fsjoW0nzuCkimZFOD1xE2k6jlAoiI50eVCCL/40KzPnSvSn02sDBPLeTc+8rNhHpIyIJIpKQkpJSkpe4PC4HlrzNDnddFrpbl/7ylVJBbY+pzXz33fw58ns4fdh2nIt4U+hSwH0luj6TMWaCMSbWGBMbHR1dkpe4PFs+g9TdjHB6YsL082Cl1OUZ5XSnDA4sH247ykW8abVkoG6e23WAQ/6J40dONiz9F9RqxiJ3rO00Sqkg9Yu5mtmutpD4CaQl245zAW8KfT3QQETqi0hZoDewwL+x/GDTNDh5ANq/TsF/dCillHfGON083yx7z26QfIosdGOMA/QDFgHbgdnGmCQR6SsifQFE5BoRSQZeAl4XkWQRCZyrRORkcuirN0l0NyBmUvDMbayUCkyHqMnkrHbkJEylzcBJtuOc59VB2MaYb4Bv8t03Ps/3h/HsiglMiZO5Vo7TP6cvunWulPKFsU4XekfG83zUPOC/bccBwuFM0eyzsHwYq12NWOVubDuNUipEpFCdKa4H6BaxHI7tsh0HCIdCX/8RnDnKMKcnunWulPKl8U4nMikLS962HQUI9ULPOg0rRsIN95JgbradRikVYo5ThU9ccbDtCziSZDtOiBf6mvGeGdLav2Y7iVIqRH3kPAzlKkP8ENtRQrjQM05wavFwvne1IGZM4J3RpZQKDWlUYnj6A7BjIY8MHG01S+gW+uqxVJGzDHd62k6ilApxk1xxnDCVeClqjtUcoVnoZ1JhzTgWulqx3VxnO41SKsSlU5EPnUfoELkJDq6zliM0C33lSMg5ywjdOldKlZJPXQ+QYqrA4sHWMoReoZ8+4rmYa5Ne7DElmhRSKaWKLYPyjHO6wL6lsG+5lQyhV+grhoMrG9q+YjuJUirMTHfdC5Wvhfi3PJe6LGWhVehpyWSt+ZhZOW2IeW+H7TRKqTCTRVleP/4g/LKaP79W+icbhVahLxuKYBh9biY0pZQqZZ+52pNsavJS1OelvpUeOoV+fB9snMosVwd+xcLFM5RSCsghilFOd5pH7IWd35bqskOn0Je9BxFRjHG62k6ilApzX7jasM99tefsUbe71JYbGoV+bBdsngmxf+Eo1W2nUUqFOReRjHR6wJGtsP3LUltuaBT6kncgqjzc/aLtJEopBcBX7rsg+maIfxvcrlJZZtAX+gMDx+PeOpcPMu4jZrC9M7SUUiovNxH89dcH4dhOXnj99VJZZtAX+otRc0inPB86j9iOopRSF/i3+w6S3NfxfNRccOX4fXnBXeiHNtExcj2TXB1Jo5LtNEopdQFDBMOdntSPOOL5nM/PgrvQ44dw0lzBROch20mUUqpAP7pvZ5P7Blj6LjhZfl1W8Bb6wXWwaxETnEc4TUXbaZRSqhDCMKcXpB2EDVP8uqTgLfT4t6BiTSa7HrSdRCmlLmm5uwnUuwuWD4OcDL8tx6tCF5E4EdkpIrtFZEABj4uIvJ/7+BYRud33UX/3p4FDYe8S3kx7kLOU9+eilFLKB4Q/7boXTv/Gm4P6+20pRRa6iEQCY4GOQCPgURFplG9YR6BB7lcfYJyPc/7OGF4q8zlHTDWmue7322KUUsqX1ppbWO66lb5RCyAr3S/L8GYLvSWw2xiz1xiTDcwCuuQb0wWYYjzWANVEpJaPs3rsWUyriB2McbqSRVm/LEIppfxhuNOLaDkF6yb45fWjvBhTGziY53Yy0MqLMbWB3/IOEpE+eLbgAdJFZGex0l7g/dyvgFQTOGY7RADQ9eCh6+F3Yb0uDgAC8I+Xa8LLJV0PhV5X05tClwLuyz8npDdjMMZMAPzzX1MAEZEEY0ys7Ry26Xrw0PXwO10XHv5aD97sckkG6ua5XQc4VIIxSiml/MibQl8PNBCR+iJSFugNLMg3ZgHw59yjXVoDacaY3/K/kFJKKf8pcpeLMcYRkX7AIiASmGSMSRKRvrmPjwe+AR4CdgNngSf9FzkohPxuJS/pevDQ9fA7XRceflkPYixcyFQppZTvBe+ZokoppS6gha6UUiFCC72ERGSSiBwVkW2FPN5ORNJEZFPu199LO2NpEJG6IhIvIttFJElEni9gTKlODWGDl+sh5N8TIlJeRNaJyObc9fCPAsaE/PsBvF4Xvn1PGGP0qwRfwD3A7cC2Qh5vByy0nbMU1kMt4Pbc7ysDPwON8o15CPgWz/kKrYG1tnNbWg8h/57I/R1Xyv2+DLAWaB1u74dirAufvid0C72EjDHLgOO2c9hmjPnNGLMh9/vTwHY8ZwnnVXpTQ1ji5XoIebm/43MTlZTJ/cp/5EXIvx/A63XhU1ro/nVn7p9b34pIY9th/E1EYoDb8GyJ5FXY1BAh6RLrAcLgPSEikSKyCTgKfG+MCdv3gxfrAnz4ntBC958NwHXGmGbAaGC+5Tx+JSKVgLnAC8aYU/kfLuApIXm8bBHrISzeE8YYlzGmOZ4zxluKyK35hoTN+8GLdeHT94QWup8YY06d+3PLGPMNUEZEalqO5RciUgZPiU03xnxRwJCwmBqiqPUQTu8JAGPMSWAJEJfvobB4P+RV2Lrw9XtCC91PROQaEZHc71viWdepdlP5Xu7POBHYbowZXsiwkJ8awpv1EA7vCRGJFpFqud9XAO4DduQbFvLvB/BuXfj6PeHNbIuqACIyE88n1DVFJBn4PzwfemA80yH0BP4qIg6QAfQ2uR9rh5g/AI8DW3P3FQK8CtSDsJoawpv1EA7viVrAp+K5ME4EMNsYszBMpwrxZl349D2hp/4rpVSI0F0uSikVIrTQlVIqRGihK6VUiNBCV0qpEKGFrpRSIUILXSmlQoQWulJKhQgtdKVyicgdufNzlxeRK3LnsM4/94ZSAUtPLFIqDxEZDJQHKgDJxpi3LUdSymta6ErlISJlgfVAJnCXMcZlOZJSXtNdLkpd6EqgEp6rDpW3nEWpYtEtdKXyEJEFwCygPlDLGNPPciSlvKazLSqVS0T+DDjGmBm5M+StEpEOxpjFtrMp5Q3dQldKqRCh+9CVUipEaKErpVSI0EJXSqkQoYWulFIhQgtdKaVChBa6UkqFCC10pZQKEf8flDxjWnSb1dMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Numerical approach\n",
    "N = 10000000\n",
    "A = np.random.uniform(a1, a2, N)\n",
    "bd = 1/(b2-b1)\n",
    "B = np.random.uniform(b1, b2, N)\n",
    "C = A + B\n",
    "\n",
    "plt.hist(C, bins=100, density=True, label='numerical PDF')\n",
    "\n",
    "\n",
    "# Analytical approach\n",
    "scale = 1/(a2-a1)/(b2-b1)\n",
    "\n",
    "def f_c(x, a1, a2, b1, b2):\n",
    "    if a1 + b2 < b1 + a2:\n",
    "        a1, a2, b1, b2 = b1, b2, a1, a2\n",
    "\n",
    "    y = np.zeros_like(x)\n",
    "    region1 = (b1 + a1 <= x) & (x <= b1 + a2)\n",
    "    region2 = (b1 + a2 <= x) & (x <= b2 + a1)\n",
    "    region3 = (b2 + a1 <= x) & (x <= b2 + a2)\n",
    "    y[region1] = scale * (x[region1] - b1 - a1)\n",
    "    y[region2] = scale * (a2 - a1)\n",
    "    y[region3] = scale * (a2 - (x[region3] - b2))\n",
    "    return y\n",
    "\n",
    "x = np.linspace(a1+b1, a2+b2, 1000)\n",
    "y = f_c(x, a1, a2, b1, b2)\n",
    "plt.plot(x, y, label='Analytical PDF')\n",
    "plt.xlabel('x')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Going back to the original problem of summing not just two random variables but $N$ of them, it seems reasonable to expect that increasing $N$ will increase the order of the polynomials involved, and visualizing it for $N=3$ does show some higher-order behavior:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAQcklEQVR4nO3dbYydaV3H8e+P2W1EIq7aUUgf2EZrCBowOHQlEHnQmvJkIZJYUIkPSVNCfUpQ6xsSwxs3vJAoxaYhjRqjDYk8NFAoRoOgsDpTXFZaKJlUpEM1210iZJGwdPn7Yk43Z8+emXPP9JyZOdd8P8mk577vK2euK/fm12v/53/upqqQJE2/p2z2BCRJ42GgS1IjDHRJaoSBLkmNMNAlqREGuiQ1olOgJzmU5EqSxSQnhlz/vST3934+l+SxJN8//ulKklaSUX3oSWaALwIHgSVgHnhDVV1eYfxrgN+tqpePea6SpFV02aEfABar6mpVPQqcBQ6vMv4NwN+OY3KSpO7u6DBmF3Ct73gJuGfYwCTfDRwCjo960507d9bdd9/d4ddLkm65ePHiQ1U1O+xal0DPkHMr1WleA/xLVX116BslR4GjAHv37mVhYaHDr5ck3ZLkv1a61qXksgTs6TveDVxfYewRVim3VNXpqpqrqrnZ2aF/wUiS1qlLoM8D+5PsS7KD5dA+NzgoyfcCLwE+ON4pSpK6GFlyqaqbSY4DF4AZ4ExVXUpyrHf9VG/o64CPVdU3JjZbSdKKRrYtTsrc3FxZQ5ektUlysarmhl3zm6KS1AgDXZIaYaBLUiMMdElqhIEuSY3o8k1RaSrcfeLDQ89/6Y9fNXJMV/3v1WUeXcZL42LboqbC7QbxVmPQa71Wa1t0hy5tAnfxmgRr6JLUCHfo2rJaK7OsxN26xsVA15ayXUJcmgQDXdpC3K3rdlhDl6RGuEPXprPMMpy7da2VO3RJaoSBLkmNsOSiTWGZRRo/A12aAtbT1YUlF0lqhIEuSY2w5KINY918PCy/aCXu0CWpEZ0CPcmhJFeSLCY5scKYlya5P8mlJP803mlKkkYZWXJJMgOcBA4CS8B8knNVdblvzF3Au4FDVfXlJD84qQlrulhmmSzLL+rXZYd+AFisqqtV9ShwFjg8MOaNwPuq6ssAVfXgeKcpSRqlS6DvAq71HS/1zvX7UeD7knw8ycUkbxrXBCVJ3XTpcsmQc4P/EOkdwE8CPwM8Ffh0kvuq6otPeKPkKHAUYO/evWufrSRpRV126EvAnr7j3cD1IWM+WlXfqKqHgE8Azxt8o6o6XVVzVTU3Ozu73jlLkoboskOfB/Yn2Qd8BTjCcs283weBdyW5A9gB3AP8yTgnqunhB6Gbww9INTLQq+pmkuPABWAGOFNVl5Ic610/VVWfT/JR4AHgO8B7qupzk5y4JOmJOn1TtKrOA+cHzp0aOH4H8I7xTU2StBZ+U1SSGmGgS1IjfDiX1CA/IN2eDHSNhZ0t0uaz5CJJjTDQJakRBrokNcIautbNuvl08APS7cMduiQ1wkCXpEYY6JLUCANdkhphoEtSI+xy0ZrY2TLd7Hhpmzt0SWqEgS5JjTDQJakRBrokNcJAl6RG2OUibVN2vLTHQNdItipK08GSiyQ1olOgJzmU5EqSxSQnhlx/aZKvJbm/9/O28U9VkrSakSWXJDPASeAgsATMJzlXVZcHhn6yql49gTlKkjroskM/ACxW1dWqehQ4Cxye7LQkSWvV5UPRXcC1vuMl4J4h416Y5LPAdeCtVXVpDPOTtAEGP/i262U6dQn0DDlXA8efAZ5VVY8keSXwAWD/k94oOQocBdi7d+8ap6qNZGeLNH26lFyWgD19x7tZ3oU/rqq+XlWP9F6fB+5MsnPwjarqdFXNVdXc7OzsbUxbkjSoS6DPA/uT7EuyAzgCnOsfkOQZSdJ7faD3vg+Pe7KSpJWNLLlU1c0kx4ELwAxwpqouJTnWu34KeD3w5iQ3gW8CR6pqsCwjSZqgTt8U7ZVRzg+cO9X3+l3Au8Y7NUnSWvhNUUlqhIEuSY3w4Vx6nK2K0nQz0CU9iY/WnU6WXCSpEQa6JDXCQJekRhjoktQIA12SGmGXyzZnq6JGseNlerhDl6RGGOiS1AgDXZIaYaBLUiMMdElqhF0ukjqz42VrM9C3IVsVpTZZcpGkRhjoktQIA12SGmGgS1IjDHRJakSnQE9yKMmVJItJTqwy7gVJHkvy+vFNUZLUxci2xSQzwEngILAEzCc5V1WXh4y7F7gwiYlK2lrsSd96uvShHwAWq+oqQJKzwGHg8sC43wT+DnjBWGeosbD3XGpfl5LLLuBa3/FS79zjkuwCXgecGt/UJElr0SXQM+RcDRy/E/iDqnps1TdKjiZZSLJw48aNrnOUJHXQpeSyBOzpO94NXB8YMwecTQKwE3hlkptV9YH+QVV1GjgNMDc3N/iXgiTpNnQJ9Hlgf5J9wFeAI8Ab+wdU1b5br5P8BfChwTCXJE3WyECvqptJjrPcvTIDnKmqS0mO9a5bN5e2OTtetoZOT1usqvPA+YFzQ4O8qn719qel22VXi7T9+E1RSWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1Aj/keiG2KqorcCe9M3jDl2SGmGgS1IjDHRJaoSBLkmNMNAlqRF2uUw5O1sk3WKgS5oYWxg3liUXSWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AjbFqeQveeaRrYwTp47dElqRKdAT3IoyZUki0lODLl+OMkDSe5PspDkxeOfqiRpNSNLLklmgJPAQWAJmE9yrqou9w37B+BcVVWS5wLvBZ49iQlLkobrskM/ACxW1dWqehQ4CxzuH1BVj1RV9Q6fBhSSpA3VJdB3Adf6jpd6554gyeuSfAH4MPDr45meJKmrLoGeIeeetAOvqvdX1bOB1wJvH/pGydFejX3hxo0ba5upJGlVXdoWl4A9fce7gesrDa6qTyT54SQ7q+qhgWungdMAc3NzlmXWwFZFtcQWxsnoskOfB/Yn2ZdkB3AEONc/IMmPJEnv9fOBHcDD456sJGllI3foVXUzyXHgAjADnKmqS0mO9a6fAn4BeFOSbwPfBH6x70NSSdIG6PRN0ao6D5wfOHeq7/W9wL3jnZokaS38pqgkNcJAl6RG+HAuSZvKjpfxMdC3MFsVJa2FJRdJaoSBLkmNMNAlqREGuiQ1wkCXpEbY5bLF2Nkiab0MdElbhj3pt8eSiyQ1wkCXpEYY6JLUCANdkhphoEtSI+xykbQl2fGydgb6FmDvuaRxsOQiSY0w0CWpEQa6JDXCQJekRnQK9CSHklxJspjkxJDrv5Tkgd7Pp5I8b/xTlSStZmSgJ5kBTgKvAJ4DvCHJcwaG/Sfwkqp6LvB24PS4JypJWl2XtsUDwGJVXQVIchY4DFy+NaCqPtU3/j5g9zgn2SJbFaXu7EnvpkvJZRdwre94qXduJb8BfOR2JiVJWrsuO/QMOVdDByYvYznQX7zC9aPAUYC9e/d2nKIkqYsuO/QlYE/f8W7g+uCgJM8F3gMcrqqHh71RVZ2uqrmqmpudnV3PfCVJK+gS6PPA/iT7kuwAjgDn+gck2Qu8D/iVqvri+KcpSRplZMmlqm4mOQ5cAGaAM1V1Kcmx3vVTwNuAHwDenQTgZlXNTW7akqRBqRpaDp+4ubm5WlhY2JTfvRXY5SLdvu3Y8ZLk4kobZp+2uIEMcUmT5Ff/JakRBrokNcJAl6RGGOiS1AgDXZIaYZeLpKnlQ7ueyECfMFsVJW0USy6S1AgDXZIaYaBLUiMMdElqhIEuSY2wy2UC7GyRNp4tjO7QJakZBrokNcJAl6RGGOiS1AgDXZIaYZeLpOZs144XA31MbFWUtNksuUhSIzoFepJDSa4kWUxyYsj1Zyf5dJJvJXnr+KcpSRplZMklyQxwEjgILAHzSc5V1eW+YV8Ffgt47URmKUkaqcsO/QCwWFVXq+pR4CxwuH9AVT1YVfPAtycwR0lSB10+FN0FXOs7XgLumcx0posfhEpb33bqeOmyQ8+Qc7WeX5bkaJKFJAs3btxYz1tIklbQJdCXgD19x7uB6+v5ZVV1uqrmqmpudnZ2PW8hSVpBl0CfB/Yn2ZdkB3AEODfZaUmS1mpkDb2qbiY5DlwAZoAzVXUpybHe9VNJngEsAE8HvpPkd4DnVNXXJzh3SVKfTt8UrarzwPmBc6f6Xv8Py6UYSdqyWv+A1K/+r5GdLZK2Kr/6L0mNMNAlqREGuiQ1whp6B9bNpfa0+AGpO3RJaoSBLkmNMNAlqRHW0CVte63U0w30FfhBqKRpY8lFkhphoEtSIyy59LHMImma6+nu0CWpEQa6JDXCkoskrWCwDLvVSzDbPtCtm0tqhSUXSWrEttyhuyuXtB5bvQPGHbokNcJAl6RGbJuSi2UWSeO0FcsvnXboSQ4luZJkMcmJIdeT5E971x9I8vzxT1WStJqRO/QkM8BJ4CCwBMwnOVdVl/uGvQLY3/u5B/jz3p+byl25pI2wVXbrXUouB4DFqroKkOQscBjoD/TDwF9VVQH3JbkryTOr6r/HPuMRDHFJm2kzw71LoO8CrvUdL/Hk3fewMbuAiQe6AS5pq9rocO8S6BlyrtYxhiRHgaO9w0eSXOnw+zfaTuChzZ7EBLW+Pmh/ja2vDxpcY+59wuHtrO9ZK13oEuhLwJ6+493A9XWMoapOA6c7/M5Nk2ShquY2ex6T0vr6oP01tr4+aH+Nk1pfly6XeWB/kn1JdgBHgHMDY84Bb+p1u/wU8LXNqJ9L0nY2codeVTeTHAcuADPAmaq6lORY7/op4DzwSmAR+D/g1yY3ZUnSMJ2+WFRV51kO7f5zp/peF/CW8U5t02zpktAYtL4+aH+Nra8P2l/jRNaX5SyWJE07n+UiSY3YdoGe5LuS/FuSzya5lOSPhoyZ6kcZdFzjS5N8Lcn9vZ+3bcZcb0eSmST/nuRDQ65N9T28ZcQaW7iHX0ryH735Lwy5PtX3scP6xnoPt83Dufp8C3h5VT2S5E7gn5N8pKru6xuzJR9lsAZd1gjwyap69SbMb1x+G/g88PQh16b9Ht6y2hph+u8hwMuqaqWe7Bbu42rrgzHew223Q69lj/QO7+z9DH6Q8PijDHoheFeSZ27kPG9HxzVOtSS7gVcB71lhyFTfQ+i0xu1g6u/jRtp2gQ6P/2/s/cCDwN9X1b8ODFnpUQZTo8MaAV7YK8t8JMmPbfAUb9c7gd8HvrPC9am/h4xeI0z3PYTljcbHklzsfZN80LTfx1HrgzHew20Z6FX1WFX9BMvfaD2Q5McHhnR6lMFW1mGNnwGeVVXPA/4M+MBGz3G9krwaeLCqLq42bMi5qbmHHdc4tfewz4uq6vksl1bekuSnB65P9X1k9PrGeg+3ZaDfUlX/C3wcODRwqdOjDKbBSmusqq/fKsv0vmdwZ5KdGz/DdXkR8PNJvgScBV6e5K8Hxkz7PRy5xim/hwBU1fXenw8C72f56a79pvo+jlrfuO/htgv0JLNJ7uq9firws8AXBoZN9aMMuqwxyTOSpPf6AMv/LTy80XNdj6r6w6raXVV3s/woin+sql8eGDbV97DLGqf5HgIkeVqS77n1Gvg54HMDw6b2PnZZ37jv4Xbscnkm8JdZ/oc7ngK8t6o+lLYeZdBlja8H3pzkJvBN4EhN+bfMGruHQzV2D38IeH8vz+4A/qaqPtrQfeyyvrHeQ78pKkmN2HYlF0lqlYEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1Ij/h+SMZMpZYATtwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "A = np.random.uniform(a1, a2, N)\n",
    "bd = 1/(b2-b1)\n",
    "B = np.random.uniform(b1, b2, N)\n",
    "C = np.random.uniform(1.6, 2.0, N)\n",
    "D = A + B + C\n",
    "\n",
    "plt.hist(D, bins=100, density=True, label='numerical PDF');"
   ]
  }
 ],
 "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}