{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Adaptive rejection sampling\n",
"\n",
"Rejection sampling (RS) is a useful method for sampling intractable distributions. It defines an envelope function which upper-bounds the target unnormalised probability density to be sampled. It then proceeds to sample points in the area under the envelope, rejecting those points which fall above the target and accepting the rest. The accepted points are independent and identically distributed samples from the target distribution. There are two important issues with RS. The first is that if the envelope is a very loose upper bound, then most samples will be rejected and the scheme will be slow. The second is that for rejection sampling to work, we must be certain that the envelope is an upper bound to the target, which in practice may be a challenging task.\n",
"\n",
"Adaptive rejection sampling (ARS) {cite}`gilks1992ars` is an efficient method for sampling log-concave targets, which deals with both of these issues. It is origially defined for univariate distributions, but can also be extended to multivariate distributions via Gibbs sampling.{cite}`bishop2006PRML` ARS maintains an envelope which adapts as more points are sampled, becoming a progressively tighter bound to the target, thereby avoiding the inefficiency of regular RS. Further, the way that ARS constructs this envelope guarantees that the envelope is in fact an upper bound to the target, which sidesteps the second difficulty described above."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"tags": [
"remove-cell"
]
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/var/folders/w_/5zj48w1d0xb7ycgdm6pk40v00000gn/T/ipykernel_98848/165326261.py:5: DeprecationWarning: `set_matplotlib_formats` is deprecated since IPython 7.23, directly use `matplotlib_inline.backend_inline.set_matplotlib_formats()`\n",
" set_matplotlib_formats('pdf', 'svg')\n"
]
}
],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"from IPython.display import HTML, set_matplotlib_formats\n",
"set_matplotlib_formats('pdf', 'svg')\n",
"#css_style = open('../../../_static/custom_style.css', 'r').read()\n",
"#HTML(f'')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## An adaptive envelope function\n",
"\n",
"Suppose we wish to sample from a log-concave univariate distribution with unnormalised distribution function $f$. Whereas RS defines a fixed envelope, ARS will define an envelope $g_u$ that upper bounds $f$ and adapts its shape as the sampling procedure progresses. By adapting its shape, using the infromation that $f$ is log-concave, the envelope reduces the probability of future rejections. In addition the envelope function, ARS can use an optional function $g_l$ which lower bounds $f$, called the squeezing function. The squeezing function can be used to avoid evaluating $f$ in the rejection step, which can be especially useful if $f$ is computationally expensive to evaluate.\n",
"\n",
"Given the an ordered set of points $x_1 < x_2 < ... < x_K$, ARS defines the log-envelope $\\log ~ g_u$ to be the minimum over the tangents to $h = \\log f$ at these points. The log-squeezing function is defined to be the piecewise linear function which joins the points $(x_k, h(x_k))$ inside the interval $[x_1, x_K]$ and is equal to $-\\infty$ outside this innterval. Examples of envelope and squeezing functions are shown below.\n",
"\n",
"
\n",
" \n",
"**Definition (Abscissa set, envelope function and squeezing function)** Let $f(x)$ be a univariate log-concave function, with non-zero domain $D = \\{x : f(x) > 0\\}$. An abscissa set $T_K$ is an ordered set of points in $D$ such that\n",
" \n",
"$$T_k = \\{x_1 < x_2 < ... < x_K\\}.$$\n",
" \n",
"The envelope function $g_u(x)$ defined by $T_k$ is\n",
" \n",
"$$g_u(x) = \\min_{k} g_{u, k}(x)$$\n",
" \n",
"where $g_{u, k}(x), k = 1, 2, ..., K$ are piecewise exponential functions such that\n",
" \n",
"$$g_{u, k}(x_k) = f(x_k) \\text{ and } g_{u, k}'(x_k) = \\log f'(x_k).$$\n",
" \n",
"The squeezing function $g_l(x)$ defined by $T_k$ is\n",
" \n",
"$$g_l(x) = \\begin{cases} \\min_{k} g_{l, k}(x) & \\text{ if } x_1 \\leq x \\leq x_k, \\\\ 0 & \\text{ otherwise.} \\end{cases}$$\n",
" \n",
"where $g_{u, k}(x), k = 1, 2, ..., K - 1$ are piecewise exponential functions such that\n",
" \n",
"$$g_{l, k}(x_k) = f(x_k) \\text{ and } g_{l, k}(x_{k+1}) = \\log f(x_{k+1}).$$\n",
" \n",
"
\n",
" \n",
"\n",
"Below are functions implementing the necessary calculations to determine the envelope and squeezing functions, all of which are cheap operations. The functions take in points at input locations $x$ and corresponding $h$ and $h'$ values and carry out computations in log-space, before exponentiating the result at the end."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"def g_u(x, xs, hs, dhdxs):\n",
" \n",
" z, _ = compute_points_of_intersection_and_intercepts(xs, hs, dhdxs)\n",
" i = np.searchsorted(z, x)\n",
" \n",
" return np.exp(dhdxs[i] * (x - xs[i]) + hs[i])\n",
" \n",
"\n",
"def g_l(x, xs, hs):\n",
" \n",
" if all(x < xs) or all(x > xs):\n",
" return 0.\n",
" \n",
" else:\n",
" i = np.searchsorted(xs, x)\n",
" m = (hs[i] - hs[i-1]) / (xs[i] - xs[i-1])\n",
" \n",
" return np.exp(hs[i-1] + (x - xs[i-1]) * m)\n",
"\n",
"\n",
"def compute_points_of_intersection_and_intercepts(x, h, dhdx):\n",
" \n",
" # y-intercepts c of envelope function line segments, intersection points z\n",
" c = h - dhdx * x\n",
" z = (c[1:] - c[:-1]) / (dhdx[:-1] - dhdx[1:])\n",
" \n",
" return z, c"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now let's define a log unnormalised Gaussian log density, and use this to illustrate the envelope and squeezing function defined by an abcissa set with three points."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"def log_gaussian(mean, variance):\n",
" return lambda x : (- 0.5 * (x - mean) ** 2 / variance, - (x - mean) / variance)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"tags": [
"hide-input",
"center-output"
]
},
"outputs": [
{
"data": {
"application/pdf": "JVBERi0xLjQKJazcIKu6CjEgMCBvYmoKPDwgL1R5cGUgL0NhdGFsb2cgL1BhZ2VzIDIgMCBSID4+CmVuZG9iago4IDAgb2JqCjw8IC9Gb250IDMgMCBSIC9YT2JqZWN0IDcgMCBSIC9FeHRHU3RhdGUgNCAwIFIgL1BhdHRlcm4gNSAwIFIKL1NoYWRpbmcgNiAwIFIgL1Byb2NTZXQgWyAvUERGIC9UZXh0IC9JbWFnZUIgL0ltYWdlQyAvSW1hZ2VJIF0gPj4KZW5kb2JqCjExIDAgb2JqCjw8IC9UeXBlIC9QYWdlIC9QYXJlbnQgMiAwIFIgL1Jlc291cmNlcyA4IDAgUgovTWVkaWFCb3ggWyAwIDAgMzcxLjIgMjAyLjE0MTg3NSBdIC9Db250ZW50cyA5IDAgUiAvQW5ub3RzIDEwIDAgUiA+PgplbmRvYmoKOSAwIG9iago8PCAvTGVuZ3RoIDEyIDAgUiAvRmlsdGVyIC9GbGF0ZURlY29kZSA+PgpzdHJlYW0KeJzVl09vHEUQxe/zKfqYHNzuruq/Bw4JAUtIHBIsceFmEgfLJoIIwsfnVzW73tldSwlcEJbW2n7TU1316lV17eWrt3/+cvP2zdXL8PUPy+VhdfNxyeGOz21I4Y7Pp5DDFZ/bJbF6WLTnKHy7332TJDGXPHoFSsfL98vybrl8wcsfeeNqkWkvjNjEN2Cslc363td5ljj3FvyNI2RnU1abt3h2x+eT2bf18hu+p3CRAPPMkf8pzhpuHpaX1xZCzLVJDdc/L5ff8nCE63fLs7+eh+u75Zvr5TXv26sXZkQsvt5jH372wYTKyHsbsrdx7zZqTCnN3vyQZFueffAHOceUpc2q+fHJrT8RfJoyek5Tx/pk49o739Oiapu5pEezm5N/euZ7ekwipWGl1zMza4Rw0apo2ew4WHm+ISGcpEq1xBFyww8Jv78NP4Zfg4TvQo4V7o/1syatpyhpn+Yap8xWRugt9opwSGQRfDE4DI1JW6+OzqhjzFLDmDHVNOcErS3mPOrUgBQ6lmoBbRJHGVIsx5zWNedm8Iy1d4wA11hGmsN2c7QmDicMaCiSUjXTfnjPowGTvJm0mB+c3lsbih7zjGXmoaCz8ZU4IE0UbkoVOzEnjSItK55I58hRZDg+0R7O4IqS/6F2PHhusY1apIOXWPDQAzKjJVe8BMdO69rdviKRUrrgTSHS2nJa8RZHJ5v4A529t1TE8IJzqSSDiUkIz0oHy8hICQsctSTEuOIsmmrGnTJiF03TzTSFbDKEO7yaWx/VjyW31eMFzxE5jZUGakX8AbjwlbQYa3AX56AaHSYxOrrBM8UO6RV3KgTWMbuzMwm9ZCFEw7XOvvaCVCJHG1WGY8+IBbcKJ5tJ3X6vZa5wj00THjqcyizJzFMjbpNN5j1kJzUvxTM0oRkWZhyKJorjPeI8ZDk7mudMxo6UEmuBfktKi1UpTwuLsNHi0Gnsk9wx6+pPHVA+UvFc5U4bcHcaiibNAxYU5Q4id3foriYANY1QNbSc6Sz0QegNxYBbmrV0d3PQJKWRxmByESqneQPN0ZrRMGniZaE1ufeTULyislCIuVUvP6ViLCXdqmQgxtaqBUVpUWuWW3CyWfpc4UGtGQnBtE4mioWqUqk1tmOdpJE48cpUJQ/UF/XP8Sgre9NXNNeU+MJUc0rcBsotU6EgjM4h0tRNmLDEK31kq8vsm+F2VukkrBNASXN/mxza0Pl6+WF5Hf5hq7N7IYU3V0/2PAQ9567nZVzqBb9NZtBR2lqDNC36T63qsCAz8XxQa6CE4LUmbZRyKg/rja2O1PexYHZ6tz9f/5vY7GqrCXGlPop6mCejwWLJT9VqwK7IB18rlWq1R3JRyzytVImTriteGqhr8L21tRMkmR6kctNONVGjF1pj63QmxyGr4VCx087X+yCJwIcVYju9jM5GjvOJAivHc8jD+RzCns8OL/s9mxeftvRlidlNMj7IZEh0gmvfjSLh8vsUXn143DKprM70RQOR/MTzTDfkor+wTtXoDIcdr92fyxdqE9Xj3Ic369RHYyEkSvQQUu1bYK0+9b8tTo/1jtZMB7Iuubx2T28WByjtjRkD4KVXoVNvrYI+Wu3b5c3q3gGwRsiQY39HG0G3xtdl35pJW/fuV+Dg/tbqIcYTbm5sOn25n069sGw63RXXac+gY+NDtrtuL80jwLrxCYB4kov9djfzHs27tJJgd/Dem8d5dR0b9+Ne2oyrNjQMCngzlH7YDZz0/9w2A+duWuUiQLLc8ShxO62mzbRK527jaMRcn361H0JpB9yl5XhQXfe8P0yhn+26zmCrJwwegB2DW+CMwbNfDc5iy/8li7f7mb/sfkhccIlvjs+h+74/nuLqiS7+FG1cJse0HYAdbVvgC2njzv4f0HZ/oO318jcioQ8mCmVuZHN0cmVhbQplbmRvYmoKMTIgMCBvYmoKMTQzMQplbmRvYmoKMTAgMCBvYmoKWyBdCmVuZG9iagoxOCAwIG9iago8PCAvTGVuZ3RoIDE0NyAvRmlsdGVyIC9GbGF0ZURlY29kZSA+PgpzdHJlYW0KeJw9j7sRwzAMQ3tOgRHMvzRPcrkUzv5tYEVOI74jRAL0NByowSd6oOvAU+XGj0QF1cQpzt6myEVuimLXj0ZV4iE2DcWNNqjYhHUg51K4J6M5HWqbfE6aXmSpmzSvKKfUr2jdQueGMfa49r1SaZcrmk5HBY2VEWlsUWjGZAAON3+71TrN+fOqj/+xp7zl9QW05DRICmVuZHN0cmVhbQplbmRvYmoKMTkgMCBvYmoKPDwgL0xlbmd0aCA0MDQgL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicPZJLjuQwDEP3OYUu0ID1s53zVKMxi6r7b+fRmelFIMWWKIr0zLZhGfYVY1nf03oO+/brpLXsc/VIm9veV+2mrKx62NcdVoPoue11Zd2kuyy2P1BO/0lel8/6d5YnKXN+SAaXWxh3W7cA9qHiS5iCvJ25C0gfbTfd7m3uGhYh6Eqlt/k6XJK5+0leVxxM0hwC9WWZdO4pwos77vMWllvpE9kCfGvbvBM9LCGf6nQ7EwP8ZTED2QIxxllwEz2gtk1tzeF0g31N8wqFGJqrpKcqkn9IZBfHrFaQ9KCjgKGktnaQFI1awf4dNKds0dCY+NVpwQbyLFlUkd3mc1LDqQARfTo1pjQaHWqBD9G6qZpaq8dj9/vX+Pf1B7ZuCeXPk6ka74tFhVk4BObY4CHy4oQXkuUnwkOTyELb4XNQVZgpZ6QMisK5alugW6VuXKx4cL4kO6KezR6r50FyPJSzJ7rmKHN2dvg471lvwf1Bc9uIh6xIHfMEEdtH+ok5QzqB2ppRWKX3ebhh6v/9X2jx8xcGNZaOCmVuZHN0cmVhbQplbmRvYmoKMjAgMCBvYmoKPDwgL0xlbmd0aCAyMTYgL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicPVAxkgMxCOv9Cp4ACLD9nr25uSL3/zbCzqTYFYOQjCiEqCD5yyqZpfJjI3MKoPI/YoPUa4R+sMj4PNUWhEtMoz6JEJB9RswUzM1OCda6uPMwJeHUcC4CEqn8djNop5BQzi+n/uIzeo+uvJ0qxZM6gIwbXVTFdnCDKbaI2ox5nm1xwRZOstfw+a2MK5d1BdpsE2f+ZBfqkucBkMlkPC6WVIVd5Lrlt7O3ZHDKqXUjk8kgDNvXjFwX7TDVZ+A8Q2BRU3Gcm2Hg0uOIoNPn/M/4G79vORVOrgplbmRzdHJlYW0KZW5kb2JqCjIxIDAgb2JqCjw8IC9MZW5ndGggNTQgL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicMzK3UDBQMDYFEkbmxgrmZgYKKYZchhYQZi4XTDCHy9AIpCiHyxhCwZTkcGVwpQEA0K4NdAplbmRzdHJlYW0KZW5kb2JqCjIyIDAgb2JqCjw8IC9MZW5ndGggMjEzIC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nDVQOXLEMAzr9Qp8YGcEHpL4HmcyKTb/bwM63sIGRIkgwOWBidx40ZGnkLHxxbE2jBu/g8679h5mH0a9N5s3c9gMMAnWFOrbhWswDiiVG/M8yM+NuU4LnOrbjrP6QgKbmhtYAYtsuIaVN3E5PRO+NUFerxEzpZoIk56ch5eauiOqHqefRO8Rx9GOvf4xyJ6qwkk9gWfCCq5YtJ6q8qtdnHs1NutGuVy8WUeYKG0iVZUrhdNezu4/vYtNZKhBWbMUtRcgAUp9tSrxbPoaP+P7D9CGSicKZW5kc3RyZWFtCmVuZG9iagoyMyAwIG9iago8PCAvTGVuZ3RoIDk1IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nD2MQQ7AIAgE77xiP9AEERX/0zQ92P9fu0bbC0x2YUo3KA4rnFUVxRvOJB8+kr3DWseQoplHQ5zd3BYOS40Uq1gWFp5hEaS0Ncz4vChrYEop6mln9b+75XoB/58cLAplbmRzdHJlYW0KZW5kb2JqCjE2IDAgb2JqCjw8IC9UeXBlIC9Gb250IC9CYXNlRm9udCAvR0NXWERWK0RlamFWdVNhbnMtT2JsaXF1ZSAvRmlyc3RDaGFyIDAKL0xhc3RDaGFyIDI1NSAvRm9udERlc2NyaXB0b3IgMTUgMCBSIC9TdWJ0eXBlIC9UeXBlMwovTmFtZSAvR0NXWERWK0RlamFWdVNhbnMtT2JsaXF1ZSAvRm9udEJCb3ggWyAtMTAxNiAtMzUxIDE2NjAgMTA2OCBdCi9Gb250TWF0cml4IFsgMC4wMDEgMCAwIDAuMDAxIDAgMCBdIC9DaGFyUHJvY3MgMTcgMCBSCi9FbmNvZGluZyA8PCAvVHlwZSAvRW5jb2RpbmcgL0RpZmZlcmVuY2VzIFsgMTAyIC9mIC9nIC9oIDEwOCAvbCAxMTcgL3UgMTIwIC94IF0gPj4KL1dpZHRocyAxNCAwIFIgPj4KZW5kb2JqCjE1IDAgb2JqCjw8IC9UeXBlIC9Gb250RGVzY3JpcHRvciAvRm9udE5hbWUgL0dDV1hEVitEZWphVnVTYW5zLU9ibGlxdWUgL0ZsYWdzIDk2Ci9Gb250QkJveCBbIC0xMDE2IC0zNTEgMTY2MCAxMDY4IF0gL0FzY2VudCA5MjkgL0Rlc2NlbnQgLTIzNiAvQ2FwSGVpZ2h0IDAKL1hIZWlnaHQgMCAvSXRhbGljQW5nbGUgMCAvU3RlbVYgMCAvTWF4V2lkdGggMTM1MCA+PgplbmRvYmoKMTQgMCBvYmoKWyA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMAo2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDMxOCA0MDEgNDYwIDgzOCA2MzYKOTUwIDc4MCAyNzUgMzkwIDM5MCA1MDAgODM4IDMxOCAzNjEgMzE4IDMzNyA2MzYgNjM2IDYzNiA2MzYgNjM2IDYzNiA2MzYgNjM2CjYzNiA2MzYgMzM3IDMzNyA4MzggODM4IDgzOCA1MzEgMTAwMCA2ODQgNjg2IDY5OCA3NzAgNjMyIDU3NSA3NzUgNzUyIDI5NQoyOTUgNjU2IDU1NyA4NjMgNzQ4IDc4NyA2MDMgNzg3IDY5NSA2MzUgNjExIDczMiA2ODQgOTg5IDY4NSA2MTEgNjg1IDM5MCAzMzcKMzkwIDgzOCA1MDAgNTAwIDYxMyA2MzUgNTUwIDYzNSA2MTUgMzUyIDYzNSA2MzQgMjc4IDI3OCA1NzkgMjc4IDk3NCA2MzQgNjEyCjYzNSA2MzUgNDExIDUyMSAzOTIgNjM0IDU5MiA4MTggNTkyIDU5MiA1MjUgNjM2IDMzNyA2MzYgODM4IDYwMCA2MzYgNjAwIDMxOAozNTIgNTE4IDEwMDAgNTAwIDUwMCA1MDAgMTM1MCA2MzUgNDAwIDEwNzAgNjAwIDY4NSA2MDAgNjAwIDMxOCAzMTggNTE4IDUxOAo1OTAgNTAwIDEwMDAgNTAwIDEwMDAgNTIxIDQwMCAxMDI4IDYwMCA1MjUgNjExIDMxOCA0MDEgNjM2IDYzNiA2MzYgNjM2IDMzNwo1MDAgNTAwIDEwMDAgNDcxIDYxNyA4MzggMzYxIDEwMDAgNTAwIDUwMCA4MzggNDAxIDQwMSA1MDAgNjM2IDYzNiAzMTggNTAwCjQwMSA0NzEgNjE3IDk2OSA5NjkgOTY5IDUzMSA2ODQgNjg0IDY4NCA2ODQgNjg0IDY4NCA5NzQgNjk4IDYzMiA2MzIgNjMyIDYzMgoyOTUgMjk1IDI5NSAyOTUgNzc1IDc0OCA3ODcgNzg3IDc4NyA3ODcgNzg3IDgzOCA3ODcgNzMyIDczMiA3MzIgNzMyIDYxMSA2MDgKNjMwIDYxMyA2MTMgNjEzIDYxMyA2MTMgNjEzIDk5NSA1NTAgNjE1IDYxNSA2MTUgNjE1IDI3OCAyNzggMjc4IDI3OCA2MTIgNjM0CjYxMiA2MTIgNjEyIDYxMiA2MTIgODM4IDYxMiA2MzQgNjM0IDYzNCA2MzQgNTkyIDYzNSA1OTIgXQplbmRvYmoKMTcgMCBvYmoKPDwgL2YgMTggMCBSIC9nIDE5IDAgUiAvaCAyMCAwIFIgL2wgMjEgMCBSIC91IDIyIDAgUiAveCAyMyAwIFIgPj4KZW5kb2JqCjI4IDAgb2JqCjw8IC9MZW5ndGggNzQgL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicszC2UDBQMDQwUzA0N1IwNzZSMDE1UUgx5AIJgZi5XDDBHDDLGKgsByyLYEFkM8BsI1NTqB4QC6LHEK4SwYLIZnClAQBRvhkWCmVuZHN0cmVhbQplbmRvYmoKMjkgMCBvYmoKPDwgL0xlbmd0aCAzNDEgL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicRVJLbkQxCNu/U3CBSOGXkPO0qrqY3n9bm0zVzeAJYGx4y1OmZMqwuSUjJNeUT30iQ6ym/DRyJCKm+EkJBXaVj8drS6yN7JGoFJ/a8eOx9Eam2RVa9e7Rpc2iUc3KyDnIEKGeFbqye9QO2fB6XEi675TNIRzL/1CBLGXdcgolQVvQd+wR3w8droIrgmGway6D7WUy1P/6hxZc7333YscugBas577BDgCopxO0BcgZ2u42KWgAVbqLScKj8npudqJso1Xp+RwAMw4wcsCIJVsdvtHeAJZ9XehFjYr9K0BRWUD8yNV2wd4xyUhwFuYGjr1wPMWZcEs4xgJAir3iGHrwJdjmL1euiJrwCXW6ZC+8wp7a5udCkwh3rQAOXmTDraujqJbt6TyC9mdFckaM1Is4OiGSWtI5guLSoB5a41w3seJtI7G5V9/uH+GcL1z26xdL7ITECmVuZHN0cmVhbQplbmRvYmoKMzAgMCBvYmoKPDwgL0xlbmd0aCA0NyAvRmlsdGVyIC9GbGF0ZURlY29kZSA+PgpzdHJlYW0KeJwzMrdQMFCwNAEShhYmCuZmBgophlyWEFYuF0wsB8wC0ZZwCiKewZUGALlnDScKZW5kc3RyZWFtCmVuZG9iagozMSAwIG9iago8PCAvTGVuZ3RoIDIxOCAvRmlsdGVyIC9GbGF0ZURlY29kZSA+PgpzdHJlYW0KeJw9ULmNBDEMy12FGljAeu2pZxaLS6b/9Ej59iLRFkVSKjWZkikvdZQlWVPeOnyWxA55huVuZDYlKkUvk7Al99AK8X2J5hT33dWWs0M0l2g5fgszKqobHdNLNppwKhO6oNzDM/oNbXQDVocesVsg0KRg17YgcscPGAzBmROLIgxKTQb/rnKPn16LGz7D8UMUkZIO5jX/WP3ycw2vU48nkW5vvuJenKkOAxEckpq8I11YsS4SEWk1QU3PwFotgLu3Xv4btCO6DED2icRxmlKOob9rcKXPL+UnU9gKZW5kc3RyZWFtCmVuZG9iagozMiAwIG9iago8PCAvTGVuZ3RoIDE1MCAvRmlsdGVyIC9GbGF0ZURlY29kZSA+PgpzdHJlYW0KeJw9TzkOwzAM2/0KfiCAdVi23pMi6JD+f63ooB0EEaB4yLKjYwUOMYFJxxyJl7Qf/DSNQCyDmiN6QsUwLHA2SYGHQVZJVz5bnEwhtQVeSPjWFDwbTWSCnseIHbiTyegD71JbsXXoAe0QVSRdswxjsa26cD1hBDXFehXm9TBjiZJHn1VL6wEFE/jS+X/ubu92fQFgxTBdCmVuZHN0cmVhbQplbmRvYmoKMzMgMCBvYmoKPDwgL0xlbmd0aCAxNTEgL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicNY/LDcMwDEPvmoILBNDPsjxPiqCHdP9rJacFDJgwySfZFoORjENMYOyYY+ElVE+tPiQjt7pJORCpUDcET2hMDDOcpEvglem+ZTy3eDmt1AWdkMjdWW00RBnNPIajp+wVTvovc5OolRllDsisU91OyMqCFZgX1HLfz7itcqETHrYrw6I7xYhymxlp+P3vpDddX9x4MNUKZW5kc3RyZWFtCmVuZG9iagoyNiAwIG9iago8PCAvVHlwZSAvRm9udCAvQmFzZUZvbnQgL0JNUVFEVitEZWphVnVTYW5zIC9GaXJzdENoYXIgMCAvTGFzdENoYXIgMjU1Ci9Gb250RGVzY3JpcHRvciAyNSAwIFIgL1N1YnR5cGUgL1R5cGUzIC9OYW1lIC9CTVFRRFYrRGVqYVZ1U2FucwovRm9udEJCb3ggWyAtMTAyMSAtNDYzIDE3OTQgMTIzMyBdIC9Gb250TWF0cml4IFsgMC4wMDEgMCAwIDAuMDAxIDAgMCBdCi9DaGFyUHJvY3MgMjcgMCBSCi9FbmNvZGluZyA8PCAvVHlwZSAvRW5jb2RpbmcKL0RpZmZlcmVuY2VzIFsgNDAgL3BhcmVubGVmdCAvcGFyZW5yaWdodCA2MSAvZXF1YWwgMTAzIC9nIDEwOCAvbCAxMTEgL28gXQo+PgovV2lkdGhzIDI0IDAgUiA+PgplbmRvYmoKMjUgMCBvYmoKPDwgL1R5cGUgL0ZvbnREZXNjcmlwdG9yIC9Gb250TmFtZSAvQk1RUURWK0RlamFWdVNhbnMgL0ZsYWdzIDMyCi9Gb250QkJveCBbIC0xMDIxIC00NjMgMTc5NCAxMjMzIF0gL0FzY2VudCA5MjkgL0Rlc2NlbnQgLTIzNiAvQ2FwSGVpZ2h0IDAKL1hIZWlnaHQgMCAvSXRhbGljQW5nbGUgMCAvU3RlbVYgMCAvTWF4V2lkdGggMTM0MiA+PgplbmRvYmoKMjQgMCBvYmoKWyA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMAo2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDMxOCA0MDEgNDYwIDgzOCA2MzYKOTUwIDc4MCAyNzUgMzkwIDM5MCA1MDAgODM4IDMxOCAzNjEgMzE4IDMzNyA2MzYgNjM2IDYzNiA2MzYgNjM2IDYzNiA2MzYgNjM2CjYzNiA2MzYgMzM3IDMzNyA4MzggODM4IDgzOCA1MzEgMTAwMCA2ODQgNjg2IDY5OCA3NzAgNjMyIDU3NSA3NzUgNzUyIDI5NQoyOTUgNjU2IDU1NyA4NjMgNzQ4IDc4NyA2MDMgNzg3IDY5NSA2MzUgNjExIDczMiA2ODQgOTg5IDY4NSA2MTEgNjg1IDM5MCAzMzcKMzkwIDgzOCA1MDAgNTAwIDYxMyA2MzUgNTUwIDYzNSA2MTUgMzUyIDYzNSA2MzQgMjc4IDI3OCA1NzkgMjc4IDk3NCA2MzQgNjEyCjYzNSA2MzUgNDExIDUyMSAzOTIgNjM0IDU5MiA4MTggNTkyIDU5MiA1MjUgNjM2IDMzNyA2MzYgODM4IDYwMCA2MzYgNjAwIDMxOAozNTIgNTE4IDEwMDAgNTAwIDUwMCA1MDAgMTM0MiA2MzUgNDAwIDEwNzAgNjAwIDY4NSA2MDAgNjAwIDMxOCAzMTggNTE4IDUxOAo1OTAgNTAwIDEwMDAgNTAwIDEwMDAgNTIxIDQwMCAxMDIzIDYwMCA1MjUgNjExIDMxOCA0MDEgNjM2IDYzNiA2MzYgNjM2IDMzNwo1MDAgNTAwIDEwMDAgNDcxIDYxMiA4MzggMzYxIDEwMDAgNTAwIDUwMCA4MzggNDAxIDQwMSA1MDAgNjM2IDYzNiAzMTggNTAwCjQwMSA0NzEgNjEyIDk2OSA5NjkgOTY5IDUzMSA2ODQgNjg0IDY4NCA2ODQgNjg0IDY4NCA5NzQgNjk4IDYzMiA2MzIgNjMyIDYzMgoyOTUgMjk1IDI5NSAyOTUgNzc1IDc0OCA3ODcgNzg3IDc4NyA3ODcgNzg3IDgzOCA3ODcgNzMyIDczMiA3MzIgNzMyIDYxMSA2MDUKNjMwIDYxMyA2MTMgNjEzIDYxMyA2MTMgNjEzIDk4MiA1NTAgNjE1IDYxNSA2MTUgNjE1IDI3OCAyNzggMjc4IDI3OCA2MTIgNjM0CjYxMiA2MTIgNjEyIDYxMiA2MTIgODM4IDYxMiA2MzQgNjM0IDYzNCA2MzQgNTkyIDYzNSA1OTIgXQplbmRvYmoKMjcgMCBvYmoKPDwgL2VxdWFsIDI4IDAgUiAvZyAyOSAwIFIgL2wgMzAgMCBSIC9vIDMxIDAgUiAvcGFyZW5sZWZ0IDMyIDAgUgovcGFyZW5yaWdodCAzMyAwIFIgPj4KZW5kb2JqCjMgMCBvYmoKPDwgL0YxIDE2IDAgUiAvRjIgMjYgMCBSID4+CmVuZG9iago0IDAgb2JqCjw8IC9BMSA8PCAvVHlwZSAvRXh0R1N0YXRlIC9DQSAwIC9jYSAxID4+Ci9BMiA8PCAvVHlwZSAvRXh0R1N0YXRlIC9DQSAxIC9jYSAxID4+Ci9BMyA8PCAvVHlwZSAvRXh0R1N0YXRlIC9DQSAwLjggL2NhIDAuOCA+PiA+PgplbmRvYmoKNSAwIG9iago8PCA+PgplbmRvYmoKNiAwIG9iago8PCA+PgplbmRvYmoKNyAwIG9iago8PCAvTTAgMTMgMCBSID4+CmVuZG9iagoxMyAwIG9iago8PCAvVHlwZSAvWE9iamVjdCAvU3VidHlwZSAvRm9ybSAvQkJveCBbIC04IC04IDggOCBdIC9MZW5ndGggMTMxCi9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nG2QQQ6EIAxF9z1FL/BJS0Vl69JruJlM4v23A3FATN000L48flH+kvBOpcD4JAlLTrPketOQ0rpMjBjm1bIox6BRLdbOdTioz9BwY3SLsRSm1NboeKOb6Tbekz/6sFkhRj8cDq+EexZDJlwpMQaH3wsv28P/EZ5e1MAfoo1+Y1pD/QplbmRzdHJlYW0KZW5kb2JqCjIgMCBvYmoKPDwgL1R5cGUgL1BhZ2VzIC9LaWRzIFsgMTEgMCBSIF0gL0NvdW50IDEgPj4KZW5kb2JqCjM0IDAgb2JqCjw8IC9DcmVhdG9yIChNYXRwbG90bGliIHYzLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZykKL1Byb2R1Y2VyIChNYXRwbG90bGliIHBkZiBiYWNrZW5kIHYzLjkuMSkKL0NyZWF0aW9uRGF0ZSAoRDoyMDI0MDczMTIwNDEwMy0wNCcwMCcpID4+CmVuZG9iagp4cmVmCjAgMzUKMDAwMDAwMDAwMCA2NTUzNSBmIAowMDAwMDAwMDE2IDAwMDAwIG4gCjAwMDAwMDg4MTUgMDAwMDAgbiAKMDAwMDAwODMwMiAwMDAwMCBuIAowMDAwMDA4MzQ1IDAwMDAwIG4gCjAwMDAwMDg0ODcgMDAwMDAgbiAKMDAwMDAwODUwOCAwMDAwMCBuIAowMDAwMDA4NTI5IDAwMDAwIG4gCjAwMDAwMDAwNjUgMDAwMDAgbiAKMDAwMDAwMDMzOSAwMDAwMCBuIAowMDAwMDAxODY2IDAwMDAwIG4gCjAwMDAwMDAyMDggMDAwMDAgbiAKMDAwMDAwMTg0NSAwMDAwMCBuIAowMDAwMDA4NTYxIDAwMDAwIG4gCjAwMDAwMDQwMjMgMDAwMDAgbiAKMDAwMDAwMzgwOCAwMDAwMCBuIAowMDAwMDAzNDUxIDAwMDAwIG4gCjAwMDAwMDUwNzYgMDAwMDAgbiAKMDAwMDAwMTg4NiAwMDAwMCBuIAowMDAwMDAyMTA2IDAwMDAwIG4gCjAwMDAwMDI1ODMgMDAwMDAgbiAKMDAwMDAwMjg3MiAwMDAwMCBuIAowMDAwMDAyOTk4IDAwMDAwIG4gCjAwMDAwMDMyODQgMDAwMDAgbiAKMDAwMDAwNzE0NiAwMDAwMCBuIAowMDAwMDA2OTM5IDAwMDAwIG4gCjAwMDAwMDY1NzUgMDAwMDAgbiAKMDAwMDAwODE5OSAwMDAwMCBuIAowMDAwMDA1MTU4IDAwMDAwIG4gCjAwMDAwMDUzMDQgMDAwMDAgbiAKMDAwMDAwNTcxOCAwMDAwMCBuIAowMDAwMDA1ODM3IDAwMDAwIG4gCjAwMDAwMDYxMjggMDAwMDAgbiAKMDAwMDAwNjM1MSAwMDAwMCBuIAowMDAwMDA4ODc1IDAwMDAwIG4gCnRyYWlsZXIKPDwgL1NpemUgMzUgL1Jvb3QgMSAwIFIgL0luZm8gMzQgMCBSID4+CnN0YXJ0eHJlZgo5MDMyCiUlRU9GCg==",
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# The log unnormalised density to illustrate\n",
"log_prob = log_gaussian(0., 1.)\n",
"\n",
"# Points in the abcissa set and corresponding log-probabilities and gradients\n",
"xs = np.array([-1., 0.1, 1.5])\n",
"hs, dhdxs = log_prob(xs)\n",
"\n",
"# Locations to plot the log unnorm. density and envelope/squeezing functions\n",
"x_plot = np.linspace(-2, 2, 200)\n",
"log_probs = [log_prob(x)[0] for x in x_plot]\n",
"gu = [g_u(x, xs, hs, dhdxs) for x in x_plot]\n",
"gl = [g_l(x, xs, hs) for x in x_plot]\n",
"\n",
"# Plot the log unnormalised density, the envelope and squeezing functions\n",
"plt.figure(figsize=(6, 3))\n",
"plt.scatter(xs, hs, color='k', zorder=3)\n",
"plt.plot(x_plot, log_probs, color='black', label='$\\log~f = h$')\n",
"plt.plot(x_plot, np.log(gu), color='red', label='$\\log~g_u$')\n",
"\n",
"# Handle the case of negatively infinite gl, for plotting presentation\n",
"floored_log_gl = np.log(np.maximum(np.array(gl), np.ones_like(gl) * 1e-9))\n",
"plt.plot(x_plot, floored_log_gl, color='green', label='$\\log~g_l$')\n",
"\n",
"# Plot formatting\n",
"plt.xlim([-2, 2])\n",
"plt.ylim([-3, 1])\n",
"plt.xticks([])\n",
"plt.yticks([])\n",
"plt.xlabel('$x$', fontsize=18)\n",
"plt.ylabel('$\\log~f(x)$', fontsize=18)\n",
"plt.legend()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Adaptive rejection sampling\n",
"\n",
"As observed above, by vitrue of the log-concavity of $f$, the envelope and squeezing functions defined in this way are upper and lower bounds to $f$. If we then sample a point at random from the area under $g_u$, and this point also happens to be in the area under $f$, then the point is uniformly distributed in the area under $f$, and is an exact sample from the target distribution. Further, if the point happened to lie in the area under the squeezing function $g_l$, it is certain to also lie in the area under $f$, so we need not check this latter condition explicitly. This shortcut is particularly useful if the function $f$ is expensive to evaluate, because it avoids some of the evaluations of $f$. Combining these checks, we arrive at the ARS algorithm below.\n",
"\n",
"
\n",
" \n",
"**Algorithm (Adaptive Rejection Sampling)** Given a univariate un-normalised probability density $f(x)$, perform the following initialisation, sampling and update steps:\n",
" \n",
"1. Initialise an abscissa set $T_k$, such that $f'(x_1) > 0$ and $f'(x_k) < 0$, as well as the corresponding envelope and squeezing functions $g_u$ and $g_l$. This can be efficiently achieved by starting from an initial guess and stepping out in steps of exponentially increasing size.\n",
"2. Sample \\\\[x' \\sim \\frac{g_u(x)}{\\int g_u(x') dx'} \\text{ and } z \\sim \\text{Unifrom}(0, 1),\\\\] and perform the following squeezing and rejection tests. If \\\\[z \\leq \\frac{g_l(x')}{g_u(x')}\\\\] holds, then accept $x'$ otherwise perform the following rejection test \\\\[z \\leq \\frac{h(x')}{g_u(x')}\\\\] If this holds, accept the point and otherwise reject it.\n",
"3. If $x'$ was accepted at the squeezing test, go to step 2 immediately. Otherwise insert $x'$ into $T_k$ to obtain $T_{k+1}$, update the piecewise exponential functions $g_l$ and $g_u$ accordingly and then return to step 2.\n",
" \n",
"
\n",
" \n",
"\n",
"Below are functions which implement envelope sampling, that is drawing \n",
"\n",
"$$ x' \\sim \\frac{g_u(x)}{\\int g_u(x') dx'}. $$\n",
"\n",
"The first function determines the left and right limits as well as the the unnormalised probabilities $\\int g_{u, k}(x') dx'$ of each piecewise exponential. The second samples uniformly from the area under the envelope function $g_u$."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"def envelope_limits_and_unnormalised_probabilities(xs, hs, dhdxs):\n",
" \n",
" # Compute the points of intersection of the lines making up the envelope\n",
" z, c = compute_points_of_intersection_and_intercepts(xs, hs, dhdxs)\n",
" \n",
" # Left-right endpoints for each piece in the piecewise envelope\n",
" limits = np.concatenate([[float('-inf')], z, [float('inf')]])\n",
" limits = np.stack([limits[:-1], limits[1:]], axis=-1)\n",
" \n",
" probs = (np.exp(dhdxs * limits[:, 1]) - np.exp(dhdxs * limits[:, 0])) * np.exp(c)\n",
" \n",
" # Catch any intervals where dhdx was zero\n",
" idx_nonzero = np.where(dhdxs != 0.)\n",
" probs[idx_nonzero] = probs[idx_nonzero] / dhdxs[idx_nonzero]\n",
" \n",
" idx_zero = np.where(dhdxs == 0.)\n",
" probs[idx_zero] = ((limits[:, 1] - limits[:, 0]) * np.exp(c))[idx_zero]\n",
" \n",
" return limits, probs\n",
"\n",
"\n",
"def sample_envelope(xs, hs, dhdxs):\n",
" \n",
" limits, probs = envelope_limits_and_unnormalised_probabilities(xs, hs, dhdxs)\n",
" probs = probs / np.sum(probs)\n",
" \n",
" # Randomly chosen interval in which the sample lies\n",
" i = np.random.choice(np.arange(probs.shape[0]), p=probs)\n",
" \n",
" # Sample u = Uniform(0, 1)\n",
" u = np.random.uniform()\n",
" \n",
" # Invert i^th piecewise exponential CDF to get a sample from that interval\n",
" if dhdxs[i] == 0.:\n",
" return u * (limits[i, 1] - limits[i, 0]) + limits[i, 0]\n",
" \n",
" else:\n",
" x = np.log(u * np.exp(dhdxs[i] * limits[i, 1]) \\\n",
" + (1 - u) * np.exp(dhdxs[i] * limits[i, 0]))\n",
" x = x / dhdxs[i] \n",
" \n",
" return x"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If we draw samples from the envelope defined by the three previous points without the rejection step, we obtain the following distribution."
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"tags": [
"hide-input",
"center-output"
]
},
"outputs": [
{
"data": {
"application/pdf": "JVBERi0xLjQKJazcIKu6CjEgMCBvYmoKPDwgL1R5cGUgL0NhdGFsb2cgL1BhZ2VzIDIgMCBSID4+CmVuZG9iago4IDAgb2JqCjw8IC9Gb250IDMgMCBSIC9YT2JqZWN0IDcgMCBSIC9FeHRHU3RhdGUgNCAwIFIgL1BhdHRlcm4gNSAwIFIKL1NoYWRpbmcgNiAwIFIgL1Byb2NTZXQgWyAvUERGIC9UZXh0IC9JbWFnZUIgL0ltYWdlQyAvSW1hZ2VJIF0gPj4KZW5kb2JqCjExIDAgb2JqCjw8IC9UeXBlIC9QYWdlIC9QYXJlbnQgMiAwIFIgL1Jlc291cmNlcyA4IDAgUgovTWVkaWFCb3ggWyAwIDAgMzcwLjYyMTg3NSAyMDIuMTQxODc1IF0gL0NvbnRlbnRzIDkgMCBSIC9Bbm5vdHMgMTAgMCBSID4+CmVuZG9iago5IDAgb2JqCjw8IC9MZW5ndGggMTIgMCBSIC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nL2cy64ctxGG9/0UvZQXarF451JOHAFeBLAjIEAuC0OWFQtnfBNi5/HzF3mmWeTUGMmCI+PAZ/7Tw/6KXc2qJov96o/vf/3+3fuv33y+/+Ev26v+6d2njfaP+Pmwm/0jfn7baX+Dnw+bwafL5pI5oqWcAj4+yY/W2IN8/fUJBw8f/7Vt322vXqOZT/jam22z+fza+Rsaj+7ws/wkZSr+KNdW+zGDXE/2866dwjl/5J1iPJzdf3m//3X/YX/12lasIxgq0aTsXbX+zTZK2/aS8PVCzhAN2C8Jn4opoQzYo5yDS/jtaW7l1Cv2V/sacAWwclvIkeIN+D1ds3MtuELC5BSh2+TzRHhP1y1aiq6hVHRz5FKcTTP6PV01aS26gsLoBs3i8ps4Id7TdZOWomsoQC+404r1eRxZ7sl3DFoJrpEwtz1yBoefuLtcDkcp2DyZI+Sl2AofsHOEHCIOGrAHuY+IQyMPGRA1PsaGv2aLfxN2l8th2vdGa4S8FFvhA3ZCoylnQyO2kKWTDNY8xEk0PmBHXOwU3DRU66pu4kpohYOZ4aqJbkK9kMsRPHIRO9oi1KXQCh6oAzw15ptOFbL068GYh/i1xsfY8NQYiKaIOMh9FBkaecgoovEB26NRTjWmaDjIAlu1ZiW2xgdsx6fPdkq+R7lj69asxNb4GBveGhDypggpZDn4iUYeNPhpfMC28NZA3k0RUsjlSJksOni25pRXYmt8jA1v9RiGpwgpZImtW7MUW+HjzBWNeh/KFCGFLAfAwZqHDIAaH+cjUG+Gv6vozGF94KxpsETKS3ORGzZ2jyO7lPwUHU8VbMWbhG+NdpziUte4Zbts7qAYYhrduYty6Ojff8zIoaBdtgIRj+DjqHGK6MpkcnZOGiHFhbgK2WWjcFCyGAIG3q464t8i0jpphRAXAmtsl83yyUNxYxTsKuCQZdjgZztOdSGyBgfkBBXimG90FXDR4b8yG3KqK5EVONx28E5E4TQObV2VY4Iw5DEjhQYHZPgnMp4y5tBdxSmhmRrhhCFCXYmswF02D//MBfKA3FVkFWQMuTwaItSFyBrcZQvwz+LMOF6cIvq1hGCJRjOEuhBYQQMvnLNEH8fEoqsuHM4jvY+DGUJcCaywXbYI3ywl5zHkdRUZMvI0smm0Q6gLkTU4IJfDGufMGPW6iv7MZGyeDBHqSmQF7rKlABXRd3TkrqI/Q7G5TIYIdSGyBnfZMh2WDPkx8nUVI4QJMUY3GiLUhcgaHJATVESGMfJ11SMsu2KTHwyR6kpkBQ4JnIMaSx4jX1e9O3y0vvanMESoK3M4BQ7I8E9rEHhH5FP1sT5xJDcaItSVyAoc8k4DB7Uu2zFaC9mnmmTSaEkXVyaeGhyYCR5qcU+N4VrIwaKlUHwcTZHySmyNj7HhpQ6DwUR9VQM4cw4hTbYIeSn0LR2YLdzUeTPNBwg5YCR2nvPkwRahroTW8JgajupSmCZfhFwjCUWTR2OkvBRb4QM2ArEFxzTVJeRkD+8TDxqDNVJeia3xARthwnrvpolFIQ/YujUrsTU+xoazejx6hAn7lBFVAoa7Ol8rGpHyUmyFjycK4K2BaJo0F3IuHLVdnRWQ1gh56WSBwsfY8NaABHka/bpc4hGiK96O1kh5KbbCB+wIbw3ZTAtCQualezwDWjuZM+grwTVCgCf4ayQ/rb8JmWys013trhQGSX0luEbI4PDY6PPU4adK3nK61B6wpDlSX4p9ywfqDIeN2U3LrkImXnULLjsazZHySmwNkLnhsAmp3BQpu0x4bik2tqdDaY/Ul4IrhAAvcNjESyYjeJfJ58N4ZH52MkjqK8E1wsuG/AIyL5WMk5BdZk/mG7BOhQ0GCX3lRKRGyODw2UxlKgESMjk8FXgXbwyS+lJwhZAnquG0OfipRknInJPg2fx5FlUaJPSls9UKIYMXTv3zVBUmZL4JY/G+BvXBIKEvBVcIAW7htMXaqRJPyBQMT+CUkCeDpL4SXCMEuIPTlhCn6kchc9bNc6hxNkjqK8E1QgaH05ZiZu6rStSmJZOdzJH6UuxbPlB7dzjDid6I3eWSD0/F1Ang3oZUV0JreExdIIfi7ER9yhnni2TqwoBoRMpLsRU+YIcAGY/kY8QUcgp17izTaI2UV2JrfMCOdDhkdWmKl12OXF1oU4yjNVJeia3xMXaCHBG+J+xTDvkoePY1frRGykuxFT5eIYW34jial0hPOdAR0AhNxnR16RKpQsfQ8FXrjJviZJd95GcxJIGjLVJeiq3wATvDV/EUEKYo2WUMHYjizqfJGiGvxNb4gF3gqxYuOsWaLrtcx7m6zCgakfJKbI2PseGsbppaO0XkSzynU58JpCVCXop8w3bZcKkhxmTH4ChkxMMYg2s51GmJEFeuo2t0DA0vRdrppzqhLstlMGnLY1bHVD5gEy80uxBH9xDyULQgrHlQ1YLGB2xLvD5e8hgbhex4Hg15aZmsEfJKbI2PseGrwXgz7V7rsiy5kdY8phJH5eNSHHhrcNmOsVHIsh5rtOYRZVoqH2PDW0Oyfi4hOmVn6xSriZM1Ql6KrfAB28Nbo4lxKiPqsqzflNY8pqxT5QM2MiEXvcljdBQy+HxC/C431pzySmyNj7HhrTGFKdZ0VWzLkE08ZreGRtc2wTqMvW6CFvJArdmyElrDY2p4asIwPEXILjvTstLRFqEuhVboeDMz/DTh3xQfhSxKfWUjD6oA1viAneGpmahM8VHI0kMGax7iIhofY8NVs480xcdBPneSjI08YoOJysd1n3DWnI2b4uMdWTdyaemnAsLY8NFCPkzxUchib4Ns5DFbHlS+y4bEArLPaYyPo9ydZGjkEU6i8gGb4K0l2zKV2g7yyXfHmpXYGh9jpwMZ/7T9YVA7tG7LUuhbOjBbB3XeUnBHVe1biaxgMHKBipuKJuY7smbKUmiFA9QuHLjoZaoUvifrxqzE1kCA7Qly8FO18Ch3rx4aeYhXa3yMDVelnKeK4VHu2Lo1S7EVPt5iwO9msXaqGr4n60Yu3WWggDA2vNXyzMyELeXe20MjD+ltjQ/YEd5qi5mqh0dZYKvWnNivXrv2Qp4PO+0f8fMbs/Dn7eeN8P+XBiIVwLUzkTlK2N9dts/f8uuADgrRhv3tt9urP+G4vL/9bnvxn8/2tx+3L95uX6ER/v5Lbol3VKbaRja8EZR/6+04m+mmoe9qQ5Fn0hDhDZ/JtEPs9ZB/vKjHpMPg2kQ8T6TzmIkH5DFY58URvZXP2iH2yDERujHY1A568ar9hfitN2TJemNvTvC3bvD/5wu18+3+Zev8I9R3LfFblr5+M72SSbzpyNnDJDhKrAn1kSjyhLlz3AW5LuhzDbPNGIp5Fs8YGw1nciEg9Teen8i5dB82Qoz2sMkQxd0BzWeqE36Rl2pKir5OuuJcdSEBTxHJ4tmhTirbhO5gNdvD82tSiIvqyWVXVwG5lBAXrKSdazu8rYtVxR4huODc7v3BS8pMUHAGY1IsvGgRChK5uq/O2CP6bGCl5zea8LQNq4wTrLe8LuNLei5Oo3oOj2QqeL4N2tMkvzYH33C4YgFN51JaqQ/hnki4lLEuSTkKrYSWN3ei8T2iZfhbqCe0ketm62KcPywVS9cKVdxUAYlQ5P0xQKrFT7g0aA/Re0+mctQUmtC1tkSewknu4BqA1KpfyxEjAGlPCT1rczul5ysZeAdI5l0Krq5gEzrHFRxbeLHVZjhCBfHlSCEEnKeg9ymmupjJVYZkkcrX4kPni2m1h/HwOcCInbsXow1s31ttXw6GV53J5IOQ+DfuyJX8PvudeG+xda374MMB/sflL+gz4w2wql6O4n0t/UL3ZHhrq/PhVSl+0OeKq4x7ObvqjpTiERO8CWf1XHiKhKfpOCjEwKUHuPToodQ6t/DLHZyBgVT3QmCwqyuMPPkLHyVf28cgato2Rx7wbC7Uil1w5+e21o6xjLh8ETrhTkrJtIoCh/EJrQCHHQTdHs9KA9iY0DcFPYxrWs/qMZwY9EPY8QBSkoutdZCVZHhLWqE6mLTGPa9hBI9uzf6IGWG7rnwF3peL+xOOkXmZ1KW6sgQnNjFw2XbCXRAjB0uWE/sUPJOdLkYLJ2U5opHik4MvOi6tN6WaD3/FheU5cS7dLMHX3uWVTXQ1HvS4/D6T8aUusiR+Iwy82/JWgtSuOsueqz75mRCOVrLztq0/4vyI4ZF3SwQ8KNa5PYtxHYdEfNMHTl9tnWXAQMQDSzB11Q+nc62nMEIkcg7DIW/GMdlWr+AlKoRAvvvQw46M8W0RCLeyDa7UDVIYwltZoTNtAhoDEe6xxC+jqNM0JvK9mnEBa27HYx6vEnBJJa4I8bDIg3zdmOLwgFBKrHvTcecluF6dLrYe/RDd825IS20mBS7jQ8QoynPfuFxUd4rwJG3GZTV1g2rEDdmmi7jTELVSnRs1LlMt9nUYYXC5eCsPjyMBraV9eMcdJ9Lw6Drm3ZG3v8gAzjEEV/A5fI+hQwtHl3tvzkOz6iv4LndfwYdv6Ge48yq/4XjR0u+e4dVr3/IUDpUfWp4Cc9trCeu9G2rBKcI2T/nGs52mtUlSDJTW5KvmMPiLo4rDjTapmY97t0kNXntyDTpHxNq8PFMu16MlU1ffXdG7xCMO/C1WGqFzadNwnquQZGNVatxP58fTONly64Sh595xXvj59rtpyexbuL8YAnfF9SoiRR6lVrsySbimfDHN9WLeJp249/hLqdt3JostFbymcKbmYH9+zgORM+B2v+Z4L358TiEJZ48ihXzxS/0DxkiOBTX3bPql6oiqHve8yCdffHNtCZkPhqN+iqf6B4yaKTsMYr2p7+/94VP9AydhGHjkOd5fzzEkqi++vebBPsXQTdirjEifkFjPqW/rlQ/zN1/SEUX3IcLU4/7dM9h7TyjtVaD1Qt2UwNXLi1SSE0h7LckapNMJeAb9lNyN9Pya0P/lqeTZQfo5moN0nzCtO/6+v/jih1/fP/340/v90zeXn57eo///ub/98jlj3/4LINwtRQplbmRzdHJlYW0KZW5kb2JqCjEyIDAgb2JqCjM3NjcKZW5kb2JqCjEwIDAgb2JqClsgXQplbmRvYmoKMTcgMCBvYmoKPDwgL0xlbmd0aCA4NyAvRmlsdGVyIC9GbGF0ZURlY29kZSA+PgpzdHJlYW0KeJw1jLsNwCAMBXtP8RaIhPEP9omiFGT/NiaExj6d9M6boeConDeKIGrHycTCHz605SBvHW6axFXRJMFSLQjJwKAVyh8Im94Crj43bJv+9qCbrhcelxi0CmVuZHN0cmVhbQplbmRvYmoKMTggMCBvYmoKPDwgL0xlbmd0aCAxNDcgL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicPY+7EcMwDEN7ToERzL80T3K5FM7+bWBFTiO+I0QC9DQcqMEneqDrwFPlxo9EBdXEKc7epshFbopi149GVeIhNg3FjTao2IR1IOdSuCejOR1qm3xOml5kqZs0ryin1K9o3ULnhjH2uPa9UmmXK5pORwWNlRFpbFFoxmQADjd/u9U6zfnzqo//sae85fUFtOQ0SAplbmRzdHJlYW0KZW5kb2JqCjE5IDAgb2JqCjw8IC9MZW5ndGggNDA0IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nD2SS47kMAxD9zmFLtCA9bOd81SjMYuq+2/n0ZnpRSDFliiK9My2YRn2FWNZ39N6Dvv266S17HP1SJvb3lftpqysetjXHVaD6LntdWXdpLsstj9QTv9JXpfP+neWJylzfkgGl1sYd1u3APah4kuYgryduQtIH2033e5t7hoWIehKpbf5OlySuftJXlccTNIcAvVlmXTuKcKLO+7zFpZb6RPZAnxr27wTPSwhn+p0OxMD/GUxA9kCMcZZcBM9oLZNbc3hdIN9TfMKhRiaq6SnKpJ/SGQXx6xWkPSgo4ChpLZ2kBSNWsH+HTSnbNHQmPjVacEG8ixZVJHd5nNSw6kAEX06NaY0Gh1qgQ/RuqmaWqvHY/f71/j39Qe2bgnlz5OpGu+LRYVZOATm2OAh8uKEF5LlJ8JDk8hC2+FzUFWYKWekDIrCuWpboFulblyseHC+JDuins0eq+dBcjyUsye65ihzdnb4OO9Zb8H9QXPbiIesSB3zBBHbR/qJOUM6gdqaUVil93m4Yer//V9o8fMXBjWWjgplbmRzdHJlYW0KZW5kb2JqCjIwIDAgb2JqCjw8IC9MZW5ndGggMjEzIC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nDVQOXLEMAzr9Qp8YGcEHpL4HmcyKTb/bwM63sIGRIkgwOWBidx40ZGnkLHxxbE2jBu/g8679h5mH0a9N5s3c9gMMAnWFOrbhWswDiiVG/M8yM+NuU4LnOrbjrP6QgKbmhtYAYtsuIaVN3E5PRO+NUFerxEzpZoIk56ch5eauiOqHqefRO8Rx9GOvf4xyJ6qwkk9gWfCCq5YtJ6q8qtdnHs1NutGuVy8WUeYKG0iVZUrhdNezu4/vYtNZKhBWbMUtRcgAUp9tSrxbPoaP+P7D9CGSicKZW5kc3RyZWFtCmVuZG9iagoyMSAwIG9iago8PCAvTGVuZ3RoIDk1IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nD2MQQ7AIAgE77xiP9AEERX/0zQ92P9fu0bbC0x2YUo3KA4rnFUVxRvOJB8+kr3DWseQoplHQ5zd3BYOS40Uq1gWFp5hEaS0Ncz4vChrYEop6mln9b+75XoB/58cLAplbmRzdHJlYW0KZW5kb2JqCjE1IDAgb2JqCjw8IC9UeXBlIC9Gb250IC9CYXNlRm9udCAvR0NXWERWK0RlamFWdVNhbnMtT2JsaXF1ZSAvRmlyc3RDaGFyIDAKL0xhc3RDaGFyIDI1NSAvRm9udERlc2NyaXB0b3IgMTQgMCBSIC9TdWJ0eXBlIC9UeXBlMwovTmFtZSAvR0NXWERWK0RlamFWdVNhbnMtT2JsaXF1ZSAvRm9udEJCb3ggWyAtMTAxNiAtMzUxIDE2NjAgMTA2OCBdCi9Gb250TWF0cml4IFsgMC4wMDEgMCAwIDAuMDAxIDAgMCBdIC9DaGFyUHJvY3MgMTYgMCBSCi9FbmNvZGluZyA8PCAvVHlwZSAvRW5jb2RpbmcgL0RpZmZlcmVuY2VzIFsgOTAgL1ogMTAyIC9mIC9nIDExNyAvdSAxMjAgL3ggXSA+PgovV2lkdGhzIDEzIDAgUiA+PgplbmRvYmoKMTQgMCBvYmoKPDwgL1R5cGUgL0ZvbnREZXNjcmlwdG9yIC9Gb250TmFtZSAvR0NXWERWK0RlamFWdVNhbnMtT2JsaXF1ZSAvRmxhZ3MgOTYKL0ZvbnRCQm94IFsgLTEwMTYgLTM1MSAxNjYwIDEwNjggXSAvQXNjZW50IDkyOSAvRGVzY2VudCAtMjM2IC9DYXBIZWlnaHQgMAovWEhlaWdodCAwIC9JdGFsaWNBbmdsZSAwIC9TdGVtViAwIC9NYXhXaWR0aCAxMzUwID4+CmVuZG9iagoxMyAwIG9iagpbIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwCjYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgMzE4IDQwMSA0NjAgODM4IDYzNgo5NTAgNzgwIDI3NSAzOTAgMzkwIDUwMCA4MzggMzE4IDM2MSAzMTggMzM3IDYzNiA2MzYgNjM2IDYzNiA2MzYgNjM2IDYzNiA2MzYKNjM2IDYzNiAzMzcgMzM3IDgzOCA4MzggODM4IDUzMSAxMDAwIDY4NCA2ODYgNjk4IDc3MCA2MzIgNTc1IDc3NSA3NTIgMjk1CjI5NSA2NTYgNTU3IDg2MyA3NDggNzg3IDYwMyA3ODcgNjk1IDYzNSA2MTEgNzMyIDY4NCA5ODkgNjg1IDYxMSA2ODUgMzkwIDMzNwozOTAgODM4IDUwMCA1MDAgNjEzIDYzNSA1NTAgNjM1IDYxNSAzNTIgNjM1IDYzNCAyNzggMjc4IDU3OSAyNzggOTc0IDYzNCA2MTIKNjM1IDYzNSA0MTEgNTIxIDM5MiA2MzQgNTkyIDgxOCA1OTIgNTkyIDUyNSA2MzYgMzM3IDYzNiA4MzggNjAwIDYzNiA2MDAgMzE4CjM1MiA1MTggMTAwMCA1MDAgNTAwIDUwMCAxMzUwIDYzNSA0MDAgMTA3MCA2MDAgNjg1IDYwMCA2MDAgMzE4IDMxOCA1MTggNTE4CjU5MCA1MDAgMTAwMCA1MDAgMTAwMCA1MjEgNDAwIDEwMjggNjAwIDUyNSA2MTEgMzE4IDQwMSA2MzYgNjM2IDYzNiA2MzYgMzM3CjUwMCA1MDAgMTAwMCA0NzEgNjE3IDgzOCAzNjEgMTAwMCA1MDAgNTAwIDgzOCA0MDEgNDAxIDUwMCA2MzYgNjM2IDMxOCA1MDAKNDAxIDQ3MSA2MTcgOTY5IDk2OSA5NjkgNTMxIDY4NCA2ODQgNjg0IDY4NCA2ODQgNjg0IDk3NCA2OTggNjMyIDYzMiA2MzIgNjMyCjI5NSAyOTUgMjk1IDI5NSA3NzUgNzQ4IDc4NyA3ODcgNzg3IDc4NyA3ODcgODM4IDc4NyA3MzIgNzMyIDczMiA3MzIgNjExIDYwOAo2MzAgNjEzIDYxMyA2MTMgNjEzIDYxMyA2MTMgOTk1IDU1MCA2MTUgNjE1IDYxNSA2MTUgMjc4IDI3OCAyNzggMjc4IDYxMiA2MzQKNjEyIDYxMiA2MTIgNjEyIDYxMiA4MzggNjEyIDYzNCA2MzQgNjM0IDYzNCA1OTIgNjM1IDU5MiBdCmVuZG9iagoxNiAwIG9iago8PCAvWiAxNyAwIFIgL2YgMTggMCBSIC9nIDE5IDAgUiAvdSAyMCAwIFIgL3ggMjEgMCBSID4+CmVuZG9iagoyNiAwIG9iago8PCAvTGVuZ3RoIDgxIC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nE3Nuw3AIAwE0J4pPALg/z5RlCLZv40NEaGxn3QnnWCHCm5xWAy0Oxyt+NRTmH3oHhKSUHPdRFgzJdqEpF/6yzDDmFjItq83V65yvhbcHIsKZW5kc3RyZWFtCmVuZG9iagoyNyAwIG9iago8PCAvTGVuZ3RoIDc3IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nDWNwQ3AMAgD/0zBCDiFUPapqj7S/b8tRHzsMwjserJwpEwT9hF8gf6c9NI4ULTITBlo2rO+2CS5g5cjlCea0qti9edFD90fyZ4YDAplbmRzdHJlYW0KZW5kb2JqCjI4IDAgb2JqCjw8IC9MZW5ndGggMzA3IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nD2SS24DMQxD9z6FLhDA+tme86Qoupjef9snJemKHNkWRWqWukxZUx6QNJOEf+nwcLGd8jtsz2Zm4Fqil4nllOfQFWLuonzZzEZdWSfF6oRmOrfoUTkXBzZNqp+rLKXdLngO1yaeW/YRP7zQoB7UNS4JN3RXo2UpNGOq+3/Se/yMMuBqTF1sUqt7HzxeRFXo6AdHiSJjlxfn40EJ6UrCaFqIlXdFA0Hu8rTKewnu295qyLIHqZjOOylmsOt0Ui5uF4chHsjyqPDlo9hrQs/4sCsl9EjYhjNyJ+5oxubUyOKQ/t6NBEuPrmgh8+CvbtYuYLxTOkViZE5yrGmLVU73UBTTucO9DBD1bEVDKXOR1epfw84La5ZsFnhK+gUeo90mSw5W2duoTu+tPNnQ9x9a13QfCmVuZHN0cmVhbQplbmRvYmoKMjkgMCBvYmoKPDwgL0xlbmd0aCAyMzEgL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicNU85kgQhDMt5hT4wVRjbQL+np7Y22Pl/upKZTpDwIcnTEx2ZeJkjI7Bmx9taZCBm4FNMxb/2tA8TqvfgHiKUiwthhpFw1qzjbp6OF/92lc9YB+82+IpZXhDYwkzWVxZnLtsFY2mcxDnJboxdE7GNda2nU1hHMKEMhHS2w5Qgc1Sk9MmOMuboOJEnnovv9tssdjl+DusLNo0hFef4KnqCNoOi7HnvAhpyQf9d3fgeRbvoJSAbCRbWUWLunOWEX712dB61KBJzQppBLhMhzekqphCaUKyzo6BSUXCpPqforJ9/5V9cLQplbmRzdHJlYW0KZW5kb2JqCjMwIDAgb2JqCjw8IC9MZW5ndGggMjQ5IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nD1QO45EIQzrOYUv8CTyI3AeRqstZu/frgOaKVBMfrYzJNARgUcMMZSv4yWtoK6Bv4tC8W7i64PCIKtDUiDOeg+IdOymNpETOh2cMz9hN2OOwEUxBpzpdKY9ByY5+8IKhHMbZexWSCeJqiKO6jOOKZ4qe594FiztyDZbJ5I95CDhUlKJyaWflMo/bcqUCjpm0QQsErngZBNNOMu7SVKMGZQy6h6mdiJ9rDzIozroZE3OrCOZ2dNP25n4HHC3X9pkTpXHdB7M+Jy0zoM5Fbr344k2B02N2ujs9xNpKi9Sux1anX51EpXdGOcYEpdnfxnfZP/5B/6HWiIKZW5kc3RyZWFtCmVuZG9iagozMSAwIG9iago8PCAvTGVuZ3RoIDcyIC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nDMyt1AwULA0ARKGFiYK5mYGCimGXEC+qYm5Qi4XSAzEygGzDIC0JZyCiGeAmCBtEMUgFkSxmYkZRB2cAZHL4EoDACXbFskKZW5kc3RyZWFtCmVuZG9iagozMiAwIG9iago8PCAvTGVuZ3RoIDQ3IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nDMyt1AwULA0ARKGFiYK5mYGCimGXJYQVi4XTCwHzALRlnAKIp7BlQYAuWcNJwplbmRzdHJlYW0KZW5kb2JqCjMzIDAgb2JqCjw8IC9MZW5ndGggMjU4IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nEWRS3IEIAhE956CI4D85DyTSmUxuf82Dc5kNnaXqP2ESiOmEiznFHkwfcnyzWS26Xc5VjsbBRRFKJjJVeixAqs7U8SZa4lq62Nl5LjTOwbFG85dOalkcaOMdVR1KnBMz5X1Ud35dlmUfUcOZQrYrHMcbODKbcMYJ0abre4O94kgTydTR8XtINnwByeNfZWrK3CdbPbRSzAOBP1CE5jki0DrDIHGzVP05BLs4+N254Fgb3kRSNkQyJEhGB2Cdp1c/+LW+b3/cYY7z7UZrhzv4neY1nbHX2KSFXMBi9wpqOdrLlrXGTrekzPH5Kb7hs65YJe7g0zv+T/Wz/r+Ax4pZvoKZW5kc3RyZWFtCmVuZG9iagozNCAwIG9iago8PCAvTGVuZ3RoIDE2MyAvRmlsdGVyIC9GbGF0ZURlY29kZSA+PgpzdHJlYW0KeJxFkDsSAyEMQ3tOoSP4IwM+z2YyKTb3b2PYbFLA01ggg7sTgtTagonogoe2Jd0F760EZ2P86TZuNRLkBHWAVqTjaJRSfbnFaZV08Wg2cysLrRMdZg56lKMZoBA6Fd7touRypu7O+UNw9V/1v2LdOZuJgcnKHQjN6lPc+TY7orq6yf6kx9ys134r7FVhaVlLywm3nbtmQAncUznaqz0/Hwo69gplbmRzdHJlYW0KZW5kb2JqCjM1IDAgb2JqCjw8IC9MZW5ndGggMjE4IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nD1QuY0EMQzLXYUaWMB67alnFotLpv/0SPn2ItEWRVIqNZmSKS91lCVZU946fJbEDnmG5W5kNiUqRS+TsCX30ArxfYnmFPfd1ZazQzSXaDl+CzMqqhsd00s2mnAqE7qg3MMz+g1tdANWhx6xWyDQpGDXtiByxw8YDMGZE4siDEpNBv+uco+fXosbPsPxQxSRkg7mNf9Y/fJzDa9TjyeRbm++4l6cqQ4DERySmrwjXVixLhIRaTVBTc/AWi2Au7de/hu0I7oMQPaJxHGaUo6hv2twpc8v5SdT2AplbmRzdHJlYW0KZW5kb2JqCjM2IDAgb2JqCjw8IC9MZW5ndGggMjM5IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nE1QyW0EMQz7uwo1MMDoHLseB4s8sv1/Q8oJkpdoS+Kh8pRblspl9yM5b8m65UOHTpVp8m7Qza+x/qMMAnb/UFQQrSWxSsxc0m6xNEkv2cM4jZdrtY7nqXuEWaN48OPY0ymB6T0ywWazvTkwqz3ODpBOuMav6tM7lSQDibqQ80KlCuse1CWijyvbmFKdTi3lGJef6Ht8jgA9xd6N3NHHyxeMRrUtqNFqlTgPMBNT0ZVxq5GBlBMGQ2dHVzQLpcjKekI1wo05oZm9w3BgA8uzhKSlrVK8D2UB6AJd2jrjNEqCjgDC3yiM9foGqvxeNwplbmRzdHJlYW0KZW5kb2JqCjM3IDAgb2JqCjw8IC9MZW5ndGggMTUwIC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nD1POQ7DMAzb/Qp+IIB1WLbekyLokP5/reigHQQRoHjIsqNjBQ4xgUnHHImXtB/8NI1ALIOaI3pCxTAscDZJgYdBVklXPlucTCG1BV5I+NYUPBtNZIKex4gduJPJ6APvUluxdegB7RBVJF2zDGOxrbpwPWEENcV6Feb1MGOJkkefVUvrAQUT+NL5f+5u73Z9AWDFMF0KZW5kc3RyZWFtCmVuZG9iagozOCAwIG9iago8PCAvTGVuZ3RoIDE1MSAvRmlsdGVyIC9GbGF0ZURlY29kZSA+PgpzdHJlYW0KeJw1j8sNwzAMQ++aggsE0M+yPE+KoId0/2slpwUMmDDJJ9kWg5GMQ0xg7Jhj4SVUT60+JCO3ukk5EKlQNwRPaEwMM5ykS+CV6b5lPLd4Oa3UBZ2QyN1ZbTREGc08hqOn7BVO+i9zk6iVGWUOyKxT3U7IyoIVmBfUct/PuK1yoRMetivDojvFiHKbGWn4/e+kN11f3Hgw1QplbmRzdHJlYW0KZW5kb2JqCjM5IDAgb2JqCjw8IC9MZW5ndGggMTYwIC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nEWQORIDMQgEc72CJ0hcgvesy7XB+v+pB9ZHoukCNBy6Fk3KehRoPumxRqG60GvoLEqSRMEWkh1Qp2OIOyhITEhjkki2HoMjmlizXZiZVCqzUuG0acXCv9la1chEjXCN/InpBlT8T+pclPBNg6+SMfoYVLw7g4xJ+F5F3Fox7f5EMLEZ9glvRSYFhImxqdm+z2CGzPcK1zjH8w1MgjfrCmVuZHN0cmVhbQplbmRvYmoKNDAgMCBvYmoKPDwgL0xlbmd0aCAzMzQgL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicLVJLcsUgDNtzCl2gM/gH5DzpdLp4vf+2kpNFRg5g9DHlholKfFkgt6PWxLeNzECF4a+rzIXPSNvIOojLkIu4ki2Fe0Qs5DHEPMSC76vxHh75rMzJswfGL9l3Dyv21IRlIePFGdphFcdhFeRYsHUhqnt4U6TDqSTY44v/PsVzLQQtfEbQgF/kn6+O4PmSFmn3mG3TrnqwTDuqpLAcbE9zXiZfWme5Oh7PB8n2rtgRUrsCFIW5M85z4SjTVka0FnY2SGpcbG+O/VhK0IVuXEaKI5CfqSI8oKTJzCYK4o+cHnIqA2Hqmq50chtVcaeezDWbi7czSWbrvkixmcJ5XTiz/gxTZrV5J89yotSpCO+xZ0vQ0Dmunr2WWWh0mxO8pITPxk5PTr5XM+shORUJqWJaV8FpFJliCdsSX1NRU5p6Gf778u7xO37+ASxzfHMKZW5kc3RyZWFtCmVuZG9iago0MSAwIG9iago8PCAvTGVuZ3RoIDU0IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nDM2NlcwAEJdSyMFYyDb3MhSIcWQy8jUBMzM5YIJ5nBZGINV5XAZQGmYohyuDK40APuEDh8KZW5kc3RyZWFtCmVuZG9iago0MiAwIG9iago8PCAvTGVuZ3RoIDE4IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nDM2tFAwgMMUQ640AB3mA1IKZW5kc3RyZWFtCmVuZG9iago0MyAwIG9iago8PCAvTGVuZ3RoIDc1IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nDO1NFIwUDA2ABKmZkYKpibmCimGXEA+iJXLZWhkCmblcBlZmilYWAAZJmbmUCGYhhwuY1NzoAFARcamYBqqP4crgysNAJWQEu8KZW5kc3RyZWFtCmVuZG9iagoyNCAwIG9iago8PCAvVHlwZSAvRm9udCAvQmFzZUZvbnQgL0JNUVFEVitEZWphVnVTYW5zIC9GaXJzdENoYXIgMCAvTGFzdENoYXIgMjU1Ci9Gb250RGVzY3JpcHRvciAyMyAwIFIgL1N1YnR5cGUgL1R5cGUzIC9OYW1lIC9CTVFRRFYrRGVqYVZ1U2FucwovRm9udEJCb3ggWyAtMTAyMSAtNDYzIDE3OTQgMTIzMyBdIC9Gb250TWF0cml4IFsgMC4wMDEgMCAwIDAuMDAxIDAgMCBdCi9DaGFyUHJvY3MgMjUgMCBSCi9FbmNvZGluZyA8PCAvVHlwZSAvRW5jb2RpbmcKL0RpZmZlcmVuY2VzIFsgMzIgL3NwYWNlIDQwIC9wYXJlbmxlZnQgL3BhcmVucmlnaHQgNDcgL3NsYXNoIDY5IC9FIDc4IC9OIDk3IC9hIDEwMCAvZAovZSAxMDUgL2kgMTA4IC9sIC9tIC9uIC9vIC9wIDExNCAvciAvcyAxMTggL3YgXQo+PgovV2lkdGhzIDIyIDAgUiA+PgplbmRvYmoKMjMgMCBvYmoKPDwgL1R5cGUgL0ZvbnREZXNjcmlwdG9yIC9Gb250TmFtZSAvQk1RUURWK0RlamFWdVNhbnMgL0ZsYWdzIDMyCi9Gb250QkJveCBbIC0xMDIxIC00NjMgMTc5NCAxMjMzIF0gL0FzY2VudCA5MjkgL0Rlc2NlbnQgLTIzNiAvQ2FwSGVpZ2h0IDAKL1hIZWlnaHQgMCAvSXRhbGljQW5nbGUgMCAvU3RlbVYgMCAvTWF4V2lkdGggMTM0MiA+PgplbmRvYmoKMjIgMCBvYmoKWyA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMAo2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDMxOCA0MDEgNDYwIDgzOCA2MzYKOTUwIDc4MCAyNzUgMzkwIDM5MCA1MDAgODM4IDMxOCAzNjEgMzE4IDMzNyA2MzYgNjM2IDYzNiA2MzYgNjM2IDYzNiA2MzYgNjM2CjYzNiA2MzYgMzM3IDMzNyA4MzggODM4IDgzOCA1MzEgMTAwMCA2ODQgNjg2IDY5OCA3NzAgNjMyIDU3NSA3NzUgNzUyIDI5NQoyOTUgNjU2IDU1NyA4NjMgNzQ4IDc4NyA2MDMgNzg3IDY5NSA2MzUgNjExIDczMiA2ODQgOTg5IDY4NSA2MTEgNjg1IDM5MCAzMzcKMzkwIDgzOCA1MDAgNTAwIDYxMyA2MzUgNTUwIDYzNSA2MTUgMzUyIDYzNSA2MzQgMjc4IDI3OCA1NzkgMjc4IDk3NCA2MzQgNjEyCjYzNSA2MzUgNDExIDUyMSAzOTIgNjM0IDU5MiA4MTggNTkyIDU5MiA1MjUgNjM2IDMzNyA2MzYgODM4IDYwMCA2MzYgNjAwIDMxOAozNTIgNTE4IDEwMDAgNTAwIDUwMCA1MDAgMTM0MiA2MzUgNDAwIDEwNzAgNjAwIDY4NSA2MDAgNjAwIDMxOCAzMTggNTE4IDUxOAo1OTAgNTAwIDEwMDAgNTAwIDEwMDAgNTIxIDQwMCAxMDIzIDYwMCA1MjUgNjExIDMxOCA0MDEgNjM2IDYzNiA2MzYgNjM2IDMzNwo1MDAgNTAwIDEwMDAgNDcxIDYxMiA4MzggMzYxIDEwMDAgNTAwIDUwMCA4MzggNDAxIDQwMSA1MDAgNjM2IDYzNiAzMTggNTAwCjQwMSA0NzEgNjEyIDk2OSA5NjkgOTY5IDUzMSA2ODQgNjg0IDY4NCA2ODQgNjg0IDY4NCA5NzQgNjk4IDYzMiA2MzIgNjMyIDYzMgoyOTUgMjk1IDI5NSAyOTUgNzc1IDc0OCA3ODcgNzg3IDc4NyA3ODcgNzg3IDgzOCA3ODcgNzMyIDczMiA3MzIgNzMyIDYxMSA2MDUKNjMwIDYxMyA2MTMgNjEzIDYxMyA2MTMgNjEzIDk4MiA1NTAgNjE1IDYxNSA2MTUgNjE1IDI3OCAyNzggMjc4IDI3OCA2MTIgNjM0CjYxMiA2MTIgNjEyIDYxMiA2MTIgODM4IDYxMiA2MzQgNjM0IDYzNCA2MzQgNTkyIDYzNSA1OTIgXQplbmRvYmoKMjUgMCBvYmoKPDwgL0UgMjYgMCBSIC9OIDI3IDAgUiAvYSAyOCAwIFIgL2QgMjkgMCBSIC9lIDMwIDAgUiAvaSAzMSAwIFIgL2wgMzIgMCBSCi9tIDMzIDAgUiAvbiAzNCAwIFIgL28gMzUgMCBSIC9wIDM2IDAgUiAvcGFyZW5sZWZ0IDM3IDAgUgovcGFyZW5yaWdodCAzOCAwIFIgL3IgMzkgMCBSIC9zIDQwIDAgUiAvc2xhc2ggNDEgMCBSIC9zcGFjZSA0MiAwIFIKL3YgNDMgMCBSID4+CmVuZG9iagozIDAgb2JqCjw8IC9GMSAxNSAwIFIgL0YyIDI0IDAgUiA+PgplbmRvYmoKNCAwIG9iago8PCAvQTEgPDwgL1R5cGUgL0V4dEdTdGF0ZSAvQ0EgMCAvY2EgMSA+PgovQTIgPDwgL1R5cGUgL0V4dEdTdGF0ZSAvQ0EgMC41IC9jYSAwLjUgPj4KL0EzIDw8IC9UeXBlIC9FeHRHU3RhdGUgL0NBIDEgL2NhIDEgPj4KL0E0IDw8IC9UeXBlIC9FeHRHU3RhdGUgL0NBIDAuOCAvY2EgMC44ID4+ID4+CmVuZG9iago1IDAgb2JqCjw8ID4+CmVuZG9iago2IDAgb2JqCjw8ID4+CmVuZG9iago3IDAgb2JqCjw8ID4+CmVuZG9iagoyIDAgb2JqCjw8IC9UeXBlIC9QYWdlcyAvS2lkcyBbIDExIDAgUiBdIC9Db3VudCAxID4+CmVuZG9iago0NCAwIG9iago8PCAvQ3JlYXRvciAoTWF0cGxvdGxpYiB2My45LjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcpCi9Qcm9kdWNlciAoTWF0cGxvdGxpYiBwZGYgYmFja2VuZCB2My45LjEpCi9DcmVhdGlvbkRhdGUgKEQ6MjAyNDA3MzEyMDQxNTUtMDQnMDAnKSA+PgplbmRvYmoKeHJlZgowIDQ1CjAwMDAwMDAwMDAgNjU1MzUgZiAKMDAwMDAwMDAxNiAwMDAwMCBuIAowMDAwMDEzNjIyIDAwMDAwIG4gCjAwMDAwMTMzMzEgMDAwMDAgbiAKMDAwMDAxMzM3NCAwMDAwMCBuIAowMDAwMDEzNTU5IDAwMDAwIG4gCjAwMDAwMTM1ODAgMDAwMDAgbiAKMDAwMDAxMzYwMSAwMDAwMCBuIAowMDAwMDAwMDY1IDAwMDAwIG4gCjAwMDAwMDAzNDQgMDAwMDAgbiAKMDAwMDAwNDIwNyAwMDAwMCBuIAowMDAwMDAwMjA4IDAwMDAwIG4gCjAwMDAwMDQxODYgMDAwMDAgbiAKMDAwMDAwNjEwNCAwMDAwMCBuIAowMDAwMDA1ODg5IDAwMDAwIG4gCjAwMDAwMDU1MzYgMDAwMDAgbiAKMDAwMDAwNzE1NyAwMDAwMCBuIAowMDAwMDA0MjI3IDAwMDAwIG4gCjAwMDAwMDQzODYgMDAwMDAgbiAKMDAwMDAwNDYwNiAwMDAwMCBuIAowMDAwMDA1MDgzIDAwMDAwIG4gCjAwMDAwMDUzNjkgMDAwMDAgbiAKMDAwMDAxMjA1MSAwMDAwMCBuIAowMDAwMDExODQ0IDAwMDAwIG4gCjAwMDAwMTE0MjAgMDAwMDAgbiAKMDAwMDAxMzEwNCAwMDAwMCBuIAowMDAwMDA3MjI5IDAwMDAwIG4gCjAwMDAwMDczODIgMDAwMDAgbiAKMDAwMDAwNzUzMSAwMDAwMCBuIAowMDAwMDA3OTExIDAwMDAwIG4gCjAwMDAwMDgyMTUgMDAwMDAgbiAKMDAwMDAwODUzNyAwMDAwMCBuIAowMDAwMDA4NjgxIDAwMDAwIG4gCjAwMDAwMDg4MDAgMDAwMDAgbiAKMDAwMDAwOTEzMSAwMDAwMCBuIAowMDAwMDA5MzY3IDAwMDAwIG4gCjAwMDAwMDk2NTggMDAwMDAgbiAKMDAwMDAwOTk3MCAwMDAwMCBuIAowMDAwMDEwMTkzIDAwMDAwIG4gCjAwMDAwMTA0MTcgMDAwMDAgbiAKMDAwMDAxMDY1MCAwMDAwMCBuIAowMDAwMDExMDU3IDAwMDAwIG4gCjAwMDAwMTExODMgMDAwMDAgbiAKMDAwMDAxMTI3MyAwMDAwMCBuIAowMDAwMDEzNjgyIDAwMDAwIG4gCnRyYWlsZXIKPDwgL1NpemUgNDUgL1Jvb3QgMSAwIFIgL0luZm8gNDQgMCBSID4+CnN0YXJ0eHJlZgoxMzgzOQolJUVPRgo=",
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"x_plot = np.linspace(-4., 4., 200)\n",
"_, probs = envelope_limits_and_unnormalised_probabilities(xs, hs, dhdxs)\n",
"\n",
"samples = [sample_envelope(xs, hs, dhdxs) for i in range(10000)]\n",
"gu = [g_u(x, xs, hs, dhdxs) / np.sum(probs) for x in x_plot]\n",
"\n",
"# Plot samples and envelope\n",
"plt.figure(figsize=(6, 3))\n",
"plt.plot(x_plot,\n",
" gu,\n",
" color='red',\n",
" label='Normalised $g_u$')\n",
"\n",
"plt.hist(samples,\n",
" density=True,\n",
" bins=100,\n",
" color='gray',\n",
" alpha=0.5,\n",
" label='Envelope samples')\n",
"\n",
"# Plot formatting\n",
"plt.title('', fontsize=20)\n",
"plt.xlim([-4, 4])\n",
"plt.ylim([0, 0.5])\n",
"plt.xticks([])\n",
"plt.yticks([])\n",
"plt.xlabel('$x$', fontsize=18)\n",
"plt.ylabel('$f(x)~/~Z$', fontsize=18)\n",
"plt.legend()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We still need to add the initialisation of the abcissa set, the (optional) squeezing test and the rejection test. For the initialisation step, we can start from an initial point, and search to the left and to the right in exponentially increasing step sizes, until we find a point on the left side with positive $h'$ and a point on the right with negative $h'$, and use these as end-points of the abscissa set. The following function `adaptive_rejection_sampling` implements this initialisation step together with the squeezing and rejection tests."
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"def initialise_abcissa(x0, log_unnorm_prob):\n",
" \n",
" # Expand to the left/right until the abcissa is correctly initialised\n",
" xs = np.array([x0])\n",
" hs, dhdxs = log_unnorm_prob(xs)\n",
" \n",
" dx = -1.\n",
" \n",
" while True:\n",
" \n",
" if dx < 0. and dhdxs[0] > 0.:\n",
" dx = 1.\n",
" \n",
" elif dx > 0. and dhdxs[-1] < 0.:\n",
" break\n",
" \n",
" insert_idx = 0 if dx < 0 else len(xs)\n",
" \n",
" x = xs[0 if dx < 0 else -1] + dx\n",
" \n",
" h, dhdx = log_unnorm_prob(x)\n",
" \n",
" xs = np.insert(xs, insert_idx, x)\n",
" hs = np.insert(hs, insert_idx, h)\n",
" dhdxs = np.insert(dhdxs, insert_idx, dhdx)\n",
" \n",
" dx = dx * 2\n",
" \n",
" return xs, hs, dhdxs"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This approach to initialising the abcissa set does not have any tunable parameters, except `x0`. Any initialisation method which guarantees $h'(x_1) > 0$ and $x'(x_K) < 0$ will give a valid abcissa set and this method is only a specific choice. Changing the `x0` value used to the does not significantly affect the efficiency of ARS, since the initialisation method will terminate quickly because of the exponentially increasing step sizes. This implementation assumes that the domain $D$ of $f$, that is the set of points where $f$ is non-zero, is all of $\\mathbb{R}$. If this is not the case, then this initialisation function will fail. A more robust initialisation method could use boolean comparisons of $h'$ values instead of $h$ values, setting $h = -\\infty$ outside $D$, but for the purposes of exposition, this illustration assumes that $D = \\mathbb{R}$ and does not bother further with this technical point. Putting the initialisation step together with the squeezing and rejection steps, we arrive at the complete ARS algorithm below."
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"def adaptive_rejection_sampling(x0, log_unnorm_prob, num_samples):\n",
" \n",
" xs, hs, dhdxs = initialise_abcissa(x0=x0, log_unnorm_prob=log_unnorm_prob)\n",
" \n",
" samples = []\n",
"\n",
" while len(samples) < num_samples:\n",
" \n",
" x = sample_envelope(xs, hs, dhdxs)\n",
" \n",
" gl = g_l(x, xs, hs)\n",
" gu = g_u(x, xs, hs, dhdxs)\n",
"\n",
" # Squeezing test\n",
" u = np.random.rand()\n",
"\n",
" if u * gu <= gl:\n",
" samples.append(x)\n",
" \n",
" h, dhdx = log_unnorm_prob(x)\n",
"\n",
" # Rejection test\n",
" if u * gu <= np.exp(h):\n",
" samples.append(x)\n",
"\n",
" i = np.searchsorted(xs, x)\n",
"\n",
" xs = np.insert(xs, i, x)\n",
" hs = np.insert(hs, i, h)\n",
" dhdxs = np.insert(dhdxs, i, dhdx)\n",
" \n",
" return samples"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we can finally use this function to sample from the example standard Gaussian distribution."
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"np.random.seed(0)\n",
"\n",
"target_mean = 0.\n",
"target_variance = 1.\n",
"\n",
"x0 = 1.\n",
"num_samples = 10000\n",
"\n",
"log_unnorm_prob = log_gaussian(mean=target_mean, variance=target_variance)\n",
"\n",
"samples = adaptive_rejection_sampling(x0=x0, log_unnorm_prob=log_unnorm_prob, num_samples=num_samples)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"tags": [
"hide-input",
"center-output"
]
},
"outputs": [
{
"data": {
"application/pdf": "JVBERi0xLjQKJazcIKu6CjEgMCBvYmoKPDwgL1R5cGUgL0NhdGFsb2cgL1BhZ2VzIDIgMCBSID4+CmVuZG9iago4IDAgb2JqCjw8IC9Gb250IDMgMCBSIC9YT2JqZWN0IDcgMCBSIC9FeHRHU3RhdGUgNCAwIFIgL1BhdHRlcm4gNSAwIFIKL1NoYWRpbmcgNiAwIFIgL1Byb2NTZXQgWyAvUERGIC9UZXh0IC9JbWFnZUIgL0ltYWdlQyAvSW1hZ2VJIF0gPj4KZW5kb2JqCjExIDAgb2JqCjw8IC9UeXBlIC9QYWdlIC9QYXJlbnQgMiAwIFIgL1Jlc291cmNlcyA4IDAgUgovTWVkaWFCb3ggWyAwIDAgMzcwLjYyMTg3NSAyMDIuMTQxODc1IF0gL0NvbnRlbnRzIDkgMCBSIC9Bbm5vdHMgMTAgMCBSID4+CmVuZG9iago5IDAgb2JqCjw8IC9MZW5ndGggMTIgMCBSIC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nL2bTY8ctxGG7/0r+rg6mMsqklXkUcqHAB8COBYQIB8HQ5YVCztJbAGxf37eImemSS43SA6tw0I77/awn2IXi1XF1uNvP/z7x/cf/vj2zf6bb7fH49P7zxvtn/Dzcff7J/z8stP+Fj8fN49Ply2od8KUNeHjU/+RPTuK9dcnXDx8/Pu2/bA9vsYwn/G1t9vG+f61+28YXIKLs/zUy1SiK7dRj2sGud7sp311ixCiyzuJuMD7zx/2P+3/2B9fc8VyyVMRrzmGav3bbZS2LWbns0YZoVN0IYfseWDu1OK8RM8RajdCp1bgb/YzkFdwl03w0LKnOBAfYnHKylFnM+7qicALNPCKKypRdQQ+1IKnnwPzZMZdPBN4wXbZNDhMVqbRkQ81eBeKr5N5DNCJJwKv0ABcXJJCKYzAdxVsmEuhMptxV89EXsBdtpxclhTLQHwXu8kcrPgCM7wAu2wF4VAocx5wDzWoK169p8GITjwReMUG4uxiyiSjDx9qQFTOKRGPdnTqmcgLuMtGPjpNMfnRjTs5FCeaUpHRlE49EXqJB2oi55PPYzw+1OgdkWqcTOnlM6EXdMasLkQlzRP0XU7JEb7mZbClV0+FXuCBmoOTGBLJSH3I2PBYBISTMZ18JvaKz7Cxl4WSp+jcyRIdvlYCj9b08qnYCz5gB3EcsDvQiH3IKg7rj0ZbDvFM5BUbkCO7FDhNYbqTMyId5CyjJb18JvaKz7CzQ/DNMvn1IRf4BIJFktGaXj4Ve8EHbIsG4JjD9SGTT64ULTVP7kYZ9DPBV4QAF7KKJMXJsw+ZiJ0ET23CO4N6/UzwFaGBq7OApmUCv8vwYodciTNPBvX6qeALQoBrdMjqeY7ch0whI4Ig9Q+TQb1+JviKEOAZOSfcNU27+yFjP3SW7SWdDOr1M8FXhAYuTrzkOXofMiW4RKbo42xQp58KviAEOJLR4lH4TT5+yIQiMXMMmSaDev1M8BWhgReHnzRH8UMmCViEJXKYDer0U8EXhJeNEY3hp6isBvBOpqSOs9ZCsbenk0/EXvIBG1kpthCeongnI94h/oWUwmxOp58JviI0cOsj5TRF8U4mq2g8MkGZDOr1U8EXhABn6yfFQmMnp5Mn8LVBZ4KvCAEerK2EPXzy8EOmiMDNrDwb1Otngq8IDdzaS5oKT+B3mby4GL1Xngzq9VPBF4QAj9ZlCmWK4p1sOWBKSZ8b1Olngq8IDby4pJ5l8vFDLkiqvPraIukG6eVTsRd8wE7WdBKZY/gh5+hUo2oZrenlM7FXfMAW6zwxtr4R+5A1oXSvMWSwppfPxF7xGbb1nspUZh5qwhaJKiHKaEsvnwr9nA7Maq2nJHPwPuQEGSl30dGWTj0TeoUH6mytJyppbGN2cqzdVsmjLZ16JvSKzqCt85R5jtqHjBKevEaabenkU7EXfMC2dkKMMpbzh4oYjWos+smWXj4TekFnzHZO5Mscrg+5b8J3tnyZ1vwS77IFb50nDX46KT3kEJyPIUiYjOnkE7GXfMAm6zwFiaOHdHJ/fDBa8yVOFZZ8hm2Np1LyRH1Xu6PGfogvcwK5ogMzW89J8McRupMLaq/kS55tuatnQq/wQB2s4cQyRete7ud6MOaLzPWKz7Ct3ZRLmfy6k5Eveewr4Zk1d/lU7AUfsK02IXjrSL1W1yaeCb3gAHOyPhPGnPz6BXlpy5nQKw6jtiZT9l4m6kPuV2M3yBdajis+eyHHWkwxzPG6kwfspTV37MfXob0A9HGn/RN+fjEW+7z9tBH+/cpDpAK4difyrqT9/WV7885eP3KUBGX0u++3x9/jury/+2F7+PXV/u7T9rt32zcYxL7/lY3E3mWtY2TvNNffjnECW4k7DfRDHQjTF6RQ9HYn3y7h2yV/fajXqBX0UXwJer9m4iFr3TImSBejvGqXMBJMJS4+sbaLHh7bXwizW4iJo+dnN/jzYfD/5wt18nn/uk2+S236p7e/pjezRFCawhGiItAiRatenVMWsrzHunfMuXopSWxnA5KchFjqCziSlIp1gZUdtsSU23sZHiPYpaouEmXMNCKip1Dbl5mtbkBQNLFIwK5kqlqg5FIzGMbM1hPkUh8zZmqHxaykYkVGwe8hlYhhk0sUkdWbqnBOr1L22ipIOdXur4ePIRKDMlidwjXFsFP2EHzElEd7kjm3jqvHGCVxrPl3wRKtF9vhNoYWsaIhMavUhjhFRPkIk6NVQom0XYwwmotS3BM5RSZZzbPTFUoxhXpuHkFVW4vwABhQsN6s9IsJQPU8idXlbO+77RLwm6dc74gJQ2Hu8WzEbMG2Um+JuREk2CHtiqtjyK3vbRmsosYpuxak4niQ7cwQGSLjSWHpZIfJS1mvZ29eA26xYyKjELeuv8DTinKCd9uax+OSSggA1K2Y7p2CbdLet+MagnfhEWCBuYRlVt+4sDOPUjA/1kXDDOYUmz2aHWIKPtWzENS9hdsRBGaOwY+bJljkPaa86sgSfVDYYSduEbUFtQ5/ch4uC+dBBDGPKe18P4tTyVjI0LOd3yK+VN0eGHm2BjA8W+FvdQ7gboS1hd3ETgow66WNX8hljACHMx1On66NbvAUuDfuK/CRHGs8RKHmKMNdSpXh39d2PtAEa4MmOTs4RixhGAQTimWgCB7jPdnbMZ0P1nnsGdkHTDir9c16m6ypzfhubXZ3c8DIvQsGhRsOc8YengWflWmOmSysYsw8PhOsJ+d9tA6unQoZQR2eEgJCEbJHDmfWVG4dbMWejWVgfoiRtBUwjFHwZaBVt/VFc+uzAkEQlszNGKs+lEjXNmYisnTQ3oYgSjUEAA8LxJsfVNyI4FKbWQhZKXjzbax7xHWqV8PnUODVHleylSBtZmC2RrU3cFRcSlzqWS7brOJmiOVwcsZ08a0RFVVh/C5WzURtnW1ESSzggjnCQ4VrwsurjDXkse0jsqIiDnWNYcoc4qhgeWBmU4Bj1iE0IQShGMXVeNDwEGntGCxmjwiBGIVHCD9uhxcIq9b4woJDSoQV0V7aYywBD4txMdZUklijNWeMgacfY30zCulIaw4Xxnrx1hSxZBDlZHsIcFxsp6HUAp/hdjlfi04EI4Yn2nYkar9ZUYd9LyXEVovkWPMWl661HnzdXsayu0Rb9VaXYC6xS/q6Gwhcrr68FwhbGzwk1JcmFfcurR4geBUedn3TDzsLNpIqI3ZhP+X6FmsKMV0zWjw+H2zXsrdx4YylZueYtmSt5NpozCy5BimMC9cst/aBYuNq8vDKctsvX5S3b/tEyPZibGDXNGjch1fb+uWlN54x7PLV6cuLr07jG+s7vPAK9nB9N9J/vcPj69jyPUs5PrZ8D+a218mx6uT+RUsk9Jomj5pllWo9OQSemxYQqLqrsAhEJrXYde+3XkMGcYcbdLu2Dt/fKZfb1T3Tob4f+A+9RRM7tIrc6xzj4o6Dqv3YcbLladTuVvd3a7PzfF7fWwr+ZnupGrhm4c+OCY4psn041pT52rl+pj8NOsLloce1fv0vAf9LRYCdZHHLls7XvNpbYuxbpvyX/eHb7y7/evrw+dX+t/3d1y1TrmlvTXmnZdbRIT7nu5VppY9W9jq83dzc39z8ZQOwXAYDrB5p1catSqjWPPzhWmrE7KPkWxnx8M9rlYLghtzrqFIefq5/QGi3ba+O3/RL1Qu2LOztXcny8N1tJOwYCOHHLZ7qHxCENQekDcdQP770h8/1D0gksEf3ZdHDh9s9hlro4ftbqRVV0mHCXmVs14p0Y66ufFekYRdE3gJXpaFI832RFpAlUJBWBM7D/HoFRtLZWd4P8qqrsrb/AOqSVSsKZW5kc3RyZWFtCmVuZG9iagoxMiAwIG9iagoyODY1CmVuZG9iagoxMCAwIG9iagpbIF0KZW5kb2JqCjE3IDAgb2JqCjw8IC9MZW5ndGggODcgL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicNYy7DcAgDAV7T/EWiITxD/aJohRk/zYmhMY+nfTOm6HgqJw3iiBqx8nEwh8+tOUgbx1umsRV0STBUi0IycCgFcofCJveAq4+N2yb/vagm64XHpcYtAplbmRzdHJlYW0KZW5kb2JqCjE4IDAgb2JqCjw8IC9MZW5ndGggMTQ3IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nD2PuxHDMAxDe06BEcy/NE9yuRTO/m1gRU4jviNEAvQ0HKjBJ3qg68BT5caPRAXVxCnO3qbIRW6KYtePRlXiITYNxY02qNiEdSDnUrgnozkdapt8TppeZKmbNK8op9SvaN1C54Yx9rj2vVJplyuaTkcFjZURaWxRaMZkAA43f7vVOs3586qP/7GnvOX1BbTkNEgKZW5kc3RyZWFtCmVuZG9iagoxOSAwIG9iago8PCAvTGVuZ3RoIDk1IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nD2MQQ7AIAgE77xiP9AEERX/0zQ92P9fu0bbC0x2YUo3KA4rnFUVxRvOJB8+kr3DWseQoplHQ5zd3BYOS40Uq1gWFp5hEaS0Ncz4vChrYEop6mln9b+75XoB/58cLAplbmRzdHJlYW0KZW5kb2JqCjE1IDAgb2JqCjw8IC9UeXBlIC9Gb250IC9CYXNlRm9udCAvR0NXWERWK0RlamFWdVNhbnMtT2JsaXF1ZSAvRmlyc3RDaGFyIDAKL0xhc3RDaGFyIDI1NSAvRm9udERlc2NyaXB0b3IgMTQgMCBSIC9TdWJ0eXBlIC9UeXBlMwovTmFtZSAvR0NXWERWK0RlamFWdVNhbnMtT2JsaXF1ZSAvRm9udEJCb3ggWyAtMTAxNiAtMzUxIDE2NjAgMTA2OCBdCi9Gb250TWF0cml4IFsgMC4wMDEgMCAwIDAuMDAxIDAgMCBdIC9DaGFyUHJvY3MgMTYgMCBSCi9FbmNvZGluZyA8PCAvVHlwZSAvRW5jb2RpbmcgL0RpZmZlcmVuY2VzIFsgOTAgL1ogMTAyIC9mIDEyMCAveCBdID4+Ci9XaWR0aHMgMTMgMCBSID4+CmVuZG9iagoxNCAwIG9iago8PCAvVHlwZSAvRm9udERlc2NyaXB0b3IgL0ZvbnROYW1lIC9HQ1dYRFYrRGVqYVZ1U2Fucy1PYmxpcXVlIC9GbGFncyA5NgovRm9udEJCb3ggWyAtMTAxNiAtMzUxIDE2NjAgMTA2OCBdIC9Bc2NlbnQgOTI5IC9EZXNjZW50IC0yMzYgL0NhcEhlaWdodCAwCi9YSGVpZ2h0IDAgL0l0YWxpY0FuZ2xlIDAgL1N0ZW1WIDAgL01heFdpZHRoIDEzNTAgPj4KZW5kb2JqCjEzIDAgb2JqClsgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAKNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCAzMTggNDAxIDQ2MCA4MzggNjM2Cjk1MCA3ODAgMjc1IDM5MCAzOTAgNTAwIDgzOCAzMTggMzYxIDMxOCAzMzcgNjM2IDYzNiA2MzYgNjM2IDYzNiA2MzYgNjM2IDYzNgo2MzYgNjM2IDMzNyAzMzcgODM4IDgzOCA4MzggNTMxIDEwMDAgNjg0IDY4NiA2OTggNzcwIDYzMiA1NzUgNzc1IDc1MiAyOTUKMjk1IDY1NiA1NTcgODYzIDc0OCA3ODcgNjAzIDc4NyA2OTUgNjM1IDYxMSA3MzIgNjg0IDk4OSA2ODUgNjExIDY4NSAzOTAgMzM3CjM5MCA4MzggNTAwIDUwMCA2MTMgNjM1IDU1MCA2MzUgNjE1IDM1MiA2MzUgNjM0IDI3OCAyNzggNTc5IDI3OCA5NzQgNjM0IDYxMgo2MzUgNjM1IDQxMSA1MjEgMzkyIDYzNCA1OTIgODE4IDU5MiA1OTIgNTI1IDYzNiAzMzcgNjM2IDgzOCA2MDAgNjM2IDYwMCAzMTgKMzUyIDUxOCAxMDAwIDUwMCA1MDAgNTAwIDEzNTAgNjM1IDQwMCAxMDcwIDYwMCA2ODUgNjAwIDYwMCAzMTggMzE4IDUxOCA1MTgKNTkwIDUwMCAxMDAwIDUwMCAxMDAwIDUyMSA0MDAgMTAyOCA2MDAgNTI1IDYxMSAzMTggNDAxIDYzNiA2MzYgNjM2IDYzNiAzMzcKNTAwIDUwMCAxMDAwIDQ3MSA2MTcgODM4IDM2MSAxMDAwIDUwMCA1MDAgODM4IDQwMSA0MDEgNTAwIDYzNiA2MzYgMzE4IDUwMAo0MDEgNDcxIDYxNyA5NjkgOTY5IDk2OSA1MzEgNjg0IDY4NCA2ODQgNjg0IDY4NCA2ODQgOTc0IDY5OCA2MzIgNjMyIDYzMiA2MzIKMjk1IDI5NSAyOTUgMjk1IDc3NSA3NDggNzg3IDc4NyA3ODcgNzg3IDc4NyA4MzggNzg3IDczMiA3MzIgNzMyIDczMiA2MTEgNjA4CjYzMCA2MTMgNjEzIDYxMyA2MTMgNjEzIDYxMyA5OTUgNTUwIDYxNSA2MTUgNjE1IDYxNSAyNzggMjc4IDI3OCAyNzggNjEyIDYzNAo2MTIgNjEyIDYxMiA2MTIgNjEyIDgzOCA2MTIgNjM0IDYzNCA2MzQgNjM0IDU5MiA2MzUgNTkyIF0KZW5kb2JqCjE2IDAgb2JqCjw8IC9aIDE3IDAgUiAvZiAxOCAwIFIgL3ggMTkgMCBSID4+CmVuZG9iagoyNCAwIG9iago8PCAvTGVuZ3RoIDc3IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nDWNwQ3AMAgD/0zBCDiFUPapqj7S/b8tRHzsMwjserJwpEwT9hF8gf6c9NI4ULTITBlo2rO+2CS5g5cjlCea0qti9edFD90fyZ4YDAplbmRzdHJlYW0KZW5kb2JqCjI1IDAgb2JqCjw8IC9MZW5ndGggMzQxIC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nDVSO9KbQQjrv1PoAp5Z3st5nMmk+HP/NgI7FSywQgLSAgeZeIkhqlGu+CVPMF4n8He9PI2fx7uQWvBUpB+4Nm3j/VizJgqWRiyF2ce+HyXkeGr8GwI9F2nCjExGDiQDcb/W5896kymH34A0bU4fJUkPogW7W8OOLwsySHpSw5Kd/LCuBVYXoQlzY00kI6dWpub52DNcxhNjJKiaBSTpE/epghFpxmPnrCUPMhxP9eLFr7fxWuYx9bKqQMY2wRxsJzPhFEUE4heUJDdxF00dxdHMWHO70FBS5L67h5OTXveXk6jAKyGcxVrCMUNPWeZkp0EJVK2cADOs174wTtNGCXdqur0r9vXzzCSM2xx2VkqmwTkO7mWTOYJkrzsmbMLjEPPePYKRmDe/iy2CK5c512T6sR9FG+mD4vqcqymzFSX8Q5U8seIa/5/f+/nz/P4HjCh+IwplbmRzdHJlYW0KZW5kb2JqCjI2IDAgb2JqCjw8IC9MZW5ndGggMzA3IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nD2SS24DMQxD9z6FLhDA+tme86Qoupjef9snJemKHNkWRWqWukxZUx6QNJOEf+nwcLGd8jtsz2Zm4Fqil4nllOfQFWLuonzZzEZdWSfF6oRmOrfoUTkXBzZNqp+rLKXdLngO1yaeW/YRP7zQoB7UNS4JN3RXo2UpNGOq+3/Se/yMMuBqTF1sUqt7HzxeRFXo6AdHiSJjlxfn40EJ6UrCaFqIlXdFA0Hu8rTKewnu295qyLIHqZjOOylmsOt0Ui5uF4chHsjyqPDlo9hrQs/4sCsl9EjYhjNyJ+5oxubUyOKQ/t6NBEuPrmgh8+CvbtYuYLxTOkViZE5yrGmLVU73UBTTucO9DBD1bEVDKXOR1epfw84La5ZsFnhK+gUeo90mSw5W2duoTu+tPNnQ9x9a13QfCmVuZHN0cmVhbQplbmRvYmoKMjcgMCBvYmoKPDwgL0xlbmd0aCAyMzEgL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicNU85kgQhDMt5hT4wVRjbQL+np7Y22Pl/upKZTpDwIcnTEx2ZeJkjI7Bmx9taZCBm4FNMxb/2tA8TqvfgHiKUiwthhpFw1qzjbp6OF/92lc9YB+82+IpZXhDYwkzWVxZnLtsFY2mcxDnJboxdE7GNda2nU1hHMKEMhHS2w5Qgc1Sk9MmOMuboOJEnnovv9tssdjl+DusLNo0hFef4KnqCNoOi7HnvAhpyQf9d3fgeRbvoJSAbCRbWUWLunOWEX712dB61KBJzQppBLhMhzekqphCaUKyzo6BSUXCpPqforJ9/5V9cLQplbmRzdHJlYW0KZW5kb2JqCjI4IDAgb2JqCjw8IC9MZW5ndGggMjQ5IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nD1QO45EIQzrOYUv8CTyI3AeRqstZu/frgOaKVBMfrYzJNARgUcMMZSv4yWtoK6Bv4tC8W7i64PCIKtDUiDOeg+IdOymNpETOh2cMz9hN2OOwEUxBpzpdKY9ByY5+8IKhHMbZexWSCeJqiKO6jOOKZ4qe594FiztyDZbJ5I95CDhUlKJyaWflMo/bcqUCjpm0QQsErngZBNNOMu7SVKMGZQy6h6mdiJ9rDzIozroZE3OrCOZ2dNP25n4HHC3X9pkTpXHdB7M+Jy0zoM5Fbr344k2B02N2ujs9xNpKi9Sux1anX51EpXdGOcYEpdnfxnfZP/5B/6HWiIKZW5kc3RyZWFtCmVuZG9iagoyOSAwIG9iago8PCAvTGVuZ3RoIDcyIC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nDMyt1AwULA0ARKGFiYK5mYGCimGXEC+qYm5Qi4XSAzEygGzDIC0JZyCiGeAmCBtEMUgFkSxmYkZRB2cAZHL4EoDACXbFskKZW5kc3RyZWFtCmVuZG9iagozMCAwIG9iago8PCAvTGVuZ3RoIDQ3IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nDMyt1AwULA0ARKGFiYK5mYGCimGXJYQVi4XTCwHzALRlnAKIp7BlQYAuWcNJwplbmRzdHJlYW0KZW5kb2JqCjMxIDAgb2JqCjw8IC9MZW5ndGggMjU4IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nEWRS3IEIAhE956CI4D85DyTSmUxuf82Dc5kNnaXqP2ESiOmEiznFHkwfcnyzWS26Xc5VjsbBRRFKJjJVeixAqs7U8SZa4lq62Nl5LjTOwbFG85dOalkcaOMdVR1KnBMz5X1Ud35dlmUfUcOZQrYrHMcbODKbcMYJ0abre4O94kgTydTR8XtINnwByeNfZWrK3CdbPbRSzAOBP1CE5jki0DrDIHGzVP05BLs4+N254Fgb3kRSNkQyJEhGB2Cdp1c/+LW+b3/cYY7z7UZrhzv4neY1nbHX2KSFXMBi9wpqOdrLlrXGTrekzPH5Kb7hs65YJe7g0zv+T/Wz/r+Ax4pZvoKZW5kc3RyZWFtCmVuZG9iagozMiAwIG9iago8PCAvTGVuZ3RoIDIxOCAvRmlsdGVyIC9GbGF0ZURlY29kZSA+PgpzdHJlYW0KeJw9ULmNBDEMy12FGljAeu2pZxaLS6b/9Ej59iLRFkVSKjWZkikvdZQlWVPeOnyWxA55huVuZDYlKkUvk7Al99AK8X2J5hT33dWWs0M0l2g5fgszKqobHdNLNppwKhO6oNzDM/oNbXQDVocesVsg0KRg17YgcscPGAzBmROLIgxKTQb/rnKPn16LGz7D8UMUkZIO5jX/WP3ycw2vU48nkW5vvuJenKkOAxEckpq8I11YsS4SEWk1QU3PwFotgLu3Xv4btCO6DED2icRxmlKOob9rcKXPL+UnU9gKZW5kc3RyZWFtCmVuZG9iagozMyAwIG9iago8PCAvTGVuZ3RoIDIzOSAvRmlsdGVyIC9GbGF0ZURlY29kZSA+PgpzdHJlYW0KeJxNUMltBDEM+7sKNTDA6By7HgeLPLL9f0PKCZKXaEviofKUW5bKZfcjOW/JuuVDh06VafJu0M2vsf6jDAJ2/1BUEK0lsUrMXNJusTRJL9nDOI2Xa7WO56l7hFmjePDj2NMpgek9MsFms705MKs9zg6QTrjGr+rTO5UkA4m6kPNCpQrrHtQloo8r25hSnU4t5RiXn+h7fI4APcXejdzRx8sXjEa1LajRapU4DzATU9GVcauRgZQTBkNnR1c0C6XIynpCNcKNOaGZvcNwYAPLs4Skpa1SvA9lAegCXdo64zRKgo4Awt8ojPX6Bqr8XjcKZW5kc3RyZWFtCmVuZG9iagozNCAwIG9iago8PCAvTGVuZ3RoIDE1MCAvRmlsdGVyIC9GbGF0ZURlY29kZSA+PgpzdHJlYW0KeJw9TzkOwzAM2/0KfiCAdVi23pMi6JD+f63ooB0EEaB4yLKjYwUOMYFJxxyJl7Qf/DSNQCyDmiN6QsUwLHA2SYGHQVZJVz5bnEwhtQVeSPjWFDwbTWSCnseIHbiTyegD71JbsXXoAe0QVSRdswxjsa26cD1hBDXFehXm9TBjiZJHn1VL6wEFE/jS+X/ubu92fQFgxTBdCmVuZHN0cmVhbQplbmRvYmoKMzUgMCBvYmoKPDwgL0xlbmd0aCAxNTEgL0ZpbHRlciAvRmxhdGVEZWNvZGUgPj4Kc3RyZWFtCnicNY/LDcMwDEPvmoILBNDPsjxPiqCHdP9rJacFDJgwySfZFoORjENMYOyYY+ElVE+tPiQjt7pJORCpUDcET2hMDDOcpEvglem+ZTy3eDmt1AWdkMjdWW00RBnNPIajp+wVTvovc5OolRllDsisU91OyMqCFZgX1HLfz7itcqETHrYrw6I7xYhymxlp+P3vpDddX9x4MNUKZW5kc3RyZWFtCmVuZG9iagozNiAwIG9iago8PCAvTGVuZ3RoIDE2MCAvRmlsdGVyIC9GbGF0ZURlY29kZSA+PgpzdHJlYW0KeJxFkDkSAzEIBHO9gidIXIL3rMu1wfr/qQfWR6LpAjQcuhZNynoUaD7psUahutBr6CxKkkTBFpIdUKdjiDsoSExIY5JIth6DI5pYs12YmVQqs1LhtGnFwr/ZWtXIRI1wjfyJ6QZU/E/qXJTwTYOvkjH6GFS8O4OMSfheRdxaMe3+RDCxGfYJb0UmBYSJsanZvs9ghsz3Ctc4x/MNTII36wplbmRzdHJlYW0KZW5kb2JqCjM3IDAgb2JqCjw8IC9MZW5ndGggMzM0IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlID4+CnN0cmVhbQp4nC1SS3LFIAzbcwpdoDP4B+Q86XS6eL3/tpKTRUYOYPQx5YaJSnxZILej1sS3jcxAheGvq8yFz0jbyDqIy5CLuJIthXtELOQxxDzEgu+r8R4e+azMybMHxi/Zdw8r9tSEZSHjxRnaYRXHYRXkWLB1Iap7eFOkw6kk2OOL/z7Fcy0ELXxG0IBf5J+vjuD5khZp95ht0656sEw7qqSwHGxPc14mX1pnuToezwfJ9q7YEVK7AhSFuTPOc+Eo01ZGtBZ2NkhqXGxvjv1YStCFblxGiiOQn6kiPKCkycwmCuKPnB5yKgNh6pqudHIbVXGnnsw1m4u3M0lm675IsZnCeV04s/4MU2a1eSfPcqLUqQjvsWdL0NA5rp69lllodJsTvKSEz8ZOT06+VzPrITkVCaliWlfBaRSZYgnbEl9TUVOaehn++/Lu8Tt+/gEsc3xzCmVuZHN0cmVhbQplbmRvYmoKMzggMCBvYmoKPDwgL0xlbmd0aCA1NCAvRmlsdGVyIC9GbGF0ZURlY29kZSA+PgpzdHJlYW0KeJwzNjZXMABCXUsjBWMg29zIUiHFkMvI1ATMzOWCCeZwWRiDVeVwGUBpmKIcrgyuNAD7hA4fCmVuZHN0cmVhbQplbmRvYmoKMzkgMCBvYmoKPDwgL0xlbmd0aCAxOCAvRmlsdGVyIC9GbGF0ZURlY29kZSA+PgpzdHJlYW0KeJwzNrRQMIDDFEOuNAAd5gNSCmVuZHN0cmVhbQplbmRvYmoKMjIgMCBvYmoKPDwgL1R5cGUgL0ZvbnQgL0Jhc2VGb250IC9CTVFRRFYrRGVqYVZ1U2FucyAvRmlyc3RDaGFyIDAgL0xhc3RDaGFyIDI1NQovRm9udERlc2NyaXB0b3IgMjEgMCBSIC9TdWJ0eXBlIC9UeXBlMyAvTmFtZSAvQk1RUURWK0RlamFWdVNhbnMKL0ZvbnRCQm94IFsgLTEwMjEgLTQ2MyAxNzk0IDEyMzMgXSAvRm9udE1hdHJpeCBbIDAuMDAxIDAgMCAwLjAwMSAwIDAgXQovQ2hhclByb2NzIDIzIDAgUgovRW5jb2RpbmcgPDwgL1R5cGUgL0VuY29kaW5nCi9EaWZmZXJlbmNlcyBbIDMyIC9zcGFjZSA0MCAvcGFyZW5sZWZ0IC9wYXJlbnJpZ2h0IDQ3IC9zbGFzaCA3OCAvTiA4MyAvUyA5NyAvYSAxMDAgL2QKL2UgMTA1IC9pIDEwOCAvbCAvbSAxMTEgL28gL3AgMTE0IC9yIC9zIF0KPj4KL1dpZHRocyAyMCAwIFIgPj4KZW5kb2JqCjIxIDAgb2JqCjw8IC9UeXBlIC9Gb250RGVzY3JpcHRvciAvRm9udE5hbWUgL0JNUVFEVitEZWphVnVTYW5zIC9GbGFncyAzMgovRm9udEJCb3ggWyAtMTAyMSAtNDYzIDE3OTQgMTIzMyBdIC9Bc2NlbnQgOTI5IC9EZXNjZW50IC0yMzYgL0NhcEhlaWdodCAwCi9YSGVpZ2h0IDAgL0l0YWxpY0FuZ2xlIDAgL1N0ZW1WIDAgL01heFdpZHRoIDEzNDIgPj4KZW5kb2JqCjIwIDAgb2JqClsgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAKNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCAzMTggNDAxIDQ2MCA4MzggNjM2Cjk1MCA3ODAgMjc1IDM5MCAzOTAgNTAwIDgzOCAzMTggMzYxIDMxOCAzMzcgNjM2IDYzNiA2MzYgNjM2IDYzNiA2MzYgNjM2IDYzNgo2MzYgNjM2IDMzNyAzMzcgODM4IDgzOCA4MzggNTMxIDEwMDAgNjg0IDY4NiA2OTggNzcwIDYzMiA1NzUgNzc1IDc1MiAyOTUKMjk1IDY1NiA1NTcgODYzIDc0OCA3ODcgNjAzIDc4NyA2OTUgNjM1IDYxMSA3MzIgNjg0IDk4OSA2ODUgNjExIDY4NSAzOTAgMzM3CjM5MCA4MzggNTAwIDUwMCA2MTMgNjM1IDU1MCA2MzUgNjE1IDM1MiA2MzUgNjM0IDI3OCAyNzggNTc5IDI3OCA5NzQgNjM0IDYxMgo2MzUgNjM1IDQxMSA1MjEgMzkyIDYzNCA1OTIgODE4IDU5MiA1OTIgNTI1IDYzNiAzMzcgNjM2IDgzOCA2MDAgNjM2IDYwMCAzMTgKMzUyIDUxOCAxMDAwIDUwMCA1MDAgNTAwIDEzNDIgNjM1IDQwMCAxMDcwIDYwMCA2ODUgNjAwIDYwMCAzMTggMzE4IDUxOCA1MTgKNTkwIDUwMCAxMDAwIDUwMCAxMDAwIDUyMSA0MDAgMTAyMyA2MDAgNTI1IDYxMSAzMTggNDAxIDYzNiA2MzYgNjM2IDYzNiAzMzcKNTAwIDUwMCAxMDAwIDQ3MSA2MTIgODM4IDM2MSAxMDAwIDUwMCA1MDAgODM4IDQwMSA0MDEgNTAwIDYzNiA2MzYgMzE4IDUwMAo0MDEgNDcxIDYxMiA5NjkgOTY5IDk2OSA1MzEgNjg0IDY4NCA2ODQgNjg0IDY4NCA2ODQgOTc0IDY5OCA2MzIgNjMyIDYzMiA2MzIKMjk1IDI5NSAyOTUgMjk1IDc3NSA3NDggNzg3IDc4NyA3ODcgNzg3IDc4NyA4MzggNzg3IDczMiA3MzIgNzMyIDczMiA2MTEgNjA1CjYzMCA2MTMgNjEzIDYxMyA2MTMgNjEzIDYxMyA5ODIgNTUwIDYxNSA2MTUgNjE1IDYxNSAyNzggMjc4IDI3OCAyNzggNjEyIDYzNAo2MTIgNjEyIDYxMiA2MTIgNjEyIDgzOCA2MTIgNjM0IDYzNCA2MzQgNjM0IDU5MiA2MzUgNTkyIF0KZW5kb2JqCjIzIDAgb2JqCjw8IC9OIDI0IDAgUiAvUyAyNSAwIFIgL2EgMjYgMCBSIC9kIDI3IDAgUiAvZSAyOCAwIFIgL2kgMjkgMCBSIC9sIDMwIDAgUgovbSAzMSAwIFIgL28gMzIgMCBSIC9wIDMzIDAgUiAvcGFyZW5sZWZ0IDM0IDAgUiAvcGFyZW5yaWdodCAzNSAwIFIKL3IgMzYgMCBSIC9zIDM3IDAgUiAvc2xhc2ggMzggMCBSIC9zcGFjZSAzOSAwIFIgPj4KZW5kb2JqCjMgMCBvYmoKPDwgL0YxIDE1IDAgUiAvRjIgMjIgMCBSID4+CmVuZG9iago0IDAgb2JqCjw8IC9BMSA8PCAvVHlwZSAvRXh0R1N0YXRlIC9DQSAwIC9jYSAxID4+Ci9BMiA8PCAvVHlwZSAvRXh0R1N0YXRlIC9DQSAwLjUgL2NhIDAuNSA+PgovQTMgPDwgL1R5cGUgL0V4dEdTdGF0ZSAvQ0EgMSAvY2EgMSA+PgovQTQgPDwgL1R5cGUgL0V4dEdTdGF0ZSAvQ0EgMC44IC9jYSAwLjggPj4gPj4KZW5kb2JqCjUgMCBvYmoKPDwgPj4KZW5kb2JqCjYgMCBvYmoKPDwgPj4KZW5kb2JqCjcgMCBvYmoKPDwgPj4KZW5kb2JqCjIgMCBvYmoKPDwgL1R5cGUgL1BhZ2VzIC9LaWRzIFsgMTEgMCBSIF0gL0NvdW50IDEgPj4KZW5kb2JqCjQwIDAgb2JqCjw8IC9DcmVhdG9yIChNYXRwbG90bGliIHYzLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZykKL1Byb2R1Y2VyIChNYXRwbG90bGliIHBkZiBiYWNrZW5kIHYzLjkuMSkKL0NyZWF0aW9uRGF0ZSAoRDoyMDI0MDczMTE3NDY1Ni0wNCcwMCcpID4+CmVuZG9iagp4cmVmCjAgNDEKMDAwMDAwMDAwMCA2NTUzNSBmIAowMDAwMDAwMDE2IDAwMDAwIG4gCjAwMDAwMTE3NzkgMDAwMDAgbiAKMDAwMDAxMTQ4OCAwMDAwMCBuIAowMDAwMDExNTMxIDAwMDAwIG4gCjAwMDAwMTE3MTYgMDAwMDAgbiAKMDAwMDAxMTczNyAwMDAwMCBuIAowMDAwMDExNzU4IDAwMDAwIG4gCjAwMDAwMDAwNjUgMDAwMDAgbiAKMDAwMDAwMDM0NCAwMDAwMCBuIAowMDAwMDAzMzA1IDAwMDAwIG4gCjAwMDAwMDAyMDggMDAwMDAgbiAKMDAwMDAwMzI4NCAwMDAwMCBuIAowMDAwMDA0NDI5IDAwMDAwIG4gCjAwMDAwMDQyMTQgMDAwMDAgbiAKMDAwMDAwMzg3MSAwMDAwMCBuIAowMDAwMDA1NDgyIDAwMDAwIG4gCjAwMDAwMDMzMjUgMDAwMDAgbiAKMDAwMDAwMzQ4NCAwMDAwMCBuIAowMDAwMDAzNzA0IDAwMDAwIG4gCjAwMDAwMTAyMjggMDAwMDAgbiAKMDAwMDAxMDAyMSAwMDAwMCBuIAowMDAwMDA5NjAzIDAwMDAwIG4gCjAwMDAwMTEyODEgMDAwMDAgbiAKMDAwMDAwNTUzNCAwMDAwMCBuIAowMDAwMDA1NjgzIDAwMDAwIG4gCjAwMDAwMDYwOTcgMDAwMDAgbiAKMDAwMDAwNjQ3NyAwMDAwMCBuIAowMDAwMDA2NzgxIDAwMDAwIG4gCjAwMDAwMDcxMDMgMDAwMDAgbiAKMDAwMDAwNzI0NyAwMDAwMCBuIAowMDAwMDA3MzY2IDAwMDAwIG4gCjAwMDAwMDc2OTcgMDAwMDAgbiAKMDAwMDAwNzk4OCAwMDAwMCBuIAowMDAwMDA4MzAwIDAwMDAwIG4gCjAwMDAwMDg1MjMgMDAwMDAgbiAKMDAwMDAwODc0NyAwMDAwMCBuIAowMDAwMDA4OTgwIDAwMDAwIG4gCjAwMDAwMDkzODcgMDAwMDAgbiAKMDAwMDAwOTUxMyAwMDAwMCBuIAowMDAwMDExODM5IDAwMDAwIG4gCnRyYWlsZXIKPDwgL1NpemUgNDEgL1Jvb3QgMSAwIFIgL0luZm8gNDAgMCBSID4+CnN0YXJ0eHJlZgoxMTk5NgolJUVPRgo=",
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Log probabilites for plotting the target\n",
"x_plot = np.linspace(-4, 4, 200)\n",
"log_probs = [np.exp(log_prob(x)[0]) / (2 * np.pi) ** 0.5 for x in x_plot]\n",
"\n",
"# Plot samples and target\n",
"plt.figure(figsize=(6, 3))\n",
"\n",
"plt.hist(samples,\n",
" density=True,\n",
" bins=50,\n",
" color='gray',\n",
" alpha=0.5,\n",
" label='Samples')\n",
"plt.plot(x_plot,\n",
" log_probs,\n",
" color='black',\n",
" label='Normalised $f(x)$')\n",
"\n",
"# Plot formatting\n",
"plt.xlim([-4, 4])\n",
"plt.ylim([0, 0.5])\n",
"plt.xticks([])\n",
"plt.yticks([])\n",
"plt.xlabel('$x$', fontsize=18)\n",
"plt.ylabel('$f(x)~/~Z$', fontsize=18)\n",
"plt.legend()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Conclusions\n",
"\n",
"ARS is an efficient method for sampling log-concave univariate distributions. Although very effective for log-concave one-dimensional differentiable distributions, the algorithm presented here has two shortcomings. First, this algorithm requires gradients of the objective with respect to the input variable. These may be expensive to compute or perhaps even may not exist if $f$ is nowhere differentiable. For this, there exists a modified version of ARS{cite}`gilks1992adaptive` which builds the envelope in a way that does not require gradients. The present page presented the gradient-based method because this is nicer for illustrative purposes. Second, although many distributions of practical interest are log-concave, there are many others which are not. In this case, the ARS algorithm does not apply since the envelope is not guaranteed to entirely contain the probability distribution. For this, there exists an extension of ARS for non-log-concave distributions called the adaptive rejection Metropolis algorithm{cite}`gilks1995adaptive` (ARMS). ARMS also builds an envelope and uses it to propose samples, which are then accepted or rejected using a Metropolis-Hastings step to ensure that the samples are distributed according to the target. Note that in general, ARMS does not produce independent samples from the target, due to the Metropolis-Hastings accept/reject step. For log-concave functions, ARMS reduces to ARS."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## References\n",
"\n",
"```{bibliography} ./ars.bib\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"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.9.6"
}
},
"nbformat": 4,
"nbformat_minor": 4
}