Skip to content

Instantly share code, notes, and snippets.

@narrowlyapplicable

narrowlyapplicable/tfp-replica_exchange.ipynb

Last active October 5, 2019 14:52
Show Gist options
  • Save narrowlyapplicable/126ef083668ce6ea60681bd2025b047b to your computer and use it in GitHub Desktop.
Save narrowlyapplicable/126ef083668ce6ea60681bd2025b047b to your computer and use it in GitHub Desktop.
Download ZIP
Tensorflow Probability の ReplicaExchangeMC で <http://statmodeling.hatenablog.com/entry/stan-parallel-tempering> を再現しようとした例。局所解を抜け出すことはできず、元記事の結果を再現できなかった。（ただしより簡単な初期値では成功した"tfp-replica_excahnge_easy.ipynb"）
Display the source blob
Display the rendered blob
Raw
tfp-replica_exchange.ipynb
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Display the source blob
Display the rendered blob
Raw
tfp-replica_exchange_easy.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Tensorflow Probabilityによるレプリカ交換モンテカルロ試行 - より近い初期値の場合\n",
"- \"tfp-replica_exchange.ipynb\"では、`tfp.mcmc.ReplicaExchangeMC`を使用したが、元にした記事の結果を再現できなかった。\n",
"- MCMCの初期値を変え、少し簡単な問題においてMCMCの設定による違いを確認した。\n",
" - Stepsizeの自動調整なしのHMCでは、真値と全く異なるところに事後分布が構成された。\n",
" - ただしサポート付近にMCMCサンプルが集中しており、bijectorを設定していればもう少し改善した可能性はある。\n",
" - Stepsize調整なしでも、レプリカ交換では真値近くに到達した。"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import warnings\n",
"warnings.simplefilter('ignore')"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"import tensorflow as tf\n",
"import tensorflow_probability as tfp\n",
"from tensorflow_probability import edward2 as ed"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"tfd = tfp.distributions\n",
"tf.compat.v1.disable_eager_execution()\n",
"\n",
"np.random.seed(123)\n",
"plt.style.use('ggplot')"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"init_g = tf.global_variables_initializer()\n",
"init_l = tf.local_variables_initializer()"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"num_results = 10000\n",
"num_burnin_steps = 2000"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 1. Data\n",
"Statmodeling Memorandum様のレプリカ交換モンテカルロに関する記事の例題を取り上げる。\n",
"> <http://statmodeling.hatenablog.com/entry/stan-parallel-tempering> "
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"N = 50\n",
"b = 0.6\n",
"s_y = 0.4\n",
"X_data = np.linspace(0.1, 4*np.pi, N)\n",
"Y_data = np.sin(b*X_data) + np.random.normal(loc=0.0, scale=s_y, size=N)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3df3BU1fk/8Pe6KyCSxM2uJIQkaxJlWqpSaVKYTBExW8axdpq0jo0VBsUOYkwBUSNR0Uz5IJloyg+BAmNMGaQTtWNw7NBidygwlVgDIVTBQmIiqyYhJpvwQ0GT7P3+4TdbQ3ZDdnP3nnNP3q8ZZ9y7d7PPkyX73HvOc8+1aJqmgYiISDJXiA6AiIgoGBYoIiKSEgsUERFJiQWKiIikxAJFRERSYoEiIiIp2UQHAAAdHR3YtGkTuru7YbFY4Ha7ceeddw7YR9M0VFZW4siRIxg7diwKCgqQnp4uKGIiIoo2KQqU1WrF/PnzkZ6ejgsXLmDFihW4+eabkZycHNjnyJEjaGtrw4YNG9DQ0ICXX34Zzz///LB+fktLS9gxOZ1OdHR0hP06WamWD8CczEK1nFTLBxCfU1JSUtDtUgzx2e32wNnQVVddhcmTJ8Pn8w3Y59ChQ7j11lthsVgwZcoUfPnll+jq6hIRLhERGUCKAvVd7e3taG5uxvXXXz9gu8/ng9PpDDx2OByDihgREalDiiG+fhcvXkR5eTnuv/9+jB8/fsBzwVZkslgsQX+Ox+OBx+MBAJSWlg4obMNls9kiep2sVMsHYE5moVpOquUDyJuTNAWqt7cX5eXlmDVrFmbMmDHoeYfDMWCMtLOzE3a7PejPcrvdcLvdgceRjK2KHpPVm2r5AMzJLFTLSbV8APE5ST0HpWkatmzZgsmTJ+Ouu+4Kuk9mZiYOHDgATdNw8uRJjB8/PmSBIiIi85PiDOrEiRM4cOAAUlNT8cQTTwAA7r333kBFnzt3Lm655RbU1dVhyZIlGDNmDAoKCkSGTEREUSZFgfre976H119/fch9LBYLfvvb3xoUERERiSZFgSIiInPxer0oKytDW1sbEhMTUVRUhNTUVF3fgwWKiIjC4vV6kZ+fj1OnTgW21dXVoaqqStciJUWTBBERmUdZWdmA4gQAp06dQllZma7vwwJFRERhaWtrC7r99OnTur4PCxQREYUlMTEx6PaEhARd34cFioiIwlJUVASXyzVgm8vlQlFRka7vwyYJIiIKS2pqKqqqqlBWVobTp08jISGBXXxERCSH1NRUbNy4MarvwSE+IiKSEs+giIhGASMurNUbCxQRkeKMurBWbxziIyJSnFEX1uqNBYqISHFGXVirNxYoIiLFGXVhrd5YoIiIFGfUhbV6Y5MEEZHijLqwVm8sUBQ1ZmxrJVKVERfW6o0FiqLCrG2tRCQPaQrU5s2bUVdXh7i4OJSXlw96/tixYygrK8PEiRMBADNmzMDdd99tdJg0TEO1tZrtKI6IxJCmQN1222244447sGnTppD7fP/738eKFSsMjIoiZXRbK4cTidQjTYGaOnUq2tvbRYehPKO+yI1sa+VwIpGapClQw3Hy5Ek88cQTsNvtmD9/PlJSUoLu5/F44PF4AAClpaVwOp1hv5fNZovodbKy2Ww4d+4c7rvvPjQ1NQW2Hz16FLt370ZaWpqu77dmzRocPXp0wHulp6djzZo1uv1e+z+jxx57LOhw4vr167F9+3Zd3ssoqv27A9TLSbV8AHlzMk2BSktLw+bNmzFu3DjU1dXhhRdewIYNG4Lu63a74Xa7A487OjrCfj+n0xnR62TldDpRXFw8oGAAQFNTE4qLi3WfF4qJicHOnTsHtbXGxMTo9nvt/4wuLU79vF6v6T5D1f7dAerlpFo+gPickpKSgm43TYEaP3584P+nT5+OiooKnD17FrGxsQKjMhej54WMams161XyRDQ006wk0d3dDU3TAACNjY3w+/2IiYkRHFX0eb1eFBYW4u6770ZhYSG8Xm/EP0vVL3KzXiVPREOT5gxq3bp1OH78OM6dO4fFixfjnnvuQW9vLwBg7ty5eO+99/DOO+/AarVizJgxWLZsGSwWi+Coo0vvyf+ioiLU1dUN+HkqfJGb9Sp5IhqaRes/LVFYS0tL2K8RPSYLAIWFhaiurh60PS8vL+yhs/58+rv4VPgiH85nZLb2cxn+3elNtZxUywcQn5Pp56BGo2jMGZlxuZNIsf2cyNxMMwc1Gqk6Z2QUs96kjYi+xQIlMU7+j4xZb9JGJAs9m7QiwSE+iXHyf2R4BkoUORmGyFmgJDea5oz0pmrXIpERZFjwmQWKlMUzUKLIyTBEzgJFSuMZKFFkZBgiZ5MEERENIkOTFs+giIhoEBmGyFmgiIgoKNFD5BziIyIiKbFAERGRlFigiIhISixQREQkJRYoIiKSEgsUERFJiQWKiIikJM11UJs3b0ZdXR3i4uJQXl4+6HlN01BZWYkjR45g7NixKCgoQHp6uoBIiYjICNKcQd1222146qmnQj5/5MgRtLW1YcOGDVi0aBFefvllA6MjIlKT1+vFggULhN3zaSjSnEFNnToV7e3tIZ8/dOgQbr31VlgsFkyZMgVffvklurq6YLfbDYySiEgdMtzzaSjSFKjL8fl8cDqdgccOhwM+ny9ogfJ4PPB4PACA0tLSAa8bLpvNFtHrZBWtfJqbm1FSUoLW1lZMmjQJJSUlSEtL0/19glHtMwKYkxmolM9jjz0W9J5P69evx/bt2wVF9T+mKVCapg3aZrFYgu7rdrvhdrsDjzs6OsJ+P6fTGdHrZBWNfIIdfdXU1Bh29KXaZwQwJzNQKZ9Li1M/r9draI5JSUlBt0szB3U5DodjwC+ss7OTw3uCDXXHTSKSnwz3fBqKaQpUZmYmDhw4AE3TcPLkSYwfP54FSjAZ7rhJRJGT4Z5PQ5FmiG/dunU4fvw4zp07h8WLF+Oee+5Bb28vAGDu3Lm45ZZbUFdXhyVLlmDMmDEoKCgQHDHJfvRFREPrv+fT+vXr4fV6hdzzaSgWLdjkjmJaWlrCfo1K48yAcXNQLpeLc1AjwJzkp1o+gPicQs1BSXMGReYjwx03iUhdLFA0IqLvuElE6jJNkwQREY0uLFBERCQlFigiIpIS56AU5PV6UVZWhra2NiQmJqKoqEiZpVmIaPRggVJMqMUf9+zZg5iYmIh/5qUFT4VOPVXzIlIFC5RiQi0/VFJSEvQ+W5cj+2rHkVI1LyKVcA5KMaGWH2ptbY3o56m63p6qeRGphAVKMaGWH5o0aVJEP0/V9fZUzYtIJSxQigm1+GNJScllX+v1elFYWDjgzpqqrrenal5EKuEclGJCLT+UlpY25FpboeZk/vCHP6Curm7QensjWe14qOYEoxoXioqKdM+LiPTFAhUB2bu/Ill+KNSczKuvvqrrentDNScAMKxxgesIklnJ/v2jJxaoMKna/TXUnIye6+1drjkh1HPRWO+P6wiS2Vzu+0e14sUCFaahvmDN/GVn1JzMUIUw1J1f2LhAo1GwYjPU909RUZFyB88sUGFStfvLqDmZSAohGxdotAl1puRwOILuf/r0aSUPnlmgwqRq95dRczKXK4RDPXfpEeWaNWsiXh2DSGahik1fX1/Q/RMSEpQ8eJamQNXX16OyshJ+vx85OTnIzc0d8Py+ffuwY8cOxMfHAwDuuOMO5OTkGB6nyt1fRszJXK4Qhnou2BHl0aNHsXPnTtMOXxCFEqrYTJw4EVarNej3T6iLzM188CxFgfL7/aioqMAzzzwDh8OB4uJiZGZmIjk5ecB+2dnZePDBBwVF+S12f43cUIUw1HPBjiibmppMPXxBFEqokRqXy4VNmzYF/f5R8eBZigLV2NiIxMTEQKXPzs5GbW3toAIli6G+YFXropGFGYYv+NmTXoYqNqG+f1Q8eJaiQPl8vgGTfw6HAw0NDYP2+/e//42PPvoIkyZNwoIFC6S7hYSqLegykH3uj5896SnSYqPapRMWLVRvr4Fqampw9OhRLF68GABw4MABNDY2YuHChYF9zp07h3HjxuHKK6/EO++8g5qaGjz33HNBf57H44HH4wEAlJaW4ptvvgk7JpvNht7e3rBes2DBgsAFp9+Vn5+P7du3hx2DniLJRybNzc2488470dTUFNiWnp6O3bt3Iy0tTWBk39Lrszf75xSMajmplg8gPqcxY8YE3S7FGZTD4UBnZ2fgcWdnJ+x2+4B9vtut5Xa7sXPnzpA/z+12w+12Bx4PtcRPKE6nM+zXXTpH0s/r9UYUg54iyUcmMTEx2Llz54Ajyv4uPhny0uuzN/vnFIxqOamWDyA+p6SkpKDbpShQGRkZaG1tRXt7O+Lj43Hw4EEsWbJkwD5dXV2BonXo0CEp56ciHYbi3MXwXDp8IfqP6rtkH4IkMiMpCpTVasXChQuxevVq+P1+zJkzBykpKXjttdeQkZGBzMxM/O1vf8OhQ4dgtVoxYcIEFBQUiA57kEi6aDh3oQYVO6iIRJNiDiraWlpawn5NpEfn/WdDw53YLCwsRHV19aDteXl5uk52ynS2oRfZcgr3sw9Gtpz0oFpOquUDiM9J6iE+lYTbRWOG9unRKJJhV9U6qIhEY4ESjHMX8rncLUE4X0hkDBYowTh3IZ9Q66A999xzOHHiBOcLiQzCW74L1n9BXl5eHrKzs5GXl8cvPMFCDbteeiABDLyXFRHpi2dQEuDchVxCDbuGwvlCoujgGRTRJYqKiuByuQZsc7lcmD59etD9OV9IFB08gyK6RKh10AAMmoPifCFR9LBAEQURathVtdWiiWTGAkUUBs4XEhmHc1BERCQlFigiIpISCxQREUmJBYqIiKTEAkVERFJigSIiIimxzZxIsO/e2sPlcmHp0qW8tooILFBEQl16a4+amhrU1NRwweBRIpL7jo0mLFBEAoW6tUdZWRkvCFbcUPcdY5H6ljQFqr6+HpWVlfD7/cjJyUFubu6A53t6erBx40Y0NTUhJiYGy5Ytw8SJEwVFS6QP3lF59OLByeVJ0STh9/tRUVGBp556CmvXrsW7776Lzz77bMA+e/fuxdVXX42XXnoJP/vZz7Bz505B0RLph3dUHr14cHJ5UhSoxsZGJCYmIiEhATabDdnZ2aitrR2wz6FDh3DbbbcBAGbOnIkPP/wQmqYJiJZIP6Fu7cEV0tXHg5PLk2KIz+fzweFwBB47HA40NDSE3MdqtWL8+PE4d+4cYmNjB/08j8cDj8cDACgtLYXT6Qw7JpvNFtHrZKVaPoAaOTmdTuzZswclJSVobW3F5MmT8eyzzyItLU10aLpR4XP6Lr3yWbNmDY4ePYqmpqbAtvT0dKxZs8bw35esn5EUBSrYmZDFYgl7n35utxtutzvwuKOjI+yYnE5nRK+TlWr5AOrkFBMTg/LycgD/y0mFvPqp8jn10yufmJgY7Ny5c9DtW2JiYgz/fYn+jJKSkoJul6JAORwOdHZ2Bh53dnbCbrcH3cfhcKCvrw9fffUVJkyYYHSoRES64e1bhibFHFRGRgZaW1vR3t6O3t5eHDx4EJmZmQP2+dGPfoR9+/YBAN577z384Ac/CHkGJSOv14vCwkLcfffdKCwshNfrFR0SEZHUhjyD6u7uxjXXXBP1IKxWKxYuXIjVq1fD7/djzpw5SElJwWuvvYaMjAxkZmbi9ttvx8aNG/G73/0OEyZMwLJly6Iel154vQMRUfgs2hCtcA8++CAWLFiAW2+91ciYdNfS0hL2a/Qcky0sLER1dfWg7Xl5eSM6vQ/nKnTRY8zRwJzMQaWcvF4v1q9fj1OnTim18oPozyiiOajHHnsMW7duxbvvvouHHnoI8fHxUQlOddG43oFnZUTG4t+c8Yacg5o6dSpefPFFuFwuPPHEE/j73/+ODz/8cMB/dHnRuN5hqKvQiUh//Jsz3mW7+K688kr86le/wmeffYY///nPiImJCTxnsVjYgTIMRUVFqKurG/CPe6QXY/IqdCJj8W/OeJctUB988AG2bduGtLQ0vPTSS4iLizMiLqWkpqaiqqpq0PUOIxkW4FXoRMbi35zxhixQf/zjH1FfX48HHngAM2fONComJel9vUM0zsqIKDT+zRlvyALV09OD8vJyXhAroWiclRGNFpHch6n/b279+vXwer38mzPAkG3mqhDdZi4D1fIBmJNZyJZTsG48l8s17G482fLRg+icQrWZS7GSBBGRUdiNZx4sUEQ0qrAbzzxYoIhoVGE3nnmwQBHRqMKbRJqHFLfbkFFzczOKi4vD6vIhIvmxA9Y8WKCC8Hq9uO+++wbc6ZJrbhGpg/dhMgcO8QVRVlY2oDgB7PIh8+C9x0gVPIMKgl0+ZFZccVucSC7+paGxQAXBLh8yq6Gu8eGQVvTwwCA6OMQXRFFREdLT0wdsY5cPmQHP/sXgxb/RIfwM6vz581i7di2++OILXHvttXj00UeDrv3361//OnAk4nQ68eSTT0YtptTUVOzevRvFxcXs8iFTUfXsX/bhMx4YRIfwArVr1y7cdNNNyM3Nxa5du7Br1y7Mmzdv0H5jxozBCy+8YFhcaWlpHBIh01FxxW0zDJ+pemAgmvAhvtraWsyePRsAMHv2bNTW1gqOiMi8+q/xycvLQ3Z2NvLy8qT6Io+EGYbPePFvdAg/gzpz5gzsdjsAwG634+zZs0H36+npwYoVK2C1WvGLX/wCP/7xj0P+TI/HA4/HAwAoLS2F0+kMOy6bzRbR62SlWj4AcwrF6XSiqqpKp4hGbqQ5+Xy+kNtFfP7B8nE6ndizZw9KSkrQ2tqKSZMmoaSkBGlpaYbHFwlZ/5YMKVCrVq1Cd3f3oO35+fnD/hmbN29GfHw8Tp8+jd///vdITU0NeVrtdrvhdrsDjyNZRl708vN6Uy0fgDmZxUhzio+PD7ldxO8qVD4xMTEoLy8fsM0sn6Xof3ehbrdhSIFauXJlyOfi4uLQ1dUFu92Orq4uxMbGBt2v/x9pQkICpk6dik8++SRkgSIidag4r0bDI3wOKjMzE/v37wcA7N+/H1lZWYP2OX/+PHp6egAAZ8+exYkTJ5CcnGxonEQkhorzajQ8wuegcnNzsXbtWuzduxdOpxPLly8HAHz88cf4xz/+gcWLF+Pzzz/Htm3bcMUVV8Dv9yM3N5cFimgU4dp5oxNv+R6C6DFZvamWD2CunIZ7HY+Zchou1XJSLR9AfE5C56CIRjMzXMejKhku8JUhBrNigSKKMq6PJ4YMBwYyxGBmwpskiFTHZXDEkOECXxliMDMWKKIo4zI4YshwYCBDDGbGAkUUZVwGRwwZDgxkiMHMWKCIoozX8YgR6YGBnnck5sHJyLBJgsgAvI7HeP0HBmVlZcO+bY7eTQ2RxED/wwJFRMoK98AgGh2XPDiJHIf4iIj+PzY1yIVnUERkanpeCMumBrmwQBGRaek9Z8SV0+XCIT4iMi29L4Rlx6VceAZFRKYVjTkjNjXIg2dQRGRanDNSGwsUkQnpeTGpmfFCWLVxiI/IZLhC9v/wQli1sUARmQxv3zEQ54zUJbxA1dTU4I033sDnn3+O559/HhkZGUH3q6+vR2VlJfx+P3JycpCbm2twpERyiEZjAG+qRzISXqBSUlLw+OOPY9u2bSH38fv9qKiowDPPPAOHw4Hi4mJkZmYiOTnZwEiJ5KB3YwCHDElWwpskkpOTQ96Pvl9jYyMSExORkJAAm82G7Oxs1NbWGhQhkVz0bgzgTfVIVsLPoIbD5/PB4XAEHjscDjQ0NAiMiCj6Qg276d0YwPXnSFaGFKhVq1ahu7t70Pb8/HxkZWVd9vWapg3aZrFYQu7v8Xjg8XgAAKWlpXA6nWFE+y2bzRbR62SlWj6A2jk1NzfjvvvuQ1NTU+C5o0ePYvfu3UhLS4PT6URVVZUu7+lyuVBTUzNoe2pqqi6/X9U+J9XyAeTNyZACtXLlyhG93uFwoLOzM/C4s7MTdrs95P5utxtutzvwuKOjI+z3dDqdEb1OVqrlA6idU3Fx8YDiBABNTU0oLi7WvWNt6dKlqKmpGbT+3NKlS3X5/ar2OamWDyA+p1DTPMLnoIYjIyMDra2taG9vR29vLw4ePIjMzEzRYRFFjZHDbpGuP8eLhSnahM9Bvf/++3jllVdw9uxZlJaW4rrrrsPTTz8Nn8+HrVu3ori4GFarFQsXLsTq1avh9/sxZ84cpKSkiA6dKGqMXsIn3GuJ2PlHRrBowSZ4FNPS0hL2a0Sf8upNtXwAtXMKVgBcLpc0BaCwsBDV1dWDtufl5Q0qdKp9TqrlA4jPKdQQn/AzKCIaTPYlfNj5R0ZggSKSlMxL+HAVcTKCKZokiEguXEWcjMAzKCIKW7SGILkmIH0XCxQRRUTvIUh2BtKlOMRHRFLgmoB0KRYoIpICOwPpUixQRCQFdgbSpVigiEgK7AykS7FJgoiGZFRnnewXJ5PxWKCIKCSjO+tkvjiZjMchPiIKiZ11JBILFBGFxM46EokFiohCYmcdicQCRUQhsbOORGKTBNEoEUk3HjvrSCQWKKJRYCTdeOysI1E4xEc0CrAbj8xI+BlUTU0N3njjDXz++ed4/vnnkZGREXS/Rx55BOPGjcMVV1wBq9WK0tJSgyMlMi9245EZCS9QKSkpePzxx7Ft27bL7vvcc88hNjbWgKiI1MJuPDIj4UN8ycnJSEpKEh0GkdLYjUdmJPwMKhyrV68GAPz0pz+F2+0OuZ/H44HH4wEAlJaWwul0hv1eNpstotfJSrV8AOYUDqfTiT179qCkpAStra2YNGkSSkpKkJaWpvt7XUq1z0m1fAB5c7JomqZF+01WrVqF7u7uQdvz8/ORlZUFACgpKcH8+fNDzkH5fD7Ex8fjzJkz+L//+z888MADmDp16rDev6WlJeyYnU4nOjo6wn6drFTLB2BOZqFaTqrlA4jPKdQomiFnUCtXrhzxz4iPjwcAxMXFISsrC42NjcMuUEREZD7C56CG4+LFi7hw4ULg///zn//wQkEiIsUJn4N6//338corr+Ds2bMoLS3Fddddh6effho+nw9bt25FcXExzpw5gxdffBEA0NfXh5/85Cf44Q9/KDhyIiKKJkPmoETjHJR6+QDMySxUy0m1fADxOYWagzLFEB8REY0+wof4iGh0MeoW8mR+LFBEZBijbyFP5sYhPiIyDBetpXCwQBGRYbhoLYWDBYqIDMNFaykcLFBEZBguWkvhYJMEERmGt5CncLBAEZGheAt5Gi4O8RERkZRYoIiISEosUEREJCUWKCIikhILFBERSYkFioiIpMQCRUREUmKBIiIiKQm/UHfHjh04fPgwbDYbEhISUFBQgKuvvnrQfvX19aisrITf70dOTg5yc3MFREtEREYRfgZ18803o7y8HC+++CImTZqE6urqQfv4/X5UVFTgqaeewtq1a/Huu+/is88+ExAtEREZRXiBmjZtGqxWKwBgypQp8Pl8g/ZpbGxEYmIiEhISYLPZkJ2djdraWqNDJSIiAwkf4vuuvXv3Ijs7e9B2n88Hh8MReOxwONDQ0BDy53g8Hng8HgBAaWkpnE5n2LHYbLaIXicr1fIBmJNZqJaTavkA8uZkSIFatWoVuru7B23Pz89HVlYWAODNN9+E1WrFrFmzBu2nadqgbRaLJeT7ud1uuN3uwOOOjo6wY3Y6nRG9Tlaq5QMwJ7NQLSfV8gHE55SUlBR0uyEFauXKlUM+v2/fPhw+fBjPPvts0MLjcDjQ2dkZeNzZ2Qm73a57nEREJA/hc1D19fV466238OSTT2Ls2LFB98nIyEBrayva29vR29uLgwcPIjMz0+BIiYjISMLnoCoqKtDb24tVq1YBAG644QYsWrQIPp8PW7duRXFxMaxWKxYuXIjVq1fD7/djzpw5SElJERw5EYXi9XpRVlaGtrY2uFwuLF26lDclpLBZtGATPIppaWkJ+zWix2T1plo+AHOSldfrRX5+Pk6dOhXY5nK5UFVVpUSRUuEzupTonELNQQkf4iMitZSVlQ0oTgBw6tQplJWVCYqIzIoFioh01dbWFnT76dOnDY6EzI4Fioh0lZiYGHR7QkKCwZGQ2bFAEZGuioqK4HK5BmxzuVwoKioSFBGZlfAuPiJSS2pqKqqqqlBWVobTp08jNTWVXXwUERYoItJdamoqNm7cCEB8hxiZF4f4iIhISixQREQkJRYoIiKSEgsUERFJiQWKiIikxAJFRERSGhWLxRIRkfnwDCqEFStWiA5BV6rlAzAns1AtJ9XyAeTNiQWKiIikxAJFRERSspaUlJSIDkJW6enpokPQlWr5AMzJLFTLSbV8ADlzYpMEERFJiUN8REQkJa5mfon6+npUVlbC7/cjJycHubm5okMakY6ODmzatAnd3d2wWCxwu9248847RYc1Yn6/HytWrEB8fLy0HUjh+PLLL7FlyxZ8+umnsFgsePjhhzFlyhTRYY3IX//6V+zduxcWiwUpKSkoKCjAmDFjRIcVls2bN6Ourg5xcXEoLy8HAJw/fx5r167FF198gWuvvRaPPvooJkyYIDjS4QuW044dO3D48GHYbDYkJCSgoKAAV199teBIAWgU0NfXpxUWFmptbW1aT0+P9vjjj2uffvqp6LBGxOfzaR9//LGmaZr21VdfaUuWLDF9TpqmaW+//ba2bt06bc2aNaJD0cVLL72keTweTdM0raenRzt//rzgiEams7NTKygo0L7++mtN0zStvLxc++c//yk2qAgcO3ZM+/jjj7Xly5cHtu3YsUOrrq7WNE3TqqurtR07dogKLyLBcqqvr9d6e3s1Tfs2P1ly4hDfdzQ2NiIxMREJCQmw2WzIzs5GbW2t6LBGxG63ByY/r7rqKkyePBk+n09wVCPT2dmJuro65OTkiA5FF1999RU++ugj3H777QAAm80mx9HrCPn9fnzzzTfo6+vDN998A7vdLjqksE2dOnXQ2VFtbS1mz54NAJg9e7bpviOC5TRt2jRYrVYAwJQpU6T5juAQ33f4fD44HI7AY4fDgYaGBoER6au9vR3Nzc24/vrrRYcyIn/6058wb948XLhwQXQoumhvb0dsbCw2b96MU6dOIT09Hffffz/GjRsnOrSIxcfH4+c//wjvaBoAAAPtSURBVDkefvhhjBkzBtOmTcO0adNEh6WLM2fOBIqt3W7H2bNnBUekr7179yI7O1t0GADYJDGAFqSh0WKxCIhEfxcvXkR5eTnuv/9+jB8/XnQ4ETt8+DDi4uKkbImNVF9fH5qbmzF37lyUlZVh7Nix2LVrl+iwRuT8+fOora3Fpk2bsHXrVly8eBEHDhwQHRZdxptvvgmr1YpZs2aJDgUAC9QADocDnZ2dgcednZ2mHJa4VG9vL8rLyzFr1izMmDFDdDgjcuLECRw6dAiPPPII1q1bhw8//BAbNmwQHdaIOBwOOBwO3HDDDQCAmTNnorm5WXBUI/PBBx9g4sSJiI2Nhc1mw4wZM3Dy5EnRYekiLi4OXV1dAICuri7ExsYKjkgf+/btw+HDh7FkyRJpDsxZoL4jIyMDra2taG9vR29vLw4ePIjMzEzRYY2IpmnYsmULJk+ejLvuukt0OCP2m9/8Blu2bMGmTZuwbNky3HjjjViyZInosEbkmmuugcPhQEtLC4Bvv9yTk5MFRzUyTqcTDQ0N+Prrr6FpGj744ANMnjxZdFi6yMzMxP79+wEA+/fvR1ZWluCIRq6+vh5vvfUWnnzySYwdO1Z0OAG8UPcSdXV12L59O/x+P+bMmYNf/vKXokMakf/+97949tlnkZqaGjgquvfeezF9+nTBkY3csWPH8PbbbyvRZv7JJ59gy5Yt6O3txcSJE1FQUGCq1uVgXn/9dRw8eBBWqxXXXXcdFi9ejCuvvFJ0WGFZt24djh8/jnPnziEuLg733HMPsrKysHbtWnR0dMDpdGL58uWm+qyC5VRdXY3e3t5AHjfccAMWLVokOFIWKCIikhSH+IiISEosUEREJCUWKCIikhILFBERSYkFioiIpMQCRUREUmKBIpLcxYsX8cgjj+Bf//pXYNuFCxfw8MMP47333hMYGVF0sUARSW7cuHFYtGgRKisrAwuTvvrqq8jIyMDMmTMFR0cUPSxQRCYwbdo0TJ8+Ha+88gqOHTuGmpoaPPjgg6LDIooqriRBZBLnz5/H8uXL0dfXh3nz5mHOnDmiQyKKKp5BEZnEhAkTkJKSgq+//tr0q9ITDQcLFJFJHDhwAO3t7bjpppvw6quvig6HKOpYoIhM4MyZM9i+fTseeughLFq0CDU1NTh+/LjosIiiigWKyAQqKiqQlZWFG2+8EXa7HfPmzcPWrVvR09MjOjSiqGGBIpLc+++/jxMnTmD+/PmBbTk5OXA4HPjLX/4iMDKi6GIXHxERSYlnUEREJCUWKCIikhILFBERSYkFioiIpMQCRUREUmKBIiIiKbFAERGRlFigiIhISixQREQkpf8H8ERoeoMpfXsAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, ax = plt.subplots()\n",
"ax.scatter(X_data, Y_data, marker=\"o\", color=\"k\")\n",
"ax.set_xlabel('X')\n",
"ax.set_ylabel('Y')\n",
"fig.tight_layout()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 2. Model\n",
"元記事でのモデルは、Stanを用いて下記のように定義されている。\n",
"\n",
"```\n",
"data {\n",
" int<lower=1> N;\n",
" vector[N] Y;\n",
" vector[N] X;\n",
" real<lower=0> Inv_T;\n",
"}\n",
"\n",
"parameters {\n",
" real<lower=0> b;\n",
" real<lower=0> s_y;\n",
"}\n",
"\n",
"transformed parameters {\n",
" real E;\n",
" {\n",
" vector[N] mu;\n",
" for (n in 1:N)\n",
" mu[n] <- sin(b * X[n]);\n",
" E <- 0;\n",
" E <- E - normal_log(b, 0, 50);\n",
" E <- E - student_t_log(s_y, 4, 0, 5);\n",
" E <- E - normal_log(Y, mu, s_y);\n",
" }\n",
"}\n",
"\n",
"model {\n",
" increment_log_prob(-Inv_T * E);\n",
"}\n",
"```\n"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"def joint_log_prob(x_data, y_data, b, s_y):\n",
" rv_b = tfd.TruncatedNormal(loc=0.0, scale=50.0, low=0.0, high=1e100, name=\"rv_b\")\n",
" #rv_sy = tfd.Uniform(low=0.0, high=1e100, name=\"rv_sy\")\n",
" rv_sy = tfd.TruncatedNormal(loc=0.0, scale=5.0, low=0.0, high=1e100, name=\"rv_sy\")\n",
" #StudentT(df=4, loc=0, scale=5, name=\"rv_sy\")\n",
" \n",
" mu = tf.sin(b*x_data)\n",
" rv_obs = tfd.Normal(loc=mu, scale=s_y, name=\"rv_obs\") \n",
" \n",
" return(\n",
" rv_b.log_prob(b)\n",
" + rv_sy.log_prob(s_y)\n",
" + tf.reduce_sum(rv_obs.log_prob(y_data))\n",
" )"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- TF-PにはStanの`<lower=0>`のような値域制限が（私の知る限りでは）ない。そのため事前分布でサポートに制限を加えている。\n",
" - bに関しては、一様事前分布の範囲を指定した。\n",
" - s_yに関しては切断t分布がtfdに無いため、切断正規分布で代用した。\n",
" - 事前分布によらず値域を制限する法は不明。`tf.gather(b, 条件)`だろうか？\n",
"- 逆温度指定の部分は`tfp.mcmc.ReplicaExchangeMC`に任せるため、モデルには明記しない。"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"def unnormalized_log_posterior(b, s_y):\n",
" return joint_log_prob(X_data, Y_data, b, s_y)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 3. HMC\n",
"### 3.1. Define the Kernel\n",
"まず最も単純なHMCカーネルを定義する。\n",
"- step_size調整無し\n",
"- bijectorによる効率化無し。"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"kernel1 = tfp.mcmc.HamiltonianMonteCarlo(\n",
" target_log_prob_fn=unnormalized_log_posterior,\n",
" step_size=0.01,\n",
" num_leapfrog_steps=3 #10\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 3.2. Inference"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- MCMCの初期値を変更"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"initial_state = [0.2, 1.0]\n",
"#[0.5, 2.0] # 真値から少しずらした初期値. 優しめの設定 -> 成功\n",
"#[24.17, 0.4048] #b=24.17, s_y=0.4048. かなり厳しめの局所解 -> 失敗\n",
"#[0.6, 0.4] #b=0.6, s_y=0.4. 真値"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### sample chain"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"WARNING:tensorflow:From /anaconda3/lib/python3.7/site-packages/tensorflow_probability/python/internal/special_math.py:154: add_dispatch_support.<locals>.wrapper (from tensorflow.python.ops.array_ops) is deprecated and will be removed in a future version.\n",
"Instructions for updating:\n",
"Use tf.where in 2.0, which has the same broadcast rule as np.where\n"
]
}
],
"source": [
"states, kernel_results = tfp.mcmc.sample_chain(\n",
" num_results=num_results, \n",
" num_burnin_steps=num_burnin_steps,\n",
" kernel=kernel1,\n",
" current_state=initial_state,\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### run"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"with tf.Session() as sess:\n",
" sess.run(init_g)\n",
" sess.run(init_l)\n",
" [states_, kernel_results_ ]= sess.run([states, kernel_results])"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"acceptance rate : 48.8 %\n"
]
}
],
"source": [
"try:\n",
" print(f'acceptance rate : {kernel_results_.inner_results.is_accepted.mean()*100} %')\n",
"except AttributeError:\n",
" print(f'acceptance rate : {kernel_results_.is_accepted.mean()*100} %')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 3.3. Result"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"#states_"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAfAAAADQCAYAAAD4dzNkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de1RU5foH8O8eruIAwQxKoBaIZpJ5Q0VNIR3NI2XmScvS8kLmojQkPV6Ol1ZmUkp4w4VmWZ08K3VlmJ2VnkOoaKZioCEeFdQ8WRjCgBeUuL2/P/gxcRmYGRhmZs98P2u5lvPO3pvnnZm9n/2+797vloQQAkRERCQrCmsHQERERKZjAiciIpIhJnAiIiIZYgInIiKSISZwIiIiGWICJyIikiEmcCIiIhliAicisnPTpk2DRqOxagxLlixBx44dIUkSPvnkE6vGYi+YwKkROe3siYmJePjhhy0XGJEFaDQaTJs2zWzbW79+PXbv3m227ZnqxIkTWL16NbZu3Yr8/Hw899xzVovFnjhbOwBqPY1Gg06dOpntrHb9+vWorq42y7ZaonZnT0lJwaBBg+Dt7d3kspmZmejTp48FoyOSj/Lycri6uja7D5m6rZbIzc2FQqHA008/3eo46E9sgZNOeXk5AMDb2xs+Pj5m2VZL1N3Z/f390a5duyaXzczMRN++fVv8t4hMFRkZiRkzZmDRokVQq9Xw8vJCdHQ07t27p1umoqICixYtQmBgIFxdXdGzZ0/885//1L1/9OhRDB06FJ6envD09ETv3r1x4MABADU9YN999x0+/fRTSJIESZJw6NAh3bobN25Ejx494O7ujm7dumHVqlWorKzUxTZz5kwsW7YM999/PwIDA3XbrNurZii+5ralT3PbmzZtGqZOnYrq6mpdfZqTnZ2NcePG4f7779ctX/vvo48+anZdhyOoTUVERIjp06eLhQsXCpVKJTw9PcXMmTPF3bt3dcuUl5eLhQsXioCAAOHi4iIefvhhsWPHDt37R44cEUOGDBFKpVIolUrx6KOPiv379wshhHj55ZcFgHr/Dh48qFt3w4YN4qGHHhJubm4iJCREvPPOO6KiokIX24wZM8TSpUuFv7+/UKvVum2OHDnS6Pia25Y+zW1PX32aUlpaKpycnERycrJ4+umnRfv27UVQUFCj2IjMKSIiQnh6eoro6Ghx7tw58fXXXws/Pz8xZ84c3TLz588Xvr6+YteuXeLChQti1apVQpIkkZqaKiorK4WPj4+YN2+euHjxorh48aLYs2ePSE9PF0IIUVJSIoYNGyYmTZok8vPzRX5+vvjjjz+EEEKsWLFCdOnSRezZs0dcvnxZ/Otf/xKdO3cWS5cu1cWmVCrFq6++KnJycsRPP/0khGi8TzcXX9166tuWPs1tr6SkRKxbt044OTnp6tOU3Nxc4eXlJV544QVx+vRpcerUKdG3b1/h5+cn/vGPfzS7riNiAm9j3NkbM9fOfuzYMQFAdOvWTXzzzTciLy9PzJkzRzg7O4uLFy+24NsiMiwiIkI88MADorKyUle2ZcsW4erqKu7cuSNKS0uFq6urSEpKqrfe+PHjxeOPPy60Wm2jE+2GRo4cKV5++eV6ZaWlpaJdu3bi22+/rVf+6aefCm9vb11s3bp1E1VVVfWWqbtPG4qvbj31bashY7a3fft24eTk1Ox2hBAiKipKREREiOrqal3ZF198IQCIoqIig+s7GibwNsadvT5z7uxJSUlCkiSRmZmpKysvLxft2rUT27ZtM7g+UUtERESIiRMn1is7e/asACDOnDkjzpw5IwCI7OzsesskJiaKDh06CCGEiI6OFq6urmLMmDFi9erV4vz58/WW1bdPnzx5UgAQHh4eon379rp/7u7uAoAoKCgQERER4vnnn28Uc9192pj4auupb1sNGbM9Y/ZprVYrXFxcxJdfflmvfN++fQKA0Gq1BmNxNBwDt4CBAwfCyclJ93ro0KEoLy/HpUuXkJeXh/LycgwfPrzeOhEREcjJyYGPjw+io6PxxBNP4C9/+Qvi4+Nx4cIFg38zJycH9+7dw1//+lcolUrdv1dffRU3b97EjRs3AAD9+/eHQtH0z8BQfHUZ2pap2zMkMzMTo0aNqjcG7uLiAjc3N6tehEeOR+h5KnPDsV4hhK7sww8/xI8//ohRo0bh8OHDeOSRR7Bly5Zm/0btb3r37t04ffq07l92djZyc3Ph6+sLAGjfvr1RMTcXXy1jt2Xs9pqTlZWFiooK9OvXr175yZMnERISAh8fH2RnZ2PQoEG6986ePYshQ4bo/fwdARO4FXBnb/3ODtQk8P79+9cry83NRUlJSaNyInPKyMhAVVWV7vUPP/wAV1dXdO3aFSEhIXBzc8Phw4frrZOeno7Q0FDd60ceeQRxcXH49ttvMXPmTGzdulX3nqura73tA0BoaCjc3d1x+fJlhISENPpXt5HQHGPjM5a5tldb39LSUl3ZzZs38fHHH+tuqQsNDcWVK1dQUVEBAIiLi8OaNWtMPnbYC95GZgG1O3vtDlZ3Z5ckSffjr/tj17ez1+7ws2fPxtatW/Hqq68CMLyzjx07tsWx1905m4vP0tsrLy9HTk4ORo0aVa88ISEB/fv3b3QWT2RORUVFeO211/DGG2/g8uXLWLZsGV555RXdSezcuXOxbNky+Pn5oU+fPti9ezf27t2L//znP8jLy8OHH36Ip556Cp07d8Zvv/2GI0eO1PvNBgUF4eDBg7h06RK8vb3h7e0NpVKJJUuWYMmSJQCAUaNGobKyEtnZ2cjKysJ7771nVOweHh7Nxmcqc21v4MCB8PLywsKFCxEfH4+SkhL87W9/Q6dOnbBgwQIAgEKhQK9evZCdnY1r167Bx8cHQ4cONTlme8EEbgHc2c2/vbNnz6K8vBx79uzBmDFjEBwcjOTkZHz22Wc4evSoyXERmeLZZ5+Fp6cnHnvsMZSXl2PixIl4//33de+vWrUKCoUCsbGxuHHjBkJCQvD5559j5MiRyM/PR25uLp5//nncuHEDKpUKUVFRWLt2rW79N998E9nZ2ejduzdKS0tx8OBBREZGYtmyZQgICMDGjRsxf/58tGvXDt27dzd50pfm4msJc2zP29sbe/bswbx589CvXz+o1WpMnDgRK1eurHf/eXh4OI4dO4atW7fi66+/blG8dsNqo+8OovY2storr5VKpZg+fbooLS3VLdPcbVW//fabeOaZZ0RgYKBwdXUV999/v4iOjhYlJSW69S9duiSGDRsm2rdv3+iCt23btonevXsLNzc3cd9994mBAweKzZs362KbOXNmo5hbehuZvm3pY2h7xlzwsm3bNhEQECD+/e9/i65duwo3NzcRERHR6EIaInMz5bdO5vf1118LX19fsXDhQmuHYnWSEA46+m8hkZGRCAkJwbZt26wdChGZAfdp68rLy8PQoUORm5sLLy8va4djVbyIjYiIZGPDhg149913HT55AxwDJyIySd1pTclyLl26hKioKAwZMgQzZsywdjg2gV3oREREMsQudCIiIhliFzqRAysvL8eKFStQWVmJqqoqhIeHY9KkSSgoKMC6detw584dBAUFYc6cOXB2dkZFRQU2bdqEy5cvw9PTE7GxsejQoYO1q0HkkNiFTuTAhBD4448/4O7ujsrKSixfvhzTpk3DN998g0GDBmHo0KHYunUrHnzwQYwePRoHDhzA1atXMWvWLHz//fc4efIk5s2bZ+1qEDkkm2mB//bbb82+r1arUVhYaKFozE/O8cs5dsD+4g8ICDDbtiVJgru7O4CaqSyrqqogSRJycnLwxhtvAKi5bWr37t0YPXo0Tp06hYkTJwKomVDj448/NmoaXEP7tzXI/XfRFNZLXlqzf9tMAici66iursbChQtx/fp1PPHEE+jYsSM8PDx0U//6+vpCq9UCALRaLVQqFQDAyckJHh4euH37dqNbelJTU5GamgoAiI+Ph1qttmCNjOPs7GyTcbUW6yUvrakXEziRg1MoFFizZg1KS0uxdu1a/Prrr00uq2/ETV/rW6PRQKPR6F7bYsvJUVp09sJR6mVKC5xXoRMRgJqnyfXs2RO5ubm4e/eu7gE5Wq1W9/Q6lUqFoqIiADVd7nfv3oVSqbRazESOjAmcyIHdunVL9/jG8vJyZGdnIzAwEKGhoTh+/DiAmolLwsLCANQ88712IpPjx48jNDTUYR/lSGRt7EIncmDFxcVISkpCdXU1hBAYPHgw+vfvj06dOmHdunX44osvEBQUhBEjRgAARowYgU2bNmHOnDlQKpWIjY21cg2IHJdDJvB9O0sAAE89d5+VIyGyrgceeKDeYzBrdezYEatXr25U7urqiri4OEuERkbgscyxsQudiIhIhpjAiYiIZIgJnIiISIaYwImIiGTIIRL4vp0luos9iIiI7IFDJHAiIiJ7wwROREQkQ3aTwNlNTkREjsRuEjgREZEjkU0C356UxxY2ERHR/5NNAiciIqI/MYETEckIeyKpFhM4ERGRDDGBExERyZDdJnB2MxERkT2z2wROROSoOC+GY2ACJyIikiFnawdARNZTWFiIpKQklJSUQJIkaDQajB07Frt27cJ3330HLy8vAMDkyZPRr18/AMBXX32FtLQ0KBQKTJ8+HX369LFmFaiFalvoTz13n5UjoZZiAidyYE5OTpg6dSqCg4Nx7949LFq0CI8++igAICoqCuPGjau3/LVr13Ds2DF88MEHKC4uxsqVK7F+/XooFOzMs2X7dpYwUdsh7nVN4BgSOQIfHx8EBwcDANq1a4fAwEBotdoml8/IyMCQIUPg4uKCDh06wN/fH3l5eZYKl4jqYAuciAAABQUFuHLlCkJCQnD+/HkcOHAA6enpCA4OxksvvQSlUgmtVotu3brp1vH19dWb8FNTU5GamgoAiI+Ph1qttlg9jOXs7GyTcRlWUifumkZG3Xo4O/95WK+7XOO6Nl7Xlsn3+2pea+rFBE5EKCsrQ0JCAqZNmwYPDw+MHj0azz77LABg586d+OyzzxATEwMhhFHb02g00Gg0uteFhYVtEndrqNVqm4zLGA3jrvu6bjKoW95UXeXyGcj5+2pOw3oFBAQYvS670IkcXGVlJRISEjBs2DAMGjQIAHDfffdBoVBAoVBg5MiRuHTpEgBApVKhqKhIt65Wq4Wvr69V4iZydEzgRA5MCIHk5GQEBgbiySef1JUXFxfr/n/y5El07twZABAWFoZjx46hoqICBQUFyM/PR0hIiMXjthctvc6G1+gQYAdd6PwRE7XchQsXkJ6eji5dumDBggUAam4Z+/777/Hzzz9DkiT4+flh1qxZAIDOnTtj8ODBiIuLg0KhwMyZM3kFOpGVyD6BE1HL9ejRA7t27WpUXnvPtz4TJkzAhAkT2jIsIjKCwQReXl6OFStWoLKyElVVVQgPD8ekSZNQUFCAdevW4c6dOwgKCsKcOXPg7OyMiooKbNq0CZcvX4anpydiY2PRoUMHS9TFZJzIgIiI5Mpg35eLiwtWrFiBNWvW4P3338fp06dx8eJFfP7554iKisKGDRvQvn17pKWlAQDS0tLQvn17bNy4EVFRUdixY0ebBN4WY0Dsjiei1jB0XGqrsWuOiTsmgwlckiS4u7sDAKqqqlBVVQVJkpCTk4Pw8HAAQGRkJDIyMgAAp06dQmRkJAAgPDwcZ8+eNfrWE3PgD5mIiByBUWPg1dXVWLhwIa5fv44nnngCHTt2hIeHB5ycnADUn8xBq9VCpVIBqJmm0cPDA7dv39bNqVzL9IkeGk46YDhJNz3ZQUmTyzU3QUJryHkSAjnHDjB+IrJPRiVwhUKBNWvWoLS0FGvXrsWvv/7a5LL6WtuSJDUqa+lED6bcyN/cZAdNlRu7jqnkPAmBnGMH7C9+UyZ6IPnhtTlkLJPu/2jfvj169uyJ3Nxc3L17F1VVVQDqT+ZQd6KHqqoq3L17F0ql0sxhExE5Bg4LUlMMJvBbt26htLQUQM0V6dnZ2QgMDERoaCiOHz8OADh06BDCwsIAAP3798ehQ4cAAMePH0doaKjeFjgRERmPiZwaMtiFXlxcjKSkJFRXV0MIgcGDB6N///7o1KkT1q1bhy+++AJBQUEYMWIEAGDEiBHYtGkT5syZA6VSidjY2DavBBERNcaEb98MJvAHHngA77//fqPyjh07YvXq1Y3KXV1dERcXZ57oiIiISC/ZzcTGM0oiIiIZJnAiIrljQ4TMgU8hICIikiEmcCIiIhlyqC50dlsREZG9YAuciIhIhuy6Bc4WNxER2Su2wImIiGTIrlvgRNS8wsJCJCUloaSkBJIkQaPRYOzYsbhz5w4SExNx48YN+Pn5Yd68eVAqlRBCYPv27cjKyoKbmxtiYmIQHBxs7WrYHD6QhCyBLXAiB+bk5ISpU6ciMTERq1atwoEDB3Dt2jWkpKSgV69e2LBhA3r16oWUlBQAQFZWFq5fv44NGzZg1qxZ2LZtm5VrYDs4ZEeWxgRuIj5QgOyJj4+PrgXdrl07BAYGQqvVIiMjAxEREQCAiIgIZGRkAABOnTqF4cOHQ5IkdO/eHaWlpSguLrZa/JbU1vu+pY4tPIbZD3aht9K+nSXsJiO7UFBQgCtXriAkJAQ3b96Ej48PgJokf+vWLQA1jw5Wq9W6dVQqFbRarW7ZWqmpqUhNTQUAxMfH11vHVjg7O5sYV03Sa3qdxkmxNlFOfy2kyWX+3F6JUeXNUavVcHbWf1g3/Hdsm+nflzy0pl5M4ESEsrIyJCQkYNq0afDw8GhyOSFEozJ9jwvWaDTQaDS614WFheYJ1IzUanWL4jL3Ok2919K/01QyMOffsYaWfl+2rmG9AgICjF6XCRwcuyLHVllZiYSEBAwbNgyDBg0CAHh7e6O4uBg+Pj4oLi6Gl5cXgJoWd92DTVFRUaPWNxFZBsfAiRyYEALJyckIDAzEk08+qSsPCwvD4cOHAQCHDx/GgAEDdOXp6ekQQuDixYvw8PBgApcJjn3bH7bAiRzYhQsXkJ6eji5dumDBggUAgMmTJ2P8+PFITExEWloa1Go14uLiAAB9+/ZFZmYm5s6dC1dXV8TExFgzfJtgTFKUa+Lk7XC2jQmcyIH16NEDu3bt0vve8uXLG5VJkoTo6Oi2DotsFBO6bWEXOhERkQwxgRMRObCmxsbl2u3vSNiF3gB/tETkyHgMlA8mcCIiE9jyOHBNbEzAjoJd6ERERDLEFrgBrTnbtuUzdSKyDeyyppZiAiciaoG2TrxM7GQIEzgREfGEQYY4Bk5ERCRDTOBEREQyxAROREQkQ0zgREREMsQETkREJENM4ERERDJk8DaywsJCJCUloaSkBJIkQaPRYOzYsbhz5w4SExNx48YN+Pn5Yd68eVAqlRBCYPv27cjKyoKbmxtiYmIQHBxsibpYHSduISIiSzGYwJ2cnDB16lQEBwfj3r17WLRoER599FEcOnQIvXr1wvjx45GSkoKUlBRMmTIFWVlZuH79OjZs2IDc3Fxs27YN7777riXq0qYa3iPZ3D2TvJ+SiIjamsEudB8fH10Lul27dggMDIRWq0VGRgYiIiIAABEREcjIyAAAnDp1CsOHD4ckSejevTtKS0tRXFzchlUgIiJyPCbNxFZQUIArV64gJCQEN2/ehI+PD4CaJH/r1i0AgFarhVqt1q2jUqmg1Wp1y9ZKTU1FamoqACA+Pr7eOvrZbqv2z9j1x6hWq+Hs7GxEHW2TnGMHGH9zNm/ejMzMTHh7eyMhIQEAsGvXLnz33Xfw8vICAEyePBn9+vUDAHz11VdIS0uDQqHA9OnT0adPnzaJi4gMMzqBl5WVISEhAdOmTYOHh0eTywkhGpVJktSoTKPRQKPR6F4XFhYaG4rNMRR7YWEh1Gq1bOso59gB+4s/ICDAbNuOjIzEmDFjkJSUVK88KioK48aNq1d27do1HDt2DB988AGKi4uxcuVKrF+/HgoFr4Ulsgaj9rzKykokJCRg2LBhGDRoEADA29tb1zVeXFysO1tXqVT1DjZFRUWNWt9EZBt69uwJpVJp1LIZGRkYMmQIXFxc0KFDB/j7+yMvL6+NIySiphhsgQshkJycjMDAQDz55JO68rCwMBw+fBjjx4/H4cOHMWDAAF35/v37MXToUOTm5sLDw4MJ/P/xKnWSiwMHDiA9PR3BwcF46aWXoFQqodVq0a1bN90yvr6+0Gq1etc3fYjM8kwfmrDdYby21nCY0Brfp9yHwprSmnoZTOAXLlxAeno6unTpggULFgCoGRMbP348EhMTkZaWBrVajbi4OABA3759kZmZiblz58LV1RUxMTEtCoyIrGP06NF49tlnAQA7d+7EZ599hpiYGL3DY02RwxCZ3IdWLKnh52SNz81ev6/WDJEZTOA9evTArl279L63fPnyRmWSJCE6OtroAIjIttx33589RCNHjsR7770HoGZ4rKioSPeeVquFr6+vxeMjohq8+oSI6ql72+fJkyfRuXNnADXDY8eOHUNFRQUKCgqQn5+PkJAQa4XZ5vbtLOGcDmTTTLqNjIjsy7p163Du3Dncvn0bs2fPxqRJk5CTk4Off/4ZkiTBz88Ps2bNAgB07twZgwcPRlxcHBQKBWbOnMkr0ImsiAncCvbtLOGFbGQTYmNjG5WNGDGiyeUnTJiACRMmtGVIRGQkJnAzYDcbERFZGvu/iIjIJLw+wDYwgRMREckQu9AtoOZMtfmzVU7yQkREpmALnIiIWoRd6dbFBE5ERCRDTOBWYujMlWe2RETUHI6BExFRs9iYsE1sgRMREckQEzgRUQNscZIcsAvdxvDAQURExmALnIiIzIINEMtiAiciIpIhJnAZ4a1lRERUiwmciIhIhpjAiYiIZIhXoRMR/T99Q1QctiJbxQRO5MA2b96MzMxMeHt7IyEhAQBw584dJCYm4saNG/Dz88O8efOgVCohhMD27duRlZUFNzc3xMTEIDg42Mo1IHJc7EIncmCRkZFYsmRJvbKUlBT06tULGzZsQK9evZCSkgIAyMrKwvXr17FhwwbMmjUL27Zts0bIZIN4ga11MIETObCePXtCqVTWK8vIyEBERAQAICIiAhkZGQCAU6dOYfjw4ZAkCd27d0dpaSmKi4stHjMR1WAXOhHVc/PmTfj4+AAAfHx8cOvWLQCAVquFWq3WLadSqaDVanXL1pWamorU1FQAQHx8fL31bIWzs7OeuNiKbK22+q71f1/y15p6MYETkVGEEI3KJEnSu6xGo4FGo9G9LiwsbLO4WkqtVttkXHLXVp+pvX5fDesVEBBg9LrsQrdxHFsiS/P29tZ1jRcXF8PLywtATYu77oGmqKhIb+ubiCyDCZyI6gkLC8Phw4cBAIcPH8aAAQN05enp6RBC4OLFi/Dw8GACJ7IidqETObB169bh3LlzuH37NmbPno1JkyZh/PjxSExMRFpaGtRqNeLi4gAAffv2RWZmJubOnQtXV1fExMRYOXoix8YETuTAYmNj9ZYvX768UZkkSYiOjm7rkIjISEzgROTweJ0JyRHHwImIiGTIYAucUy22LZ75E5E9qT2mPfXcfXpfk/kYbIFzqkUiIiLbYzCBc6pF28N7w4mIqEVj4KZOtUiWxQRP1Lx9O0uwPSnP2mEQtYpZr0I3ZapF0+dKduyEVPP5lOgp06fEwPumkfscxIyfiOxRixJ47VSLPj4+LZ5qUQ5zJdsSfZ+Poc/MXJ+p3Ocgtrf4TZkrmcha2AvY9lrUhc6pFm0Lu8yJiByPwRY4p1okIiKyPQYTOKdaJCIisj2ciY2IiEiGmMDtGMfFiYjsFx9mYkeYsImIHAdb4HaOV6gTEdknJnAiIrI4Ni5ajwmciIhIhjgGLhM8UyVLe+211+Du7g6FQgEnJyfEx8c3+ShhIkOMeaxow+McH0HaPCZwImrSihUrdFMlA38+Snj8+PFISUlBSkoKpkyZYsUIiRwXu9CJyGhNPUqYiCyPLXAiatKqVasAAKNGjYJGo2nyUcINmf60QUur6arl0JTl/flb0PfUxKafuGivT+VrTb2YwB1UU+NR+3aWcNyJAAArV66Er68vbt68iXfeecekp6DxaYPUlIa/heZ+G3Xfk/tTBZvSmqcNMoE7CGMuICGqy9fXF0DN44MHDBiAvLy8Jh8lTGQs9nqYD8fAZYz3UVJbKSsrw71793T//+mnn9ClS5cmHyVMRJbHFjgRNXLz5k2sXbsWAFBVVYXHHnsMffr0QdeuXfU+StiW1R0W4gkv2RMmcCJqpGPHjlizZk2jck9PT72PEiZqDZ5YtQy70ImISLYceSiRLXBqhBe8kb1x1AM82Te2wImIyCY1bF07cmtbH7bACQBbKGR/+Jsme8cE7mB4UCMiW9LSYxKPZexCJyIiO+QI3e1M4GTQ9qQ8u98RiIjkhl3oDo6JmewF754gwLGOaUzg1KSGOwIPkERkDTXHnrZLzA2PbXI51rELnczKEcadiIhsARM4mUzfvZlERNZkTOPB3hoYTOBks+xtZyMi8zDluGDPxxGOgZNZGHPmC9QfU5LLOBMRkS1iAqcWa8lZLZM2tTV7bW0RNcQudLIKfQdZe+7qIvPh74RawlLHHEv+NtkCpzZh6QMsW/ZEZA3WPPYwgZMsteUJwr6dJTwRkAGetJG5NdVKB2zzd9YmCfz06dPYvn07qqurMXLkSIwfP74t/gzZqaaSsymJte5Otz0pT/d/aj1r7N+2fBAl+TG2AaDvdllb+g2aPYFXV1fjo48+wtKlS6FSqbB48WKEhYWhU6dO5v5TZCeMvYLdmGVN3TaZxtL7N78/sjXNNTCaKmurpG/2BJ6Xlwd/f3907NgRADBkyBBkZGQwgZNFmHp/aF36bnFrqOFUi8Ys09zO21wMtqit9u/WfA5M8mRrLDUNtSSEEObc4PHjx3H69GnMnj0bAJCeno7c3FzMnDmz3nKpqalITU0FAMTHx5szBCJqI9y/iWyH2W8j03c+IElSozKNRoP4+Hijd+5Fixa1OjZrknP8co4dYPzm1Fb7tzXY0udqTqyXvLSmXmZP4CqVCkVFRbrXRUVF8PHxMfefISIr4P5NZDvMnsC7du2K/Px8FBQUoLsnMiwAAAeKSURBVLKyEseOHUNYWJi5/wwRWQH3byLb4fTWW2+9Zc4NKhQK+Pv7Y+PGjdi/fz+GDRuG8PBws2w7ODjYLNuxFjnHL+fYAcZvLm25f1uDrXyu5sZ6yUtL62X2i9iIiIio7XEudCIiIhliAiciIpIhm5gL3dDUjBUVFdi0aRMuX74MT09PxMbGokOHDgCAr776CmlpaVAoFJg+fTr69Okji9h/+ukn7NixA5WVlXB2dsbUqVPxyCOPWDT21sRfq7CwEPPmzcPEiRMxbtw4S4ffqvivXr2KrVu34t69e5AkCatXr4arq6vNx15ZWYnk5GRcuXIF1dXVGD58OJ555hmLxS03hj7nwsJCJCUlobS0FNXV1XjhhRfQr18/K0VrnM2bNyMzMxPe3t5ISEho9L4QAtu3b0dWVhbc3NwQExMji/FjQ/U6cuQI9u7dCwBwd3dHdHQ0HnzwQQtHaTpD9aqVl5eHv//975g3b55x15YIK6uqqhKvv/66uH79uqioqBDz588Xv/zyS71l9u/fL7Zs2SKEEOLo0aPigw8+EEII8csvv4j58+eL8vJy8fvvv4vXX39dVFVVySL2y5cvi6KiIiGEEFevXhWzZs2yWNy1WhN/rTVr1oiEhASxd+9ei8VdqzXxV1ZWijfffFNcuXJFCCHErVu3ZPPbOXLkiEhMTBRCCFFWViZiYmLE77//brHY5cSYzzk5OVkcOHBACFFzTImJibFGqCbJyckRly5dEnFxcXrf//HHH8WqVatEdXW1uHDhgli8eLGFI2wZQ/U6f/68uH37thBCiMzMTLuplxA1v9W33npLvPvuu+KHH34wartW70KvOzWjs7OzbmrGuk6dOoXIyEgAQHh4OM6ePQshBDIyMjBkyBC4uLigQ4cO8Pf3R15enixiDwoKgq+vLwCgc+fOqKioQEVFhcVib238AHDy5El07NjRatPktib+M2fOoEuXLrqzd09PTygUltsdWvvZl5WVoaqqCuXl5XB2doaHh4fFYpcTYz5nSZJw9+5dAMDdu3dlcV97z549oVQqm3z/1KlTGD58OCRJQvfu3VFaWori4mILRtgyhur10EMP6d7v1q1bvTkJbJmhegHAt99+i0GDBsHLy8vo7Vo9gWu1WqhUKt1rlUoFrVbb5DJOTk7w8PDA7du3G63r6+vbaN221JrY6zpx4gSCgoLg4uLS9kE3ERtgWvxlZWXYu3cvJk6caNGYm4oNMC3+/Px8SJKEVatWYeHChbpuOTnEHh4eDnd3d8yaNQsxMTF46qmnDB4cHJUxn/PEiRNx5MgRzJ49G6tXr8aMGTMsHabZabVaqNVq3Wt99Za7tLQ09O3b19phmIVWq8XJkycxevRok9azegIXRkzN2NQy+sotqTWx1/rll1+wY8cOvPLKK+YP0IDWxL9r1y5ERUXB3d29zeIzpDXxV1VV4fz585gzZw7efvttnDx5EtnZ2W0Wa0OtiT0vLw8KhQJbtmzBpk2bsG/fPvz+++9tFqucGfM5f//994iMjERycjIWL16MjRs3orq62lIhtglj6i1nZ8+excGDB/Hiiy9aOxSz+OSTT/Diiy+a3Ato9YvYjJmasXYZlUqFqqoq3L17F0qlstG6Wq1W1y1t67HXLr927Vq89tpr8Pf3t1jcDWOrZUr8eXl5OHHiBHbs2IHS0lJIkgRXV1eMGTNGFvGrVCr07NlT113Vt29fXLlyBb169bL52I8ePYo+ffrA2dkZ3t7eeOihh3Dp0iXdE8LoT8Z8zmlpaViyZAkAoHv37qioqMDt27fh7e1t0VjNSaVSobCwUPfanqa8vXr1KrZs2YLFixfD09PT2uGYxaVLl7B+/XoAwK1bt5CVlQWFQoGBAwc2u57VW+DGTM3Yv39/HDp0CEDN05BCQ0MhSRLCwsJw7NgxVFRUoKCgAPn5+QgJCZFF7KWlpYiPj8fkyZPRo0cPi8VcV2vif/vtt5GUlISkpCSMHTsWzzzzjEWTd2vj7927N/73v//hjz/+QFVVFf773/9adCy/NbGr1WrdeHhZWRlyc3MRGBhosdjlxJjPufbzBIBr166hoqLCpHFIWxQWFob09HQIIXDx4kV4eHjYRQIvLCzE2rVr8frrryMgIMDa4ZhN7bE0KSkJ4eHhiI6ONpi8ARuZiS0zMxOffvopqqur8fjjj2PChAnYuXMnunbtirCwMJSXl2PTpk24cuUKlEolYmNjda2NPXv24ODBg1AoFJg2bZrFx0RaGvuXX36JlJSUei3vpUuXWvysvzWffa1du3bB3d3dKreRtSb+9PR0pKSkQJIk9O3bF1OmTJFF7GVlZdi8eTOuXbsGIQQef/xxq3z2cmHoc7527Rq2bNmCsrIyAMCUKVPQu3dvK0fdvHXr1uHcuXO6noJJkyahsrISADB69GgIIfDRRx/hzJkzcHV1RUxMDLp27WrlqA0zVK/k5GScOHFCN77v5ORk00+8q2WoXnUlJSWhf//+Rt1GZhMJnIiIiExj9S50IiIiMh0TOBERkQwxgRMREckQEzgREZEMMYETERHJEBM4ERGRDDGBExERydD/AYdFsTdD9LXdAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 504x216 with 2 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, ax = plt.subplots(1,2,figsize=(7,3))\n",
"ax[0].hist(states_[0], bins=100, color=\"C2\")\n",
"ax[0].set_title(r'posterior of $b$')\n",
"ax[1].hist(states_[1], bins=100, color=\"C2\")\n",
"ax[1].set_title(r'posterior of $\\sigma_y$')\n",
"fig.tight_layout()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- b, s_yともに、真値([0.6, 0.4])とかけ離れたサンプルを得た。\n",
" - bが０付近（すなわちサポートの淵）に集中していることから、bijectorによる探索空間の変換があれば解決したかも。\n",
"- MCMCの需要率も50%を下回った。\n",
" - 通常HMCでは75%が目安とされることから、著しく低い需要率と言える。"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 4. HMC + Step Size Adaptation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 4.1. Define the Kernel\n",
"- `tfp.mcmc.SimpleStepSizeAdoptation`による自動ステップサイズ調整を用いる。\n",
" - TFP 0.7.0時点では、`tfp.mcmc.TransformedTransitionKernel`（bijectorによる探索空間の変換）や、この後使用する`tfp.mcmc.ReplicaExchangeMC`と併用するとエラーが発生する。"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"kernel2 = tfp.mcmc.SimpleStepSizeAdaptation(\n",
" inner_kernel=kernel1,\n",
" num_adaptation_steps=int(num_burnin_steps * 0.8),\n",
" adaptation_rate=0.001\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 4.2. Inference"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"states, kernel_results = tfp.mcmc.sample_chain(\n",
" num_results=num_results, \n",
" num_burnin_steps=num_burnin_steps,\n",
" kernel=kernel2,\n",
" current_state=initial_state,\n",
")"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"with tf.Session() as sess:\n",
" sess.run(init_g)\n",
" sess.run(init_l)\n",
" [states_, kernel_results_ ]= sess.run([states, kernel_results])"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"acceptance rate : 64.75999999999999 %\n"
]
}
],
"source": [
"try:\n",
" print(f'acceptance rate : {kernel_results_.inner_results.is_accepted.mean()*100} %')\n",
"except AttributeError:\n",
" print(f'acceptance rate : {kernel_results_.is_accepted.mean()*100} %')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 4.3. Result"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAfAAAADQCAYAAAD4dzNkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de1hUdf4H8PcMCAiDODB4Aa0V0S0NRcUbplCM5UZrrKlbm5a3zNU0YTXNx9uumeQlFKNVoyzLfVKfNTbbTVvEoGRNXLyhpqD2PBm4CAMoqD9u398fPEwMtxnmes7M+/U8Po9z5pwzn+9wvudzvpdzRiGEECAiIiJZUTo6ACIiIuo4JnAiIiIZYgInIiKSISZwIiIiGWICJyIikiEmcCIiIhliAiciIpIhJnAiIhcyY8YMaLVah8awYsUKdO/eHQqFAh999JFDY5EzJnBqlZwqeVJSEh5++GH7BUZkR1qtFjNmzLDa/rZt24YDBw5YbX8d9f3332PDhg3YtWsXioqK8Pvf/95hscidu6MDIOvQarXo1auX1a5mt23bhvr6eqvsyxyNlTwtLQ0jR46En59fm+vm5uYiPDzcjtERyU91dTU8PDzarUsd3Zc58vPzoVQq8cwzz1gch6tjC5wMVFdXAwD8/PygVqutsi9zNK3kPXr0QOfOndtcNzc3F0OGDDH7s4jMFR0djVmzZmH58uXQaDTo0qUL5syZg3v37unXqampwfLlyxEcHAwPDw8MGDAAf/vb3/Tvf/fddxgzZgx8fX3h6+uLwYMH48iRIwAaesKOHj2Kjz/+GAqFAgqFAt98841+2+3bt+Ohhx6Cl5cX+vXrh/Xr16O2tlYf2+zZs7Fq1Sr07NkTwcHB+n027V0zFl97+2pNe/ubMWMGpk+fjvr6en152nP+/HlMnDgRPXv21K/f+O+DDz5od1uXIMjmoqKixMyZM8WyZctEQECA8PX1FbNnzxZ3797Vr1NdXS2WLVsmgoKCRKdOncTDDz8s9u7dq3//22+/FZGRkUKlUgmVSiUGDRokDh8+LIQQ4qWXXhIADP4dO3ZMv21ycrL49a9/LTw9PUVoaKh48803RU1NjT62WbNmiZUrV4oePXoIjUaj32dMTIzJ8bW3r9a0t7/WytOWqqoq4ebmJnbs2CGeeeYZ4ePjI/r06dMiNiJbiIqKEr6+vmLOnDni4sWL4osvvhCBgYFi4cKF+nWWLFki/P39xf79+8Xly5fF+vXrhUKhEOnp6aK2tlao1WoRHx8vrly5Iq5cuSIOHjwosrKyhBBClJeXi7Fjx4qpU6eKoqIiUVRUJP7v//5PCCHEmjVrxAMPPCAOHjworl27Jv75z3+K3r17i5UrV+pjU6lU4pVXXhEXLlwQ586dE0K0rNvtxde0nK3tqzXt7a+8vFxs3bpVuLm56cvTlvz8fNGlSxfxhz/8QZw5c0acOnVKDBkyRAQGBopPPvmk3W1dBRO4HbCSt2StSp6dnS0AiH79+okvv/xSFBQUiIULFwp3d3dx5coVM/5aRKaLiooSDz74oKitrdUv27lzp/Dw8BCVlZWiqqpKeHh4iJSUFIPt4uLixGOPPSZ0Ol2LC+7mYmJixEsvvWSwrKqqSnTu3Fl89dVXBss//vhj4efnp4+tX79+oq6uzmCdpnXbWHxNy9navpozZX+7d+8Wbm5u7e5HCCFiY2NFVFSUqK+v1y/77LPPBABRWlpqdHtXwARuB6zkhqxZyVNSUoRCoRC5ubn6ZdXV1aJz584iNTXV6PZEloiKihJTpkwxWJaXlycAiLNnz4qzZ88KAOL8+fMG6yQlJYlu3boJIYSYM2eO8PDwEBMmTBAbNmwQP/zwg8G6rdXtkydPCgDC29tb+Pj46P95eXkJAKK4uFhERUWJ5557rkXMTeu2KfE1lrO1fTVnyv5Mqds6nU506tRJ/P3vfzdYfujQIQFA6HQ6o7G4Ao6B28mIESPg5uamfz1mzBhUV1fj6tWrKCgoQHV1NcaNG2ewTVRUFC5cuAC1Wo05c+bgySefxG9+8xskJibi8uXLRj/zwoULuHfvHp599lmoVCr9v1deeQUVFRW4desWAGDYsGFQKts+FIzF15SxfXV0f8bk5uZi/PjxBmPgnTp1gqenp0Mn4ZHrEq38QnPzsV4hhH7Z+++/j//+978YP348MjMz8cgjj2Dnzp3tfkbjsX3gwAGcOXNG/+/8+fPIz8+Hv78/AMDHx8ekmNuLr5Gp+zJ1f+05ffo0ampqMHToUIPlJ0+eRGhoKNRqNc6fP4+RI0fq38vLy0NkZGSr37+zYgJ3EFZyyys50JDAhw0bZrAsPz8f5eXlLZYT2UJOTg7q6ur0r//zn//Aw8MDffv2RWhoKDw9PZGZmWmwTVZWFgYOHKh//cgjjyAhIQFfffUVZs+ejV27dunf8/DwMNg/AAwcOBBeXl64du0aQkNDW/xr2lhoj6nxmcpa+2ssb1VVlX5ZRUUFPvzwQ/0tdQMHDsT169dRU1MDAEhISMCmTZs6fA6RM95GZieNlbyxYjWt5AqFQn/QNz3IW6vkjRV93rx52LVrF1555RUAxiv5U089ZXbsTStle/HZe3/V1dW4cOECxo8fb7B8y5YtGDZsWIurdyJbKC0txYIFC/Daa6/h2rVrWLVqFV5++WX9xeyiRYuwatUqBAYGIjw8HAcOHMA//vEP/Pvf/0ZBQQHef/99/Pa3v0Xv3r1RWFiIb7/91uDY7dOnD44dO4arV6/Cz88Pfn5+UKlUWLFiBVasWAEAGD9+PGpra3H+/HmcPn0ab7/9tkmxe3t7txtfR1lrfyNGjECXLl2wbNkyJCYmory8HK+//jp69eqFpUuXAgCUSiXCwsJw/vx53LhxA2q1GmPGjOlwzHLGBG4nrOTW319eXh6qq6tx8OBBTJgwASEhIdixYwf27NmD7777rsNxEZlj8uTJ8PX1xaOPPorq6mpMmTIFGzdu1L+/fv16KJVKLF68GLdu3UJoaCg+/fRTxMTEoKioCPn5+Xjuuedw69YtBAQEIDY2Fps3b9Zv/6c//Qnnz5/H4MGDUVVVhWPHjiE6OhqrVq1CUFAQtm/fjiVLlqBz587o379/hx/60l585rDG/vz8/HDw4EHEx8dj6NCh0Gg0mDJlCtatW2dw//moUaOQnZ2NXbt24YsvvjArXllz2Oi7C2m8jaxx5rVKpRIzZ84UVVVV+nXau62qsLBQ/O53vxPBwcHCw8ND9OzZU8yZM0eUl5frt7969aoYO3as8PHxaTHhLTU1VQwePFh4enqKrl27ihEjRoj33ntPH9vs2bNbxGzubWSt7as1xvZnykSX1NRUERQUJL7++mvRt29f4enpKaKiolpMoCGylY4c82R9X3zxhfD39xfLli1zdCgOoRDChUb8HSQ6OhqhoaFITU11dChEZEWs245VUFCAMWPGID8/H126dHF0OHbHSWxERCRLycnJeOutt1wyeQMcAyciMlvTx5qS/Vy9ehWxsbGIjIzErFmzHB2Ow7ALnYiISIbYhU5ERCRDTOBEREQyxAROREQkQ5KZxFZYWGjWdhqNBiUlJVaOxj4Yu/3JNW6gZexBQUEOjMZ85tZ1S8j5794alke6bFGWtuo6W+BEREQyxAROREQkQ0a70Kurq7FmzRrU1tairq4Oo0aNwtSpU1FcXIytW7eisrISffr0wcKFC+Hu7o6amhq8++67uHbtGnx9fbF48WJ069bNHmUhIiJyGUZb4J06dcKaNWuwadMmbNy4EWfOnMGVK1fw6aefIjY2FsnJyfDx8UFGRgYAICMjAz4+Pti+fTtiY2Oxd+9emxeCiIjI1RhN4AqFAl5eXgAafqO1rq4OCoUCFy5cwKhRowA0PA84JycHAHDq1ClER0cDaPilmLy8PJf6gXUiIiJ7MGkWen19PZYtW4abN2/iySefRPfu3eHt7a3/bWt/f3/odDoAgE6nQ0BAAADAzc0N3t7euHPnTotn1aanpyM9PR0AkJiYCI1GY14B3N3N3tbR7BX77pQCAMDMBaFW26dcv3e5xg3IO3Yy3aF95fjt77s6OgySAZMSuFKpxKZNm1BVVYXNmzfj559/bnPd1lrbCoWixTKtVgutVqt/be60eznffmDv2K35WXL93uUaN+A8t5ERkXV0aBa6j48PBgwYgPz8fNy9exd1dXUAGlrd/v7+AICAgACUlpYCaOhyv3v3LlQqlZXDJksc2leOQ/vKHR0GkUvrSD1knaXWGE3gt2/fRlVVFYCGGennz59HcHAwBg4ciBMnTgBo+EWeiIgIAMCwYcP0v9Bz4sQJDBw4sNUWOMkDTxxERNJktAu9rKwMKSkpqK+vhxACo0ePxrBhw9CrVy9s3boVn332Gfr06YPHH38cAPD444/j3XffxcKFC6FSqbB48WKbF4KILMdbRonkxWgCf/DBB7Fx48YWy7t3744NGza0WO7h4YGEhATrREdEdtN4y6iXlxdqa2uxevVqhIeH48svv0RsbCzGjBmDXbt2ISMjA0888YTBLaPHjx/H3r17ER8f7+hiSB57tMha+CQ2IgLAW0alhENXZArJ/JgJETmelG8ZtYS0bsEzTMwt4ypv5b1yg9fSKo/lnKk89iwLEzgR6Un5llFLSPn2wfbiav5e42spl8cczlQee/4aGRM4EbXQ2i2jbm5urd4yGhAQwFtGTcAucbI2joG7MHPG2Tg257x4y6j0sL5Re9gCp1YZO2k0Pp6Vj3x0HrxllEhemMCJCABvGSWSG3ahExERyRATOBGRTHBMnJpiFzoZ4MmByDFY96ij2AInPZ5AiORhd0oB6ysxgRMREckREzhZBVsDRET2xTFwMgkTNBGRtLAFTkQkU5yV7trYAieL8ORBROQYRhN4SUkJUlJSUF5eDoVCAa1Wi6eeegr79+/H0aNH9T8d+Pzzz2Po0KEAgM8//xwZGRlQKpWYOXMmwsPDbVsKIiIiF2M0gbu5uWH69OkICQnBvXv3sHz5cgwaNAgAEBsbi4kTJxqsf+PGDWRnZ+Odd95BWVkZ1q1bh23btkGpZG89EZEtHdpXzt8ncCFGE7harYZarQYAdO7cGcHBwdDpdG2un5OTg8jISHTq1AndunVDjx49UFBQgP79+1svajJLW93d7AYnIpKfDo2BFxcX4/r16wgNDcUPP/yAI0eOICsrCyEhIXjxxRehUqmg0+nQr18//Tb+/v6tJvz09HSkp6cDABITE6HRaMwrgLu72ds6mv1it0+ClsPfgccLSQ0voMlcJifw+/fvY8uWLZgxYwa8vb3xxBNPYPLkyQCAffv2Yc+ePZg/fz6EECbtT6vVQqvV6l+XlJR0MPQGGo3G7G0dTc6xt0YOPzEq5++8eexBQUEOjIaIHM2kgena2lps2bIFY8eOxciRIwEAXbt2hVKphFKpRExMDK5evQoACAgIQGlpqX5bnU4Hf39/G4RORCR9vNWLbMVoC1wIgR07diA4OBhPP/20fnlZWZl+bPzkyZPo3bs3ACAiIgLJycl4+umnUVZWhqKiIoSGhtoofCKyBt5tQiQ/RhP45cuXkZWVhQceeABLly4F0FCJjx8/jh9//BEKhQKBgYGYO3cuAKB3794YPXo0EhISoFQqMXv2bM5AJ5I43m1CJD9GE/hDDz2E/fv3t1jeeBXemkmTJmHSpEmWRUZEdsO7TYjkh09iIyID1rzbBLDeHSeWcOwMftuPf/9StnJZ3qngTHdY2LMsTOBEpGftu00A691xYgk5331giqZlk2M5nenvY4uytHXHCQesiAgA7zYhkhsmcCJq926TRs3vNsnOzkZNTQ2Ki4t5twmRA7ALnayu8Z5XKT/QhQzxbhMi+WECJyLebUIkQ7xkJiIikiG2wJ0YH99IROS82AInm+OzoMlV8Fgne2ICJyIikiEmcCIiIhliAiciIpIhJnAiIiIZ4ix0J8RJNESOxTpI9sAE7iT49DMi19XeBQPPDc6LCZxshq0QIiLbMZrAS0pKkJKSgvLycigUCmi1Wjz11FOorKxEUlISbt26hcDAQMTHx0OlUkEIgd27d+P06dPw9PTE/PnzERISYo+yEBFRB7GFLl9GJ7G5ublh+vTpSEpKwvr163HkyBHcuHEDaWlpCAsLQ3JyMsLCwpCWlgYAOH36NG7evInk5GTMnTsXqampNi8EERE14MNkXIfRBK5Wq/Ut6M6dOyM4OBg6nQ45OTmIiooCAERFRSEnJwcAcOrUKYwbNw4KhQL9+/dHVVWVwU8SEhERkeU6NAZeXFyM69evIzQ0FBUVFVCr1QAakvzt27cBADqdDhqNRr9NQEAAdDqdft1G6enpSE9PBwAkJiYabNOhAri7m72to1kj9t0pBQavpXzl3RjbzAWO+91oVz9eiMh5mJzA79+/jy1btmDGjBnw9vZucz0hRItlCoWixTKtVgutVqt/XVJSYmooBjQajdnbOpqcY7dE40WHI8bc5PydN489KCjIgdFQU1K+cCbnZVICr62txZYtWzB27FiMHDkSAODn54eysjKo1WqUlZWhS5cuABpa3E1PMqWlpS1a30QkLZysSiQ/RsfAhRDYsWMHgoOD8fTTT+uXR0REIDMzEwCQmZmJ4cOH65dnZWVBCIErV67A29ubCZxI4jhZ1TycMEaOZDSBX758GVlZWcjLy8PSpUuxdOlS5ObmIi4uDufOncOiRYtw7tw5xMXFAQCGDBmCbt26YdGiRdi5cyfmzJlj80K4Ep4wyBY4WdX58Fzh/Ix2oT/00EPYv39/q++tXr26xTKFQsGkTSRj1pysClhvwqolbDcBUD4Jsu3ylxt53/acaYKmPcvCJ7ERkZ61J6sC1puwagk5T160FmPld+T340x/H1uUpa0Jq/w1MiIC0P5kVQCcrOok2LXuPJjAiYiTVYlkiF3oMuUMV9B8BrN0NE5WfeCBB7B06VIAwPPPP4+4uDgkJSUhIyMDGo0GCQkJABomq+bm5mLRokXw8PDA/PnzHRk+kUtiApcJZ0jYJF2crNoxcqyPcoyZ2scudCIiIhliC5yIyAWxRS5/bIETERHJEBM4ERGRDDGBkyTxXlUiovYxgRMREckQEzhJAlvcREQdwwROREQkQ0zgREREMsQETkREJENM4ORwHPsmKeF8DJILo09ie++995Cbmws/Pz9s2bIFALB//34cPXpU/9OCzz//PIYOHQoA+Pzzz5GRkQGlUomZM2ciPDzchuETERG5JqMJPDo6GhMmTEBKSorB8tjYWEycONFg2Y0bN5CdnY133nkHZWVlWLduHbZt2walkg19Mg1bPkREpjGaWQcMGACVSmXSznJychAZGYlOnTqhW7du6NGjBwoKCiwOkoiIiAyZ/WMmR44cQVZWFkJCQvDiiy9CpVJBp9OhX79++nX8/f2h0+msEqgr4e9kE5Ej8RwkD2Yl8CeeeAKTJ08GAOzbtw979uzB/PnzIYQweR/p6elIT08HACQmJkKj0ZgTCtzd3c3e1tHajr2h8hi+55pdy9b+2zrn8WIdnO9i6NC+ciYwkjSzEnjXrr8c1DExMXj77bcBAAEBASgtLdW/p9Pp4O/v3+o+tFottFqt/nVJSYk5oUCj0Zi9raMZi12u5bIma38HznS8BAUFWXX/nO/SNrZISYrMqm1lZWX6/588eRK9e/cGAERERCA7Oxs1NTUoLi5GUVERQkNDrRMpEdkU57sQyYvRFvjWrVtx8eJF3LlzB/PmzcPUqVNx4cIF/Pjjj1AoFAgMDMTcuXMBAL1798bo0aORkJAApVKJ2bNnO+0Vub1wVjY5mqXzXaw1XGaJjg0//FLnmte/X/bhfPVSo9Fgd0pBi2X2IOehrebsWRajCXzx4sUtlj3++ONtrj9p0iRMmjTJsqiISBKsMd/FWsNllrDW0Ilch19M0VrZ7FVeOQ9tNWeLsrQ1XMbmMUkan4rlWF27doVSqYRSqURMTAyuXr0KoGPzXYjINpjAiahNnO9CJF1m3wdORM6F812MY28QSQkTOBEB4HwXIrlx/ktmIiIiJ8QETkRE7eJkUmliAiciIpIhJnAiIiIZYgInIiKSISZwIiIiGWICJyIiTlKTISZwkgXOgiUiMsQETkREJENM4ERERDLER6lKSPMuYnYZExFRW9gCJyIikiGjLfD33nsPubm58PPzw5YtWwAAlZWVSEpKwq1btxAYGIj4+HioVCoIIbB7926cPn0anp6emD9/PkJCQmxeCCIiIldjtAUeHR2NFStWGCxLS0tDWFgYkpOTERYWhrS0NADA6dOncfPmTSQnJ2Pu3LlITU21TdTksjgbnayFxxLJndEEPmDAAKhUKoNlOTk5iIqKAgBERUUhJycHAHDq1CmMGzcOCoUC/fv3R1VVFcrKymwQtnPYnVLAEwgREZnFrElsFRUVUKvVAAC1Wo3bt28DAHQ6HTQajX69gIAA6HQ6/bpEJF0cLiOSF6vOQhdCtFimUChaXTc9PR3p6ekAgMTERIPE3xHu7u5mb+sou1MKDF7/Ej9b46ZypeOlka1jj46OxoQJE5CSkqJf1jhcFhcXh7S0NKSlpWHatGkGw2X5+flITU3FW2+9ZbPYyDHYQyhtZiVwPz8/lJWVQa1Wo6ysDF26dAHQ0OIuKSnRr1daWtpm61ur1UKr1epfN92uIzQajdnbSoXc43cEVzxemsceFBRk1f0PGDAAxcXFBstycnKwdu1aAA3DZWvXrsW0adPaHC5jbxuR/ZiVwCMiIpCZmYm4uDhkZmZi+PDh+uWHDx/GmDFjkJ+fD29vb1ZosqlD+8rx2993dXQYTssaw2XW6m2zROu9Fw2ty7aWU0u2+tvJuWesOXuWxWgC37p1Ky5evIg7d+5g3rx5mDp1KuLi4pCUlISMjAxoNBokJCQAAIYMGYLc3FwsWrQIHh4emD9/vs0LQK6JXXuO1ZHhMmv1tlmivZ4XufbIOIKtvis594w1Z4uytNXbZjSBL168uNXlq1evbrFMoVBgzpw5HQyNiKTKGsNlRGQbfBIbEbWpcbgMQIvhsqysLAghcOXKFQ6XETkAn4VORAA4XEbGNQ5dcd6JNDCBExEADpcRyQ270ImIiGSILXAiIvDOBpIftsCJiIhkiC1wCeCVPxERdRRb4ERERDLEBE5ETq/5Dwg1xd8F7zh+Z9LABE5ERCRDTOBERGSWpi1xtsrtjwmciIhIhpjAiYiIZIi3kdkBnx9MJB3s5iVnwRY4OQ2OwRE5Huuh/bAFbmNND2Qe1EREZC1sgZPTYQuAiFyBRS3wBQsWwMvLC0qlEm5ubkhMTERlZSWSkpJw69YtBAYGIj4+HiqVylrxEhEREazQhb5mzRp06dJF/zotLQ1hYWGIi4tDWloa0tLSMG3aNEs/hqhNbG3bnjNcrPM4IWdj9S70nJwcREVFAQCioqKQk5Nj7Y8gMgm70q1rzZo12LRpExITEwH8crGenJyMsLAwpKWlOThCItdicQt8/fr1AIDx48dDq9WioqICarUaAKBWq3H79u1Wt0tPT0d6ejoAIDExERqNxqzPd3d3N3tbW2rv2ctkX02PD6keL6aQWuw5OTlYu3YtgIaL9bVr17K3jfSaXzw33kbL22qtx6IEvm7dOvj7+6OiogJvvvkmgoKCTN5Wq9VCq9XqX5eUlJgVg0ajMXtbcg1Njw85Hy/NY+9IfbMGR1+sW4Y9MbZkSk9X83Wc5cK6OXuWxaIE7u/vDwDw8/PD8OHDUVBQAD8/P5SVlUGtVqOsrMxgfNwVsMtWenjFbzkpXKyTc3GWC+vmbFGWtuqb2WPg9+/fx7179/T/P3fuHB544AFEREQgMzMTAJCZmYnhw4eb+xFEJBHtXawDcMmLdbIM56hYzuwWeEVFBTZv3gwAqKurw6OPPorw8HD07dsXSUlJyMjIgEajQUJCgtWCJSL7u3//PoQQ6Ny5s/5iffLkyfqL9bi4OF6sEzmA2Qm8e/fu2LRpU4vlvr6+WL16tUVBSRm7Y8nV8GKdSJr4KFUiaperXqyTfRzaV46ZC5xjApu9MYFbCcdyiIjInvgsdCIicqjdKQUtfviJjSLjmMDJZRzaV27wgB2eIIhIztiFTi6HiZuInAFb4OTS2FUnH+b8nfj3JWfGFriFeHIgIrKO5udT3rbbPiZwEzFRExGRlLALnYiISIbYAm+FKd02bJETSRfrp3NhV3rr2AJvByfAEI8BIpIqJnAiIpIVXlg3cOku9ObdMm0dEDxQnF/zv3HzrrrWuvDYrUdkX5aci52xvrp0AidqS0cu5pzxxCBXvNh2Pa5c/5jAiUj2mLhdkykX1M58bDCBE9mYK7cQrI3fJbXFmRN1W2yWwM+cOYPdu3ejvr4eMTExiIuLs9VHmYyVn+ypvROKsTF3uXB0PXfFkzaZxpRhsLZa6XKpjzZJ4PX19fjggw+wcuVKBAQE4I033kBERAR69epli49rlykVnCcBsiZLJkPK6SLTVvWcz2EgqWirO14q9dMmCbygoAA9evRA9+7dAQCRkZHIycmxWsUGjM8SZgUne3O1Y85W9ZzIXkyts209o711bb9n7QsBhRBCWLSHVpw4cQJnzpzBvHnzAABZWVnIz8/H7Nmz9eukp6cjPT0dAJCYmGjtEIjIxkyp5wDrOpGt2ORBLq1dEygUCoPXWq0WiYmJFlfo5cuXW7S9IzF2+5Nr3ID0YjelngPWq+uWkNp3ZymWR7rsWRabJPCAgACUlpbqX5eWlkKtVtvio4jIQVjPiRzLJgm8b9++KCoqQnFxMWpra5GdnY2IiAhbfBQROQjrOZFjua1du3attXeqVCrRo0cPbN++HYcPH8bYsWMxatQoa3+MXkhIiM32bWuM3f7kGjcgrdjtXc8tJaXvzhpYHumyV1lsMomNiIiIbIu/RkZERCRDTOBEREQyJOlnoRt7TOM333yDTz75BP7+/gCACRMmICYmBgBQUlKCHTt26GfJvvHGG+jWrZssYv/000+Rm5sLIQTCwsIwc+bMVm/PcUTcAJCdnY0DBw5AoVDgwQcfxGuvvaYv08GDBwEAkyZNQnR0tF1ibmRu7D/++CPef/993Lt3D0qlEi/qQGAAAAXsSURBVJMmTUJkZKQsYm909+5dxMfHY8SIES3uw3Ylpj7a9cSJE3jnnXewYcMG9O3b185Rms6S84jUWHqMS42x8nz00Ue4cOECAKC6uhoVFRX46KOPrBuEkKi6ujrx6quvips3b4qamhqxZMkS8dNPPxmsc+zYMZGamtrq9mvWrBFnz54VQghx7949cf/+fZvH3MiS2H/44QexcuVKUVdXJ+rq6sSKFStEXl6eZOIuLCwUS5cuFXfu3BFCCFFeXi6EEOLOnTtiwYIF4s6dOwb/txdLYv/5559FYWGhEEKI0tJS8fLLL4vKykpZxN7oww8/FFu3bm2zPrgCU75HIYS4e/euWL16tVixYoUoKChwQKSmsfQcKCXWOMalxNRjrdG//vUvkZKSYvU4JNuF3vQxje7u7vrHNJrixo0bqKurw6BBgwAAXl5e8PT0tGW4BiyJXaFQoLq6GrW1taipqUFdXR38/PxsHHEDU+I+evQonnzySahUKgDQx3bmzBkMGjQIKpUKKpUKgwYNwpkzZ+wSt6WxBwUFoWfPngAAf39/+Pn54fbt27KIHQCuXbuGiooKDB482G4xS5Gp9W7fvn2YOHEiOnXq5IAoTWfJeURqLD3Gpaajf5vjx4/j0UcftXocku1C1+l0CAgI0L8OCAhAfn5+i/W+//57XLp0CT179sRLL70EjUaDwsJC+Pj4YPPmzSguLkZYWBheeOEFKJX2uV6xJPb+/ftj4MCBmDt3LoQQmDBhgt2eLW1K3IWFhQCAVatWob6+HlOmTEF4eHiLbf39/aHT6ewSN2BZ7E0VFBSgtrZW/3xve7Ak9vr6euzZswevvvoq8vLy7BazFJnyPV6/fh0lJSUYNmwYDh06ZO8QO8SS84jUWKt+SoWpfxsAuHXrFoqLi/HII49YPQ7JtsCFCY9pHDZsGFJSUrB582aEhYUhJSUFQMOvJF26dAnTp0/Hhg0b8L///Q/ffPONPcIGYFnsN2/exM8//4wdO3Zg586dyMvLw8WLFyUTd319PYqKirBmzRq89tpr2LFjB6qqqlrdn73G7QHrxF5WVobt27fjj3/8o90u9gDLYv/6668xZMgQSZ607c3Y91hfX4+PP/4YL774oj3DMpsl5xGpsfa5xdFMKU+j48ePY9SoUTY5p0g2gZvymEZfX199N5hWq8W1a9cANLT++vTpg+7du8PNzQ0jRozQvyf12E+ePIl+/frBy8sLXl5eGDJkSJtXdo6I29/fH8OHD4e7uzu6deuGoKAgFBUVwd/f32BbnU5n18dqWhI70DAJLDExEc899xz69+9vt7gtjf3KlSs4fPgwFixYgE8++QRZWVnYu3evXeOXCmPf4/379/HTTz/hz3/+MxYsWID8/Hxs3LgRV69edUS4RllyHpEaS+un1HTkMcLZ2dkYM2aMTeKQbAI35TGNZWVl+v+fOnVK39UcGhqKqqoq/ThmXl6eXX/i0JLYNRoNLl26hLq6OtTW1uLixYsIDg6WTNwjRozQd9Xevn0bRUVF6N69O8LDw3H27FlUVlaisrISZ8+etWv3lyWx19bWYvPmzRg3bhxGjx5tt5itEfuiRYvw17/+FSkpKZg+fTrGjRuHF154we5lkAJj36O3tzc++OADpKSkICUlBf369cPrr78u2VnolpxHpMaSY1yKTH2McGFhIaqqqmzWKJDsGLibmxtmzZqF9evXo76+Ho899hh69+6Nffv2oW/fvoiIiMBXX32FU6dOwc3NDSqVCvPnzwfQ8IjH6dOn4y9/+QuEEAgJCYFWq5VF7KNGjUJeXh6WLFkCAAgPD7fb86VNiXvw4ME4e/Ys4uPjoVQqMW3aNPj6+gIAnn32WbzxxhsAgMmTJ+sno0g99qysLFy6dAl37tzRD7UsWLAAv/rVryQfO/3ClO9RTiw5j0iNsx3jph5r3333HSIjI202nMhHqRIREcmQZLvQiYiIqG1M4ERERDLEBE5ERCRDTOBEREQyxAROREQkQ0zgREREMsQETkREJEP/D3C9RQBoZG1sAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 504x216 with 2 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, ax = plt.subplots(1,2,figsize=(7,3))\n",
"ax[0].hist(states_[0], bins=100, color=\"C2\")\n",
"ax[0].set_title(r'posterior of $b$')\n",
"ax[1].hist(states_[1], bins=100, color=\"C2\")\n",
"ax[1].set_title(r'posterior of $\\sigma_y$')\n",
"fig.tight_layout()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- 真値([0.6, 0.4])とかなり近い結果を得た。\n",
"- ややacceptance rateも65%近くまで改善している。"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 5. Replica Exchange Monte Carlo"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 5.1. Define the Kernel\n",
"- `tfp.mcmc.ReplicaExchangeMC`を使う。\n",
"- 前述の通り、TFP 0.7.0時点では`SimpleStepSizeAdaptation`と併用できないため、ステップサイズ調整は無し。"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [],
"source": [
"def make_kernel_fn(target_log_prob_fn, seed):\n",
" return tfp.mcmc.HamiltonianMonteCarlo(\n",
" target_log_prob_fn=unnormalized_log_posterior,\n",
" step_size=0.01, #0.1, #0.001\n",
" num_leapfrog_steps=3, #100, #2, #10,\n",
" seed=seed\n",
" )\n",
"\n",
"\n",
"kernel3 = tfp.mcmc.ReplicaExchangeMC(\n",
" target_log_prob_fn=unnormalized_log_posterior,\n",
" inverse_temperatures=(0.5**np.linspace(0, -np.log(0.002)/np.log(2), 10).astype(np.float32)).tolist(),\n",
" make_kernel_fn=make_kernel_fn,\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- 逆温度は元記事どおり、1.0から0.02まで等比数列で10段階を指定。"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 5.2. Inference"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"WARNING:tensorflow:From /anaconda3/lib/python3.7/site-packages/tensorflow_probability/python/mcmc/replica_exchange_mc.py:385: setdiff1d (from tensorflow.python.ops.array_ops) is deprecated and will be removed after 2018-11-30.\n",
"Instructions for updating:\n",
"This op will be removed after the deprecation date. Please switch to tf.sets.difference().\n"
]
}
],
"source": [
"states, kernel_results = tfp.mcmc.sample_chain(\n",
" num_results=num_results, \n",
" num_burnin_steps=num_burnin_steps,\n",
" kernel=kernel3,\n",
" current_state=initial_state,\n",
")"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [],
"source": [
"with tf.Session() as sess:\n",
" sess.run(init_g)\n",
" sess.run(init_l)\n",
" [states_, kernel_results_ ]= sess.run([states, kernel_results])"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[array([0.5987548, 0.5949783, 0.6007985, ..., 0.5862992, 0.606583 ,\n",
" 0.5895981], dtype=float32),\n",
" array([0.47895777, 0.47838116, 0.5159441 , ..., 0.46003112, 0.51681954,\n",
" 0.44236308], dtype=float32)]"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"states_"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [],
"source": [
"# try:\n",
"# print(f'acceptance rate : {kernel_results_.inner_results.is_accepted.mean()*100} %')\n",
"# except AttributeError:\n",
"# print(f'acceptance rate : {kernel_results_.is_accepted.mean()*100} %')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 5.3. Result"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAfAAAADQCAYAAAD4dzNkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de1hUdf4H8PcwCAgDBIyiYLkqsqXhJcxIUijHcnMrulhtaZnSZTEvmOVlve2WyZouho2LStrNfTKfddmszVryQulasIgLajFI+zwZuAgDCqM1wHx/f/BjYmAYZuDM5TDv1/PwPMzhXD5nzvmez/l+z/d7UAghBIiIiEhWfNwdABERETmOCZyIiEiGmMCJiIhkiAmciIhIhpjAiYiIZIgJnIiISIaYwImIiGSICZyIyEvMmTMHGo3GrTGsXLkSkZGRUCgUeOutt9wai9wxgZNNcirwmZmZuOGGG1wXGJGTaTQazJkzR7L1vf7669i3b59k63PUV199hQ0bNmDHjh2oqqrCI4884rZY+gJfdwdA0tJoNBgyZIhkd7avv/46TCaTJOvqibYCn5ubi1tuuQWhoaFdzltUVIRx48a5MDoieTAajfDz87NZfhxdV0/odDr4+Pjgvvvu63UcxBo4dcFoNAIAQkNDERYWJsm6eqJ9gR80aBD69+/f5bxFRUUYP358j7dF5Ijk5GTMnTsXy5cvh1qtRkhICFJTU3H16lXzPE1NTVi+fDmio6Ph5+eHUaNG4S9/+Yv5719++SUSExMRHByM4OBgjB07Fp9++imA1tavzz//HG+//TYUCgUUCgWOHDliXnbr1q24/vrrERAQgJEjR2L9+vVobm42xzZv3jysXr0agwcPRnR0tHmd7VvUuovP1rqssbW+OXPmYPbs2TCZTOb9saWkpAT33nsvBg8ebJ6/7efNN9+0uazXEOQySUlJ4qmnnhLLli0TERERIjg4WMybN09cuXLFPI/RaBTLli0TUVFRol+/fuKGG24Qe/bsMf/9iy++EJMmTRIqlUqoVCoxZswYcfDgQSGEEE8++aQAYPFz+PBh87JZWVnil7/8pfD39xcxMTHilVdeEU1NTebY5s6dK1atWiUGDRok1Gq1eZ1Tp061Oz5b67LG1vqs7U9XDAaDUCqVIjs7W9x3330iKChIDBs2rFNsRFJJSkoSwcHBIjU1VZw5c0Z8+OGHYsCAAWLBggXmeZYuXSrCw8PFBx98IL799luxfv16oVAoRF5enmhubhZhYWEiPT1dlJWVibKyMrF//36Rn58vhBCivr5eTJ48WTz88MOiqqpKVFVViZ9++kkIIcTatWvFddddJ/bv3y8qKirExx9/LK699lqxatUqc2wqlUo8++yz4vTp0+I///mPEKJzebYVX/v9tLYua2ytr76+XmzZskUolUrz/nRFp9OJkJAQ8dhjj4ni4mJRWFgoxo8fLwYMGCDeffddm8t6EyZwF2KB70yqAn/8+HEBQIwcOVJ89NFHory8XCxYsED4+vqKsrKyHhwtItuSkpLE0KFDRXNzs3na9u3bhZ+fn2hsbBQGg0H4+fkJrVZrsVxKSoq4/fbbhV6v73ST3dHUqVPFk08+aTHNYDCI/v37i08++cRi+ttvvy1CQ0PNsY0cOVK0tLRYzNO+PHcXX/v9tLaujuxZ3+7du4VSqbS5HiGEmDFjhkhKShImk8k87f333xcARG1tbbfLewsmcBdigbckZYHXarVCoVCIoqIi8zSj0Sj69+8vcnJyul2eyFFJSUli5syZFtNKS0sFAHHq1Clx6tQpAUCUlJRYzJOZmSkGDhwohBAiNTVV+Pn5ienTp4sNGzaIb775xmJea+X566+/FgBEYGCgCAoKMv8EBAQIAKK6ulokJSWJRx99tFPM7cuzPfG17ae1dXVkz/rsKc96vV7069dP/PWvf7WYfuDAAQFA6PX6bmPxFnwG7mITJ06EUqk0f05MTITRaMS5c+dQXl4Oo9GIKVOmWCyTlJSE06dPIywsDKmpqbjrrrvwq1/9ChkZGfj222+73ebp06dx9epVPPjgg1CpVOafZ599FpcuXcLFixcBAPHx8fDx6fqU6C6+9rpbl6Pr605RURGmTZtm8Qy8X79+8Pf3d2snPPIuwsp/Z+74rFcIYZ62c+dO/Pvf/8a0adNw9OhR3Hjjjdi+fbvNbbSdz/v27UNxcbH5p6SkBDqdDuHh4QCAoKAgu2K2FV8be9dl7/psOXnyJJqamnDTTTdZTP/6668RExODsLAwlJSU4JZbbjH/rbS0FJMmTbL6/fdlTOBuxgLf+wIPtCbw+Ph4i2k6nQ719fWdphNJpaCgAC0tLebP//rXv+Dn54cRI0YgJiYG/v7+OHr0qMUy+fn5GD16tPnzjTfeiCVLluCTTz7BvHnzsGPHDvPf/Pz8LNYPAKNHj0ZAQAAqKioQExPT6ad9BcEWe+Ozl1Tra9tfg8Fgnnbp0iXs2rXLPKRu9OjR+O6779DU1AQAWLJkCV577TWHrxtyx2FkLtZW4NsKWfsCr1AozAWg/QlvrcC3FfrnnnsOO3bswLPPPgug+wJ/99139zj29gXUVnyuXp/RaMTp06cxbdo0i+mbN29GfHx8pzt5IqnU1tZi/vz5WLRoESoqKrB69Wo8/fTT5hvYhQsXYvXq1RgwYADGjRuHffv24e9//zv++c9/ory8HDt37sQ999yDa6+9FpWVlfjiiy8sztdhw4bh8OHDOHfuHEJDQxEaGgqVSoWVK1di5cqVAIBp06ahubkZJSUlOHnyJP74xz/aFXtgYKDN+Bwl1fomTpyIkJAQLFu2DBkZGaivr8dLL72EIUOG4MUXXwQA+Pj4IC4uDiUlJTh//jzCwsKQmJjocMxyxwTuYizw0q+vtLQURqMR+/fvx/Tp0zF8+HBkZ2fjnXfewZdffulwXET2euihhxAcHIzbbrsNRqMRM2fOxMaNG81/X79+PXx8fLB48WJcvHgRMTExeO+99zB16lRUVVVBp9Ph0UcfxcWLFxEREYEZM2Zg06ZN5uVfeOEFlJSUYOzYsTAYDDh8+DCSk5OxevVqREVFYevWrVi6dCn69++P2NhYh1/6Yiu+npBifaGhodi/fz/S09Nx0003Qa1WY+bMmXj55Zctxp8nJCTg+PHj2LFjBz788MMexSt7bnv67oXahpG19bxWqVTiqaeeEgaDwTyPrWFVlZWV4v777xfR0dHCz89PDB48WKSmpor6+nrz8ufOnROTJ08WQUFBnTq85eTkiLFjxwp/f39xzTXXiIkTJ4pt27aZY5s3b16nmHs6jMzauqzpbn32dHrJyckRUVFR4rPPPhMjRowQ/v7+IikpqVNnGiIpOXKek/Q+/PBDER4eLpYtW+buUNxGIYSXPfV3o+TkZMTExCAnJ8fdoRBRL7E8u1d5eTkSExOh0+kQEhLi7nDcgp3YiIhIdrKysvDqq696bfIG+AyciKhH2r/WlFzn3LlzmDFjBiZNmoS5c+e6Oxy3YhM6ERGRDLEJnYiISIaYwImIiGSICZyIiEiGPKYTW2VlpVPXr1arUVNT49RtOMoTYwIYlyM8JaaoqCh3h9AlbyzbHTFGaXhjjLbKNmvgREREMsQETkREJENM4ERERDLEBE5ERCRDTOBEREQy5DG90Mm1DuytBwDc88g1bo6EqG9rLWv1LGskOSZwAvBzQm/Fiw1RT/HmmFyFTeh93IG99R2SMxFJhWWL3Ik1cCIiCTCZk6uxBk5ERCRDrIETEfUCa97kLkzgXoYXGyKivoEJ3MsxoRMRyROfgRMREckQEzgREZEMMYGTVRw/TkTk2ZjAiYiIZIid2IiIutGT16OyBYucjQmciMhBTM7kCZjAiYjsxMRNnoTPwImIiGSICbwPY22BiKjvYgInIiKSIT4D9xKsjRMR9S2sgZNNfKELEZFnYgInIvp/rrxh5c0x9RYTOBERkQwxgRMRdcDaMclBt53YjEYj1q5di+bmZrS0tCAhIQEPP/wwqqursWXLFjQ2NmLYsGFYsGABfH190dTUhDfeeAMVFRUIDg7G4sWLMXDgQFfsCznRgb31Dr1Gkogs9eR1rES2dFsD79evH9auXYvXXnsNGzduRHFxMcrKyvDee+9hxowZyMrKQlBQEA4dOgQAOHToEIKCgrB161bMmDEDe/bscfpOEBEReZtuE7hCoUBAQAAAoKWlBS0tLVAoFDh9+jQSEhIAAMnJySgoKAAAFBYWIjk5GQCQkJCA0tJSCCGcFD4REZF3smscuMlkwrJly3DhwgXcddddiIyMRGBgIJRKJQAgPDwcer0eAKDX6xEREQEAUCqVCAwMRENDA0JCQizWmZeXh7y8PABARkYG1Gq1ZDtlja+vr9O34SjnxyTtMzx3f3/eeQypr+GzdZKKXQncx8cHr732GgwGAzZt2oQffvihy3mt1bYVCkWnaRqNBhqNxvy5pqbGnlB6TK1WO30bjvLEmGxxd6ye+H15SkxRUVHuDoGIXMyhXuhBQUEYNWoUdDodrly5gpaWFgCtte7w8HAAQEREBGprawG0NrlfuXIFKpVK4rCJiPoG9ninnuq2Bn758mUolUoEBQXBaDSipKQE9913H0aPHo0TJ04gMTERR44cwYQJEwAA8fHxOHLkCGJjY3HixAmMHj3aag2cnIcXA7IHR5gQyVu3Cbyurg5arRYmkwlCCNx6662Ij4/HkCFDsGXLFrz//vsYNmwY7rjjDgDAHXfcgTfeeAMLFiyASqXC4sWLnb4TROS4thEmAQEBaG5uxpo1azBu3Dh89NFHmDFjBhITE7Fjxw4cOnQId955p8UIk2PHjmHPnj1IT093924Qea1uE/jQoUOxcePGTtMjIyOxYcOGTtP9/PywZMkSaaIjIqexNcJk0aJFAFpHmOzbtw933nknCgsLMXPmTACtI0x27doFIQRb2IjchP+NjMiLcYRJR5aPn35er/MfS7l7NIMcRlQwxg7bcslWiMgjcYSJba4cYeDu0QyeMqLCFm+M0dYIEybwPoCd1qi3rI0wUSqVVkeYREREeM0IE5Yt8mT8ZyZkt47DXTj8Rd4uX74Mg8EAAOYRJtHR0eYRJgCsjjABwBEmRB6ANXAiL8URJkTyxgRO5KU4woRI3tiETkREJENM4ERERDLEBE5ERCRDTOBEREQyxAROREQkQ+yFTr3WNhb8nkeucXMkRD3D9xmQHLEGTkREJENM4ERERDLEBC5DfIUpERHxGbiMuSuJ8+aBSHoH9tazHwk5hDVwkgxbBoiIXIcJnIiISIaYwImIiGSICZyIiEiGmMCJiIhkiAmciIhIhpjAiYiIZIgJnIiISIaYwImIiGSICZyIiEiGmMBJcnwjGxGR8/Fd6DLCpEhERG1YAyci8jBsxSJ7MIETERHJEBM4ERGRDDGBExERyVC3ndhqamqg1WpRX18PhUIBjUaDu+++G42NjcjMzMTFixcxYMAApKenQ6VSQQiB3bt34+TJk/D390daWhqGDx/uin0hIiLyGt3WwJVKJWbPno3MzEysX78en376Kc6fP4/c3FzExcUhKysLcXFxyM3NBQCcPHkSFy5cQFZWFp555hnk5OQ4fSeIiIi8TbcJPCwszFyD7t+/P6Kjo6HX61FQUICkpCQAQFJSEgoKCgAAhYWFmDJlChQKBWJjY2EwGFBXV+fEXSAi6hvY+5wc4dA48Orqanz33XeIiYnBpUuXEBYWBqA1yV++fBkAoNfroVarzctERERAr9eb5yXv03ZBuueRa9wcCZElJkuSM7sT+I8//ojNmzdjzpw5CAwM7HI+IUSnaQqFotO0vLw85OXlAQAyMjIskr4z+Pr6On0bjnI8JnldbH7et/oOn3umbxxDz+HN/VuYuKkvsCuBNzc3Y/PmzZg8eTJuueUWAEBoaCjq6uoQFhaGuro6hISEAGitcdfU1JiXra2ttVr71mg00Gg05s/tl3EGtVrt9G04yhNjklLHfevtvnri9+UpMUVFRTm8TFv/luHDh+Pq1atYvnw5xowZgyNHjiAuLg4pKSnIzc1Fbm4uZs2aZdG/RafTIScnB6+++qoT9oaI7NHtM3AhBLKzsxEdHY1f//rX5ukTJkzA0aNHAQBHjx7FzTffbJ6en58PIQTKysoQGBjI5nMiD+St/VtY+6a+otsa+Lfffov8/Hxcd911ePHFFwEAv/nNb5CSkoLMzEwcOnQIarUaS5YsAQCMHz8eRUVFWLhwIfz8/JCWlubcPSCiXmP/FiL56TaBX3/99fjggw+s/m3NmjWdpikUCqSmpvY+MpI91nTkwfv6t8jnvHRl/wo59OdgjB225ZKtEJFHYv8Wz+bK/ZLD9+iNMdrq38JXqRJ5KfZvIZI31sA9GMdPkzOxf4vn4zWAbGECJ/JS7N9CJG9sQiciIpIhJnAiIiIZYhO6B+LwKyJqj8/CyRrWwImIiGSINXByGdYiyJ3YskV9DRO4DPDCQ0REHbEJnYiISIaYwImIiGSICZyIiEiGmMCJiIhkiAmciIhIhtgLnYhIJjqOSOGQTO/GGjgREZEMMYF7iAN76znem4iI7MYETkREJENM4ORybGkgV2LrFvVVTOAehhcaIiKyBxM4ERGRDDGBExHJFB8PeDcmcCIiIhliAiciIpIhJnByCzb9ERH1DhM4EZHM8YbYOzGBk0fgBYiIyDH8ZyZuxqRFREQ9wRo4ERGRDDGBExERyRATOBERkQwxgRMREclQt53Ytm3bhqKiIoSGhmLz5s0AgMbGRmRmZuLixYsYMGAA0tPToVKpIITA7t27cfLkSfj7+yMtLQ3Dhw93+k4QEdHPnWLveeQaN0dCrtBtDTw5ORkrV660mJabm4u4uDhkZWUhLi4Oubm5AICTJ0/iwoULyMrKwjPPPIOcnBznRN0HHNhbj93acneH4XYcPkZE1DPdJvBRo0ZBpVJZTCsoKEBSUhIAICkpCQUFBQCAwsJCTJkyBQqFArGxsTAYDKirq3NC2NRXMaG71rZt25CamooXXnjBPK2xsREvv/wyFi5ciJdffhmNjY0AACEEdu3ahQULFmDp0qWoqKhwV9jUDZYh79CjceCXLl1CWFgYACAsLAyXL18GAOj1eqjVavN8ERER0Ov15nnby8vLQ15eHgAgIyPDYjln8PX1dfo2HMMCZou1Y+V5x9AzY3JEcnIypk+fDq1Wa57W1sKWkpKC3Nxc5ObmYtasWRYtbDqdDjk5OXj11VfdGL1tTGLU10n6IhchRKdpCoXC6rwajQYajcb8uaamRspQOlGr1U7fBknH2rHyxGPoKTFFRUX1aLlRo0ahurraYlpBQQHWrVsHoLWFbd26dZg1a1aXLWzWbtCJyPl6lMBDQ0PNBbeurg4hISEAWmvc7S9mtbW1LNxEMtPbFjbPaV3z7hq4o9+7HFqTGGOHbfVkoQkTJuDo0aNISUnB0aNHcfPNN5unHzx4EImJidDpdAgMDGQCJ+oj7G1hY+uaZ2j7TuztmS6H79EbY7TVutZtAt+yZQvOnDmDhoYGPPfcc3j44YeRkpKCzMxMHDp0CGq1GkuWLAEAjB8/HkVFRVi4cCH8/PyQlpYm2U70FXwuR56OLWx9A681fV+3CXzx4sVWp69Zs6bTNIVCgdTU1N5HRURuwxY2InngfyMj8mJsYSOSLyZwIi/WF1vY2HRM3oLvQiciIpIhJnAiIiIZYgInIiKSISZwIiIiGWICJ4/Gf25CRGQde6ETUZ/AGz3yNkzg5JF4MSZyDntfrUqej03oREReiI+n5I81cBdhQemd1u+vnrUGoh7iNajvYQInIlnbrS13dwhEbsEmdCIiIhliAncSPl8iIiJnYgInIvJirGzIFxM4ERGRDLETm8Q63snyzlZaHMNKbVi2yNuxBk6yxGY/IvJ2rIETkazwxs012Nrl+ZjAJcKLChERuRKb0EnWeONEJI3d2nKr5YllzHMxgRMREckQm9AlwDtUIuoreD2TD9bAiYiIZIg1cCIisqljrZw90z0DEzjJHoe7eAc27RJZYgLvBV5QiMib8ebZvZjAe4CJ2zPxYtJ3scx5Fh4Pz8BObNTn8DWrRO7BsudaTOBEROQUTOjOxQROfR4vIkSegeVQWnwGbgeedETuwbLXN/A4OgcTuA086Yjcg2VPXjoeLx4/13BKAi8uLsbu3bthMpkwdepUpKSkOGMzDtmtLQfwcw/lA3vrO/VWZi/mvqWri0hPjjPPjZ95Yvkm+bC3LHWcr/VzfYfPljr+ra+XV8kTuMlkwptvvolVq1YhIiICK1aswIQJEzBkyBCpN2UXW3eCXR1k3j16B2sXCJ4LtjmrfPN79j5dvd2NtXn7KYQQQsoVlpWVYd++ffjd734HAPjb3/4GALj//vttLldZWdntuttfYK3fnRG5RlfnnbXpPa0FOFKLiIqK6tE2HNWT8m1v2SZyl3seucbhc7CrcqlWq1FTU2NzWanKtuQJ/MSJEyguLsZzzz0HAMjPz4dOp8O8efMs5svLy0NeXh4AICMjQ8oQiMhJ7CnfLNtEriH5MDJr9wMKhaLTNI1Gg4yMDJcV8OXLl7tkO47wxJgAxuUIT4zJmewp3yzbnTFGaTBGS5In8IiICNTW1po/19bWIiwsTOrNEJEbsHwTeQ7JE/iIESNQVVWF6upqNDc34/jx45gwYYLUmyEiN2D5JvIcynXr1q2TcoU+Pj4YNGgQtm7dioMHD2Ly5MlISEiQchM9Nnz4cHeH0IknxgQwLkd4YkzO4qnlWw7HgDFKgzH+TPJObEREROR8fBc6ERGRDDGBExERyZDs34Xe3Wsdjxw5gnfffRfh4eEAgOnTp2Pq1KkAgJqaGmRnZ5t71a5YsQIDBw50e1zvvfceioqKIIRAXFwcnnrqKatD8ZwRFwAcP34c+/btg0KhwNChQ7Fo0SJzzPv37wcAPPDAA0hOTnZrTP/973+xc+dOXL16FT4+PnjggQcwadIkSWLqTVxtrly5gvT0dEycOLHTexDIcfa+wvXEiRP405/+hA0bNmDEiBEeFaOtcu8pMQK2z2tPiPGtt97C6dOnAQBGoxGXLl3CW2+95VEx1tTUQKvVwmAwwGQy4bHHHsNNN90kbRBCxlpaWsTzzz8vLly4IJqamsTSpUvF999/bzHP4cOHRU5OjtXl165dK06dOiWEEOLq1avixx9/dHtc33zzjVi1apVoaWkRLS0tYuXKlaK0tNRlcVVWVooXX3xRNDQ0CCGEqK+vF0II0dDQIObPny8aGhosfndnTD/88IOorKwUQghRW1srnn76adHY2NjrmHobV5tdu3aJLVu2dHn+kf3sOR5CCHHlyhWxZs0asXLlSlFeXu5xMdq6HrmCFOe1J8TY3j/+8Q+h1WpdGKF9MWZnZ4tPP/1UCCHE999/L9LS0iSPQ9ZN6OXl5Rg0aBAiIyPh6+uLSZMmoaCgwK5lz58/j5aWFowZMwYAEBAQAH9/f7fHpVAoYDQa0dzcjKamJrS0tCA0NNRlcX3++ee46667oFKpAMC87eLiYowZMwYqlQoqlQpjxoxBcXGxW2OKiorC4MGDAQDh4eEIDQ3F5cuXex1Tb+MCgIqKCly6dAljx46VJB5vZ2+Z2rt3L+69917069fPY2N0p96e154SY3vHjh3Dbbfd5sII7YtRoVDgypUrAFpb45zxvgRZN6Hr9XpERESYP0dERECn03Wa76uvvsLZs2cxePBgPPnkk1Cr1aisrERQUBA2bdqE6upqxMXF4fHHH4ePT+/vaXoTV2xsLEaPHo1nnnkGQghMnz5dsn8EY09cbe+tXr16NUwmE2bOnIlx48Z1WjY8PBx6vd6tMbVXXl6O5uZmREZG9jqm3sZlMpnwzjvv4Pnnn0dpaakk8Xg7e47Hd999h5qaGsTHx+PAgQOuDrFX5d6TYrSnvLk7xjYXL15EdXU1brzxRleFB8C+GGfOnIlXXnkFBw8exE8//YTVq1dLHoesa+DCjtc6xsfHQ6vVYtOmTYiLi4NWqwXQ+l+Vzp49i9mzZ2PDhg343//+hyNHjrg9rgsXLuCHH35AdnY2tm/fjtLSUpw5c8ZlcZlMJlRVVWHt2rVYtGgRsrOzYTAYrK5PiufyUsRUV1eHrVu34re//a0kN2C9jeuzzz7D+PHjXXph7uu6Ox4mkwlvv/02nnjiCVeGZaE35d5VpL4GOIM9MbY5duwYEhISJCv39rInxmPHjiE5ORnZ2dlYsWIFtm7dCpPJJGkcsk7g9rzWMTg42NycptFoUFFRAaC1Bjls2DBERkZCqVRi4sSJ5r+5M66vv/4aI0eOREBAAAICAjB+/Pgu7z6dEVd4eDhuvvlm+Pr6YuDAgYiKikJVVRXCw8MtltXr9ZI0CfUmJqC1aSojIwOPPvooYmNjex2PFHGVlZXh4MGDmD9/Pt59913k5+djz549ksXmjbo7Hj/++CO+//57/P73v8f8+fOh0+mwceNGnDt3zmNiBLou954Uo63y5ikxtjl+/DgSExNdFZqZPTEeOnQIt956KwAgNjYWTU1NaGhokDQOWSdwe17rWFdXZ/69sLDQ3BwdExMDg8FgfmZaWloqWVN1b+JSq9U4e/YsWlpa0NzcjDNnziA6OtplcU2cONHc7Hv58mVUVVUhMjIS48aNw6lTp9DY2IjGxkacOnVKkma13sTU3NyMTZs2YcqUKeaCIpXexLVw4UL8+c9/hlarxezZszFlyhQ8/vjjksbnbbo7HoGBgXjzzTeh1Wqh1WoxcuRIvPTSSy7thd6bcu9JMXZ1XntSjEBrU7/BYJD0xl3KGNVqtfl7PH/+PJqamhASEiJpHLJ+Bq5UKjF37lysX78eJpMJt99+O6699lrs3bsXI0aMwIQJE/DJJ5+gsLAQSqUSKpUKaWlpAFpfCTl79mz84Q9/gBACw4cPh0ajcXtcCQkJKC0txdKlSwEA48aNk+xd0/bENXbsWJw6dQrp6enw8fHBrFmzEBwcDAB48MEHsWLFCgDAQw89ZO7k4q6Y8vPzcfbsWTQ0NJgff8yfPx+/+MUv3BoXSc+e4+FuvSn3nhSju89re4/1l19+iUmTJkk2xFbqGK7SJEoAAABPSURBVJ944gls374dH3/8MQAgLS1N8lj5KlUiIiIZknUTOhERkbdiAiciIpIhJnAiIiIZYgInIiKSISZwIiIiGWICJyIikiEmcCIiIhn6P6ANAiifUTpPAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 504x216 with 2 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, ax = plt.subplots(1,2,figsize=(7,3))\n",
"ax[0].hist(states_[0], bins=100, color=\"C2\")\n",
"ax[0].set_title(r'posterior of $b$')\n",
"ax[1].hist(states_[1], bins=100, color=\"C2\")\n",
"ax[1].set_title(r'posterior of $\\sigma_y$')\n",
"fig.tight_layout()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- こちらも真値([0.6, 0.4])に近いサンプルを得ることができた。\n",
" - stepsize調整はないが、レプリカ交換により局所解を脱したと考えられる。"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"'0.7.0'"
]
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"tfp.__version__"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 補足. ステップサイズ調整との併用試行\n",
"- `SimpleStepSizeAdaptation`と`ReplicaExchangeMC`を併用するとエラーが生じることの確認。"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [],
"source": [
"def make_kernel_fn(target_log_prob_fn, seed):\n",
" return tfp.mcmc.SimpleStepSizeAdaptation(\n",
" inner_kernel=tfp.mcmc.HamiltonianMonteCarlo(\n",
" target_log_prob_fn=target_log_prob_fn,\n",
" step_size=0.01, #0.1, #0.001\n",
" num_leapfrog_steps=3, #100, #2, #10,\n",
" seed=seed\n",
" ),\n",
" num_adaptation_steps=int(num_burnin_steps * 0.8),\n",
" adaptation_rate=0.001\n",
" )\n",
"\n",
"\n",
"kernel4 = tfp.mcmc.ReplicaExchangeMC(\n",
" target_log_prob_fn=unnormalized_log_posterior,\n",
" inverse_temperatures=(0.5**np.linspace(0, -np.log(0.002)/np.log(2), 10).astype(np.float32)).tolist(),\n",
" make_kernel_fn=make_kernel_fn,\n",
")"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"ename": "TypeError",
"evalue": "Cannot extract target_log_prob from SimpleStepSizeAdaptationResults(inner_results=MetropolisHastingsKernelResults(accepted_results=UncalibratedHamiltonianMonteCarloKernelResults(log_acceptance_correction=<tf.Tensor 'mcmc_sample_chain_3/trace_scan/while/smart_for_loop/while/remc_one_step/simple_step_size_adaptation___init___1/_one_step/mh_one_step/choose_inner_results/Select:0' shape=() dtype=float32>, target_log_prob=<tf.Tensor 'mcmc_sample_chain_3/trace_scan/while/smart_for_loop/while/remc_one_step/simple_step_size_adaptation___init___1/_one_step/mh_one_step/choose_inner_results_1/Select:0' shape=() dtype=float32>, grads_target_log_prob=[<tf.Tensor 'mcmc_sample_chain_3/trace_scan/while/smart_for_loop/while/remc_one_step/simple_step_size_adaptation___init___1/_one_step/mh_one_step/choose_inner_results_2/Select:0' shape=() dtype=float32>, <tf.Tensor 'mcmc_sample_chain_3/trace_scan/while/smart_for_loop/while/remc_one_step/simple_step_size_adaptation___init___1/_one_step/mh_one_step/choose_inner_results_2/Select_1:0' shape=() dtype=float32>], step_size=<tf.Tensor 'mcmc_sample_chain_3/trace_scan/while/smart_for_loop/while/Identity_42:0' shape=() dtype=float32>, num_leapfrog_steps=<tf.Tensor 'mcmc_sample_chain_3/trace_scan/while/smart_for_loop/while/Identity_28:0' shape=() dtype=int32>), is_accepted=<tf.Tensor 'mcmc_sample_chain_3/trace_scan/while/smart_for_loop/while/remc_one_step/simple_step_size_adaptation___init___1/_one_step/mh_one_step/Less:0' shape=() dtype=bool>, log_accept_ratio=<tf.Tensor 'mcmc_sample_chain_3/trace_scan/while/smart_for_loop/while/remc_one_step/simple_step_size_adaptation___init___1/_one_step/mh_one_step/compute_log_accept_ratio/Select:0' shape=() dtype=float32>, proposed_state=[<tf.Tensor 'mcmc_sample_chain_3/trace_scan/while/smart_for_loop/while/remc_one_step/simple_step_size_adaptation___init___1/_one_step/mh_one_step/hmc_kernel_one_step/leapfrog_integrate/while/Exit_3:0' shape=() dtype=float32>, <tf.Tensor 'mcmc_sample_chain_3/trace_scan/while/smart_for_loop/while/remc_one_step/simple_step_size_adaptation___init___1/_one_step/mh_one_step/hmc_kernel_one_step/leapfrog_integrate/while/Exit_4:0' shape=() dtype=float32>], proposed_results=UncalibratedHamiltonianMonteCarloKernelResults(log_acceptance_correction=<tf.Tensor 'mcmc_sample_chain_3/trace_scan/while/smart_for_loop/while/remc_one_step/simple_step_size_adaptation___init___1/_one_step/mh_one_step/hmc_kernel_one_step/compute_log_acceptance_correction/safe_sum/Select:0' shape=() dtype=float32>, target_log_prob=<tf.Tensor 'mcmc_sample_chain_3/trace_scan/while/smart_for_loop/while/remc_one_step/simple_step_size_adaptation___init___1/_one_step/mh_one_step/hmc_kernel_one_step/leapfrog_integrate/while/Exit_5:0' shape=() dtype=float32>, grads_target_log_prob=[<tf.Tensor 'mcmc_sample_chain_3/trace_scan/while/smart_for_loop/while/remc_one_step/simple_step_size_adaptation___init___1/_one_step/mh_one_step/hmc_kernel_one_step/leapfrog_integrate/while/Exit_6:0' shape=() dtype=float32>, <tf.Tensor 'mcmc_sample_chain_3/trace_scan/while/smart_for_loop/while/remc_one_step/simple_step_size_adaptation___init___1/_one_step/mh_one_step/hmc_kernel_one_step/leapfrog_integrate/while/Exit_7:0' shape=() dtype=float32>], step_size=<tf.Tensor 'mcmc_sample_chain_3/trace_scan/while/smart_for_loop/while/Identity_42:0' shape=() dtype=float32>, num_leapfrog_steps=<tf.Tensor 'mcmc_sample_chain_3/trace_scan/while/smart_for_loop/while/Identity_28:0' shape=() dtype=int32>), extra=[]), target_accept_prob=<tf.Tensor 'mcmc_sample_chain_3/trace_scan/while/smart_for_loop/while/Identity_39:0' shape=() dtype=float32>, adaptation_rate=<tf.Tensor 'mcmc_sample_chain_3/trace_scan/while/smart_for_loop/while/Identity_40:0' shape=() dtype=float32>, step=<tf.Tensor 'mcmc_sample_chain_3/trace_scan/while/smart_for_loop/while/remc_one_step/simple_step_size_adaptation___init___1/_one_step/add_2:0' shape=() dtype=int32>, new_step_size=<tf.Tensor 'mcmc_sample_chain_3/trace_scan/while/smart_for_loop/while/remc_one_step/simple_step_size_adaptation___init___1/_one_step/Select:0' shape=() dtype=float32>)",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mTypeError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m<ipython-input-30-e3866d51e9fa>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[0mnum_burnin_steps\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mnum_burnin_steps\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 4\u001b[0m \u001b[0mkernel\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mkernel4\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 5\u001b[0;31m \u001b[0mcurrent_state\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0minitial_state\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 6\u001b[0m )\n",
"\u001b[0;32m/anaconda3/lib/python3.7/site-packages/tensorflow_probability/python/mcmc/sample.py\u001b[0m in \u001b[0;36msample_chain\u001b[0;34m(num_results, current_state, previous_kernel_results, kernel, num_burnin_steps, num_steps_between_results, trace_fn, return_final_kernel_results, parallel_iterations, name)\u001b[0m\n\u001b[1;32m 359\u001b[0m trace_fn(*state_and_results)),\n\u001b[1;32m 360\u001b[0m \u001b[0;31m# pylint: enable=g-long-lambda\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 361\u001b[0;31m parallel_iterations=parallel_iterations)\n\u001b[0m\u001b[1;32m 362\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 363\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mreturn_final_kernel_results\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/anaconda3/lib/python3.7/site-packages/tensorflow_probability/python/mcmc/internal/util.py\u001b[0m in \u001b[0;36mtrace_scan\u001b[0;34m(loop_fn, initial_state, elems, trace_fn, parallel_iterations, name)\u001b[0m\n\u001b[1;32m 368\u001b[0m \u001b[0mbody\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0m_body\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 369\u001b[0m \u001b[0mloop_vars\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0minitial_state\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtrace_arrays\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 370\u001b[0;31m parallel_iterations=parallel_iterations)\n\u001b[0m\u001b[1;32m 371\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 372\u001b[0m \u001b[0mstacked_trace\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtf\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnest\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmap_structure\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;32mlambda\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstack\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtrace_arrays\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/anaconda3/lib/python3.7/site-packages/tensorflow/python/ops/control_flow_ops.py\u001b[0m in \u001b[0;36mwhile_loop\u001b[0;34m(cond, body, loop_vars, shape_invariants, parallel_iterations, back_prop, swap_memory, name, maximum_iterations, return_same_structure)\u001b[0m\n\u001b[1;32m 3499\u001b[0m \u001b[0mops\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0madd_to_collection\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mops\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mGraphKeys\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mWHILE_CONTEXT\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mloop_context\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3500\u001b[0m result = loop_context.BuildLoop(cond, body, loop_vars, shape_invariants,\n\u001b[0;32m-> 3501\u001b[0;31m return_same_structure)\n\u001b[0m\u001b[1;32m 3502\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mmaximum_iterations\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3503\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mresult\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/anaconda3/lib/python3.7/site-packages/tensorflow/python/ops/control_flow_ops.py\u001b[0m in \u001b[0;36mBuildLoop\u001b[0;34m(self, pred, body, loop_vars, shape_invariants, return_same_structure)\u001b[0m\n\u001b[1;32m 3010\u001b[0m \u001b[0;32mwith\u001b[0m \u001b[0mops\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_default_graph\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_mutation_lock\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0;31m# pylint: disable=protected-access\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3011\u001b[0m original_body_result, exit_vars = self._BuildLoop(\n\u001b[0;32m-> 3012\u001b[0;31m pred, body, original_loop_vars, loop_vars, shape_invariants)\n\u001b[0m\u001b[1;32m 3013\u001b[0m \u001b[0;32mfinally\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3014\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mExit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/anaconda3/lib/python3.7/site-packages/tensorflow/python/ops/control_flow_ops.py\u001b[0m in \u001b[0;36m_BuildLoop\u001b[0;34m(self, pred, body, original_loop_vars, loop_vars, shape_invariants)\u001b[0m\n\u001b[1;32m 2935\u001b[0m expand_composites=True)\n\u001b[1;32m 2936\u001b[0m \u001b[0mpre_summaries\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mops\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_collection\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mops\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mGraphKeys\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_SUMMARY_COLLECTION\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;31m# pylint: disable=protected-access\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2937\u001b[0;31m \u001b[0mbody_result\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mbody\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mpacked_vars_for_body\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2938\u001b[0m \u001b[0mpost_summaries\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mops\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_collection\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mops\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mGraphKeys\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_SUMMARY_COLLECTION\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;31m# pylint: disable=protected-access\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2939\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mnest\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mis_sequence_or_composite\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mbody_result\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/anaconda3/lib/python3.7/site-packages/tensorflow_probability/python/mcmc/internal/util.py\u001b[0m in \u001b[0;36m_body\u001b[0;34m(i, state, trace_arrays)\u001b[0m\n\u001b[1;32m 357\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 358\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m_body\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstate\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtrace_arrays\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 359\u001b[0;31m \u001b[0mstate\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mloop_fn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mstate\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0melems_array\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mread\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 360\u001b[0m trace_arrays = tf.nest.pack_sequence_as(trace_arrays, [\n\u001b[1;32m 361\u001b[0m a.write(i, v) for a, v in zip(\n",
"\u001b[0;32m/anaconda3/lib/python3.7/site-packages/tensorflow_probability/python/mcmc/sample.py\u001b[0m in \u001b[0;36m_trace_scan_fn\u001b[0;34m(state_and_results, num_steps)\u001b[0m\n\u001b[1;32m 343\u001b[0m \u001b[0mbody_fn\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mkernel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mone_step\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 344\u001b[0m \u001b[0minitial_loop_vars\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mlist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mstate_and_results\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 345\u001b[0;31m parallel_iterations=parallel_iterations)\n\u001b[0m\u001b[1;32m 346\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mnext_state\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcurrent_kernel_results\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 347\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/anaconda3/lib/python3.7/site-packages/tensorflow_probability/python/mcmc/internal/util.py\u001b[0m in \u001b[0;36msmart_for_loop\u001b[0;34m(loop_num_iter, body_fn, initial_loop_vars, parallel_iterations, name)\u001b[0m\n\u001b[1;32m 284\u001b[0m \u001b[0mbody\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mlambda\u001b[0m \u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mlist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mbody_fn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 285\u001b[0m \u001b[0mloop_vars\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mint32\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0minitial_loop_vars\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 286\u001b[0;31m \u001b[0mparallel_iterations\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mparallel_iterations\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 287\u001b[0m )[1:]\n\u001b[1;32m 288\u001b[0m \u001b[0mresult\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0minitial_loop_vars\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/anaconda3/lib/python3.7/site-packages/tensorflow/python/ops/control_flow_ops.py\u001b[0m in \u001b[0;36mwhile_loop\u001b[0;34m(cond, body, loop_vars, shape_invariants, parallel_iterations, back_prop, swap_memory, name, maximum_iterations, return_same_structure)\u001b[0m\n\u001b[1;32m 3499\u001b[0m \u001b[0mops\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0madd_to_collection\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mops\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mGraphKeys\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mWHILE_CONTEXT\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mloop_context\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3500\u001b[0m result = loop_context.BuildLoop(cond, body, loop_vars, shape_invariants,\n\u001b[0;32m-> 3501\u001b[0;31m return_same_structure)\n\u001b[0m\u001b[1;32m 3502\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mmaximum_iterations\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3503\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mresult\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/anaconda3/lib/python3.7/site-packages/tensorflow/python/ops/control_flow_ops.py\u001b[0m in \u001b[0;36mBuildLoop\u001b[0;34m(self, pred, body, loop_vars, shape_invariants, return_same_structure)\u001b[0m\n\u001b[1;32m 3010\u001b[0m \u001b[0;32mwith\u001b[0m \u001b[0mops\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_default_graph\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_mutation_lock\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0;31m# pylint: disable=protected-access\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3011\u001b[0m original_body_result, exit_vars = self._BuildLoop(\n\u001b[0;32m-> 3012\u001b[0;31m pred, body, original_loop_vars, loop_vars, shape_invariants)\n\u001b[0m\u001b[1;32m 3013\u001b[0m \u001b[0;32mfinally\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3014\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mExit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/anaconda3/lib/python3.7/site-packages/tensorflow/python/ops/control_flow_ops.py\u001b[0m in \u001b[0;36m_BuildLoop\u001b[0;34m(self, pred, body, original_loop_vars, loop_vars, shape_invariants)\u001b[0m\n\u001b[1;32m 2935\u001b[0m expand_composites=True)\n\u001b[1;32m 2936\u001b[0m \u001b[0mpre_summaries\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mops\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_collection\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mops\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mGraphKeys\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_SUMMARY_COLLECTION\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;31m# pylint: disable=protected-access\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2937\u001b[0;31m \u001b[0mbody_result\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mbody\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mpacked_vars_for_body\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2938\u001b[0m \u001b[0mpost_summaries\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mops\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_collection\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mops\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mGraphKeys\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_SUMMARY_COLLECTION\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;31m# pylint: disable=protected-access\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2939\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mnest\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mis_sequence_or_composite\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mbody_result\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/anaconda3/lib/python3.7/site-packages/tensorflow_probability/python/mcmc/internal/util.py\u001b[0m in \u001b[0;36m<lambda>\u001b[0;34m(i, *args)\u001b[0m\n\u001b[1;32m 282\u001b[0m return tf.while_loop(\n\u001b[1;32m 283\u001b[0m \u001b[0mcond\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mlambda\u001b[0m \u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mi\u001b[0m \u001b[0;34m<\u001b[0m \u001b[0mloop_num_iter\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 284\u001b[0;31m \u001b[0mbody\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mlambda\u001b[0m \u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mlist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mbody_fn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 285\u001b[0m \u001b[0mloop_vars\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mint32\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0minitial_loop_vars\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 286\u001b[0m \u001b[0mparallel_iterations\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mparallel_iterations\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/anaconda3/lib/python3.7/site-packages/tensorflow_probability/python/mcmc/replica_exchange_mc.py\u001b[0m in \u001b[0;36mone_step\u001b[0;34m(self, current_state, previous_kernel_results)\u001b[0m\n\u001b[1;32m 380\u001b[0m exchanged_states = self._get_exchanged_states(\n\u001b[1;32m 381\u001b[0m \u001b[0mold_states\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mexchange_proposed\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mexchange_proposed_n\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 382\u001b[0;31m sampled_replica_states, sampled_replica_results)\n\u001b[0m\u001b[1;32m 383\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 384\u001b[0m no_exchange_proposed, _ = tf.compat.v1.setdiff1d(\n",
"\u001b[0;32m/anaconda3/lib/python3.7/site-packages/tensorflow_probability/python/mcmc/replica_exchange_mc.py\u001b[0m in \u001b[0;36m_get_exchanged_states\u001b[0;34m(self, old_states, exchange_proposed, exchange_proposed_n, sampled_replica_states, sampled_replica_results)\u001b[0m\n\u001b[1;32m 426\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mreplica\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnum_replica\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 427\u001b[0m replica_log_prob = _get_field(sampled_replica_results[replica],\n\u001b[0;32m--> 428\u001b[0;31m 'target_log_prob')\n\u001b[0m\u001b[1;32m 429\u001b[0m \u001b[0minverse_temp\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minverse_temperatures\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mreplica\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 430\u001b[0m \u001b[0mtarget_log_probs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mreplica_log_prob\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0minverse_temp\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/anaconda3/lib/python3.7/site-packages/tensorflow_probability/python/mcmc/replica_exchange_mc.py\u001b[0m in \u001b[0;36m_get_field\u001b[0;34m(kernel_results, field_name)\u001b[0m\n\u001b[1;32m 573\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mhasattr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkernel_results\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'accepted_results'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 574\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mgetattr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkernel_results\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0maccepted_results\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfield_name\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 575\u001b[0;31m \u001b[0;32mraise\u001b[0m \u001b[0mTypeError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'Cannot extract %s from %s'\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mfield_name\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mkernel_results\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
"\u001b[0;31mTypeError\u001b[0m: Cannot extract target_log_prob from SimpleStepSizeAdaptationResults(inner_results=MetropolisHastingsKernelResults(accepted_results=UncalibratedHamiltonianMonteCarloKernelResults(log_acceptance_correction=<tf.Tensor 'mcmc_sample_chain_3/trace_scan/while/smart_for_loop/while/remc_one_step/simple_step_size_adaptation___init___1/_one_step/mh_one_step/choose_inner_results/Select:0' shape=() dtype=float32>, target_log_prob=<tf.Tensor 'mcmc_sample_chain_3/trace_scan/while/smart_for_loop/while/remc_one_step/simple_step_size_adaptation___init___1/_one_step/mh_one_step/choose_inner_results_1/Select:0' shape=() dtype=float32>, grads_target_log_prob=[<tf.Tensor 'mcmc_sample_chain_3/trace_scan/while/smart_for_loop/while/remc_one_step/simple_step_size_adaptation___init___1/_one_step/mh_one_step/choose_inner_results_2/Select:0' shape=() dtype=float32>, <tf.Tensor 'mcmc_sample_chain_3/trace_scan/while/smart_for_loop/while/remc_one_step/simple_step_size_adaptation___init___1/_one_step/mh_one_step/choose_inner_results_2/Select_1:0' shape=() dtype=float32>], step_size=<tf.Tensor 'mcmc_sample_chain_3/trace_scan/while/smart_for_loop/while/Identity_42:0' shape=() dtype=float32>, num_leapfrog_steps=<tf.Tensor 'mcmc_sample_chain_3/trace_scan/while/smart_for_loop/while/Identity_28:0' shape=() dtype=int32>), is_accepted=<tf.Tensor 'mcmc_sample_chain_3/trace_scan/while/smart_for_loop/while/remc_one_step/simple_step_size_adaptation___init___1/_one_step/mh_one_step/Less:0' shape=() dtype=bool>, log_accept_ratio=<tf.Tensor 'mcmc_sample_chain_3/trace_scan/while/smart_for_loop/while/remc_one_step/simple_step_size_adaptation___init___1/_one_step/mh_one_step/compute_log_accept_ratio/Select:0' shape=() dtype=float32>, proposed_state=[<tf.Tensor 'mcmc_sample_chain_3/trace_scan/while/smart_for_loop/while/remc_one_step/simple_step_size_adaptation___init___1/_one_step/mh_one_step/hmc_kernel_one_step/leapfrog_integrate/while/Exit_3:0' shape=() dtype=float32>, <tf.Tensor 'mcmc_sample_chain_3/trace_scan/while/smart_for_loop/while/remc_one_step/simple_step_size_adaptation___init___1/_one_step/mh_one_step/hmc_kernel_one_step/leapfrog_integrate/while/Exit_4:0' shape=() dtype=float32>], proposed_results=UncalibratedHamiltonianMonteCarloKernelResults(log_acceptance_correction=<tf.Tensor 'mcmc_sample_chain_3/trace_scan/while/smart_for_loop/while/remc_one_step/simple_step_size_adaptation___init___1/_one_step/mh_one_step/hmc_kernel_one_step/compute_log_acceptance_correction/safe_sum/Select:0' shape=() dtype=float32>, target_log_prob=<tf.Tensor 'mcmc_sample_chain_3/trace_scan/while/smart_for_loop/while/remc_one_step/simple_step_size_adaptation___init___1/_one_step/mh_one_step/hmc_kernel_one_step/leapfrog_integrate/while/Exit_5:0' shape=() dtype=float32>, grads_target_log_prob=[<tf.Tensor 'mcmc_sample_chain_3/trace_scan/while/smart_for_loop/while/remc_one_step/simple_step_size_adaptation___init___1/_one_step/mh_one_step/hmc_kernel_one_step/leapfrog_integrate/while/Exit_6:0' shape=() dtype=float32>, <tf.Tensor 'mcmc_sample_chain_3/trace_scan/while/smart_for_loop/while/remc_one_step/simple_step_size_adaptation___init___1/_one_step/mh_one_step/hmc_kernel_one_step/leapfrog_integrate/while/Exit_7:0' shape=() dtype=float32>], step_size=<tf.Tensor 'mcmc_sample_chain_3/trace_scan/while/smart_for_loop/while/Identity_42:0' shape=() dtype=float32>, num_leapfrog_steps=<tf.Tensor 'mcmc_sample_chain_3/trace_scan/while/smart_for_loop/while/Identity_28:0' shape=() dtype=int32>), extra=[]), target_accept_prob=<tf.Tensor 'mcmc_sample_chain_3/trace_scan/while/smart_for_loop/while/Identity_39:0' shape=() dtype=float32>, adaptation_rate=<tf.Tensor 'mcmc_sample_chain_3/trace_scan/while/smart_for_loop/while/Identity_40:0' shape=() dtype=float32>, step=<tf.Tensor 'mcmc_sample_chain_3/trace_scan/while/smart_for_loop/while/remc_one_step/simple_step_size_adaptation___init___1/_one_step/add_2:0' shape=() dtype=int32>, new_step_size=<tf.Tensor 'mcmc_sample_chain_3/trace_scan/while/smart_for_loop/while/remc_one_step/simple_step_size_adaptation___init___1/_one_step/Select:0' shape=() dtype=float32>)"
]
}
],
"source": [
"states, kernel_results = tfp.mcmc.sample_chain(\n",
" num_results=num_results, \n",
" num_burnin_steps=num_burnin_steps,\n",
" kernel=kernel4,\n",
" current_state=initial_state,\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- ダメでした。\n",
" - この辺りの改善は今後のTFPに期待。"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment