narrowlyapplicable/tfp-replica_exchange.ipynb

## tfp-replica_exchange.ipynb
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tensorflow Probabilityによるレプリカ交換モンテカルロ試行\n",
    "- `tfp.mcmc.ReplicaExchangeMC`を使用したが、例題を再現できなかった。"
   ]
  },
  {
   "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": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "initial_state = [24.17, 0.4048] #[0.6, 0.4] #b=24.17, s_y=0.4048"
   ]
  },
  {
   "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 : 96.66 %\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+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deVxU9f4/8NcMiwIDODC4gEsiel0iN9wwBRXL1IrrVctSc0sNr6ZcTfOb2uOhJFclFKNQKcush8vjGteuNy1EISUTwwWxFFB7ZGIswyZobJ/fH/6YC7INMMs5zOv5ePh4yOGcw2tmzue855zzOZ+jEEIIEBERkawozR2AiIiImo4FnIiISIZYwImIiGSIBZyIiEiGWMCJiIhkiAWciIhIhljAiYiIZIgFnIjIAsyZMwcBAQFmzbB27Vp06NABCoUCn376qVmztAYs4NQgOTX68PBw9OnTx3TBiIwoICAAc+bMMdj6duzYgcOHDxtsfU31448/YvPmzdi9ezcyMzPx0ksvmS1La2Ft7gBkWAEBAejcubPBvt3u2LEDlZWVBllXc1Q1+piYGAwbNgzOzs71zpucnIwBAwaYMB2R9JWWlsLW1rbBttPUdTVHWloalEolXnzxxRbnoEd4BE51Ki0tBQA4OztDrVYbZF3NUb3Rd+zYEXZ2dvXOm5ycjIEDBzb7bxHpy9/fH/PmzcOaNWug0Wjg5OSEBQsW4MGDB7p5ysrKsGbNGnh4eMDW1hZ9+/bFl19+qfv9mTNnMHLkSDg6OsLR0RH9+/fHiRMnADw683Xy5El89tlnUCgUUCgUOH36tG7ZnTt3onfv3mjbti169uyJkJAQlJeX67LNnz8f69atQ6dOneDh4aFbZ/WzaY3la2hddWlofXPmzMGsWbNQWVmpez0NSUlJwQsvvIBOnTrp5q/69/HHHze4rEURZDJ+fn5i7ty5YvXq1cLV1VU4OjqK+fPni5KSEt08paWlYvXq1cLd3V3Y2NiIPn36iC+++EL3+++//174+voKlUolVCqVeOqpp8Tx48eFEEK89tprAkCNf6dOndItGxERIf7yl7+INm3aCC8vL7Fp0yZRVlamyzZv3jzxzjvviI4dOwqNRqNb57hx4/TO19C66tLQ+up6PfUpLi4WVlZWIioqSrz44ovCwcFBdO/evVY2IkPw8/MTjo6OYsGCBeLatWvi6NGjws3NTSxdulQ3z8qVK4WLi4s4dOiQuH79uggJCREKhULExsaK8vJyoVarxYoVK8SNGzfEjRs3xJEjR0RCQoIQQoj8/HwxatQoMX36dJGZmSkyMzPFn3/+KYQQYsOGDaJr167iyJEj4ubNm+LYsWOiS5cu4p133tFlU6lUYtGiRSI1NVVcuXJFCFG7LTeUr/rrrGtddWloffn5+WL79u3CyspK93rqk5aWJpycnMQrr7wiLl26JC5cuCAGDhwo3NzcxOeff97gspaGBdyE2OhrM1SjT0xMFABEz549xX/+8x+Rnp4uli5dKqytrcWNGzea8WkR1c/Pz09069ZNlJeX66bt2rVL2Nraivv374vi4mJha2srIiMjaywXGBgoxowZI7Raba0v2I8bN26ceO2112pMKy4uFnZ2duKbb76pMf2zzz4Tzs7Oumw9e/YUFRUVNeap3pYby1f9dda1rsfps769e/cKKyurBtcjhBCTJk0Sfn5+orKyUjftwIEDAoDIzc1tdHlLwgJuQmz0NRmy0UdGRgqFQiGSk5N100pLS4WdnZ2Ijo5udHmipvDz8xPTpk2rMe3q1asCgLh8+bK4fPmyACBSUlJqzBMeHi7at28vhBBiwYIFwtbWVkyYMEFs3rxZ/PLLLzXmrastnz9/XgAQ9vb2wsHBQfevbdu2AoDIysoSfn5+4uWXX66VuXpb1idf1eusa12P02d9+rRlrVYrbGxsxL/+9a8a07/++msBQGi12kazWBJeAzexoUOHwsrKSvfzyJEjUVpaioyMDKSnp6O0tBSjR4+usYyfnx9SU1OhVquxYMECPPvss3juuecQGhqK69evN/o3U1NT8eDBA/ztb3+DSqXS/Vu0aBEKCgqQnZ0NABg8eDCUyvo3icbyVdfYupq6vsYkJydj/PjxNa6B29jYoE2bNmbthEeWQ9TxZObHr/UKIXTT9uzZg59++gnjx49HfHw8nnzySezatavBv1G1LR8+fBiXLl3S/UtJSUFaWhpcXFwAAA4ODnplbihfFX3Xpe/6GnLx4kWUlZVh0KBBNaafP38eXl5eUKvVSElJwbBhw3S/u3r1Knx9fet8/1s7FnAzY6NveaMHHhXwwYMH15iWlpaG/Pz8WtOJDCEpKQkVFRW6n3/44QfY2tqiR48e8PLyQps2bRAfH19jmYSEBPTr10/385NPPong4GB88803mD9/Pnbv3q37na2tbY31A0C/fv3Qtm1b3Lx5E15eXrX+VT84aIi++fRlqPVVvd7i4mLdtIKCAnzyySe6W+r69euHW7duoaysDAAQHByMrVu3Nnmf0RrwNjITq2r0VQ2teqNXKBS6RlB9o6+r0Vc1/MWLF2P37t1YtGgRgMYb/cSJE5udvXojbSifqddXWlqK1NRUjB8/vsb0sLAwDB48uNa3eSJDyM3NxZIlS/Dmm2/i5s2bWLduHV5//XXdl9dly5Zh3bp1cHNzw4ABA3D48GH8+9//xnfffYf09HTs2bMHzz//PLp06YK7d+/i+++/r7Gtdu/eHadOnUJGRgacnZ3h7OwMlUqFtWvXYu3atQCA8ePHo7y8HCkpKbh48SL++c9/6pXd3t6+wXxNZaj1DR06FE5OTli9ejVCQ0ORn5+Pt956C507d8aqVasAAEqlEt7e3khJScGdO3egVqsxcuTIJmduDVjATYyN3vDru3r1KkpLS3HkyBFMmDABnp6eiIqKwr59+3DmzJkm5yLSx9SpU+Ho6Iinn34apaWlmDZtGrZs2aL7fUhICJRKJZYvX47s7Gx4eXlh//79GDduHDIzM5GWloaXX34Z2dnZcHV1xaRJk7Bt2zbd8v/4xz+QkpKC/v37o7i4GKdOnYK/vz/WrVsHd3d37Ny5EytXroSdnR169erV5EFfGsrXHIZYn7OzM44cOYIVK1Zg0KBB0Gg0mDZtGjZu3Fjj/vPhw4cjMTERu3fvxtGjR5uVt1Uw29V3C1R1G1lVz2uVSiXmzp0riouLdfM0dFvV3bt3xV//+lfh4eEhbG1tRadOncSCBQtEfn6+bvmMjAwxatQo4eDgUKvDW3R0tOjfv79o06aNaNeunRg6dKj48MMPddnmz59fK3NzbyOra111aWx9+nR8iY6OFu7u7uLbb78VPXr0EG3atBF+fn61OtQQGUpTtnEyvKNHjwoXFxexevVqc0cxK4UQFnjl30z8/f3h5eWF6Ohoc0chohZgWzav9PR0jBw5EmlpaXBycjJ3HLNhJzYiIpKViIgIvPfeexZdvAFeAyciarLqw5qS6WRkZGDSpEnw9fXFvHnzzB3H7HgKnYiISIZ4Cp2IiEiGWMCJiIhkiAWciIhIhiTTie3u3btNXkaj0SAnJ8cIaVpGirmYSX9SzKVPJnd3dxOlaZrmtG1Dk+JnCkgzlxQzAdLMZapM9bVtHoETERHJEAs4ERGRDLGAExERyRALOBERkQyxgBMREckQCzgZxNcH8/H1wXxzxyAym8baANsHGRoLODULCzYRkXmxgBMREckQCzgREZEMsYBTg3iqnIhImiQzlCq1PlWF//mX2pk5CZH58AswGQuPwImIiGSIR+DUJDyaICKSBh6BExERyRCPwImIDKg5Z6nYX4Sag0fgREREMsQCTkREJEONnkIvLS3Fhg0bUF5ejoqKCgwfPhzTp09HVlYWtm/fjvv376N79+5YunQprK2tUVZWhg8++AA3b96Eo6Mjli9fjvbt25vitZCEsLMbUeP2RqabOwLJWKNH4DY2NtiwYQO2bt2KLVu24NKlS7hx4wb279+PSZMmISIiAg4ODoiLiwMAxMXFwcHBATt37sSkSZPwxRdfGP1FEBERWZpGC7hCoUDbtm0BABUVFaioqIBCoUBqaiqGDx8OAPD390dSUhIA4MKFC/D39wcADB8+HFevXoUQwkjxSWo4chsRkWno1Qu9srISq1evxr179/Dss8+iQ4cOsLe3h5WVFQDAxcUFWq0WAKDVauHq6goAsLKygr29PYqKiuDk5FRjnbGxsYiNjQUAhIaGQqPRND28tXWzljM2KeZqfqZHxfh/y9YszvVNr3seQ2UyLinmkmImIjIvvQq4UqnE1q1bUVxcjG3btuH333+vd966jrYVCkWtaQEBAQgICND9nJOTo0+UGjQaTbOWMzYp5mpppvqu1emzzvrmkeL7BEgzlz6Z3N3dTZSGWopnqcgQmnQfuIODA/r27Yu0tDSUlJSgoqICVlZW0Gq1cHFxAQC4uroiNzcXrq6uqKioQElJCVQqlVHCExHJCQs3GVKj18ALCwtRXFwM4FGP9JSUFHh4eKBfv344d+4cAOD06dPw8fEBAAwePBinT58GAJw7dw79+vWr8wiciIiImq/RI/C8vDxERkaisrISQgiMGDECgwcPRufOnbF9+3YcOHAA3bt3x9ixYwEAY8eOxQcffIClS5dCpVJh+fLlRn8RRERElqbRAt6tWzds2bKl1vQOHTpg8+bNtabb2toiODjYMOmIiCTm8WFPeVqczIVjoRNZKA7SRCRvHEqVyEJxkCYieWMBJ7JQHKSJSN54Cp3IghljkCZLwWvfZG4s4EQWzBiDNBlilEVDM+xIdsYr3K3vvTIcKeYydyYWcCIy6CBNhhhl0dCkOLpeXaSQUarvlRRzmSpTfaMs8ho4GR0fcCJNHKSJSN54BE4twsIsXxykiUjeWMCJLBQHaSKSN55CJyIikiEegRMR6cEUl4seH6aVqCE8AieTY6c2IqKW4xE4mQyLNhGR4bCAUw0sskRE8tBoAc/JyUFkZCTy8/OhUCgQEBCAiRMn4tChQzh58qRuGMUZM2Zg0KBBAICvvvoKcXFxUCqVmDt3LgYMGGDcV0FERGRhGi3gVlZWmDVrFjw9PfHgwQOsWbMGTz31FABg0qRJeOGFF2rMf+fOHSQmJuL9999HXl4eNm7ciB07dkCp5OV2IiIiQ2m0qqrVanh6egIA7Ozs4OHhoXu4QV2SkpLg6+sLGxsbtG/fHh07dkR6errhEhMRtXLs6En6aNI18KysLNy6dQteXl745ZdfcOLECSQkJMDT0xOzZ8+GSqWCVqtFz549dctUf5pRdYZ44IG5B5KvjxRz6Z/JdDsNKb5PgDRzSTETEZmX3gX84cOHCAsLw5w5c2Bvb49nnnkGU6dOBQAcPHgQ+/btQ1BQkN7PBzbEAw+kOLg9IM1cUsxUXl4uuUyANN8rfTLV98ADImqd9LowXV5ejrCwMIwaNQrDhg0DALRr1w5KpRJKpRLjxo1DRkYGgP89sahK9acZETWEpw2JiPTX6BG4EAJRUVHw8PDA5MmTddPz8vKgVqsBAOfPn0eXLl0AAD4+PoiIiMDkyZORl5eHzMxMeHl5GSk+EZFxSOHLJEdmo4Y0WsCvX7+OhIQEdO3aFatWrQLw6Jaxs2fP4vbt21AoFHBzc8PChQsBAF26dMGIESMQHBwMpVKJ+fPnswe6DEhhZ0VERPprtID37t0bhw4dqjW96p7vukyZMgVTpkxpWTJq9fZGPro7gUcXRERNx0NjIiIiGWIBJ54+JyKSIY6FbsFYuInqxrZBcsAjcCIiIhliASciIpIhnkIns+PpSiKipuMROBGRxHGUQqoLCzgREZEMsYATERHJEAs4EZFMPH4qnafVLRsLuAWR23U0ueUlIjIlFnAiIiIZYgEnIiKSoUbvA8/JyUFkZCTy8/OhUCgQEBCAiRMn4v79+wgPD0d2djbc3NywYsUKqFQqCCGwd+9eXLx4EW3atEFQUBA8PT1N8VqIiCwCLy0RoMcRuJWVFWbNmoXw8HCEhITgxIkTuHPnDmJiYuDt7Y2IiAh4e3sjJiYGAHDx4kXcu3cPERERWLhwIaKjo43+IoiIiCxNowVcrVbrjqDt7Ozg4eEBrVaLpKQk+Pn5AQD8/PyQlJQEALhw4QJGjx4NhUKBXr16obi4GHl5eUZ8CdRUcuscJre8RESm0KShVLOysnDr1i14eXmhoKAAarUawKMiX1hYCADQarXQaDS6ZVxdXaHVanXzVomNjUVsbCwAIDQ0tMYyeoe3tm7WcsYmxVzW1vIZNbe+Ym2q91Sqn5/UMpE0VLWX519qZ+YkZGp679UfPnyIsLAwzJkzB/b29vXOJ4SoNU2hUNSaFhAQgICAAN3POTk5+kbR0Wg0zVrO2KSYqzXs/E31nkr182ssk7u7u4nSEJEU6NULvby8HGFhYRg1ahSGDRsGAHB2dtadGs/Ly4OTkxOAR0fc1Xc0ubm5tY6+iYiIqGUaPQIXQiAqKgoeHh6YPHmybrqPjw/i4+MRGBiI+Ph4DBkyRDf9+PHjGDlyJNLS0mBvb88CTiRBvMOESN4aLeDXr19HQkICunbtilWrVgEAZsyYgcDAQISHhyMuLg4ajQbBwcEAgIEDByI5ORnLli2Dra0tgoKCjPsKqF41ryWzExjVVHWHiaenJx48eIA1a9bgqaeewunTp+Ht7Y3AwEDExMQgJiYGM2fOrHGHSVpaGqKjo/Hee++Z+2UQWaxGC3jv3r1x6NChOn+3fv36WtMUCgUWLFjQ8mREZFRqtVp3duzxO0zeffddAI/uMHn33Xcxc+bMeu8w4Rk2IvPgSGxE1KI7TIjIPORzbxERGYWh7zAxxC2ihtb02/Dkd8nJUO+zVG9ZlGIuc2diASeyYA3dYaJWq5t1h4khbhE1NCneGmhohnp9Un2vpJjLVJnqu0WUp9CJLFRjd5gAqHWHSUJCAoQQuHHjRqu8w4Sj/pGc8AicyELxDhMieWMBb4V4BEH64B0mRPLGU+hEREQyxAJOREQkQyzgJBu8NEBE9D8s4ERERDLEAk5ERCRD7IVORBavNVyeqXoNz7/UzsxJyFR4BE5ERCRDLOBEREQy1Ogp9A8//BDJyclwdnZGWFgYAODQoUM4efKkbozkGTNmYNCgQQCAr776CnFxcVAqlZg7dy4GDBhgxPhERESWqdEC7u/vjwkTJiAyMrLG9EmTJuGFF16oMe3OnTtITEzE+++/j7y8PGzcuBE7duyAUskDfWPitS8iIsvTaGXt27cvVCqVXitLSkqCr68vbGxs0L59e3Ts2BHp6ektDklEREQ1NbsX+okTJ5CQkABPT0/Mnj0bKpUKWq0WPXv21M3j4uICrVZb5/KGeGawuZ/FWh/T53p0BN4aetI2xhTvqxS3KylmIiLzalYBf+aZZzB16lQAwMGDB7Fv3z4EBQVBCKH3OgzxzGApPh8WkG6u1sAU76sUPz99MtX3zGCyLLykZjmadXG6Xbt2UCqVUCqVGDduHDIyMgAArq6uyM3N1c2n1Wrh4uJimKRERESk06wCnpeXp/v/+fPn0aVLFwCAj48PEhMTUVZWhqysLGRmZsLLy8swSYmq+fpgvkVcMiDj4nZEctboKfTt27fj2rVrKCoqwuLFizF9+nSkpqbi9u3bUCgUcHNzw8KFCwEAXbp0wYgRIxAcHAylUon58+ezBzoREZERNFrAly9fXmva2LFj651/ypQpmDJlSstSEdWDR0tkCNyOqDXg4TEREZEMsYATERHJEJ9GJmM8DUhEZLl4BE6yxl7ERGSpWMCJiIhkiAWciIhIhljAiYhaIV5eav1YwImIWjEW8taLBZyIiEiGeBuZDPHbNBE1FZ9S1vrwCJxaBZ4mJCJLwyNwIrIY/JJHrQmPwImIiGSo0SPwDz/8EMnJyXB2dkZYWBgA4P79+wgPD0d2djbc3NywYsUKqFQqCCGwd+9eXLx4EW3atEFQUBA8PT2N/iJaM163IiKiujR6BO7v74+1a9fWmBYTEwNvb29ERETA29sbMTExAICLFy/i3r17iIiIwMKFCxEdHW2c1ERERBau0QLet29fqFSqGtOSkpLg5+cHAPDz80NSUhIA4MKFCxg9ejQUCgV69eqF4uJi5OXlGSE2ERGRZWvWNfCCggKo1WoAgFqtRmFhIQBAq9VCo9Ho5nN1dYVWqzVATGIvayIiqs6gvdCFELWmKRSKOueNjY1FbGwsACA0NLRG4deXtbV1s5YzNsPmqlm0H62Xhbw+VV9y5i7xavY6pLhdGSMT+7cQyVuzCrizszPy8vKgVquRl5cHJycnAI+OuHNycnTz5ebm6o7UHxcQEICAgADdz9WX05dGo2nWcsZmzFxSfL1S1JL3SYrblT6Z3N3dm7ROf39/TJgwAZGRkbppVf1bAgMDERMTg5iYGMycObNG/5a0tDRER0fjvffea9ZrMYe9kenmjiAZ7BjbejTrFLqPjw/i4+MBAPHx8RgyZIhuekJCAoQQuHHjBuzt7est4ERkXuzfYtl4WU7+Gj0C3759O65du4aioiIsXrwY06dPR2BgIMLDwxEXFweNRoPg4GAAwMCBA5GcnIxly5bB1tYWQUFBRn8BRGQ4Te3fUtcXdENcHjM8Fqr6PP75SPESEiDNXObO1GgBX758eZ3T169fX2uaQqHAggULWp6KiCSlKf1bDHF5jEzn8c9HipeQAGnmMlWm+i6PcSQ2apV4erB5qvq3AGh2/xYiMg0WcCLSYf8WIvngw0yILJQl9G/hWRhqzVjAiSwU+7cQwNvK5IwFnIhaHR55kyXgNXCZ4Y6pedipjYhaGxZwsigs5ETUWrCASwwLDBER6YMFnIiISIbYiY2IWg2evWq+/z3JT1rDlVL9WMAlijsiIiJqCE+hExGRDh+9Kh8s4ERERDLEAk5ERCRDLboGvmTJErRt2xZKpRJWVlYIDQ3F/fv3ER4ejuzsbLi5uWHFihVQqVSGyttq8Zq3cfB9JWq6x9sNh1mVphZ3YtuwYYPukYMAEBMTA29vbwQGBiImJgYxMTGYOXNmS/8MkVF8fTCfOycikiWDn0JPSkqCn58fAMDPzw9JSUmG/hNELVZ9wBwOnkNEctTiI/CQkBAAwPjx4xEQEICCggLdc4LVajUKCwvrXC42NhaxsbEAgNDQUGg0Tb/30NraulnLGVtTcrHHp3RUfWZS3K6kmImIzKtFBXzjxo1wcXFBQUEBNm3aBHd3d72XDQgIQEBAgO7nnJycJv99jUbTrOWMTaq5qGFVn5kUPz99MjWl/RGR/LXoFLqLiwsAwNnZGUOGDEF6ejqcnZ2Rl5cHAMjLy6txfZyIiIgMo9kF/OHDh3jw4IHu/1euXEHXrl3h4+OD+Ph4AEB8fDyGDBlimKRERESk0+xT6AUFBdi2bRsAoKKiAk8//TQGDBiAHj16IDw8HHFxcdBoNAgODjZYWDmr6iRV1eOZnaaIiKglml3AO3TogK1bt9aa7ujoiPXr17coVGvGwk1kHGxbZGn4MBOi/49PY5InFm7je/wMIkkDh1IlIiKSIR6BExGRXuobYpVH6ObBAm5E3KjljZ8fUcN4+cK8WMBNgBu5vHB0PCKSAxZwI2ABaJ14RE5EUsICTkSyxDNbZOnYC52IiAyCT/YzLRZwIpIVFgn5aOhz4ufYcjyF3gK8JmqZ+LmbFnfy8sb2Yjws4AbADbR1YwEhahq2GdPgKXSiZnr8FCB3WkRkSjwCb4b6dtTcgRMR1a2x/SbPYDYdC/hjvj6Yzw2JmqSuHRN3Ss3HL8KWrb62UzW+xvMvtTN4+5Lrft9oBfzSpUvYu3cvKisrMW7cOAQGBhrrTxkNd8LUVI8Xn/rGjpar1tCuSd74Be9/jFLAKysr8fHHH+Odd96Bq6sr3n77bfj4+KBz587G+HNN8nhRrqtIN7YTJmopOX45NFa7luN7QYYnlf2unLZHoxTw9PR0dOzYER06dAAA+Pr6IikpySQNvb4P/fFlWJTJHOrbSdW3TUvp1J6x2nUVtklqKlNsM/rWlOass6VtWyGEEC1aQx3OnTuHS5cuYfHixQCAhIQEpKWlYf78+bp5YmNjERsbCwAIDQ01dAQiMjB92jXAtk1kKka5jayu7wQKhaLGzwEBAQgNDW1RA1+zZk2zlzUmKeZiJv1JMZcUMunTrgHDtG1Dk8L7Vxcp5pJiJkCaucydySgF3NXVFbm5ubqfc3NzoVarjfGniMhE2K6JpMUoBbxHjx7IzMxEVlYWysvLkZiYCB8fH2P8KSIyEbZrImmxevfdd9819EqVSiU6duyInTt34vjx4xg1ahSGDx9u6D8DAPD09DTKeltKirmYSX9SzGXuTKZs18Zg7vevPlLMJcVMgDRzmTOTUTqxERERkXFxLHQiIiIZYgEnIiKSIcmMhZ6Tk4PIyEjk5+dDoVAgICAAEydO1P3+6NGj2L9/P6Kjo+Hk5FRr+ZCQEKSlpaF37941uvZHREQgIyMD1tbW6NGjBxYuXAhra/1etrEyVfnkk09w6tQpfP7553rlMWYmIQQOHDiAc+fOQalUYvz48TXWa65cKSkp2L9/PyorK9G2bVssWbIEHTt2NHqm27dvY8+ePXjw4AGUSiWmTJkCX19fAEBWVha2b9+O+/fvo3v37li6dKlJtqmGMrVkO29NGhvq9fTp0/j888/h4uICAJgwYQLGjRtn1kwAkJiYiMOHD0OhUKBbt2548803jZpJn1yffvopUlNTAQClpaUoKCjAp59+atZMVe2nuLgYlZWVeOWVVzBo0CCjZtInV3Z2Nj766CMUFhZCpVJh6dKlcHV1NXouCInQarUiIyNDCCFESUmJWLZsmfjtt9+EEEJkZ2eLTZs2iTfeeEMUFBTUufyVK1dEUlKS2Lx5c43pP/30k6isrBSVlZUiPDxcnDhxwuyZhBAiPT1dREREiJkzZ+qdx5iZ4uLixM6dO0VFRYUQQoj8/HxJ5Kq+nuPHj4sPPvjAJJl+//13cffuXSGEELm5ueL1118X9+/fF0IIERYWJs6cOSOEEGLXrl0m26YaytSS7by1qKioEH//+9/FvXv3RFlZmVi5cqXuva1y6tQpESJex/kAAAabSURBVB0dLalMd+/eFatWrRJFRUVCiKa3PWPlqu6///2viIyMNHumqKgo3bb922+/iaCgIKNm0jdXWFiYOHXqlBBCiJSUFBEREWH0XEIIIZlT6Gq1Wtebz87ODh4eHtBqtQCAzz77DK+++mqdg0ZU8fb2hp2dXa3pgwYNgkKhgEKhgJeXV437WM2VqbKyEvv378fMmTP1zmLsTN9++y2mTp0KpfLRJuHs7CyJXADw4MEDAEBJSUmT7jtuSSZ3d3d06tQJAODi4gJnZ2cUFhZCCIHU1FRd72t/f38kJSWZNRPQsu28tag+1Ku1tbVuqFepZzp58iSeffZZqFQqAE1ve8bKVd3Zs2fx9NNPmz2TQqFASUkJgKbvD4yZ686dO/D29gYA9OvXDxcuXDB6LkBCp9Cry8rKwq1bt+Dl5YULFy7AxcUFTzzxRIvWWV5eju+//x5z5swxe6bjx49j8ODBLd74DJnpjz/+QGJiIs6fPw8nJyfMnTtXVyzMmWvx4sXYvHkzbG1tYWdnh5CQEJNnSk9PR3l5OTp06ICioiLY29vDysoKwKNCWlWAzZWpupZu53Km1WprnLZ0dXVFWlparfl+/PFH/Pzzz+jUqRNee+01aDQas2a6e/cuAGDdunWorKzEtGnTMGDAAKNl0jdXlezsbGRlZeHJJ580e6Zp06Zh06ZNOH78OP7880+sW7fOqJn0zdWtWzf8+OOPmDhxIs6fP48HDx6gqKgIjo6ORs0mmSPwKg8fPkRYWBjmzJkDKysrHDlyBC+99FKL1xsdHY0+ffqgT58+Zs2k1Wrxww8/4LnnnmvW8sbIBABlZWWwsbFBaGgoxo0bh48++kgSuY4dO4a3334bUVFRGDNmDPbt22fSTHl5edi5cyfeeOMN3dkJQzBWppZs53In9BjqdfDgwYiMjMS2bdvg7e2NyMhIs2eqrKxEZmYmNmzYgDfffBNRUVEoLi42e64qZ8+exfDhww26/Tc309mzZ+Hv74+oqCi8/fbb2LlzJyorK82ea9asWbh27RreeustXLt2DS4uLrov+cYkqQJeXl6OsLAwjBo1CsOGDcMff/yBrKwsrFq1CkuWLEFubi5Wr16N/PymPYHm8OHDKCwsxOzZs82e6fbt27h37x6WLVuGJUuWoLS0FEuXLjVrJuDRt8phw4YBAIYOHYpff/21SZmMkauwsBC//vorevbsCeDR06+uX79uskwlJSUIDQ3Fyy+/jF69egEAHB0dUVJSgoqKCgCPvpBVdYgyV6YqLdnOWwN9hnp1dHSEjY0NgEdjtt+8edPsmVxcXDBkyBBYW1ujffv2cHd3R2ZmptlzVUlMTMTIkSONmkffTHFxcRgxYgQAoFevXigrK0NRUZHZc7m4uGDlypXYsmULZsyYAQCwt7c3ai5AQqfQhRCIioqCh4cHJk+eDADo2rUroqOjdfMsWbIEmzdvrrMXc31OnjyJy5cvY/369U3+BmmMTIMGDcKePXt0P8+aNQs7d+40ayYAGDJkCK5evYqxY8fi2rVrcHd313tZY+VycHBASUkJ7t69C3d3d1y5cgUeHh4myVReXo5t27Zh9OjRuh0G8Oibd79+/XDu3DmMHDkSp0+fbtJwosbIBLRsO28tqg/16uLigsTERCxbtqzGPHl5ebqd74ULFwz2KNSWZBo6dCjOnDkDf39/FBYWIjMzs9alEXPkAh6d3i8uLq71ZdFcmTQaDa5evQp/f3/cuXMHZWVlTdrPGStXVe9zpVKJr776CmPGjDFqpiqSKeDXr19HQkICunbtilWrVgEAZsyYUe8tAhkZGfjuu+90jzZcv349fv/9dzx8+BCLFy/G4sWLMWDAAOzZswdubm74v//7PwDAsGHDMHXqVLNmagljZQoMDERERASOHTuGtm3bYtGiRZLItWjRIoSFhUGpVMLBwQFvvPGGSTIlJibi559/RlFREU6fPg3gUWF94okn8Oqrr2L79u04cOAAunfvjrFjx5o9U0u289bCysoK8+bNQ0hICCorKzFmzBh06dIFBw8eRI8ePeDj44NvvvkGFy5cgJWVFVQqFYKCgsyeqX///rh8+TJWrFgBpVKJmTNnGv3aqT65AODMmTPw9fVtsAOqKTPNnj0bu3btwrFjxwAAQUFBRs+mT65r167hyy+/hEKhQJ8+fWo9YtdYOJQqERGRDFnmuTYiIiKZYwEnIiKSIRZwIiIiGWIBJyIikiEWcCIiIhliASciIpIhFnAiIiIZ+n89I3ELRP/MAQAAAABJRU5ErkJggg==\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": [
    "- 局所解から抜け出せていない。\n",
    "  - 真値付近から開始(`initial_state = [0.6, 0.4]`)した場合には、その周辺の探索は可能。"
   ]
  },
  {
   "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 : 60.74 %\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+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dfVRUdf4H8PcAgsAADYypYJqKbkbkAz4lKah4Mq2N3c0etXygNFxMXU1zUztnM1mTg0EUKmUP1rE867K2blSkSEYmhBpiqwNWpwIDGRAcLB7m+/uDH7M8DPMAM3Pnzrxf53iOc+fOnfcM9zuf+/3e79xRCCEEiIiISFY8pA5ARERE1mMBJyIikiEWcCIiIhliASciIpIhFnAiIiIZYgEnIiKSIRZwIiIiGWIBJyJyQ4sXL0ZcXJykGTZt2oSBAwdCoVDgzTfflDSLHLGAUydyatSpqakYM2aM44IRSSguLg6LFy+22fZefvllHDx40Gbbs9ZXX32F7du3Y8+ePaisrMSDDz4oWRa58pI6APVNXFwchgwZYrOj15dffhl6vd4m2+qN9kadnZ2NKVOmICgoqMd1i4uLMW7cOAemI5K/pqYmeHt7m2xb1m6rNzQaDTw8PHDffff1OYe7Yg+cALQ1RAAICgqCSqWyybZ6o2OjHjRoEHx9fXtct7i4GOPHj+/1cxHZSmxsLJYuXYqNGzdCrVYjMDAQCQkJuH79umGd5uZmbNy4EWFhYfD29satt96K9957z3D/iRMnEB0djYCAAAQEBGDs2LH4+OOPAbSNjH322Wd46623oFAooFAokJeXZ3hseno6brnlFvTv3x+jRo3Ctm3b0NLSYsi2bNkybN68GYMHD0ZYWJhhmx1H28zlM7UtY0xtb/HixVi0aBH0er3h9ZhSUlKC3//+9xg8eLBh/fZ/r7/+usnHujRBdhMTEyOWLFkiNmzYIEJCQkRAQIBYtmyZaGxsNKzT1NQkNmzYIEJDQ0W/fv3EmDFjxLvvvmu4//PPPxfTpk0TSqVSKJVKcfvtt4ucnBwhhBCPP/64ANDp37FjxwyPTUtLE7/73e+Ej4+PCA8PFy+88IJobm42ZFu6dKl47rnnxKBBg4RarTZsc/bs2RbnM7UtY0xtz9jr6YlOpxOenp4iMzNT3HfffcLf318MHz68WzYiR4iJiREBAQEiISFBnD9/Xhw+fFgMGDBAJCUlGdZZt26dCA4OFh988IG4cOGC2LZtm1AoFCI3N1e0tLQIlUol1qxZIy5evCguXrwoDh06JPLz84UQQtTV1Ynp06eLBx54QFRWVorKykrx22+/CSGE2Lp1qxg6dKg4dOiQuHTpkjhy5Ii46aabxHPPPWfIplQqxfLly0Vpaan45ptvhBDd27qpfB1fp7FtGWNqe3V1dWLXrl3C09PT8Hp6otFoRGBgoHjkkUfEmTNnRFFRkRg/frwYMGCAeOedd0w+1tWxgNsRG3V3tmrUBQUFAoAYNWqU+Pe//y3KyspEUlKS8PLyEhcvXuzFX4uo92JiYsSwYcNES0uLYdnu3buFt7e3uHbtmtDpdMLb21tkZGR0elx8fLyYOXOm0Gq13Q7Au5o9e7Z4/PHHOy3T6XTC19dXfPTRR52Wv/XWWyIoKMiQbdSoUaK1tbXTOh3burl8HV+nsW11Zcn29u3bJzw9PU1uRwgh5s+fL2JiYoRerzcsO3DggAAgampqzD7elbGA2xEbdWe2bNQZGRlCoVCI4uJiw7Kmpibh6+srsrKyzD6eyJZiYmLEggULOi07d+6cACDOnj0rzp49KwCIkpKSTuukpqaKG2+8UQghREJCgvD29hZz584V27dvF//97387rWusrZ86dUoAEH5+fsLf39/wr3///gKAqKqqEjExMeKhhx7qlrljW7ckX/vrNLatrizZniVtXavVin79+ol//OMfnZZ/+OGHAoDQarVms7gyngO3s8mTJ8PT09NwOzo6Gk1NTSgvL0dZWRmampowY8aMTo+JiYlBaWkpVCoVEhIScNddd+Huu+9GcnIyLly4YPY5S0tLcf36dfzpT3+CUqk0/Fu+fDmuXr2K6upqAEBUVBQ8PHreBczl68jctqzdnjnFxcWYM2dOp3Pg/fr1g4+Pj6ST8IjaCSO/1Nz1XK8QwrBs7969+PrrrzFnzhwcP34ct912G3bv3m3yOdr39YMHD+LMmTOGfyUlJdBoNAgODgYA+Pv7W5TZVL52lm7L0u2Zcvr0aTQ3N2PChAmdlp86dQrh4eFQqVQoKSnBlClTDPedO3cO06ZNM/r+uxoWcAdjo+57owbaCnhUVFSnZRqNBnV1dd2WEzlCYWEhWltbDbe//PJLeHt7Y+TIkQgPD4ePjw+OHz/e6TH5+fmIiIgw3L7tttuwdu1afPTRR1i2bBn27NljuM/b27vT9gEgIiIC/fv3x6VLlxAeHt7tX8fOgymW5rOUrbbX/np1Op1h2dWrV/HGG28YvlIXERGB7777Ds3NzQCAtWvX4qWXXrL6M0WO+DUyO2tv1O0NqWOjVigUhp28405trFG3N+wVK1Zgz549WL58OQDzjXrevHm9zt6xEZrK5+jtNTU1obS0FHPmzOm0PCUlBVFRUd2O1okcoaamBitXrsTTTz+NS5cuYfPmzXjiiScMB7erVq3C5s2bMWDAAIwbNw4HDx7Ev/71L3z66acoKyvD3r17ce+99+Kmm25CRUUFPv/880778vDhw3Hs2DGUl5cjKCgIQUFBUCqV2LRpEzZt2gQAmDNnDlpaWlBSUoLTp0/j73//u0XZ/fz8TOazlq22N3nyZAQGBmLDhg1ITk5GXV0dnnnmGQwZMgTr168HAHh4eCAyMhIlJSX46aefoFKpEB0dbXVmOWIBtzM2attv79y5c2hqasKhQ4cwd+5cjBgxApmZmXj77bdx4sQJq3MR2cL999+PgIAA3HnnnWhqasKCBQuwY8cOw/3btm2Dh4cHVq9ejerqaoSHh2P//v2YPXs2KisrodFo8NBDD6G6uhohISGYP38+du7caXj8X/7yF5SUlGDs2LHQ6XQ4duwYYmNjsXnzZoSGhiI9PR3r1q2Dr68vRo8ebfVFX0zl6w1bbC8oKAiHDh3CmjVrMGHCBKjVaixYsAB/+9vfOn3/fOrUqSgoKMCePXtw+PDhXuWVJcnOvruB9q+Rtc+8ViqVYsmSJUKn0xnWMfW1qoqKCvGHP/xBhIWFCW9vbzF48GCRkJAg6urqDI8vLy8X06dPF/7+/t0mvGVlZYmxY8cKHx8fccMNN4jJkyeLV1991ZBt2bJl3TL39mtkxrZljLntWTKxJSsrS4SGhopPPvlEjBw5Uvj4+IiYmJhuE2aIHMWaNkC2d/jwYREcHCw2bNggdRSHUgjhBmf6JRIbG4vw8HBkZWVJHYWI7IhtXVplZWWIjo6GRqNBYGCg1HEchpPYiIhI1tLS0vDiiy+6VfEGeA6ciKjPOl7WlBynvLwc8+fPx7Rp07B06VKp4zgch9CJiIhkiEPoREREMsQCTkREJEMs4ERERDLkNJPYKioqLFpPrVbjypUrdk5jOeYxjXlMsyZPaGiondNIw9K2bw/Otj+0Yy7ruHIuU+2ePXAiIiIZYgEnIiKSIRZwIiIiGWIBJyIikiEWcCIiIhliASejPny/Dh++Xyd1DCKXwfZEtsYCTkRkBzwIJntjASciIpIhFnAiIiIZcporsZF8tA8L3vvgDRInIZK3rkPsbFNkDfbAiYiIZIgFnIjIyezLKJM6AskACzgREZEMsYBTr/FrMkRE0uEkNiIiO+JBLtkLe+BEREQyxAJOREQkQyzgREREMsQCTkREJEOcxEZERjU1NWHr1q1oaWlBa2srpk6digceeABVVVXYtWsXrl27huHDhyMpKQleXl5obm7GK6+8gkuXLiEgIACrV6/GjTfeKPXLIHJZ7IETkVH9+vXD1q1b8dJLL2HHjh04c+YMLl68iP3792P+/PlIS0uDv78/jh49CgA4evQo/P39kZ6ejvnz5+Pdd9+V+BUQuTYWcCIySqFQoH///gCA1tZWtLa2QqFQoLS0FFOnTgUAxMbGorCwEABQVFSE2NhYAMDUqVNx7tw5CCEkyU7kDjiETkQ90uv12LBhAy5fvoy77roLAwcOhJ+fHzw9PQEAwcHB0Gq1AACtVouQkBAAgKenJ/z8/NDQ0IDAwEDJ8ssNvzNO1mABJ6IeeXh44KWXXoJOp8POnTvx888/97iusd62QqHotiw3Nxe5ubkAgOTkZKjVatsFtpKXl5cdn797Me7+XD0XbCnfl57Y9/3qPXfNxQJORGb5+/vj1ltvhUajQWNjI1pbW+Hp6QmtVovg4GAAQEhICGpqahASEoLW1lY0NjZCqVR221ZcXBzi4uIMt69cueKw19GVWq222fNb8jO77T9SYsnPhkr5vvTElu+XLblyrtDQ0B7v4zlwIjKqvr4eOp0OQNuM9JKSEoSFhSEiIgInT54EAOTl5WHixIkAgKioKOTl5QEATp48iYiICKM9cOLvCJBtsAdOREbV1tYiIyMDer0eQgjccccdiIqKwpAhQ7Br1y4cOHAAw4cPx6xZswAAs2bNwiuvvIKkpCQolUqsXr1a4ldA5NpYwInIqGHDhmHHjh3dlg8cOBDbt2/vttzb2xtr1651RDSnwp40SYVD6ATANkN6/CAjInIcsz1wXo2JiEh6lkySI/ditgfOqzG5l5564uxdExE5F7MFnFdjIiIicj4WTWKzx9WYensxB2f7wr6c8+zLKMOSleH/f8t4D/t/26rr9H/z61ifxxGYh4hchUUF3B5XY+rtxRyc7Qv7cs9jbt2O9/e0rql15P7+2Js1eUxd0IGI3I9VXyOz5dWYyHXw/DiR/bB9UU/MngPn1ZiIiByPV2sjc8z2wHk1JiIiIudjtoDzakzyx++PEhG5Hl6JjYiISIZ4LXQ3x3NsRL3H9kNSYgEni/HDikh6PCVG7VjAiYhsiAe65Cg8B05ERCRD7IGTSexNEHXHdkHOgD1wIiIiGWIBJyIikiEWcCIiIhliASciIpIhFnAiIiIZYgEnu+AvKRER2RcLOBERkQyxgJNDsWdOZBtsS8QCTkREJEMs4ERERDLES6m6IQ67ERHJHwu4C2OhJrINufyEp1xykm1wCJ2IiEiGWMDJrvZllHEkgMjBOEPdPXAInWyKHxpERI7BAu5GWFyJiFwHCzgRkczx4Nw98Rw4ERGRDLGAExERyRALOBGRjHH43H3xHDgRGXXlyhVkZGSgrq4OCoUCcXFxmDdvHq5du4bU1FRUV1djwIABWLNmDZRKJYQQ2LdvH06fPg0fHx8kJiZixIgRUr8MIpfFAk5ERnl6emLRokUYMWIErl+/jo0bN+L2229HXl4eIiMjER8fj+zsbGRnZ2PhwoU4ffo0Ll++jLS0NGg0GmRlZeHFF1+U+mXYFHu75Ew4hE6S4Aeh81OpVIYetK+vL8LCwqDValFYWIiYmBgAQExMDAoLCwEARUVFmDFjBhQKBUaPHg2dTofa2lrJ8hO5OvbAicisqqoqfPfddwgPD8fVq1ehUqkAtBX5+vp6AIBWq4VarTY8JiQkBFqt1rBuu9zcXOTm5gIAkpOTOz3G0by8vCx8fnkecNr6vbX8/XIsd83FAk4OwR63fP36669ISUnB4sWL4efn1+N6QohuyxQKRbdlcXFxiIuLM9y+cuWKbYL2glqtlvT57c3Wr81Z3y9XzhUaGtrjfRxCJ6IetbS0ICUlBdOnT8eUKVMAAEFBQYah8draWgQGBgJo63F3/LCqqanp1vsmItsx2wPnTFQi9ySEQGZmJsLCwnDPPfcYlk+cOBHHjx9HfHw8jh8/jkmTJhmW5+TkIDo6GhqNBn5+fizgRHZktoBzJiqRe7pw4QLy8/MxdOhQrF+/HgDw8MMPIz4+HqmpqTh69CjUajXWrl0LABg/fjyKi4uxatUqeHt7IzExUcr4RC7PbAFXqVSGo+iuM1Gff/55AG0zUZ9//nksXLiwx5moPBInkpdbbrkFH3zwgdH7tmzZ0m2ZQqFAQkKCvWMR0f+zahKbM8xEdbbZhs6dx7knjknxvjn334uIyHIWF3BnmYnqbLMNpc7TPrv73gdvcIo81pAip7O9P9bkMTUblYjcj0Wz0DkTlezpw/fr+DUzIiIrmS3g5maiAug2EzU/Px9CCFy8eJEzUYlItnhwSc7M7BA6Z6LKy76MMsNwOhERuS6zBZwzUYmIiJwPr8RGREQkQyzgRERugufzXQt/zMRFdGyYcm+kXb8aR0RE3bEHTkTkhjjDXv5YwImIiGSIBZyIyEWxl+3aWMCJiIhkiAWciIhIhljAiYiIZIgFnIiISIZYwMlpcQIOEVHPWMCJiIhkiAWcZIG9cSKizngpVZJM14LMAk1EZDkWcCIiF8eDY9fEAk5E1AULHskBz4ETERHJEAu4THASFxERdcQCTkREJEM8B05E5EY4kuc62AMnIiKSIRZwcnrsMRARdccCTkREJEMs4ERERDLEAk5ERCRDnIUuUzwvTNQ3bW3of+3o3gdvkC4MUS+wB05ERCRDLOBEREQyxCF0meHQORERASzgRNSDV199FcXFxQgKCkJKSgoA4Nq1a0hNTUV1dTUGDBiANWvWQKlUQgiBffv24fTp0/Dx8UFiYiJGjBgh8SuwTvvBMc+Fk1xwCJ2IjIqNjcWmTZs6LcvOzkZkZCTS0tIQGRmJ7OxsAMDp06dx+fJlpKWl4cknn0RWVpYUkYncCgs4ERl16623QqlUdlpWWFiImJgYAEBMTAwKCwsBAEVFRZgxYwYUCgVGjx4NnU6H2tpah2cmcidmh9DdbRjNGfG8NzmLq1evQqVSAQBUKhXq6+sBAFqtFmq12rBeSEgItFqtYd2OcnNzkZubCwBITk7u9DjHMt6u2vK4T5uz5v338vKS8O/VM3fNZbaAx8bGYu7cucjIyDAsax9Gi4+PR3Z2NrKzs7Fw4cJOw2gajQZZWVl48cUX7RaeiJyDEKLbMoVCYXTduLg4xMXFGW5fuXLFbrl6Y19GmdQRHMqa91+tVjvd3wtw7VyhoaE93md2CJ3DaNL58P069r57wPdGGkFBQYY2XVtbi8DAQABtPe6OH1Q1NTVGe9/k/Ni25KNXs9ClHEZztqES++ZhI+rqf+91XZfblnGv/cf2Jk6ciOPHjyM+Ph7Hjx/HpEmTDMtzcnIQHR0NjUYDPz8/FnCZ+/D9Os7Id3I2/RqZI4bRnG2oxNnyuLqu77W1772z/b2syWNqKM0edu3ahfPnz6OhoQErVqzAAw88gPj4eKSmpuLo0aNQq9VYu3YtAGD8+PEoLi7GqlWr4O3tjcTERIdmpb5jr1t+elXA24fRVCoVh9GIXNTq1auNLt+yZUu3ZQqFAgkJCfaORHbAwi1fvfoaWfswGoBuw2j5+fkQQuDixYscRiMip8NzvOQqzPbAOYxGRHLFq6uRKzNbwDmMRs6MH9BkKfa6+679K3Zsb86BV2IjIiKSIf6YCckKe1FERG3YAyciIqM44c+5sQdORC6PRahv+P45JxZwJ8TGQkRE5rCAOxEWbiIishTPgRMREckQCzi5BE62IXI8tjtpsYATERHJEAs4ERH1CXvh0mABJ5fCDxIichechS4hXsebiOSIB8rOgQWcXA4PjIgFhtwBh9CJiIhkiAWciIhIhjiE7gQ43GcfHEonIlfGHrhEWLSJyJXwoi6OxwLuINy5iYjIljiE7mAs4kTkynjqynHYAyciIpvjqKP9sYATERHJEAs4uTz2BIjIFbGA2wmLhvPj34eI5IyT2GyIBYHIcT58v44TpcitsYCT22g7wOKHPhG5Bg6hExGR3fG0ou2xB25n3GGdD/8mRI5jrL11XcZRsd5hAScil8GDM3InLOB90H5OlYiILGOqR961J86JiqbxHLgVeA6HiIicBXvgFmDRdl2mjvyNLSci+zHW7tgWe8YC3kXXnYXF2z309Hfmh4dz49/HNfFz1zIcQieyAD9QiKTFU5jd2aUHfubMGezbtw96vR6zZ89GfHy8PZ7GrrijUFf86ot5UrV9tlf3YXyU1D0nu9m8gOv1erz++ut47rnnEBISgmeffRYTJ07EkCFDbP1UNmXNBwA/LIi6s1fb5zA5GdP1c9jUjHVTB99y3r9sXsDLysowaNAgDBw4EAAwbdo0FBYW2qwRL1mp7nS7nbE331Sh5TlushVj+1BP+1f35a7Tc7BX2yeylKWf55ZcXMYUa9qssc8AWx00KIQQok9b6OLkyZM4c+YMVqxYAQDIz8+HRqPBsmXLOq2Xm5uL3NxcAEBycrItIxCRBNj2iRzL5pPYjB0PKBSKbsvi4uKQnJxsdQPeuHFjr7PZA/OYxjymOVuevrB327cHZ33/mcs67prL5gU8JCQENTU1hts1NTVQqVS2fhoicjJs+0SOZfMCPnLkSFRWVqKqqgotLS0oKCjAxIkTbf00RORk2PaJHMvz+eeff96WG/Tw8MCgQYOQnp6OnJwcTJ8+HVOnTrXlU2DEiBE23V5fMY9pzGOas+XpLUe0fXtw1vefuazjjrlsPomNiIiI7I9XYiMiIpIhFnAiIiIZkuzHTK5cuYKMjAzU1dVBoVAgLi4O8+bNM9x/+PBh7N+/H1lZWQgMDOz2+G3btkGj0eCWW27pNFU/LS0N5eXl8PLywsiRI/Hkk0/Cy8v8y7RXnnZvvPEGjh07hnfeecdsFnvmEULgwIEDOHnyJDw8PDBnzpxO25UiU0lJCfbv3w+9Xo/+/ftj5cqVGDRokF3zfP/999i7dy+uX78ODw8P/PGPf8S0adMAAFVVVdi1axeuXbuG4cOHIykpye77kKk8vd2n6X/MXeI1Ly8P77zzDoKDgwEAc+fOxezZsyXPBQAFBQU4ePAgFAoFhg0bhqefflryXG+++SZKS0sBAE1NTbh69SrefPNNyXO1t0GdTge9Xo9HHnkEEyZMkDxXdXU1XnvtNdTX10OpVCIpKQkhISF9f2IhEa1WK8rLy4UQQjQ2NopVq1aJH3/8UQghRHV1tXjhhRfEU089Ja5evWr08d98840oLCwU27dv77T866+/Fnq9Xuj1epGamio+/vhjSfMIIURZWZlIS0sTCxcutCiLPfMcPXpUpKeni9bWViGEEHV1dZJn6ridnJwc8corr9g9z88//ywqKiqEEELU1NSIJ554Qly7dk0IIURKSoo4ceKEEEKI3bt3O2QfMpWnt/s0tWltbRV//vOfxeXLl0Vzc7NYt26d4e/S7tixYyIrK8vpclVUVIj169eLhoYGIYR17dWeuTr6z3/+IzIyMpwiV2ZmpqF9/PjjjyIxMdEpcqWkpIhjx44JIYQoKSkRaWlpNnluyYbQVSqVYXaer68vwsLCoNVqAQBvvfUWHn30UaMXgWgXGRkJX1/fbssnTJgAhUIBhUKB8PDwTt9LlSKPXq/H/v37sXDhQoty2DvPJ598gvvvvx8eHm1/+qCgIMkzAcD169cBAI2NjRZ/d7gveUJDQzF48GAAQHBwMIKCglBfXw8hBEpLSw2zp2NjY1FYWChZHqD3+zS16XiJVy8vL8MlXqVmSa7PPvsMd911F5RKJQDr2qs9c3X0xRdf4M4773SKXAqFAo2NjQCs+yyxd66ffvoJkZGRAICIiAgUFRXZ5LmdYhyuqqoK3333HcLDw1FUVITg4GDcfPPNfdpmS0sLPv/8cyxevFjSPDk5OYiKiurTjmTLPL/88gsKCgpw6tQpBAYGYsmSJYbCIVWmFStWYPv27fD29oavry+2bdvm0DxlZWVoaWnBwIED0dDQAD8/P3h6egJoK6btRViKPB31ZZ92Z1qtttNwZUhICDQaTbf1vvrqK3z77bcYPHgwHn/8cajVaslzVVRUAAA2b94MvV6PBQsWYNy4cZLnalddXY2qqircdtttds1kaa4FCxbghRdeQE5ODn777Tds3rzZKXINGzYMX331FebNm4dTp07h+vXraGhoQEBAQJ+eW/JJbL/++itSUlKwePFieHp64tChQ3jwwQf7vN2srCyMGTMGY8aMkSyPVqvFl19+ibvvvrtXj7d1HgBobm5Gv379kJycjNmzZ+O1116TPNORI0fw7LPPIjMzEzNnzsTbb7/tsDy1tbVIT0/HU089ZRiV6Ct75entPu3uhAWXeI2KikJGRgZ27tyJyMhIZGRkOEUuvV6PyspKbN26FU8//TQyMzOh0+kkz9Xuiy++wNSpU23WdkyxJNcXX3yB2NhYZGZm4tlnn0V6ejr0er3kuRYtWoTz58/jmWeewfnz5xEcHGzoJPSFpAW8paUFKSkpmD59OqZMmYJffvkFVVVVWL9+PVauXImamhps2LABdXXW/WLYwYMHUV9fj8cee0zSPN9//z0uX76MVatWYeXKlWhqakJSUpJkeYC2o8MpU6YAACZPnowffvjB4sfaI1N9fT1++OEHjBo1CkDbL1hduHDBIXkaGxuRnJyMhx56CKNHjwYABAQEoLGxEa2trQDaDsLaJzZJkaddb/dpsuwSrwEBAejXrx+Atmu1X7p0ySlyBQcHY9KkSfDy8sKNN96I0NBQVFZWSp6rXUFBAaKjo+2ax5pcR48exR133AEAGD16NJqbm9HQ0CB5ruDgYKxbtw47duzAww8/DADw8/Pr83NLNoQuhEBmZibCwsJwzz33AACGDh2KrKwswzorV67E9u3bjc5o7slnn32Gs2fPYsuWLVYdFdojz4QJE7B3717D7UWLFiE9PV2yPAAwadIknDt3DrNmzcL58+cRGhpq8WPtkcnf3x+NjY2oqKhAaGgovvnmG4SFhdk9T0tLC3bu3IkZM2YYGjzQduQcERGBkydPIjo6Gnl5eRZfDtQeeYDe79PUpuMlXoODg1FQUIBVq1Z1Wqe2ttbwoVtUVOSQn0C1JNfkyZNx4sQJxMbGor6+HpWVld1OrUiRC2gb3tfpdN0ONqXMpVarce7cOcTGxuKnn35Cc3OzVZ+P9srVPvvcw8MD//znPzFz5kybPLdkBfzChQvIz8/H0KFDsX79egDAww8/3OOU//Lycg/t/OEAAAFQSURBVHz66aeGnyrcsmULfv75Z/z6669YsWIFVqxYgXHjxmHv3r0YMGAA/vrXvwIApkyZgvvvv1+yPL1lrzzx8fFIS0vDkSNH0L9/fyxfvlzyTMuXL0dKSgo8PDzg7++Pp556yu55CgoK8O2336KhoQF5eXkA2orrzTffjEcffRS7du3CgQMHMHz4cMyaNUvSPL3dp6mNp6cnli5dim3btkGv12PmzJm46aab8P7772PkyJGYOHEiPvroIxQVFcHT0xNKpRKJiYlOkWvs2LE4e/Ys1qxZAw8PDyxcuLDP501tkQsATpw4gWnTppmcuOroXI899hh2796NI0eOAAASExPtns+SXOfPn8d7770HhUKBMWPGdPuJ3d7ipVSJiIhkiONxREREMsQCTkREJEMs4ERERDLEAk5ERCRDLOBEREQyxAJOREQkQyzgREREMvR/u8DQCbjNVusAAAAASUVORK5CYII=\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": [
    "- ややacceptance rateが変動するものの、結果に大きな変動は見られなかった。"
   ]
  },
  {
   "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時点では`SimpleStepSizeAdoptation`と併用できないため、ステップサイズ調整は無し。"
   ]
  },
  {
   "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([24.160742, 24.183947, 24.15525 , ..., 24.159927, 24.178778,\n",
       "        24.16116 ], dtype=float32),\n",
       " array([0.59501934, 0.53401023, 0.5231044 , ..., 0.55630195, 0.5180631 ,\n",
       "        0.5327149 ], 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+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3daVRU5xkH8P9lEAQGCDBu4BIEbQwSF1yIRkHFExuz2DbG7DFKEoPFoHWLidqexkhUAkFJUUnVZmmMp9aYptmIIi4xgYAGsJVFe44JWJYBgQHLMm8/UKYsIzPAzNy58P+dwwfu3Ln3uctzn/u+dxlJCCFAREREiuIgdwBERETUfSzgRERECsQCTkREpEAs4ERERArEAk5ERKRALOBEREQKxAJORESkQCzgRET9wNKlSxERESFrDJs2bcKQIUMgSRIOHjwoayx9AQs4dUlJSR8fH49x48bZLjAiK4qIiMDSpUstNr233noLR44csdj0uuvbb7/F9u3bsW/fPpSUlGDJkiWyxdJXOModAFlWREQEhg8fbrGz27feegt6vd4i0+qJ1qQ/duwYpk+fDk9Pz1uOm5WVhYkTJ9owOiL719DQACcnpy5zp7vT6omCggI4ODjgoYce6nUc1IItcDKqoaEBAODp6QkvLy+LTKsn2ib90KFD4eLicstxs7KyMGnSpB7Pi8hc4eHhWLZsGTZu3AiNRgMPDw9ERkaivr7eME5jYyM2btwIPz8/ODk54c4778QHH3xg+PzMmTOYOXMm3N3d4e7ujgkTJuCLL74A0NLz9fXXX+PQoUOQJAmSJCEtLc3w3d27d+OOO+7AwIEDMWbMGGzbtg1NTU2G2JYvX47Nmzdj2LBh8PPzM0yzbW+aqfi6mpYxXU1v6dKleOqpp6DX6w3L05WcnBw8+OCDGDZsmGH81r933nmny+/2K4JsJiwsTDz77LNiw4YNwsfHR7i7u4vly5eLuro6wzgNDQ1iw4YNwtfXVwwYMECMGzdOvP/++4bPT58+LWbMmCHUarVQq9XirrvuEp9//rkQQohnnnlGAGj3d/LkScN3ExMTxc9+9jPh7OwsAgMDxWuvvSYaGxsNsS1btky8+uqrYujQoUKj0RimOW/ePLPj62paxnQ1PWPLcys6nU6oVCqRnJwsHnroIeHm5ib8/f07xUZkCWFhYcLd3V1ERkaKS5cuiePHj4tBgwaJ6Ohowzhr164V3t7e4qOPPhKXL18W27ZtE5IkidTUVNHU1CS8vLzE6tWrRX5+vsjPzxdHjx4V6enpQgghqqqqxKxZs8QjjzwiSkpKRElJifjPf/4jhBBi69atYuTIkeLo0aPiypUr4tNPPxUjRowQr776qiE2tVotXnjhBZGXlyd++OEHIUTnXO4qvrbLaWxaxnQ1vaqqKpGQkCBUKpVheW6loKBAeHh4iMcff1xcuHBBZGZmikmTJolBgwaJd999t8vv9jcs4DbEpO/MUkl/7tw5AUCMGTNG/O1vfxOFhYUiOjpaODo6ivz8/B5sLaJbCwsLE6NGjRJNTU2GYXv37hVOTk6itrZW6HQ64eTkJJKSktp9b9GiRWLOnDlCq9V2OsHuaN68eeKZZ55pN0yn0wkXFxfx2WeftRt+6NAh4enpaYhtzJgxorm5ud04bXPZVHxtl9PYtDoyZ3oHDhwQKpWqy+kIIcTChQtFWFiY0Ov1hmEffvihACAqKipMfr8/YQG3ISZ9e5ZM+qSkJCFJksjKyjIMa2hoEC4uLiIlJcXk94m6IywsTCxevLjdsNzcXAFAXLx4UVy8eFEAEDk5Oe3GiY+PF4MHDxZCCBEZGSmcnJzEggULxPbt28U///nPduMay+XvvvtOABCurq7Czc3N8Ddw4EABQJSWloqwsDDx6KOPdoq5bS6bE1/rchqbVkfmTM+cXNZqtWLAgAHiL3/5S7vhn3zyiQAgtFqtyVj6E14Dt7Fp06ZBpVIZ/p85cyYaGhpQVFSEwsJCNDQ0YPbs2e2+ExYWhry8PHh5eSEyMhL33nsvfv7znyM2NhaXL182Oc+8vDzU19fjV7/6FdRqteHvhRdewI0bN1BWVgYACAkJgYPDrXcJU/G1ZWpa3Z2eKVlZWZg/f367a+ADBgyAs7OzrDfhUf8hjPwyc8drvUIIw7D9+/fj+++/x/z583Hq1CmMHz8ee/fu7XIerfvykSNHcOHCBcNfTk4OCgoK4O3tDQBwc3MzK+au4mtl7rTMnV5XsrOz0djYiMmTJ7cb/t133yEwMBBeXl7IycnB9OnTDZ/l5uZixowZRtd/X8cCLjMmfe+THmgp4CEhIe2GFRQUoKqqqtNwIkvIyMhAc3Oz4f9vvvkGTk5OCAgIQGBgIJydnXHq1Kl230lPT0dQUJDh//Hjx2PNmjX47LPPsHz5cuzbt8/wmZOTU7vpA0BQUBAGDhyIK1euIDAwsNNf28ZBV8yNz1yWml7r8up0OsOwGzdu4I9//KPhkbqgoCBcvXoVjY2NAIA1a9Zg586d3T5m9AV8jMzGWpO+NdHaJr0kSYYkaLvTG0v61sRfsWIF9u3bhxdeeAGA6aS/7777ehx72yTtKj5bT6+hoQF5eXmYP39+u+FxcXEICQnpdDZPZAkVFRVYuXIlXnrpJVy5cgWbN2/Gc889Zzh5XbVqFTZv3oxBgwZh4sSJOHLkCD7++GN89dVXKCwsxP79+/HAAw9gxIgRKC4uxunTp9vtq/7+/jh58iSKiorg6ekJT09PqNVqbNq0CZs2bQIAzJ8/H01NTcjJyUF2djbeeOMNs2J3dXXtMr7ustT0pk2bBg8PD2zYsAGxsbGoqqrC+vXrMXz4cKxbtw4A4ODggODgYOTk5ODHH3+El5cXZs6c2e2Y+wIWcBtj0lt+erm5uWhoaMDRo0exYMECjB49GsnJyfjTn/6EM2fOdDsuInM8/PDDcHd3xz333IOGhgYsXrwYO3bsMHy+bds2ODg4ICYmBmVlZQgMDMR7772HefPmoaSkBAUFBXj00UdRVlYGHx8fLFy4ELt27TJ8/ze/+Q1ycnIwYcIE6HQ6nDx5EuHh4di8eTN8fX2xe/durF27Fi4uLhg7dmy3X/rSVXw9YYnpeXp64ujRo1i9ejUmT54MjUaDxYsX4/e//327589DQ0Nx7tw57Nu3D8ePH+9RvH2CbFff+6HWx8ha77xWq9Xi2WefFTqdzjBOV49VFRcXi1/84hfCz89PODk5iWHDhonIyEhRVVVl+H5RUZGYNWuWcHNz63TDW0pKipgwYYJwdnYWt912m5g2bZp4++23DbEtX768U8w9fYzM2LSMMTU9c258SUlJEb6+vuLLL78UAQEBwtnZWYSFhXW6oYbIUrqzj5PlHT9+XHh7e4sNGzbIHYqsJCH64ZV/mYSHhyMwMBApKSlyh0JEvcBclldhYSFmzpyJgoICeHh4yB2ObHgTGxERKUpiYiJef/31fl28AV4DJyLqtravNSXbKSoqwsKFCzFjxgwsW7ZM7nBkxy50IiIiBWIXOhERkQKxgBMRESkQCzgREZEC2c1NbMXFxTabl0ajQXl5uc3m11NKiFMJMQL9I05fX18LR2MZ5ua2vW4jxtU9jKt7TMXVVV6zBU5ERKRALOBEREQKxAJORESkQCzgRERECsQCTkREpEAs4GTUJ4er8MnhKrnDIOpXmHfUHSzgRER2joWdjGEBJyIiUiAWcCIiIgViASciIlIgFnAiIiIFYgEnIiJSILv5MROyD7zTlYhIGdgCpy51fHyFj7MQEdkHtsCJiOwUT5apK2yBExERKRBb4NQrB5IKAQAPLLlN5kiIlIstbeoJFnAyCw8wRLbDfCNzsAudiIhIgVjAiYiIFIgFnIhIIfgYJ7XFAk5ERKRALOBEREQKxLvQqUfYjUdEJC+2wImIFIYn0ASwgPdr1jgItN5kwwMMEZF1sYCT1bCQE1kf86z/MnkNvKGhAVu3bkVTUxOam5sRGhqKRx55BKWlpUhISEBtbS38/f0RHR0NR0dHNDY2Ys+ePbhy5Qrc3d0RExODwYMH22JZSEY8gCgPc1vZWl9j3Ko1B/la4/7DZAt8wIAB2Lp1K3bu3IkdO3bgwoULyM/Px3vvvYeFCxciMTERbm5uOHHiBADgxIkTcHNzw+7du7Fw4UK8//77Vl8IIuo+5jaRspks4JIkYeDAgQCA5uZmNDc3Q5Ik5OXlITQ0FAAQHh6OjIwMAEBmZibCw8MBAKGhocjNzYUQwkrhE1FPMbf7Jnap9x9mPUam1+uxYcMGXL9+Hffeey+GDBkCV1dXqFQqAIC3tze0Wi0AQKvVwsfHBwCgUqng6uqKmpoaeHh4tJtmamoqUlNTAQCxsbHQaDQWWyhTHB0dbTq/nrJ+nLZJdHtY19zmxtlTbtvrNrJNXJbPQ7nWZf/ejt3Xm7jMKuAODg7YuXMndDoddu3ahZ9++umW4xo7I5ckqdOwiIgIREREGP4vLy83JxSL0Gg0Np1fTyklTlPsYRmUsi57E6evr2+3v2NPuW2v28gacdnierVc67I/bUdLMBVXV3ndrbvQ3dzccOedd6KgoAB1dXVobm4G0HJm7u3tDQDw8fFBRUUFgJZuubq6OqjV6u7MhvoYdunZP+Y2kfKYLODV1dXQ6XQAWu5azcnJgZ+fH4KCgnD+/HkAQFpaGqZMmQIACAkJQVpaGgDg/PnzCAoKMnqWTkTyYm7bHk9myZJMdqFXVlYiKSkJer0eQgjcfffdCAkJwfDhw5GQkIAPP/wQ/v7+mDt3LgBg7ty52LNnD6Kjo6FWqxETE2P1hSCi7mNuEymbJOzkNtLi4mKbzcter4V0ZO04bd0SkPP51P6wzXtyDdwWzM1te91GloyrY8615oQ1clGufOsP29GSbHYNnIiIiOwDCzgREZECsYATEREpEAs4ERGRArGAk83wERoiIsthASciIlIgFnAiIiIFMutd6NS3sBubyD4wF6k32AInIiJSILbA+xGe7RMR9R0s4P0ACzeR7bXNOzlfI0x9F7vQiYj6ID622fexgBMRESkQCzjJhi0EIqKe4zVwIiIr44kqWQNb4ERERArEAk5ERKRA7EInIrIgdpeTrbAFTkREpEAs4ERERArEAk5E1A/wsc2+hwWciIhIgXgTG9kcWwFEtsN867vYAiciIlIgFnAiIiIFYgEnIiJSIBZwIiIiBWIBJyIiUiDehd6H8e5TIqK+iy1wIiIiBWIBJ9nxDVFERN3HAk5ERKRALOBkN9gSJ7Id5pvymbyJrby8HElJSaiqqoIkSYiIiMB9992H2tpaxMfHo6ysDIMGDcLq1auhVqshhMCBAweQnZ0NZ2dnREVFYfTo0bZYFiLqBuZ2/9SxaH9yuAoPLLlNpmioN0wWcJVKhaeeegqjR49GfX09Nm7ciLvuugtpaWkIDg7GokWLcOzYMRw7dgxPPvkksrOzcf36dSQmJqKgoAApKSl4/fXXbbEsRNQNzG3LYmuWbM1kF7qXl5fhLNvFxQV+fn7QarXIyMhAWFgYACAsLAwZGRkAgMzMTMyePRuSJGHs2LHQ6XSorKy04iIQUU8wt4mUrVvPgZeWluLq1asIDAzEjRs34OXlBaDlQFBdXQ0A0Gq10Gg0hu/4+PhAq9Uaxm2VmpqK1NRUAEBsbGy771ibo6OjTefXU72PU5ktAmtsm/6zzXvGHnLbXreR+XEpM98Ay+ac8rejbfUmLrML+M2bNxEXF4elS5fC1dX1luMJIToNkySp07CIiAhEREQY/i8vLzc3lF7TaDQ2nV9PmRtna9ddX7mOZY1t09e2uTG+vr49+p695La9biN7jcuSLLl89rq+lBpXV3lt1l3oTU1NiIuLw6xZszB9+nQAgKenp6H7rLKyEh4eHgBazsrbBlNRUdHpDJ2I7ANzmwDeka5UJgu4EALJycnw8/PD/fffbxg+ZcoUnDp1CgBw6tQpTJ061TA8PT0dQgjk5+fD1dWVSW5jTEYyB3ObSNlMdqFfvnwZ6enpGDlyJNatWwcAeOyxx7Bo0SLEx8fjxIkT0Gg0WLNmDQBg0qRJyMrKwqpVq+Dk5ISoqCjrLgER9Qhzm0jZTBbwO+64Ax999JHRz7Zs2dJpmCRJiIyM7H1kRGRVzG0iZeOb2IiIiBSIBZyIiEiB+HvgRES9wBtGSS4s4H0IDyREtsN8I7mxC52IiADwEVSlYQEnIiJSIBZwIiIiBWIBJyIiUiAWcCIiIgViASciIlIgFnAiIiIFYgEnIiJSIL7IRcH66vOarcv1wJLbZI6EiMh+sQVORESkQCzgRERECsQCTkREpEAs4GS3+uo1fiIiS2ABJyIiUiDehU5E1A3sGSJ7wRY4ERGRArGAExFRO/xdcGVgASdF4AGFyPaYd/aNBZyIiEiBWMCJiIgUiAWciIhIgfgYGdk1Xn8jIjKOLXAiIiIFYgEnIurCgaRC9gT9D+9Kty8s4ERERArEAk5ERKRAvInNzrV2Vz2w5LZOw4iIqP9iAVcQFm4i+TD/yN6wC52IiEiBTLbA3377bWRlZcHT0xNxcXEAgNraWsTHx6OsrAyDBg3C6tWroVarIYTAgQMHkJ2dDWdnZ0RFRWH06NFWXwjqP4xdUqCeYW4TKZvJFnh4eDg2bdrUbtixY8cQHByMxMREBAcH49ixYwCA7OxsXL9+HYmJiXj++eeRkpJinaj7IXbfkaUxt8lcfHzMPpks4HfeeSfUanW7YRkZGQgLCwMAhIWFISMjAwCQmZmJ2bNnQ5IkjB07FjqdDpWVlVYIm4h6i7lNpGw9uontxo0b8PLyAgB4eXmhuroaAKDVaqHRaAzj+fj4QKvVGsZtKzU1FampqQCA2NjYdt+zNkdHR5vOr6ccHXmP4a10d/spaZvLGaecuS33st8aW54dtbbGn10Z2Okze92OfTEui1YIIUSnYZIkGR03IiICERERhv/Ly8stGUqXNBqNTefXU/a4s9mL7m4/JW3znsbp6+tr4Wj+zxa5rZRtRP9nbHvZ63ZUalxd5XWP7kL39PQ0dJ9VVlbCw8MDQMtZedtAKioqjJ6hE5F9Ym4TKUePCviUKVNw6tQpAMCpU6cwdepUw/D09HQIIZCfnw9XV1cmOZGCMLeJlMNkF3pCQgIuXbqEmpoarFixAo888ggWLVqE+Ph4nDhxAhqNBmvWrAEATJo0CVlZWVi1ahWcnJwQFRVl9QXoq1quMfHaG1kPc5sshY93ysNkAY+JiTE6fMuWLZ2GSZKEyMjI3kdFRFbH3CZSNt7mTIrU8ZlUnvkTUX/DAk5EZARfXEL2jgWciIh6hCc58uKPmRARESkQCzj1CXxXMxH1N+xCJyJqgyeCpBRsgRMRESkQW+B2hmf/RERkDhZwIiKyiLZvkOS7GayPXehEREQKxAJORERWwydErIcFnPokHjSIqK9jASci+h+e9JGSsIATEZHVsVfM8ngXusz4O7pE8mNhISViAbcTPIBYBtcjEfUXLOA2xgJDRP0Bj3XWx2vgRERECsQCTkREpEAs4NSn8c5XIvvE3Ow9XgOnfuWTw1W8458MWEBsr+M655M4PccCTkT9Dgs39QXsQiciIrvDLnbT2AK3Ee6I8jqQVCh3CETUBXaldx9b4ERERArEFriVseVNRGQ+HjPNxwJuJdwJ7Re76oiUibnbHrvQqd/iTTL9B7e1cnHb3Rpb4ETUb7AQ9A1sibdgAbeAjgeF/r5TKQ0PBkT2jydfnbGAWwF3tL6BhV25mIPUH7CAE/0PD/pEytLdk+y+dlLOAt4DPND3L7xEQqQMtzo2d8xZY4VcicXdKgX8woULOHDgAPR6PebNm4dFixZZYzZWc6sNycJNgDIT3VLsPbeZo9RbSspvixdwvV6Pd955B6+++ip8fHzw8ssvY8qUKRg+fLilZ9Ur3dlIPCgQ2V9u85flqJWpY3TL553H6WofUkIht3gBLywsxNChQzFkyBAAwIwZM5CRkWGRJG+7Qk2t3I6ft9/AxjekOcOIWrXdP4zvZ53dajxj+7G9HUCslduty/nsSk27/zsytu6Yo9Rb5u5D5ozXcR+19smBJIQQvZpCB+fPn8eFCxewYsUKAEB6ejoKCgqwfPnyduOlpqYiNTUVABAbG2vJEIjICpjbRPbF4m9iM3Y+IElSp2ERERGIjY2VJcE3btxo83n2hBLiVEKMAOO0BGvntr0uO+PqHsbVPb2Jy+IF3MfHBxUVFYb/Kyoq4OXlZenZEJGNMbeJ7IvFC3hAQABKSkpQWlqKpqYmnDt3DlOmTLH0bIjIxpjbRPZF9dvf/va3lpygg4MDhg4dit27d+Pzzz/HrFmzEBoaaslZWMTo0aPlDsEsSohTCTECjLO3bJHb9rrsjKt7GFf39DQui9/ERkRERNbHnxMlIiJSIBZwIiIiBVLcu9DLy8uRlJSEqqoqSJKEiIgI3HfffYbPjx8/jvfeew8pKSnw8PDo9P1t27ahoKAAd9xxR7vb9xMTE1FUVARHR0cEBATg+eefh6OjI/Ly8rBjxw4MHjwYADB9+nQ8/PDDssWZlJSES5cuwdXVFQCwcuVK3H777RBC4MCBA8jOzoazszOioqJMXlexVoxbtmxBfX09AKC6uhoBAQFYv369LOvyX//6F/bv34/6+no4ODjgl7/8JWbMmAEAKC0tRUJCAmpra+Hv74/o6Gg4OjqisbERe/bswZUrV+Du7o6YmBhDzHLEael9U27mvo71/PnzePPNN7F9+3YEBATIHldaWhreffddeHt7AwAWLFiAefPmyR4XAJw7dw5HjhyBJEkYNWoUXnrpJdnjOnjwIPLy8gAADQ0NuHHjBg4ePCh7XK15qtPpoNfr8fjjj2Py5Mmyx1VWVoY//OEPqK6uhlqtRnR0NHx8fLqeqFAYrVYrioqKhBBC1NXViVWrVolr164JIYQoKysTr732mnjxxRfFjRs3jH7/hx9+EBkZGWL79u3thn///fdCr9cLvV4v4uPjxRdffCGEECI3N7fTuHLGuWfPHvHNN990Gv/7778X27ZtE3q9Xly+fFm8/PLLssXY1s6dO0VaWpoQQp51+dNPP4ni4mIhhBAVFRXiueeeE7W1tUIIIeLi4sSZM2eEEELs3bvXsM0///xzsXfvXiGEEGfOnBFvvvmmrHFaet+UU3Nzs/j1r38trl+/LhobG8XatWsN66ituro6sWXLFrFp0yZRWFhoF3GdPHlSpKSkWD2W7sZVXFws1q1bJ2pqaoQQQlRVVdlFXG39/e9/F0lJSXYRV3JysiGHrl27JqKiouwirri4OHHy5EkhhBA5OTkiMTHR5HQV14Xu5eVlaFm6uLjAz88PWq0WAHDo0CE88cQTRl8u0So4OBguLi6dhk+ePBmSJEGSJAQGBrZ73tWe4ryVzMxMzJ49G5IkYezYsdDpdKisrJQ1xvr6euTl5WHq1KlmL4el4/T19cWwYcMAAN7e3vD09ER1dTWEEMjLyzPcRR0eHo6MjAwALesyPDwcABAaGorc3FyjLzGxRZyA5fdNObV9Haujo6PhdawdHT58GA8++CAGDBhgV3HZmjlxff3117j33nuhVqsBAJ6ennYRV1tnz57FPffcYxdxSZKEuro6AEBdXZ1N3mVgTlw//vgjgoODAQBBQUHIzMw0OV3FFfC2SktLcfXqVQQGBiIzMxPe3t64/fbbezXNpqYmnD59GhMnTjQMy8/Px7p16/D666/j2rVrssf55z//GWvXrsXBgwfR2NgIANBqtdBoNIZxfHx8DMVDjhgB4LvvvsP48eMN3f2AvOuysLAQTU1NGDJkCGpqauDq6gqVSgWgpWi2ri+tVmvoulKpVHB1dUVNTY0scbZljX3T1tquW8D4fnr16lWUl5cjJCTEruICgG+//RZr165FXFwcysvL7SKu4uJilJSUYPPmzXjllVdw4cIFu4irVVlZGUpLSzF+/Hi7iGvx4sU4ffo0VqxYge3bt2PZsmV2EdeoUaPw7bffAmg5dtbX15s87ii2gN+8eRNxcXFYunQpVCoVjh49iiVLlvR6uikpKRg3bhzGjRsHAPD398fbb7+NnTt3YsGCBdi5c6escT7++ONISEjA9u3bUVtbi48//hiA+a+5tEWMrc6ePYuZM2ca/pdzXVZWVmL37t148cUX4eDQ9W7fm3VpzTgtvW/KwdS61ev1OHToEJ5++mlbhmXWNg8JCUFSUhJ27dqF4OBgJCUl2UVcer0eJSUl2Lp1K1566SUkJydDp9PJHlers2fPIjQ01GTeWYI5cZ09exbh4eFITk7Gyy+/jN27d0Ov18se11NPPYVLly5h/fr1uHTpEry9vQ0NjFtRZAFvampCXFwcZs2ahenTp+Pf//43SktLsW7dOqxcuRIVFRXYsGEDqqq690tFR44cQXV1dbuDh6urKwYOHAigpSuzubnZ0LUpR5xeXl6QJAkDBgzAnDlzUFhYCKDljK5ti8Dc11xaa13W1NSgsLCw3c0hcq3Luro6xMbG4tFHH8XYsWMBAO7u7qirq0NzczOAljPk1puT2r4ytLm5GXV1dYbuSVvH2crS+6ZcTL2O9ebNm7h27Rp+97vfYeXKlSgoKMCOHTtQVFQka1xAyz7T2qUfERGBK1euWDUmc+Py9vbG1KlT4ejoiMGDB8PX1xclJSWyx9Xq3Llz7U7k5Y7rxIkTuPvuuwEAY8eORWNjY7d72KwRl7e3N9auXYsdO3bgscceA4B2vZfGKK6ACyGQnJwMPz8/3H///QCAkSNHIiUlBUlJSUhKSoKPjw/eeOMN3Hab+T/V9vXXX+PixYuIiYlpd6ZYVVVlOHsqLCyEXq+Hu7u7bHG2XtcWQiAjIwMjRowAAEyZMgXp6ekQQiA/Px+urq4mC7i1YgSAb775BpMnT4aTk5NhmBzrsqmpCbt27cLs2bMNSQu0nP0GBQXh/PnzAFruMG59LWhISAjS0tIAtNwJHRQUZFYL3BpxApbfN+Vk6nWsrq6ueOeddwzra8yYMVi/fr3V70I35zWxbe8pyczMtMnvoJsT17Rp05Cbmwug5amPkpKSTpdf5IgLaOne1+l0nU5I5YxLo9EY1tePP/6IxsZGo0/Z2G6euWAAAAGfSURBVDqu6upqQ0/AX//6V8yZM8fkdBX3GNnly5eRnp6OkSNHYt26dQCAxx577JaPARQVFeGrr74y/ATili1b8NNPP+HmzZtYsWIFVqxYgYkTJ2L//v0YNGgQXnnlFQD/fyTn/Pnz+PLLL6FSqeDk5ISYmBizDubWijMxMdHQyho1ahSef/55AMCkSZOQlZWFVatWwcnJCVFRUbLFCLScdXd8TEKOdXnu3Dn84x//QE1NjaEotz5698QTTyAhIQEffvgh/P39MXfuXADA3LlzsWfPHkRHR0OtViMmJsZkjNaM09L7ppxUKhWWLVuGbdu2Qa/XY86cORgxYgQOHz6MgIAA2d6tbk5cn332GTIzM6FSqaBWq83KMVvENWHCBFy8eBGrV6+Gg4MDnnzySaufyJm7Hc+cOYMZM2bYbL80J66nn34ae/fuxaeffgoAiIqKsnp85sR16dIlfPDBB5AkCePGjev0M73G8FWqRERECqS4LnQiIiJiASciIlIkFnAiIiIFYgEnIiJSIBZwIiIiBWIBJyIiUiAWcCIiIgX6L75zjDHvpYJeAAAAAElFTkSuQmCC\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": [
    "- 局所解から移動できていない！\n",
    "  - ステップサイズが小さすぎる？\n",
    "    - これ以上大きくすると、そもそものサンプリングができなくなる。\n",
    "    - 逆温度ごとに変えられれば良いが、現状のtf-pでは（筆者の知る限りにおいて）不可能\n",
    "  - リープフロッグの回数が不足？\n",
    "    - 100回に設定しても結果は変わらなかった。\n",
    "    - こちらも逆温度ごとに変える設定はできないはず。。\n",
    "  - 交換頻度が高すぎる？\n",
    "    - TF-P側での設定方法が不明\n",
    "- bijectorの使用は？\n",
    "  - TF-P 0.7時点において、`tfp.mcmc.TransformedTransitionKernel`を`ReplicaExchangeMC`に与えるとエラーを生じる。\n",
    "  - `StepSizeAdaptation`との併用と同様、他の手法との併用を考慮できていない様子。。。\n",
    "- TF-P 0.8(?)でのNUTS実装が待たれる。\n",
    "  - レプリカごとにリープフロッグ回数が調整されれば、移動距離の問題は解決するはず。"
   ]
  },
  {
   "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
}

## tfp-replica_exchange_easy.ipynb

      
Display the source blob

    
Display the rendered blob

    
    Raw
  

              tfp-replica_exchange_easy.ipynb
            
          
        Loading

      Sorry, something went wrong. Reload?
      Sorry, we cannot display this file.
      Sorry, this file is invalid so it cannot be displayed.
      
          Viewer requires iframe.